################################### ## DISEÑOS FACTORIALES GENERALES ## ## n-way AN(c)OVA ## ## Luis M. Carrascal ## ################################### ## Importar y Exportar datos ## datos <- read.table("clipboard", header=T, sep="\t", dec=".") ## importar desde portapapeles ## importar en MAC desde el portapapeles: datos <- read.table(pipe("pbpaste"), header=T, sep="\t", dec=".") ## write.table(datos, "nombre.txt", sep="\t", row.names=FALSE) ## lo guarda en el directorio por defecto de RStudio ## de Internet y ponderlo en uso en "environment" meconecto <- url("http://www.lmcarrascal.eu/cursos/nwayMANOVA.RData") ## abro una conexión con internet load(meconecto) ## cargo el archivo de la conexión close(meconecto); rm(meconecto) ## cierro la conexión ## para descargarse un archivo de internet a nuestro directorio por defecto de R; no pondremos en uso ese archivo download.file(url='http://www.lmcarrascal.eu/cursos/nwayMANOVA.RData', destfile ="para_practicar.RData") ## mirad: https://www.datacamp.com/community/tutorials/r-data-import-tutorial ## ÍNDICE (buscar con la "lupa" los caracteres "**** número ****") ## **** 1 **** Cargar paquetes de análisis ## **** 2 **** CONSTRUIMOS EL MODELO Y ALGUNOS MANEJOS DE DATOS ## **** 3 **** ANALISIS DE LOS RESIDUOS Y OTROS SUPUESTOS CANONICOS DE LOS MODELOS ## **** 4 **** TEST DE HOMOGENEIDAD DE VARIANZAS EN LOS NIVELES DE LOS FACTORES ## **** 5 **** PRIMEROS RESULTADOS ## **** 6 **** SIGNIFICACIONES GLOBALES DE EFECTOS ## **** 7 **** IMPORTANCIA DE LOS EFECTOS EN EL MODELO (EFFECT SIZES); PARTICIÓN DE LA VARIANZA ## **** 8 **** TESTS POST HOC ## **** 9 **** EXAMEN DEL EFECTO DE LA VIOLACIÓN DEL SUPUESTO DE HOMOCEDASTICIDAD ## **** 10 **** VALIDACIÓN CRUZADA DEL MODELO (CROSS-VALIDATION) ## **** 11 **** ESTIMAS ROBUSTAS ## **** 12 **** CONSTRUCCIÓN DE HIPÓTESIS NULAS PARA LAS F DE LOS EFECTOS EN EL MODELO ## **** 13 **** REDUCIR MANUALMENTE EL MODELO y MULTIMODEL INFERENCE ## **** 14 **** ESPECIFICANDO TABLAS DE CONTRASTES ## **** 15 **** TEST DE PARALELISMO EN EL ANCOVA ## **** 16 **** TRANSFORMACIÓN BOX-COX ## **** 17 **** ANÁLISIS CON EL PAQUETE "afex" ## **** 18 **** APROXIMACIONES BAYESIANAS ## **** 19 **** REFERENCIAS ## **** 1 **** ## install.packages(c("afex", "car", "DAAG", "fit.models", "heplots", "lme4", "lmerTest", "LMERConvenienceFunctions", "lmtest", "MASS", "moments", "MuMIn", "pbkrtest", "phia", "psych", "robust", "robustbase", "sandwich", "wle")) ## cargamos paquetes de análisis library(car) ## para obtener VIF's y usar Anova, bootCase, y transformación Box-Cox library(MASS) ## para usar dropterm y rlm library(MuMIn) ## para calcular AICc library(lmtest) ## para usar coeftest y correcciones por heterocedasticidad library(sandwich) ## para correcciones de heterocedasticidad usando vcovHC library(heplots) ## para la estima de los valores de partial eta2 library(moments) ## para obtener el sesgo, kurtosis y sus significaciones: kurtosis, anscombe.test, skewness, agostino.test library(psych) ## para construir tablas con describe library(DAAG) ## para hacer validaciones cruzadas del modelo usando cv.lm library(robust) ## para hacer estimas robustas usando rlm library(robustbase) ## para hacer estimas robustas usando lmrob library(phia) ## para tests a posteriori y análisis de interacciones library(emmeans) ## para tests a posteriori ## Las líneas de código siguientes funcionan mayoritariamente SIN NINGUNA MODIFICACIÓN, ## si nuestro data.frame se llama "datos" y hemos introducido la fórmula del modelo en "eqt". ## En las ocasiones en las que haga falta modificar las líneas de código con el nombre de nuestras variables, ## se indica (ANTES O A LA DERECHA): ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## **** 2 **** ## tipos de variables que tenemos str(datos) ## para conocer la naturaleza de las variables en el dataframe datos class(datos$CODZONA) ## para conocer la naturaleza de una variable ## trabajemos con factores datos$FZONA <- as.factor(datos$CODZONA) ## convertimos una variable numérica en un factor (al que se aplican contrastes "contr.sum") class(datos$DOMINANCIA) ## para conocer la naturaleza de una variable ## convertimos una variable continua (con valores 1 2 3 y 4) a otra con valores (0 1) ## a la condición "<3" se le asigna el valor 0, y a lo demás el valor 1 datos$DOMINANCIA2 <- ifelse(datos$DOMINANCIA<3, 0, 1) datos$FDOMINANCIA <- as.ordered(datos$DOMINANCIA) ## convertimos una variable numérica en un factor ordinal "ordered" (al que se aplican contrastes "contr.poly") ## establecemos una multinomial ordinal, poniendo niveles de corte a la variable fuente (DOMINANCIA) ## en <=1 (de 0 a 1), 1-2, 2-3 y 3-4 datos$FDOMINANCIA2 <- cut(datos$DOMINANCIA, c(0,1,2,3,4), labels=c("MUY SUBORDINADO", "SUBORDINADO", "DOMINANTE", "MUY DOMINANTE")) ## gestión de "missing cells" o valores perdidos ## http://www.ats.ucla.edu/stat/r/faq/missing.htm ## podemos añadir dentro del modelo lm(...) los argumentos: ## na.action = na.omit ... para omitir los valores perdidos de las unidades muestrales a través de todas las variables ## na.action = na.exclude ... para omitir los valores perdidos pero manteniendo sus posiciones para los residuos y valores predichos ## válido también para modelos glm, glm.nb, lmer, glmer v1 <- c(1, 2, 3, NA, 5, 6, NA, 7, 7, 9, NA) v2 <- c(NA, 2, NA, 4, 5, 6, 7, NA, 8, NA, 10) dd <- data.frame(v1, v2) ## si el dataframe datos tiene valores NA podemos convertirlo en otro sin ellos dd2 <- na.omit(dd) dim(dd); dim(dd2) ## para estadísticos básicos con valores NA mean(v1) mean(na.omit(v1)) cor(dd$v1, dd$v2) cor(dd2$v1, dd2$v2) plot(dd$v1, dd$v2) plot(dd2$v1, dd2$v2) summary(lm(v1~v2, data=dd)) summary(lm(v1~v2, data=na.omit(dd))) summary(lm(v1~v2, data=dd, na.action=na.omit)) ## borramos todo menos el objedo dataframe "datos" rm(list=setdiff(ls(), "datos")) ## eliminamos todos los plots dev.off() ## CONSTRUIMOS EL MODELO ### ¡¡ ojo !! poned aquí VUESTROS DATOS dentro de as.formula(...) eqt <- as.formula(MUSCULO~TARSO+ZONA*SEXO*EDAD) ## primero definimos la función de efectos (respuesta ~ predictores) ## ## antes cargamos la siguiente línea de código para obtener los mismos resultados ## que STATISTICA y SPSS utilizando type III SS ## As the treatment factors are unordered factors, R will use contr.sum in order to construct a contrast matrix of the apropriate order. options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## modelo <- lm(eqt, data=datos) ## tras haber corrido options(contrasts(...)), es igual que lo siguiente ## ... y sin haber corrido previamente potions(contrasts(...)) modelo <- lm(eqt, data=datos, contrasts=c(factor="contr.sum", ordered="contr.poly")) ## ## si quisiesemos más detalles en los contrastes, tendríamos que definirlos del siguiente modo: ### ¡¡ ojo !! poned aquí VUESTROS DATOS en list(...) modelo <- lm(MUSCULO~TARSO+ZONA*SEXO*EDAD, data=datos, contrasts=list(ZONA=contr.sum, SEXO=contr.sum, EDAD=contr.sum)) ## **** 3 **** ## EXPLORACIÓN VISUAL GLOBAL DE LOS RESIDUOS ## NORMALIDAD ## ambas formas son equivalentes residuals(modelo) modelo$residuals ## hist(residuals(modelo), main="residuos del modelo") hist(residuals(modelo), main="residuos del modelo", col="yellow", freq=F) lines(density(residuals(modelo)), col="red", lwd=2) hist(residuals(modelo), density=5, freq=FALSE, main="residuos del modelo") curve(dnorm(x, mean=mean(residuals(modelo)), sd=sd(residuals(modelo))), col="red", lwd=2, add=TRUE, yaxt="n") qqnorm(residuals(modelo), main="residuos del modelo") qqline(residuals(modelo), col="red", lwd=2) qqnorm(residuals(modelo), main="residuos del modelo", datax=TRUE, xlab="Expected Normal Value", ylab="Residuals") ## datax=TRUE para cambiar el orden de los ejes qqline(residuals(modelo), col="red", lwd=2, datax=TRUE) ## ## valoramos cómo podrían ser unos residuos Normales con una muestra similar a la nuestra ## repetid estas líneas varias veces para ver diferentes posibilidades de histogramas y qqplots ## seleccionad todo el bloque de las siguientes líneas y dad click a la flechita verde de Run media_residuos <- mean(residuals(modelo)) sd_residuos <- sd(residuals(modelo)) n_residuos <- length(residuals(modelo)) residuos_alzar_normales <- rnorm(n=n_residuos, mean=media_residuos, sd=sd_residuos) par(mfcol=c(1,1)) ## fija un sólo panel gráfico (de Plots) par(mfcol=c(1,2)) ## fija dos paneles según 1 fila y 2 columnas hist(residuos_alzar_normales) qqnorm(residuos_alzar_normales); qqline(residuos_alzar_normales) par(mfcol=c(1,1)) ## fija un sólo panel gráfico (de Plots) ## ## HOMOCEDASTICIDAD plot(modelo$fitted.values, residuals(modelo), main="¿HAY HETEROCEDASTICIDAD?") abline(h=0, col="red") ## traza una línea horizontal (h) por el Y=0 ## si hay heterocedasticidad de los residuos, que no haya una clara relación en lo siguiente ## si la relación es positiva, aumenta el error de tipo I ## si la relación es negativa, aumenta el error de tipo II ### ¡¡ ojo !! poned aquí VUESTROS DATOS medias.originales <- aggregate(MUSCULO ~ ZONA*SEXO*EDAD, data=datos, FUN=mean)[4] varianzas.residuos <- aggregate(residuals(modelo) ~ ZONA*SEXO*EDAD, data=datos, FUN=var)[4] cor.media.varianzaresids <- round(cor(medias.originales[,1], varianzas.residuos[,1]), 3) plot(medias.originales[,1], varianzas.residuos[,1], ylab="VARIANZAS DE LOS RESIDUOS DE LAS CELDAS", xlab="MEDIAS DE LAS CELDAS", main=c("CORRELACION (r) ENTRE CELDAS DE LOS FACTORES",cor.media.varianzaresids)) abline(lm((varianzas.residuos[,1]) ~ medias.originales[,1]), col="red", lwd=2) ## ## PLOTS MÚLTIPLES ## lo mismo de antes pero de una sola vez y ... más; saca los cuatro paneles par(mfcol=c(1,1)) ## fija un sólo panel gráfico par(mfcol=c(2,2)) ## fija cuatro paneles según 2 columnas y 2 filas plot(modelo, c(1:2,4,5)) par(mfcol=c(1,1)) ## volvemos al modo gráfico de un solo panel ## en otras ocasiones puede ser interesante (y en una sola línea): par(mfcol=c(1,1)); par(mfcol=c(1,2)); plot(modelo, 1:2); par(mfcol=c(1,1)) ## ## Parametrización de la NORMALIDAD DE LOS RESIDUOS ## kurtosis(residuals(modelo)) ## kurtosis de Pearson, Ho = 3 anscombe.test(residuals(modelo)) ## skewness(residuals(modelo)) ## sesgo Ho = 0 agostino.test(residuals(modelo)) ## shapiro.test(residuals(modelo)) ## test de shapiro de normalidad de residuos ## EXPLORACIÓN DE DATOS RESIDUALES DE MANERA INDIVIDUALIZADA por cada unidad muestral ## LEVERAGE: https://en.wikipedia.org/wiki/Leverage_(statistics) ## https://stats.stackexchange.com/questions/208242/hat-matrix-and-leverages-in-classical-multiple-regression) ## niveles críticos: 2*(g.l. del modelo)/(Número de datos) leverage.modelo <- hat(model.matrix(lm(formula(modelo), modelo$model))) ## otra opción es leverage.modelo <- hatvalues(modelo) ## si no vemos la línea roja es que no hay problemas con puntos con gran leverage plot(leverage.modelo) abline(h=2*(length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals), col="red") ## ## DFFITS (https://en.wikipedia.org/wiki/DFFITS) ## niveles críticos: 2*raiz((g.l. del modelo)/(Número de datos)) ## si no vemos la línea roja es que no hay problemas con puntos con grandes valores de dffits plot(dffits(modelo)) abline(h=2*((length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals))^0.5, col="red") abline(h=-2*((length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals))^0.5, col="red") identify(dffits(modelo)) ## terminad dando clik en Finish ## plot(modelo$residuals, dffits(modelo)) abline(lm(dffits(modelo) ~ modelo$residuals), col="red") ## ## DISTANCIA DE COOK (https://en.wikipedia.org/wiki/Cook%27s_distance) ## niveles críticos: preocupante > 4/(número de datos); muy preocupante > 1 plot(cooks.distance(modelo)) abline(h=4/length(modelo$residuals), col="red") identify(cooks.distance(modelo)) ## terminad dando clik en Finish ## ## PLOTS COMBINADOS plot(dffits(modelo)~leverage.modelo, col="darkgreen") abline(h=2*((length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals))^0.5, col="red") abline(h=-2*((length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals))^0.5, col="red") abline(v=2*(length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals), col="red") identify(dffits(modelo)~leverage.modelo) ## terminad dando clik en Finish ## plot(cooks.distance(modelo)~dffits(modelo), col="darkgreen") abline(h=4/length(modelo$residuals), col="red") abline(v=2*((length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals))^0.5, col="red") abline(v=-2*((length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals))^0.5, col="red") identify(cooks.distance(modelo)~dffits(modelo)) ## terminad dando clik en Finish ## plot(cooks.distance(modelo)~leverage.modelo, col="darkgreen") abline(h=4/length(modelo$residuals), col="red") abline(v=2*(length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals), col="red") identify(cooks.distance(modelo)~leverage.modelo) ## terminad dando clik en Finish ## ## RESIDUOS ESTUDENTIZADOS (https://en.wikipedia.org/wiki/Studentized_residual) ## Los residuos studentizados nos pueden dar una indicación de qué datos es probable que ## no pertenezcan a la población. El valor crítico lo define la t de Student teniendo en cuenta ## los g.l. (error) del modelo. Valores mayores o menores de 2.0 (aprox p=0.05) o 2.7 (aprox p=0.01) son peligrosos. ## podemos calcular la t de Student crítica para una alfa determinada (0.05 ó 0.01) y los df de nuestro modelo con qt(...): t_05 <- qt(0.05/2, df=modelo$df.residual, lower.tail=FALSE) t_01 <- qt(0.01/2, df=modelo$df.residual, lower.tail=FALSE) ## plot(studres(modelo) ~ modelo$fitted.values) abline(h=t_05, col="red", lwd=1); abline(h=-t_05, col="red", lwd=1) abline(h=t_01, col="red", lwd=3); abline(h=-t_01, col="red", lwd=3) identify(studres(modelo) ~ modelo$fitted.values) ## terminad dando clik en Finish (boton sup-decho panel de plots) ## ## INDEPENDENCIA ENTRE LAS PREDICTORAS - VIF (https://en.wikipedia.org/wiki/Variance_inflation_factor) ## VIF = 1 / (1 - R2), de cada predictora por las restantes ## GVIF or VIF specifically indicate the magnitude of the inflation in the standard errors ## associated with a particular beta weight that is due to multicollinearity ## la raiz cuadrada (sqrt(...)) de un valor VIF o GVIF el número de veces que se inflan los errores standard de esa predictora vif(modelo) sqrt(vif(modelo)) ## ## REPRESENTACIÓN GRÁFICA DE LAS RELACIONES ENTRE VARIABLES ## Mirad que no haya "missing cells" en los factores ## si hemos introducido la ecuación de efectos manualmente dentro de lm(...) en vez de con eqt, ## podemos extraerla con as.formula(modelo$call$formula) ## panel con valores de correlaciones en valor absoluto (pueden ser + ó -) panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...) { usr <- par("usr") on.exit(par(usr)) par(usr = c(0, 1, 0, 1)) r <- abs(cor(x, y, use="complete.obs")) txt <- format(c(r, 0.123456789), digits=digits)[1] txt <- paste(prefix, txt, sep="") if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt) text(0.5, 0.5, txt, col="blue", cex = cex.cor * (1 + r) / 1) ## cambiar el tamaño de la letra modificando "1" en el denominador de "cex = cex.cor * (1 + r) / 1" } ## ampliarla con Zoom dentro de la solapa Plots pairs(eqt, data=datos, cex.labels=2, pch="o", lower.panel = panel.cor) par(mfcol=c(1,1)) ## **** 4 **** ## TEST DE HOMOGENEIDAD DE VARIANZAS a través de los niveles cruzados de los factores del diseño ## varianzas de las diferentes celdas de interacción entre factores ### ¡¡ ojo !! poned aquí VUESTROS DATOS aggregate(MUSCULO ~ ZONA*SEXO*EDAD, data=datos, FUN=var) aggregate(MUSCULO ~ ZONA*SEXO*EDAD, data=datos, FUN=mean) aggregate(MUSCULO ~ ZONA*SEXO*EDAD, data=datos, FUN=length) ## ## disponemos de diferentes tests ## The Levene test works well in the ANOVA framework, providing there are small to moderate deviations from the normality. ## In this case, it outperfoms the Bartlett test. If the distribution are nearly normal, however, the Bartlett test is better. ## The Fligner-Killeen test is to be prefered in case of strong departure from the normality (to which the Bartlett test is sensible). ## En el test de Levene otra opción es center=mean (e.g., la del programa STATISTICA) ## consultad: http://www.cookbook-r.com/Statistical_analysis/Homogeneity_of_variance/ ## para un solo factor leveneTest(MUSCULO~ZONA, data=datos, center=median) ### ¡¡ ojo !! poned aquí VUESTROS DATOS leveneTest(MUSCULO~ZONA, data=datos, center=mean) ### ¡¡ ojo !! poned aquí VUESTROS DATOS bartlett.test(MUSCULO~ZONA, data=datos) ### ¡¡ ojo !! poned aquí VUESTROS DATOS fligner.test(MUSCULO~ZONA, data=datos) ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## para más de un factor leveneTest(MUSCULO~ZONA:SEXO:EDAD, data=datos, center=median) ### ¡¡ ojo !! poned aquí VUESTROS DATOS leveneTest(MUSCULO~ZONA:SEXO:EDAD, data=datos, center=mean) ### ¡¡ ojo !! poned aquí VUESTROS DATOS bartlett.test(MUSCULO~interaction(ZONA,SEXO,EDAD), data=datos) ### ¡¡ ojo !! poned aquí VUESTROS DATOS fligner.test(MUSCULO~interaction(ZONA,SEXO,EDAD), data=datos) ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## ## relación entre medias y varianzas ## si hay heterocedasticidad de los residuos, que no haya una clara relación en lo siguiente ## si la relación es positiva, aumenta el error de tipo I ## si la relación es negativa, aumenta el error de tipo II ### ¡¡ ojo !! poned aquí VUESTROS DATOS medias.originales <- aggregate(MUSCULO ~ ZONA*SEXO*EDAD, data=datos, FUN=mean)[4] varianzas.originales <- aggregate(MUSCULO ~ ZONA*SEXO*EDAD, data=datos, FUN=var)[4] cor.media.varianza <- round(cor(medias.originales[,1], varianzas.originales[,1]), 3) plot(medias.originales[,1], varianzas.originales[,1], ylab="VARIANZAS DE LAS CELDAS", xlab="MEDIAS DE LAS CELDAS", main=c("CELDAS DE LOS FACTORES (r)",cor.media.varianza)) abline(lm((varianzas.originales[,1]) ~ medias.originales[,1]), col="red", lwd=2) par(mfcol=c(1,1)) ## **** 5 **** ## PRIMEROS RESULTADOS summary(modelo) ## OMNIBUS TEST - TEST DE SIGNIFICACIÓN GLOBAL DEL MODELO ## si la última línea de summary(modelo) no da un resultados significativo ... NO MIRAREMOS LOS DEMÁS RESULTADOS ## otros modos de efectuar el omnibus test de todo el modelo son: ## test de Wald waldtest(modelo, test="F") ## comparación del modelo de interés con el modelo nulo modelo.nulo <- update(modelo, . ~ 1) anova(modelo.nulo, modelo) ## ## vemos los COEFICIENTES (usando los contrastes especificados): valores (Estimates), errores estándard, t y P ## da igual que coeftest(modelo) y que coeftest(modelo, vcovHC(modelo, type="const")) coeftest(modelo) coeftest(modelo, vcovHC(modelo, type="const")) ## ## COEFCIENTES DE REGRESIÓN STANDARIZADOS ## zeta-standarizamos (a MEDIA=0 Y SD=1) std.coef(modelo, partial.sd=F) ## ## VARIANZA TOTAL EXPLICADA ## en la penúltima línea de summary(modelo), ó con: summary(modelo)$r.squared summary(modelo)$adj.r.squared ## algunos parámetros globales del modelo relacionados con la varianza (sumas de cuadrados, SS) ## los datos seleccionados para correr nuestro modelo se encuentran en modelo$model (la primera columna es la Respuesta) ## suma de cuadrados (SStotal) de la variable respuesta (varianza original) SStotal <- sum((modelo$model[,1]-mean(modelo$model[,1]))^2) ## suma de cuadrados error (SSerror) o residual del modelo SSerror <- sum(residuals(modelo)^2) ## suma de cuadrados del modelo (SSmodelo) SSmodelo <- SStotal-SSerror ## R2 del modelo estimada manualmente R2 <- summary(modelo)$r.squared ## que es lo mismo que calcular SSmodelo/SStotal print(c(SStotal, SSmodelo, SSerror, R2), digits=4) print(c("R2 % =",round(R2*100, 2)), quote=FALSE) ## RELACIÓN ENTRE VALORES OBSERVADOS Y PREDICHOS titulo <- paste("R2 (%) =",round(R2*100,1)) par(mfcol=c(1,1)) plot(fitted(modelo), modelo$model[,1], xlab="PREDICHO POR EL MODELO", ylab="VARIABLE RESPUESTA", main=titulo) ## modelo$model[,1] es la variable respuesta original abline(lm(modelo$model[,1]~fitted(modelo)), col="red", lwd=2) ## ## EXPLORACIÓN VISUAL DE LOS EFECTOS DE LOS FACTORES SOBRE LA RESPUESTA ## medias de las diferentes celdas de interacción entre factores, sin tener en cuenta el efecto de las covariantes ### ¡¡ ojo !! poned aquí VUESTROS DATOS (en el nombre de la respuesta y los factores) aggregate(MUSCULO ~ ZONA*SEXO*EDAD, data=datos, FUN=mean) ## ## resultado numérico de los valores ajustados (por la(s) covariante(s)) y el standard error interactionMeans(modelo) ## otras opciones, con iguales resultados, son: ### ¡¡ ojo !! poned aquí VUESTROS DATOS interactionMeans(modelo, pairwise=c("ZONA", "SEXO", "EDAD")) interactionMeans(modelo, factors=c("ZONA", "SEXO", "EDAD")) ## y restringiendo los factores interactionMeans(modelo, factors=c("ZONA", "SEXO")) ## otra forma es usando el paquete {emmeans} ### ¡¡ ojo !! poned aquí VUESTROS DATOS emmeans(modelo, c("ZONA", "SEXO", "EDAD")) ## resultado gráfico con una figura y múltiples paneles; media +/- std_error ## si se solapan las cajas de los factores con el panel izquierdo ... ampliad el ancho de Plots plot(interactionMeans(modelo), legend.margin=0.2) ## ## VISUALIZACION DE EFECTOS PARCIALES, incluyendo covariante(s) ## "Partial residual plots" (https://en.wikipedia.org/wiki/Partial_residual_plot) ## They are most commonly used to identify the nature of the relationship between Y and Xi (given the effect of the other independent variables in the model). ## no admite interacciones; tendremos que recalcular nuestro modelo sin ellas ### ¡¡ ojo !! poned aquí VUESTROS DATOS modelo.noint <- lm(MUSCULO~TARSO+ZONA+SEXO+EDAD, data=datos) ## crPlots(modelo.noint, smooth=FALSE, pch=15, lwd=2, col.lines="blue", grid=FALSE, ylab="partial residuals", main="TITULA COMO QUIERAS") ## esta otra línea presenta los plots de uno en uno; hay que teclear cada vez crPlots(modelo.noint, smooth=FALSE, pch=15, lwd=2, col.lines="blue", grid=FALSE, ylab="partial residuals", main="TITULA COMO QUIERAS", layout=c(1, 1)) ## en esta otra versión podemos identificar posibles desvíos de los efectos lienales ## (modificado con el argumento por defecto smooth=TRUE) crPlots(modelo.noint) par(mfcol=c(1,1)) ## **** 6 **** ## SIGNIFICACIONES GLOBALES DE EFECTOS, USANDO SS DE TIPO I, II Y III ## The standard R "anova" function calculates sequential ("type-I") tests. ## These rarely test interesting hypotheses in unbalanced designs ## Type-II tests are calculated according to the principle of marginality, testing each term after all others, ## except ignoring the term's higher-order relatives; ## so-called type-III tests violate marginality, testing each term in the model after all of the others. ## para significación SS type-I ## en la mayoría de las veces esto es irrelevante y de nula utilidad ## depende del orden de entrada de las variables anova(modelo) ## ## para significación SS type-III ## varios modos convergentes (las dos siguinetes líneas producen los mismos resultados) drop1(modelo, ~., test="F") dropterm(modelo, ~., test="F", sorted=FALSE) ## similar al anterior y con muestra del intercepto y SS error de tipo III Anova(modelo, type=3, test="F") ## para significación SS type-II ## caso conveniente para modelos con muuuchas interacciones y ... poco tamaño muestral Anova(modelo, type=2, test="F") ## **** 7 **** ## IMPORTANCIA DE LOS EFECTOS EN EL MODELO (EFFECT SIZES); PARTICIÓN DE LA VARIANZA ## con SS de tipo III SSerror <- sum(residuals(modelo)^2) SStotal <- sum((modelo$model[,1]-mean(modelo$model[,1]))^2) SSmodelo <- SStotal - SSerror tabla.ss <- as.data.frame(Anova(modelo, type=3, test="F")) tabla.ss <- tabla.ss[-1,] names(tabla.ss) <- c("SS", "df", "F_value", "P_value") eta2 <- tabla.ss[,1]/SStotal partialeta2 <- tabla.ss[,1]/(tabla.ss[,1]+SSerror) tabla.ss <- data.frame(eta2, partialeta2, tabla.ss) ## relación entre eta2 y partial_eta2 plot(tabla.ss$eta2, tabla.ss$partialeta2, xlab="eta2 (% de SS)", ylab="eta2 parcial") abline(a=0, b=1, col="red") ## concomitancias y efectos exclusivos exclusivos <- sum(tabla.ss[c(1:(dim(tabla.ss)[1]-1)),1]) concomitancias <- SSmodelo/SStotal - exclusivos ## eta2 : la proporción de la varianza explicada por los efectos ## partial eta2 es otra medida relativa; estimada como SSefecto / (SSefecto + SSresidual) round(tabla.ss, 4) print(paste("Proporción de la varianza exclusiva de los efectos =", round(exclusivos, 3)), quote=F) print(paste("Proporción de la varianza compartida por todos los efectos =", round(concomitancias, 3)), quote=F) print(paste("Proporción de la varianza explicada por el Modelo (R2) =", round(SSmodelo/SStotal, 3)), quote=F) ## **** 8 **** ## TESTS POST HOC ## ## algunas figuras ### ¡¡ ojo !! poned aquí VUESTROS DATOS boxplot(MUSCULO ~ ZONA * SEXO * EDAD, data=datos) ## no acepta covariantes; si las hay, introducir la fórmula a mano, en vez de con as.formula(modelo$call$formula) dotplot(MUSCULO~ZONA|EDAD, data=datos) ## sólo es posible con dos factores que introducimos manualmente interaction.plot(datos$ZONA,datos$EDAD,datos$MUSCULO) ## en la forma de (un factor, otro factor, variable respuesta) ## ## la más potente, con el paquete {emmeans}; consultad con [vignette("models", "emmeans")] los modelos que puede usar ## sí admite covariantes; por tanto, sí controla por ellas ## modo usando "," entre factores para efectuar test de tukey, scheffe, bonferroni, none, fdr, holm, hommel, hochberg, ... ## (mirad "P-value adjustments" en el Help de summary.emmGrid) pairs(emmeans(modelo, c("ZONA", "SEXO", "EDAD")), adjust="tukey") ### ¡¡ ojo !! poned aquí VUESTROS DATOS pairs(emmeans(modelo, c("ZONA", "SEXO", "EDAD")), adjust="holm") ### ¡¡ ojo !! poned aquí VUESTROS DATOS pairs(emmeans(modelo, c("ZONA", "SEXO", "EDAD")), adjust="scheffe") ### ¡¡ ojo !! poned aquí VUESTROS DATOS pairs(emmeans(modelo, c("ZONA", "SEXO", "EDAD")), adjust="fdr") ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## ## otro forma con Tukey Honest Significant Differences; construimos antes un modelo [aov] que incluimos en el comando TukeyHSD ## sólo es para factores; quitar las covariantes si el modelo las tiene, introduciendo la fórmula tras el paréntesis de aov(, en sustitución de as.formula(...). TukeyHSD(aov(MUSCULO ~ ZONA : EDAD : SEXO, data=datos)) ### ¡¡ ojo !! poned aquí VUESTROS DATOS plot(TukeyHSD(aov(MUSCULO ~ ZONA : EDAD : SEXO, data=datos))) ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## ## otra forma con la versión de los "pairwise t-tests" con correcciones de significación (primero la variable respuesta en x, y luego los factores en g) ## no admite covariantes; por tanto, no controla por ellas ## ## el argumento adjustment puede ser "holm" (por defecto), "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none" ## (mirad p.adjust en el Help y las referencias al final) ## method="fdr": false discovery rate is designed to control the expected proportion of "discoveries" (rejected null hypotheses) that are false (incorrect rejections). It has greater power, at the cost of increased numbers of Type I errors ## method="holm": familywise error rate (the probability of making one or more false discoveries, or type I errors when performing multiple hypotheses tests) ### ¡¡ ojo !! poned aquí VUESTROS DATOS pairwise.t.test(x=datos$MUSCULO, g=datos$DOMINANCIA, p.adj="fdr") ## para un factor pairwise.t.test(x=datos$MUSCULO, g=datos$ZONA:datos$EDAD, p.adj="fdr") ## para la interacción de dos factores pairwise.t.test(x=datos$MUSCULO, g=datos$ZONA:datos$EDAD:datos$SEXO, p.adj="fdr") ## para la interacción de tres factores pairwise.t.test(x=datos$MUSCULO, g=datos$ZONA:datos$EDAD:datos$SEXO, p.adj="holm") ## produce un resultado muy parecido a TukeyHSD ## ## Opción de tests de interacciones con correcciones de significación utilizando el paquete "phia". Admite covariantes. ## el argumento adjustment puede ser "holm" (por defecto), "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none". ## mirad la ayuda de "p.adjust {stats}" ## Para una selección de tests: testInteractions(modelo, adjustment="holm", fixed="...el factor fijado...", across="...el factor de interés...") ## el resultado bajo "Value" da la diferencia entre los niveles de referencia del factor incluido en across. ### ¡¡ ojo !! poned aquí VUESTROS DATOS testInteractions(modelo, adjustment="holm", fixed="ZONA", across="SEXO") ## tests entre SEXOs (across) dentro de ZONA testInteractions(modelo, adjustment="holm", fixed=c("ZONA", "EDAD"), across="SEXO") ## tests entre SEXOs (across) dentro de las celdas ZONA:EDAD ## Para el test manual de interacciones (da lo mismo que en summary(modelo)), usamos pairwise="...factor de interés..." testInteractions(modelo, adjustment="holm") ## test de la interacción ZONA:SEXO:EDAD. Lo mismo que lo siguiente. testInteractions(modelo, adjustment="holm", pairwise=c("ZONA", "SEXO", "EDAD")) ## test de la interacción ZONA:SEXO:EDAD. Lo mismo que lo anterior. testInteractions(modelo, adjustment="holm", pairwise=c("ZONA", "SEXO")) ## test de la interacción ZONA:SEXO ## **** 9 **** ## EXAMEN DEL EFECTO DE LA VIOLACIÓN DEL SUPUESTO DE HOMOCEDASTICIDAD ## otras opciones son vcov=sandwich (que es type="HC0"), "const", "HC0", "HC1", "HC2", "HC4", "HC4m", "HC5" ## HC3 (que es igual que vcovHC) provides the best performance in small samples as it gives less weight to influential observations ## HC4 further improves small sample performance, especially in the presence of influential observations ## it is tailored to take into account the effect of leverage points in the design matrix on associated quasi-t tests ## HC4m testing inference is more reliable than that based on HC4 standard errors under both normal and nonnormal random errors ## HC5 standard errors deliver much more reliable inference in finite samples than those based on the most commonly used heteroskedasticity-consistent standard errors, in particular HC0- and HC3-based tests. coeftest(modelo, df=modelo$df.residual, vcov=vcovHC) ## es igual que coeftest(modelo, vcovHC(modelo, type="HC3")) coeftest(modelo, df=modelo$df.residual, vcovHC(modelo, type="HC4m")) coeftest(modelo, df=modelo$df.residual, vcovHC(modelo, type="HC5")) coeftest(modelo, df=modelo$df.residual, vcovHC(modelo, type="const")) ## comparadlo con la salida inicial sin efectuar correcciones ## tabla1 <- as.table(coeftest(modelo, df=modelo$df.residual, vcovHC(modelo, type="const"))) tabla2 <- as.table(coeftest(modelo, df=modelo$df.residual, vcovHC(modelo, type="HC4m"))) tabla.efecto.heterocedst <- data.frame(tabla1[,2], tabla2[,2], tabla1[,4], tabla2[,4]) colnames(tabla.efecto.heterocedst) <- c("se.original", "se.HC4m", "P.original", "P.HC4m") round(tabla.efecto.heterocedst, 5) ## ## otra opción, para la estima de la significación global de los factores es mediante el argumento Anova(..., white.adjust) ## white.adjust puede ser: "hc0", "hc1", "hc2", "hc3", "hc4" Anova(modelo, type=3, test="F", white.adjust="hc4") ## **** 10 **** ## VALIDACIÓN CRUZADA DEL MODELO (CROSS-VALIDATION) ## la tabla de Analysis of Variance Table se corresponde con SS tipo-I; por tanto, ¡¡no utilizadla!! ## los argumentos de cv.lm(...) son (data.frame, modelo del tipo lm, m número de partes o folds, valor numérico para diferentes particiones) ## cambiad el valor de seed para diferentes procesos de validación cruzada; ## aquí determinado con el valor de la user-written function salida.cv(...) salida.cv <- function(x){ cv.linmod <- cv.lm(datos, modelo, m=5, seed=x, plotit=FALSE) ## aquí en seed=x está el parámetro x que modificaremos R2cvmod <- cor(cv.linmod$cvpred, cv.linmod$Predicted)^2 R2obscv <- cor(cv.linmod$cvpred, modelo$model[,1])^2 encabezado <- paste("R2 model-crossval =", round(R2cvmod*100, 2)," %") Pobscv <- round(cor.test(cv.linmod$cvpred, modelo$model[,1], alternative="greater")$p.value,5) encabezados <- paste("R2 observed-crossval =", round(R2obscv*100, 2)," % ; P = ", Pobscv) plot(cv.linmod$cvpred, modelo$model[,1], main=encabezados, xlab="PREDICHOS POR CROSSVALIDATION", ylab="VALORES OBSERVADOS", col="blue") abline(a=0, b=1, col="red", lwd=2) } ## introducir en el paréntesis el valor semilla (seed) salida.cv(1) ## o correr este bucle con valores de seed de 1 a 10; otros podrían ser "(i in 23:32)" for (i in 1:10) { salida.cv(i) } ## otra opción más general y elaborada a partir del uso de comandos de uso frecuente ## VALIDACIÓN CRUZADA DE MODELOS LM Gausianos ## ## número de particiones de datos ## si definimos tantos N_fold como tamaño muestral hay (dim(datos)[1]), ## entonces el proceso se denomina Jackknife resampling (https://en.wikipedia.org/wiki/Jackknife_resampling) ## y sólo es posible un resultado de validación cruzada N_fold <- 10 N_fold <- dim(datos)[1] ## para jackknife resampling ## ## ejecución de un proceso de validación cruzada - INICIO ## mirad la línea #231 para cambiar de modelo glm(...) a glm.nb(...) ## preparando datos datos2 <- modelo$model datos2$azar <- runif(n=dim(datos2)[1], min=1, max=1000) datos2 <- datos2[order(datos2$azar),] datos2$azar <- NULL datos2$id <- rep(c(1:N_fold), length.out=dim(datos2)[1]) predichos <- 99999; observados <- 99999 ; id_cv <- 0 matriz_cv <- data.frame(predichos, observados, id_cv) ## ## bucle de validación cruzada for (i in 1:N_fold) { datos3 <- subset(datos2, id !=i) respuesta <- datos3[,1] datos3 <- datos3[, -1] mod.predict <- lm(respuesta~.-id, data=datos3) ## para el caso de modelos binomiales negativos; sustituir la siguiente línea por la anterior ## mod.predict <- glm.nb(respuesta~.-id, data=datos3, link=log) datos4 <- subset(datos2, id==i) predichos <- predict(mod.predict, newdata=datos4, type="response") observados <- datos4[,1] id_cv <- datos4$id datos_cv <- data.frame(predichos, observados, id_cv) matriz_cv <- rbind(matriz_cv, datos_cv) } matriz_cv <- matriz_cv[-1,] ## FINAL del proceso de validación cruzada ## ## RESULTADOS de la validación cruzada para RESPUESTAS CONTINUAS r.obs_pred <- round(cor(matriz_cv$predichos, matriz_cv$observados), 3) resultado <- paste("R2 observado - predicho =", r.obs_pred^2) plot(matriz_cv$predichos, matriz_cv$observados, xlab="VALORES PREDICHOS", ylab="VALORES OBSERVADOS", cex.lab=1) abline(lm(matriz_cv$observados ~ matriz_cv$predichos), col="red", lwd=2) title(main=resultado, cex.main=2, col.main="red") ## podremos repetir este proceso varias veces para valorar los posibles resultados ## con distintos procesos de aleatorización ## seleccionaremos todo el bloque de líneas y daremos click al botón [-> Run] ## **** 11 **** ## ESTIMAS ROBUSTAS - teniendo en cuenta el leverage y la distancia de Cook de las observaciones (casos) ## teniendo en cuenta el leverage y la distancia de Cook de las observaciones (casos) ## http://www.ats.ucla.edu/stat/r/dae/rreg.htm ## http://quantoid.net/files/reg3/lecture10_2015.pdf ## http://www.dst.unive.it/rsr/BelVenTutorial.pdf ## ## aproximación robusta que estima nuevos coeficientes y errores standard (MASS:rlm) ## proporciona una gran flexibilidad de resultados con tablas de SS para efectos y estimas asumiendo violación de homocedasticidad de residuos ## Fit a linear model by robust regression using an M estimator. ## Fitting is done by iterated re-weighted least squares (IWLS or IRwLS). ## Selecting method = "MM" selects a specific set of options which ensures that the estimator has a high breakdown point. ## Case weights are not supported for method = "MM". ## Psi functions are supplied for the Huber, Hampel and Tukey bisquare proposals as psi.huber, psi.hampel and psi.bisquare. ## Huber's corresponds to a convex optimization problem and gives a unique solution (up to collinearity). ## The other two will have multiple local minima, and a good starting point is desirable. ## Different functions have advantages and drawbacks. Huber weights can have difficulties with severe outliers ## (Huber's proposals have the drawback that large outliers are not down weighted to zero). ## This can be achieved with the Tukey's bisquare proposal, but bisquare weights can have difficulties converging, ## or may yield multiple solutions. Tukey's bisquare proposal gives extreme observations zero weights. ## In general, the standard error of these estimates are slightly smaller than in the Huber fit but the results are qualitatively similar. ## Mirad: https://goo.gl/l0WCWz https://goo.gl/Aemvfe https://goo.gl/08Bxt5 ## previamente hemos creado un modelo lm denominado "modelo" (usando eqt y datos) que ahora queremos valorar robustamente modelo.r31 <- rlm(eqt, data=datos, method="MM", maxit=100) F_r31 <- waldtest(modelo.r31)$F[2] p_r31 <- 1-pf(F_r31, df1=dim(modelo$model)[1]-modelo$df.residual-1, df2=modelo$df.residual) print(c("F del modelo robusto =", round(F_r31, 3)), quote="FALSE") print(c("P del modelo robusto =", p_r31), quote="FALSE") summary(modelo.r31) coeftest(modelo.r31, df=modelo$df.residual) R2.r <- cor(modelo.r31$model[,1],modelo.r31$fitted.values)^2 print(c("R2 del modelo robusto =", round(R2.r,3)), quote="FALSE") Anova(modelo.r31, type=3, test="F") coeftest(modelo.r31, vcovHC(modelo.r31, type="HC4m")) ## otra opción: vcovHC(modelo, type="HC3") plot(cooks.distance(modelo)~modelo.r31$w) ## relación de los pesos del modelo.r3 con la distancia de Cook del modelo original identify(cooks.distance(modelo)~modelo.r31$w) ## acabad dando clik en el botón [Finish] de la ventana gráfica de Plots (esquina sup-decha) ## modelo.r32 <- rlm(eqt, data=datos, scale.est="Huber", maxit=100) F_r32 <- waldtest(modelo.r32)$F[2] p_r32 <- 1-pf(F_r32, df1=dim(modelo$model)[1]-modelo$df.residual-1, df2=modelo$df.residual) print(c("F del modelo robusto =", round(F_r32, 3)), quote="FALSE") print(c("P del modelo robusto =", p_r32), quote="FALSE") summary(modelo.r32) coeftest(modelo.r32, df=modelo$df.residual) R2.r <- cor(modelo.r32$model[,1],modelo.r32$fitted.values)^2 print(c("R2 del modelo robusto =", round(R2.r,3)), quote="FALSE") Anova(modelo.r32, type=3, test="F") coeftest(modelo.r32, vcovHC(modelo.r32, type="HC4m")) ## otra opción: vcovHC(modelo, type="HC3") plot(cooks.distance(modelo)~modelo.r32$w) ## relación de los pesos del modelo.r3 con la distancia de Cook del modelo original identify(cooks.distance(modelo)~modelo.r32$w) ## acabad dando clik en el botón [Finish] de la ventana gráfica de Plots (esquina sup-decha) ## modelo.r33 <- rlm(eqt, data=datos, psi=psi.bisquare, maxit=100) F_r33 <- waldtest(modelo.r33)$F[2] p_r33 <- 1-pf(F_r33, df1=dim(modelo$model)[1]-modelo$df.residual-1, df2=modelo$df.residual) print(c("F del modelo robusto =", round(F_r33, 3)), quote="FALSE") print(c("P del modelo robusto =", p_r33), quote="FALSE") summary(modelo.r33) coeftest(modelo.r33, df=modelo$df.residual) R2.r <- cor(modelo.r33$model[,1],modelo.r33$fitted.values)^2 print(c("R2 del modelo robusto =", round(R2.r,3)), quote="FALSE") Anova(modelo.r33, type=3, test="F") coeftest(modelo.r33, vcovHC(modelo.r33, type="HC4m")) ## otra opción: vcovHC(modelo, type="HC3") plot(cooks.distance(modelo)~modelo.r33$w) ## relación de los pesos del modelo.r3 con la distancia de Cook del modelo original identify(cooks.distance(modelo)~modelo.r33$w) ## acabad dando clik en el botón [Finish] de la ventana gráfica de Plots (esquina sup-decha) ## aproximación robusta que estima nuevos coeficientes y errores standard (package:robustbase) ## "specifying how points (potential outliers) in x-space are downweighted" ## lmrob uses the latest of the fast-S algorithms and heteroscedasticity and autocorrelation corrected (HAC) standard errors ## Los "methods" pueden ser MM (M-regression estimate) que genera weights más pequeños para más observaciones ## o SMDM (intial S-estimate, followed by an M-estimate, a Design Adaptive Scale estimate and a final M-step) ## que principalmente penaliza a las observaciones "más outliers". ## previamente hemos creado un modelo lm denominado "modelo" (usando eqt y datos) que ahora queremos valorar robustamente modelo.r <- lmrob(eqt, data=datos, method="SMDM", setting = "KS2014") waldtest(modelo.r, test="F") summary(modelo.r) coeftest(modelo.r, df=modelo$df.residual) Anova(modelo.r, type=3, test="F") R2.r <- cor(modelo.r$model[,1],modelo.r$fitted.values)^2 print(c("R2 del modelo robusto =", round(R2.r,3)), quote="FALSE") plot(cooks.distance(modelo)~modelo.r$rweights) ## relación de los pesos del modelo.r con la distancia de Cook del modelo original identify(cooks.distance(modelo)~modelo.r$rweights) hist(modelo.r$rweights, main="PESOS-SMDM DE LA OBSERVACIONES") ## histograma de los pesos del modelo lmRob plot(modelo.r$rweights, dffits(modelo)) ## relación de los pesos del modelo.r y los valores de dffits del modelo original plot(modelo.r$rweights, hatvalues(modelo)) ## relación de los pesos del modelo.r y los valores de leverage del modelo original ## gráfico de relaciones entre weights panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...) { usr <- par("usr") on.exit(par(usr)) par(usr = c(0, 1, 0, 1)) r <- abs(cor(x, y, use="complete.obs")) txt <- format(c(r, 0.123456789), digits=digits)[1] txt <- paste(prefix, txt, sep="") if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt) text(0.5, 0.5, txt, col="blue", cex = cex.cor * (1 + r) / 1) ## cambiar el tamaño de la letra modificando "1" en el denominador de "cex = cex.cor * (1 + r) / 1" } ## ampliarla con Zoom dentro de la solapa Plots pairs(residuals(modelo)~modelo.r31$w+modelo.r32$w+modelo.r33$w+modelo.r$rweights, data=datos, cex.labels=2, pch="o", lower.panel = panel.cor) pairs(dffits(modelo)~modelo.r31$w+modelo.r32$w+modelo.r33$w+modelo.r$rweights, data=datos, cex.labels=2, pch="o", lower.panel = panel.cor) pairs(cooks.distance(modelo)~modelo.r31$w+modelo.r32$w+modelo.r33$w+modelo.r$rweights, data=datos, cex.labels=2, pch="o", lower.panel = panel.cor) pairs(hatvalues(modelo)~modelo.r31$w+modelo.r32$w+modelo.r33$w+modelo.r$rweights, data=datos, cex.labels=2, pch="o", lower.panel = panel.cor) ## **** 12 **** ## CONSTRUCCIÓN DE HIPÓTESIS NULAS PARA LAS F DE LOS EFECTOS EN EL MODELO ## valores nuevos de significación recalculados a partir de F obtenidas de modelos aleatorios ## desacoplando los valores reales de la respuesta de los de las predictoras ## ## generamos 4.000 modelos aleatorios desacoplando la respuesta de los predictores n.efectos <- length(Anova(modelo.nulo, type=3, test="F")$Df)-1 nef <- n.efectos F.nulas <- matrix(nrow=4000, ncol=nef-1) ## proceso muy lento; efectúa 4000 modelos lm con datos aleatorizados for (i in 1:4000) { analizar.nulos <- modelo$model analizar.nulos[,1] <- sample(modelo$model[,1], replace=FALSE) ## sin reemplazo modelo.nulo <- lm(eqt, data=analizar.nulos) F.esta.vez <- Anova(modelo.nulo, type=3, test="F")[c(2:nef),3] F.nulas[i,] <- F.esta.vez } F.nulas <- as.data.frame(F.nulas) colnames(F.nulas) <- rownames(Anova(modelo.nulo, type=3, test="F"))[c(2:nef)] ## ## calculamos los cuantiles de las F nulas tabla.quantiles <- matrix(nrow=(nef-1), ncol=5) colnames(tabla.quantiles) <- c("Fp_0.1", "Fp_0.05", "Fp_0.01", "Fp_0.005", "Fp_0.001") for (j in 1:(nef-1)) { tabla.quantiles[j,] <- quantile(F.nulas[,j], probs=c(0.9, 0.95, 0.99, 0.995, 0.999), na.rm=TRUE) } tabla.quantiles <- as.data.frame(tabla.quantiles) tabla.quantiles$F_observada <- Anova(modelo, type=3, test="F")[c(2:nef),3] rownames(tabla.quantiles) <- rownames(Anova(modelo, type=3, test="F"))[c(2:nef)] ## tabla con las F de los modelos nulos a diferentes niveles de P (0.1, 0.05, 0.01, 0.005, 0.001) y las F reales de los efectos ## i.e., tabla de significación para las F con nuestros propios datos round(tabla.quantiles, 3) ## ## obtención de valores críticos de F para p determinadas (0.95 = 1-p cuando p=0.05) qf(0.90, df1=1, df2=modelo$df.residual) qf(0.95, df1=1, df2=modelo$df.residual) qf(0.99, df1=1, df2=modelo$df.residual) qf(0.995, df1=1, df2=modelo$df.residual) qf(0.999, df1=1, df2=modelo$df.residual) ## ## calculamos las nuevas significaciones tabla.F <- matrix(nrow=(nef-1), ncol=3) colnames(tabla.F) <- c("F_modelo", "p_modelo", "p_nueva") tabla.F <- as.data.frame(tabla.F) tabla.F[,1] <- Anova(modelo, type=3, test="F")[c(2:nef),3] tabla.F[,2] <- Anova(modelo, type=3, test="F")[c(2:nef),4] for (j in 1:(nef-1)) { tabla.F[j,3] <- 1-ecdf(F.nulas[,j])(tabla.F[j,1]) } rownames(tabla.F) <- rownames(Anova(modelo, type=3, test="F"))[c(2:nef)] ## Tabla con las F y P originales del modelo, y las nuevas P (p_nueva) a partir de las F de modelos aleatorios round(tabla.F[-1,], 5) ## **** 13 **** ## REDUCIR MANUALMENTE EL MODELO ## podemos reducir el modelo original del siguiente modo (quitando efectos con el signo menos "-" delante del término a eliminar) ### ¡¡ ojo !! poned aquí VUESTROS DATOS modelo.reducido <- update(modelo, .~. -ZONA:SEXO:EDAD) modelo.efectos_principales <- update(modelo, .~. -ZONA:SEXO -ZONA:EDAD -SEXO:EDAD -ZONA:SEXO:EDAD) modelo.nulo <- lm(MUSCULO ~ 1, data=datos) ## tablas de Anova Anova(modelo, test="F", type=3) Anova(modelo.reducido, test="F", type=3) Anova(modelo.efectos_principales, test="F", type=3) ## comparación de modelos haciendo uso de AICc (AIC corregido para pequeños tamaños muestrales) AICc(modelo.nulo, modelo, modelo.reducido, modelo.efectos_principales) ## Otra posibilidad es ordenar los posibles modelos atendiendo a los valores de AICc, y teniendo en cuenta la combinación de diferentes términos ## esto podemos efectuarlo con el comando dredge del paquete MuMIn. ## utilizamos la versión de los coeficientes de regresión estandarizados (betas) con beta="sd" que es la más recomendada ## http://onlinelibrary.wiley.com/doi/10.1890/14-1639.1/abstract ## consultad también: http://warnercnr.colostate.edu/~kenb/pdfs/KenB/AICRelativeVariableImportanceWeights-Burnham.pdf ## no admite una fórmula previamente establecida (e.g., puesta en eqt) options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ### ¡¡ ojo !! poned aquí VUESTROS DATOS modelo <- lm(MUSCULO~TARSO+SEXO*EDAD*ZONA, data=datos) options(na.action = "na.fail") dredge.modelo <- dredge(modelo, beta="sd", rank="AICc", extra="R^2") ## resumen de modelos ordenados por AICc dredge.modelo coefTable(dredge.modelo) ## promedio de los modelos ## Std. Error son los "unconditional standard errors"; Adjusted SE es la versión revisada de los inflated standard errors ## 'conditional average' only averages over the models where the parameter appears ## 'full average' assumes that a variable is included in every model, but in some models the corresponding coefficient is set to zero summary(model.avg(dredge.modelo)) ## para todos los modelos summary(model.avg(dredge.modelo[c(1:2)])) ## para solo los dos primeros modelos de dredge.modelo summary(model.avg(dredge.modelo, subset=delta<4)) ## para los modelos con una delta <4 summary(model.avg(dredge.modelo, cumsum(weight)<=0.99)) ## seleccionando para promediar a aquellos modelos cuyo sumatorio de weight sea una cierta cantidad (de o a 1) confint(model.avg(dredge.modelo), level=0.95) ## para todos los modelos al 95% de probabilidad confint(model.avg(dredge.modelo, subset=delta<4), level=0.95) ## para los modelos con incremento_AICc<4 al 95% de probabilidad ## pesos de las variables en los modelos (suma de los valores weight de los modelos que incluyen cada variable) sw(model.avg(dredge.modelo)) ## con TODOS los modelos sw(model.avg(dredge.modelo, subset=delta<4)) ## seleccionando los modelos con delta < 4 ## versión más elaborada de tabla (e.g., con la estima 'full average') dredge.tabla <- data.frame(summary(model.avg(dredge.modelo))$coefmat.full, confint(model.avg(dredge.modelo), level=0.95)) colnames(dredge.tabla)[5] <- "signif_P" colnames(dredge.tabla)[2] <- "Uncondit.SE" colnames(dredge.tabla)[6] <- "ICinf95%" colnames(dredge.tabla)[7] <- "ICsup95%" print(dredge.tabla, digits=3) ## ## podemos efectuar una selección de modelos para promediar a aquellos cuyo sumatorio de weight sea una cierta cantidad (de o a 1) ## si sólo se selecciona un modelo en tonces produce Error: 'object' consists of only one model ## **** 14 **** ## ESPECIFICANDO TABLAS DE CONTRASTES PARTICULARES (mirar en Help "contr.treatment") ## consultad http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm ## los resultados de los contrastes planificados no afectan al "efecto global" que tiene el factor en la variable respuesta ## creamos una tabla de contraste de desvío para un factor con cuatro (4) niveles ### ¡¡ ojo !! poned aquí VUESTROS DATOS contr.sum(4) ## para ver cómo es la tabla de contrastes datos$FDOMINANCIA <- as.ordered(datos$DOMINANCIA) ## convertimos una variable numérica en un factor ordinal "ordered" matdesv <- matrix(c(1,0,0,-1, 0,1,0,-1, 0,0,1,-1), ncol=3) ## n-1 columnas (ncol=3) para n niveles del factor summary(lm(MUSCULO~FDOMINANCIA, data=datos, contrasts=list(FDOMINANCIA=matdesv))) ## que es lo mismo que lo siguiente summary(lm(MUSCULO~FDOMINANCIA, data=datos, contrasts=list(FDOMINANCIA=contr.sum))) ## que es lo mismo que la salida de STATISTICA y SPSS ## contrastes de Helmert contr.helmert(4) ## primer nivel vs siguiente, dos primeros niveles vs el siguiente, ... mathelmert <- matrix(c(-1,1,0,0, -1,-1,2,0, -1,-1,-1,3), ncol=3) ## n-1 columnas (ncol=3) para n niveles del factor summary(lm(MUSCULO~FDOMINANCIA, data=datos, contrasts=list(FDOMINANCIA=mathelmert))) summary(lm(MUSCULO~FDOMINANCIA, data=datos, contrasts=list(FDOMINANCIA=contr.helmert))) ## creamos una tabla de contrastes lineales para un factor con cuatro niveles (valores 1, 2, 3 y 4) ## It returns contrasts based on orthogonal polynomials matlin <- matrix(c(-3,-1,1,3, 1,-1,-1,1, -1,3,-3,1), ncol=3) ## manteniendo proporcionalidad y sumando CERO los coeficientes de contraste. summary(lm(MUSCULO~FDOMINANCIA, data=datos, contrasts=list(FDOMINANCIA=matlin))) ## da como una regresión polinomial (en STATISTICA) ## elegimos el tipo de contraste que queremos, de los tres que salen contr.poly(4) ## para ver cómo es la tabla de contrastes lineal (primera columna) summary(lm(MUSCULO~FDOMINANCIA, data=datos, contrasts=list(FDOMINANCIA=contr.poly))) ## introduce los efectos lineal y cuadrático; da como STATISTICA ## **** 15 **** ## TEST DE PARALELISMO EN EL ANCOVA ## hacemos el modelo con la covariante (modelo) y otro modelo que introduce la interacción de los factores y la covariante (modelo.testparal) ### ¡¡ ojo !! poned aquí VUESTROS DATOS modelo <- lm(MUSCULO ~ TARSO + ZONA*SEXO*EDAD, data=datos) modelo.testparal <- lm(MUSCULO ~ TARSO*ZONA*SEXO*EDAD, data=datos) anova(modelo, modelo.testparal) ## test de comparación de dos modelos; para examinar SÓLO la interacción de los factores y la covariante ## si el anterior test no dió significativo, no deberíamos mirar los siguientes resultados ## de los que sólo deberíamos valorar las significaciones de las interacciones COVARIANTE x FACTOR(es). Anova(modelo.testparal, test="F", type=3) ## tabla de Anova (tipo III) para exámen de significaciones drop1(modelo.testparal, ~., test="F") ## equivalente al anterior; miramos el resultado SÓLO de la línea que incluye el efecto de la interacción de los factores y la covariante ## examen de las interacciones concretas, convergentes con los resultados de la tabla de Anova ### ¡¡ ojo !! poned aquí VUESTROS DATOS testInteractions(modelo.testparal, pairwise=c("ZONA", "SEXO", "EDAD"), slope="TARSO") testInteractions(modelo.testparal, pairwise=c("ZONA", "SEXO"), slope="TARSO") ## cuantificación de las pendientes en todas las celdas del diseño factorial deseado: ### ¡¡ ojo !! poned aquí VUESTROS DATOS interactionMeans(modelo.testparal, slope="TARSO") ## es lo mismo que la siguiente línea de código interactionMeans(modelo.testparal, pairwise=c("ZONA", "SEXO", "EDAD"), slope="TARSO") ## con los tres factores del modelo ## ## si existe ausencia de paralelismo debemos utilizar un MODELO DE DIFERENTES PENDIENTES modelo.difpend <- lm(MUSCULO~ZONA*SEXO*EDAD/TARSO, data=datos) summary(modelo.difpend) ## miramos las últimas líneas que valoran el efecto de las diferentes regresiones con la covariante drop1(modelo.difpend, ~., test="F") ## para valorar el efecto global (i.e., combinado) del efecto de la covariante; SS tipo-III Anova(modelo.difpend, type=3, test="F") ## para valorar el efecto global (i.e., combinado) del efecto de la covariante; SS tipo-III ## **** 16 **** ## TRANSFORMACIÓN BOX-COX ## añadir +1 a la variable respuesta SI tiene ceros, en cuyo caso creamos una ecuación de modelo nueva ## la transformación es: Y' = (Y^lambda - 1)/lambda ,ó, Y' = ((Y+1)^lambda - 1)/lambda ## la siguiente línea de código SÓLO TRANSFORMA la respuesta de nuestra ecuación eqt ## podemos afinar y restringir aun más el rango de lambda; por ejemplo: lambda=seq(0,1, 1/1000)) print(boxCox(modelo, lambda=seq(-2,2, 1/100))) ## with(..., x[which.max(y)]) para seleccionar lambda; en ... poner un objeto de car::boxCox o toda la instrucción boxCox(....) lambda.bc <- with(boxCox(modelo, lambda=seq(-2,2, 1/1000)), x[which.max(y)]) print(c("parámetro lambda de la transformación Box-Cox =",round(lambda.bc, 3)), quote=FALSE) ## transformamos la variable a continuación ... que llamamos respuesta.bc respuesta.bc <- ((modelo$model[,1]^lambda.bc)-1)/(lambda.bc) ## para explorar la relación entre la respuesta original y su transformación Box-Cox plot(respuesta.bc, modelo$model[,1], xlab="variable respuesta transformada", ylab="variable respuesta original") ## ahora rehacemos el modelo sustituyendo la variable respuesta del modelo por datos$respuesta.bc ## Aquí abajo, formulado automáticamente: datos.bc <- modelo$model ## capturamos en un nuevo data-frame los datos de trabajo de nuestro "modelo" datos.bc[,1] <- respuesta.bc ## asignamos a la primera variable (¡la respuesta!) el valor de su transformada Box-Cox modelo.bc <- lm(eqt, data=datos.bc) ## corremos el mismo modelo pero con la respuesta transformada ## sus resultados son: summary(modelo.bc) Anova(modelo.bc, type=3, test="F") coeftest(modelo.bc, df=modelo$df.residual, vcovHC(modelo.bc, type="HC4")) Anova(modelo.bc, type=3, test="F", white.adjust="hc4") ## a este nuevo modelo le podemos aplicar todo el resto de las cosas vistas hasta este momento ## **** 17 **** ## ANÁLISIS CON EL PAQUETE "afex" library(afex) ## nuestro análisis requiere tener una variable que identifique numéricamente las observaciones diferentes; ## en este caso es un valor por cada fila del data frame (juego de datos en "datos"); podemos crearla con ndatos <- dim(datos)[1] ## dim(datos) da las dimensiones de la matriz: casos-filas [1], variables-columnas [2] id <- seq(1, ndatos, 1) ## construimos un vector id, del valor 1, al ndatos, de uno-en-uno ## construimos el modelo; en error se especifica lo que define la unidad muestral (el vedtor id); en esta ocasión todas las muestras son independientes ## factorice=FALSE es para no convertir en factor a la covariante (MUSCULO); debemos comprobar que los factores ya sean "factores" (ZONA, SEXO y EDAD) ## es="pes" se refiere a partial eta squared, para particionar la varianza de los efectos ## utiliza por defecto la tabla de contrastes de desvío (Contrasts set to contr.sum for the following variables: ...) ### ¡¡ ojo !! poned aquí VUESTROS DATOS aovcar <- aov_car(MUSCULO~TARSO+ZONA*SEXO*EDAD+Error(id), factorize=FALSE, anova_table=list(correction="none", es="pes"), type=3, data=datos) ## resultados; de diferentes modos, equivalentes, con más o menos información ## aovcar$lm contiene toda la información de un objeto lm(...) summary(aovcar) ## es equivalente a aovcar$anova_table aovcar$Anova ## Anova con A mayúsula; llama al comando Anova{car} summary(aovcar$lm) ## resultados de la R2, adj R2, F df y p del modelo globalal final de la salida. ## importancia de los efectos en el modelo (effect sizes); calculamos eta2 y partial eta2 ## con SS de tipo III; no considerar la última fila de "Residuals" SSerror <- sum(residuals(aovcar$lm)^2) SStotal <- sum((aovcar$lm$model[,1]-mean(aovcar$lm$model[,1]))^2) SSmodelo <- SStotal-SSerror R2 <- SSmodelo/SStotal print(c("R2 (%) = ", round(R2*100,2)), quote=FALSE) ## R2 del modelo tabla.ss <- as.data.frame(aovcar$Anova) tabla.ss <- tabla.ss[-1,] eta2 <- tabla.ss[,1]/SStotal ## proporción de la varianza explicada por los efectos partialeta2 <- tabla.ss[,1]/(tabla.ss[,1]+SSerror) tabla.ss <- data.frame(eta2, partialeta2,tabla.ss) print(tabla.ss, digits=4) ## no tendremos en cuenta los valores de la última fila (Residuals) ## examen de la violación del supuesto de homocedasticidad coeftest(aovcar$lm, vcovHC(aovcar$lm, type="const")) ## comparadlo con la salida inicial sin efectuar correcciones coeftest(aovcar$lm, vcov=vcovHC) ## es igual que coeftest(modelo, vcovHC(modelo, type="HC3")) coeftest(aovcar$lm, vcovHC(aovcar$lm, type="HC4m")) ## HC4m testing inference is more reliable under both normal and nonnormal random errors ## exploración visual global de los residuos hist(residuals(aovcar$lm), density=5, freq=FALSE, main="residuos del modelo") curve(dnorm(x, mean=mean(residuals(aovcar$lm)), sd=sd(residuals(aovcar$lm))), col="red", lwd=2, add=TRUE, yaxt="n") qqnorm(residuals(aovcar$lm), main="residuos del modelo") qqline(residuals(aovcar$lm), col="red") plot(aovcar$lm$fitted.values, residuals(aovcar$lm), main="¿HAY HETEROCEDASTICIDAD?") abline(h=0, col="red") ## traza una línea horizontal (h) por el Y=0 ## lo mismo de antes pero de una sola vez y ... más par(mfcol=c(1,1)) ## fija un sólo panel gráfico par(mfcol=c(2,2)) ## fija cuatro paneles según 2 columnas y 2 filas plot(aovcar$lm) par(mfcol=c(1,1)) ## volvemos al modo gráfico de un solo panel ## shapiro.test(residuals(aovcar$lm)) ## test de shapiro de normalidad de residuos ## ## Podemos aplicar todo lo aplicado en las líneas previas utilizando otros paquetes. ## Sólo tendremos que pasar el objeto aovcar$lm a modelo con: modelo <- aovcar$lm ## **** 18 **** ## APROXIMACIONES BAYESIANAS ## ## Aproximación bayesiana usando la simulación para representar la incertidumbre en los coeficientes de regresión. ## The sim function gets posterior simulations of sigma and beta from a lm object. ## sim() uses maximum a posteriori probability estimation (MAP) in all cases. ## Consultad (pag. 236): https://pdfs.semanticscholar.org/475d/df234f49d15fa13443a3c606c07d8747624b.pdf ## It uses the values for the initial parameters, the given data, the model and the prior distribution ## to determine whether the initial values for the parameters are good values or not. ## The MCMC algorithm (Markov Chain Monte Carlo) will (ideally) converge towards the joint posterior distribution of the parameters ## given the priors, the initial parameter values, the model and the data. ## Thereby fitting the model, and sampling from the posterior distribution at the same time. ## detach(package:psych) library(matrixStats) ## para cuantiles de matrices ## sim.modelo <- sim(modelo, n.sims=1000) ## library(psych) ## tabla.simulada <- describe(coef(sim.modelo))[,c(2:4,11:12)] colnames(tabla.simulada)[3] <- "Std.Error" quantiles95 <- colQuantiles(coef(sim.modelo), probs=c(0.025, 0.975)) round(quantiles95, 3) tabla.simulada <- data.frame(quantiles95, tabla.simulada, summary(modelo)$coefficients[,c(1:2)]) colnames(tabla.simulada)[1] <- "percentil_2.5" colnames(tabla.simulada)[2] <- "percentil_97.5" colnames(tabla.simulada)[8] <- "Coef. original" colnames(tabla.simulada)[9] <- "se original" print(tabla.simulada, digits=5) ## Los valores mean se refieren a la media de los coeficientes de regresión ## Tanto más parecidos son mean y Coef. original, y Std.Error y se original ... menos efecto de puntos influentes y perdidos ## Los sesgos y las kurtosis de los coeficientes simulados deberían tomar valores próximos a "cero". ## De ese modo podríamos efectuar estimas paramétricas de intervalos de confianza mediante: ## mean +/- Std.Error*1.96 qt(1-0.05/2, df=1000) ## Aproximación utilizando MCMCregress::MCMCpack ## http://www.stat.umn.edu/geyer/3701/notes/mcmc-bayes.html ## https://www.jstatsoft.org/article/view/v042i09/v42i09.pdf ## https://www.baseballprospectus.com/news/article/38289/bayesian-bagging-generate-uncertainty-intervals-catcher-framing-story/ ## https://en.wikipedia.org/wiki/Credible_interval ## http://mc-stan.org/rstanarm/reference/posterior_interval.stanreg.html library(MCMCpack) ## eqt b0.OLSpriors <- as.numeric(coefficients(lm(eqt, data=datos))) ## mcmc.modelo <- MCMCregress(eqt, data=datos, burnin=1000, mcmc=40000, marginal.likelihood="none") ## la siguiente opción es más elaborada, definiendo los priors a partir del modelo lm(eqt, data=datos) y estableciendo la precisión de ellos mcmc.modelo2 <- MCMCregress(eqt, data=datos, burnin=1000, mcmc=40000, marginal.likelihood="Chib95", b0=b0.OLSpriors, B0=abs(b0.OLSpriors/1000)) ## summary(modelo)$coefficients summary(mcmc.modelo) summary(mcmc.modelo2) ## par(mfcol=c(1,1)) plot(mcmc.modelo) ## An estimate I (the `dependence factor') of the extent to which autocorrelation inflates ## the required sample size is also provided. Values of I larger than 5 indicate strong autocorrelation ## which may be due to a poor choice of starting value, high posterior correlations or `stickiness' of ## the MCMC algorithm. raftery.diag(mcmc.modelo, q=0.025, r=0.005, s=0.95, converge.eps=0.001) raftery.diag(mcmc.modelo2, q=0.025, r=0.005, s=0.95, converge.eps=0.001) ## **** 19 **** ## REFERENCIAS ## ## Experimental Design and Data Analysis for Biologists. Quinn & Keough 2002 ## http://www.lacbiosafety.org/wp-content/uploads/2011/09/experimental-design-and-data-analysis-for-biologists1.pdf ## R: Analysis of variance (ANOVA) (con diferencias entre R y SPSS en relación con Type-I, Type-II y Type-III) ## https://egret.psychol.cam.ac.uk/statistics/R/anova.html ## Block designs ## https://onlinecourses.science.psu.edu/stat503/node/18 ## Formulae in R: ANOVA and other models ## http://conjugateprior.org/2013/01/formulae-in-r-anova/ ## Obtaining the same ANOVA results in R as in SPSS - the difficulties with Type II and Type III sums of squares ## http://myowelt.blogspot.com.es/2008/05/obtaining-same-anova-results-in-r-as-in.html ## Analysis of Variance and Covariance in R (introducción y manejo de R; múltiples tipos de modelos) ## http://www.southampton.ac.uk/~cpd/anovas/datasets/ANOVA%20in%20R.htm ## Contrast Coding Systems for categorical variables ## http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm ## vcovHC Heteroskedasticity-Consistent Covariance Matrix Estimation (páginas 17-18) ## https://cran.r-project.org/web/packages/sandwich/sandwich.pdf ## Econometric Computing with HC and HAC Covariance Matrix Estimators (página 4 para los tipos de HC0-HC4) ## https://cran.r-project.org/web/packages/sandwich/vignettes/sandwich.pdf ## Randomized Block Designs and Latin Square Design ## http://web.udl.es/Biomath/Bioestadistica/Dossiers/Temas%20especiales/ANOVA/The%20randomized%20complete%20block%20design%20(R).pdf ## FACTORIAL BETWEEN SUBJECTS ANOVA using aov (interesante por representaciones y test a posteriori TukeyHSD) ## http://ww2.coastal.edu/kingw/statistics/R-tutorials/factorial.html ## Análisis de la varianza usando AOV ## http://www.personality-project.org/r/r.anova.html ## How can I test contrasts in R? ## http://www.ats.ucla.edu/stat/r/faq/testing_contrasts.htm ## A stagewise rejective multiple test procedure based on a modified Bonferroni test ## http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/MultipleComparision/Hommel88.pdf ## A simple sequentially rejective multiple test procedure ## http://www.ime.usp.br/~abe/lista/pdf4R8xPVzCnX.pdf ## AOV Assumption Checking and Transformations ## http://goo.gl/FCeT07 ## A sharper Bonferroni procedure for multiple tests of significance ## http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/MultipleComparision/Hochberg88.pdf ## Sum of Squares Type I, II, III: the underlying hypotheses, model comparisons, and their calculation in R ## http://www.uni-kiel.de/psychologie/dwoll/r/ssTypes.php ## Analysing interactions of fitted models ## https://cran.r-project.org/web/packages/phia/vignettes/phia.pdf ## ## ## ## OTRAS COSAS de interés. ## para estimar significaciones a partir de valores de F recalculados: pf(F, df_efectos, df_error, lower.tail=FALSE) pf(5.0825, 1, 150, lower.tail=FALSE) ## para estimar significaciones a partir de valores de t de Student recalculados: pt(t, df, lower.tail=FALSE) ## esto es para una cola; para la p con dos colas hay que multiplicar por "dos" el valor de la p. pt(1.965, df=200, lower.tail=FALSE)*2 ## ## para estimar quantiles. deciles: seq(0, 1, 0.1); quintiles: seq(0, 1, 02); cuartiles: seq(0, 1, 0.25); terciles: seq(0, 1, 1/3) quantile(datos$MUSCULO, prob=seq(0, 1, 0.1)) ## ## Colores para gráficos en R ## http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf ##