################################### ## DISEÑOS FACTORIALES GENERALES ## ## ANOVA DE MEDIDAS REPETIDAS ## ## n-way within-subjects ## ## package car ## ## Luis M. Carrascal ## ################################### ## ## para importar y exportar datos ## datos <- read.table("clipboard", header=T, sep="\t", dec=".") ## importar desde portapapeles ## write.table(datos, "nombre.txt", sep="\t") ## lo guarda en el directorio por defecto de RStudio ## ## de internet y ponderlo en uso en "environment" meconecto <- url("http://www.lmcarrascal.eu/cursos/nwayREPMEAS_ANOVA.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/nwayREPMEAS_ANOVA.RData', destfile ="para_practicar.RData") ## mirad: https://www.datacamp.com/community/tutorials/r-data-import-tutorial ## ## usar el juego de datos "nwayREPMEAS_ANOVA.RData" ## names(datos) ## ## ## CARGAMOS LOS PAQUETES: library(car) ## para obtener VIF's y usar Anova, bootCase library(biotools) ## para efectuar el test de la M de Box library(MVN) ## para el test de la normalidad multivariante ## ## ## ## ## ESPECIFICAMOS LAS TABLAS DE CONTRASTES PARA LOS FACTORES ## Cargamos la siguiente línea de código para obtener contrastes definidos ## Se obtienen 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")) ## ## ORGANIZACIÓN DE DATOS DE EFECTOS DENTRO DE SUJETOS (WITHIN-SUBJECTS) ## ## combinamos con cbind los efectos within-subjects, ## sean de factores o covariantes cambiantes ## ## para la RESPUESTA, FACTOR DENTRO-DE two-way (dos factores que se ordenan según una secuencia) ## ponemos ordenadamente la secuencia de mis variables respuesta ## definir cómo organizar esas columnas: primero la FUENTE (niveles: inf y luz) y dentro de eso la POSICIÓN (niveles: inv y nor) ## aquí están los valores: respuesta2w.dentro <- cbind(datos$tasa_inf_inv, datos$tasa_inf_nor, datos$tasa_luz_inv, datos$tasa_luz_nor) ## otro modelo asumiendo que solo existe un factor repetido dentro-de con cuatro niveles (lo mismo que lo anterior pero llamado con distinto nombre) respuesta1w.dentro <- cbind(datos$tasa_inf_inv, datos$tasa_inf_nor, datos$tasa_luz_inv, datos$tasa_luz_nor) ## ## Construimos un nuevo dataframe con los niveles ordenados de los factores ## según hemos re-estructurado los valores de la respuesta en respuesta2w.dentro ## definimos los niveles del primer factor con los nombres que queramos FUENTE <- factor(c("Infra","Infra","Luz","Luz")) ## definimos los niveles del otro factor con los nombres sencillos que queramos POSICION <- factor(c("Invertido","Normal","Invertido","Normal")) ## ## combinamos los dos factores en un dataframe ## aquí están los nombres inteligibles factor2w.dentro <- data.frame(FUENTE, POSICION) ## ## para el caso de un solo factor repetido con cuatro niveles hacemos: EXPERIMENTO <- factor(c("Infra-Invertido", "Infra-Normal", "Luz-Invertido", "Luz-Normal")) factor1w.dentro <- data.frame(EXPERIMENTO) ## ## COVARIANTE CAMBIANTE DENTRO DE SUJETOS covariante2w.dentro <- cbind(datos$temp_inf_inv, datos$temp_inf_nor, datos$temp_luz_inv, datos$temp_luz_nor) ## ## para el caso de un solo factor repetido con cuatro niveles hacemos lo mismo (cambio solo el nombre) covariante1w.dentro <- cbind(datos$temp_inf_inv, datos$temp_inf_nor, datos$temp_luz_inv, datos$temp_luz_nor) ## ## ## ## CONSTRUIMOS EL MODELO AN(c)OVA ## que tiene un submodelo (modelo.entre) que se incluye en Anova(.) ## tiene la forma: respuesta-repetida ~ ## factor_fijo + covariante_fija + covariante_cambiante ## idata define cuál es la estructura del factor-dentro-de-sujetos ## idesign define qué diseño within-subjects queremos hacer (aquí two-way con interacción) ## type=3 establece un diseño de suma de cuadrados (SS) de tipo 3 (efectos parciales) ## ## two-way within-subjects factors & one between-subjects factor, with one fixed covariate plus one changing covariate submodelo <- lm(respuesta2w.dentro ~ datos$factorfijo + datos$covarfija + covariante2w.dentro) modelo <- Anova(submodelo, idata=factor2w.dentro, idesign=~FUENTE*POSICION, type="III") modelo2 <- Anova(submodelo, idata=factor2w.dentro, idesign=~FUENTE*POSICION, type="II") ## ## ## ## EXPLORACIÓN DE LOS SUPUESTOS CANÓNICOS DEL MODELO ## ## Normalidad de las variables respuesta originales ## mirad: https://cran.r-project.org/web/packages/MVN/vignettes/MVN.pdf ## descriptores de nuestras variables originales y test de Shapiro-Wilk de desvío de la normalidad uniNorm(as.matrix(respuesta2w.dentro), type = "SW", desc = TRUE) ## plots: de histogramas y Q-Q plots par(mfcol=c(1,1)) uniPlot(as.matrix(respuesta2w.dentro), type="histogram") uniPlot(as.matrix(respuesta2w.dentro), type="qqplot") par(mfcol=c(1,1)) ## ## Normalidad Multivariante ## test de Royston de normalidad multivariante; para cuando el tamaño muestral (N) sea < 50 ## usaremos el test de Henze-Zirkler cuando N > 100: hzTest(...) roystonTest(as.matrix(respuesta2w.dentro), qqplot=TRUE) hzTest(as.matrix(respuesta2w.dentro), qqplot=TRUE) ## ## efectuamos el test de la M de Box, SÓLO con cada factor "entre sujetos" del análisis ## "sacamos" las variables respuesta de nuestro objeto creado con cbind(...) ## Si no podemos efectuar la interacción completa (por falta de muestra), ## probad con efectos principales e interacciones dobles ## ¡¡poned aquí vuestros datos!! reemplazando los del ejemplo boxM(respuesta2w.dentro, datos$factorfijo) ## ## TEST DE HOMOGENEIDAD DE VARIANZAS a través de los niveles de los factores "entre sujetos" ## no admite covariantes ## otra opción es center=mean (la del programa STATISTICA) ## lo hacemos con cada una de las variables respuesta ### ¡¡ ojo !! poned aquí VUESTROS DATOS nvarresp <- dim(factor2w.dentro)[1] for (i in 1:nvarresp) { print(c("Respuesta", i), quote=FALSE) print(leveneTest(respuesta2w.dentro[,i]~datos$factorfijo, data=datos, center=median)) print("******************************************************", quote=FALSE) } ## ## ## ## RESIDUOS ## los visualizamos ## es para el conjunto de las N filas-observaciones x r variables respuesta hist(residuals(submodelo)) qqnorm(residuals(submodelo)) qqline(residuals(submodelo), col="red") ## ## test de Shapiro de normalidad de residuos shapiro.test(residuals(submodelo)) ## ## Homocedasticidad de los residuos a través de las predicciones del modelo plot(fitted(submodelo), residuals(submodelo), main="relación entre residuos y predicciones") abline(h=0, col="red") ## ## ## ## ## RESULTADOS ## tabla sintética de efectos con la aproximación within-subjects ## al final se presentan los resultados de los tests de: "Mauchly Tests for Sphericity" y "Greenhouse-Geisser and Huynh-Feldt Corrections" ## en el caso de que los factores "dentro-de" hubiesen tenido más de tres niveles summary(modelo, multivariate=FALSE) summary(modelo2, multivariate=FALSE) ## lo mismo pero añadiendo las partial eta2 para cada efecto tabla <- data.frame(summary(modelo, multivariate=FALSE)$univariate.tests[,]) partialeta2 <- tabla$SS/(tabla$SS+tabla$Error.SS) tabla.anova <- data.frame(tabla, partialeta2) round(tabla.anova, 5) ## ## resultados pormenorizados para cada variable respuesta (niveles de los factores "dentro de") summary(submodelo) ## salida multivariante al estilo MANOVA por si hubiésemos violado los supuestos de ## esfericidad y simetría compuesta; al final proporciona el modelo within-subjeects summary(modelo, multivariate=TRUE) Anova(submodelo, type=3, test="Wilks") ## otras posibilidades son "Pillai", "Hotelling-Lawley", "Roy" ## ## ## ## ## para el otro diseño con solo un factor repetido dentro-de con cuatro niveles: ## one-way within-subjects factors & one between-subjects factor, with one fixed covariate plus one changing covariate submodelo <- lm(respuesta1w.dentro ~ datos$factorfijo + datos$covarfija + covariante1w.dentro) modelo <- Anova(submodelo, idata=factor1w.dentro, idesign=~EXPERIMENTO, type="III") ## ## resultados summary(modelo, multivariate=FALSE) ## lo mismo pero añadiendo las partial eta2 para cada efecto tabla <- data.frame(summary(modelo, multivariate=FALSE)$univariate.tests[,]) partialeta2 <- tabla$SS/(tabla$SS+tabla$Error.SS) tabla.anova <- data.frame(tabla, partialeta2) round(tabla.anova, 5) ## ## resultados pormenorizados para cada variable respuesta (niveles de los factores "dentro de") summary(submodelo) ## salida multivariante al estilo MANOVA por si hubiésemos violado los supuestos de ## esfericidad y simetría compuesta; al final proporciona el modelo within-subjeects summary(modelo, multivariate=TRUE) Anova(submodelo, type=3, test="Wilks") ## otras posibilidades son "Pillai", "Hotelling-Lawley", "Roy" ##