##################################### ## Modelos de efectos aleatorios ## ## MODELOS MIXTOS GENERALES ## ##################################### ## Luis M. Carrascal ## https://lmcarrascal.eu ## revisión 27/02/2022 ## REFERENCIAS EN LINEA ## https://en.wikipedia.org/wiki/Multilevel_model ## https://www.sciencedirect.com/science/article/pii/S0169534709000196 ## https://peerj.com/articles/4794/ ## https://peerj.com/articles/9522/ ## https://docplayer.es/20740859-Modelos-mixtos-lineales-una-introduccion-para-el-usuario-temeroso.html ## http://www.john-ros.com/Rcourse/lme.html?utm_source=pocket_mylist ## http://www.bristol.ac.uk/cmm/learning/videos/random-slopes.html?utm_source=pocket_mylist ## https://yury-zablotski.netlify.app/post/mixed-effects-models-1/?utm_source=pocket_mylist ## https://yury-zablotski.netlify.app//post/mixed-models/?utm_source=pocket_mylist ## http://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf (páginas 6 y 30) ## https://rpsychologist.com/r-guide-longitudinal-lme-lmer ## https://fromthebottomoftheheap.net/2021/02/02/random-effects-in-gams/ ## https://www.crumplab.com/psyc7709_2019/book/docs/a-tutorial-for-using-the-lme-function-from-the-nlme-package-.html ## https://stat.ethz.ch/pipermail/r-help/2006-October/115572.html ## https://peerj.com/articles/12794/ ## https://cran.r-project.org/web/packages/glmmTMB/vignettes/glmmTMB.pdf ## https://backend.orbit.dtu.dk/ws/files/154739064/Publishers_version.pdf (mirad los ejemplos del Appendix A) ## https://cran.r-project.org/web/packages/glmmTMB/vignettes/model_evaluation.pdf ## https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html ## https://besjournals.onlinelibrary.wiley.com/doi/pdfdirect/10.1111/2041-210X.13434 ## https://m-clark.github.io/posts/2019-05-14-shrinkage-in-mixed-models/ ## https://stats.stackexchange.com/questions/37647/what-is-the-minimum-recommended-number-of-groups-for-a-random-effects-factor?utm_source=pocket_mylist ## IMPORTACIÓN DE DATOS ## para importar los datos fuente ## 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=".") ## para exportar ## write.table(datos, "nombre.txt", sep="\t", row.names=FALSE) ## lo guarda en el directorio por defecto de RStudio ## datos de ejemplo obtenidos de una URL meconecto <- url("http://www.lmcarrascal.eu/cursos/nwayREPEATED.RData") ## abro una conexión con internet load(meconecto) ## cargo el archivo de la conexión close(meconecto); rm(meconecto) ## cierro la conexión ## convertimos el identificador numérico de los individuos (id; factor aleatorio) en un factor datos$id <- as.factor(datos$id) ## CARGAMOS PAQUETES DE ANÁLISIS library(lme4) ## modelos generales lmer library(nlme) ## modelos generales lineales lme library(DHARMa) ## para diagnosis de modelos library(moments) ## para obtener el sesgo, kurtosis y sus significaciones: kurtosis, anscombe.test, skewness, agostino.test library(lmerTest) ## para MS, df, p ... usando type 3/type 1 hypotheses with "Satterthwaite" and "Kenward-Roger" library(pbkrtest) ## necesario para lmerTest library(LMERConvenienceFunctions) library(car) ## para Anova(modelo, type=3) para equivalente a suma de tipo III; para boxCox(modelo, lambda=seq(-2,2, 1/100)) library(MuMIn) ## para AICc library(lmtest) ## para lrtest library(lattice) ## para plots de residuos library(psych) ## para estadísticas resumen de tablas library(performance) ## para pruebas de los modelos library(sjPlot) ## para efectos parciales ## Definimos los TIPOS DE CONTRASTES del modelo ## https://rstudio-pubs-static.s3.amazonaws.com/84177_4604ecc1bae246c9926865db53b6cc29.html ## https://rpubs.com/monajhzhu/608609 ## https://www.clayford.net/statistics/tag/sum-contrasts/ ## options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## ## otra posibilidad es incluir lo siguiente en la función: ## contrasts=list(factor.mio1=contr.sum, factor.mio2=contr.sum, factor.mio3=contr.sum, ...) ## Por defecto R opera con contr.treatment que no es correcto: ## "[It] contrasts each level with the baseline level (specified by base); the baseline level is omitted. ## Note that this does not produce 'contrasts' as defined in the standard theory for linear models ## as they are not orthogonal to the intercept" ## contr.treatment: obtiene para los niveles de los factores los EFECTOS MARGINALES respecto al intercepto, que es el valor medio en el nivel de referencia ## los coeficientes de los niveles de los factores miden los cambios respecto al valor del nivel de referencia ## contr.sum: obtiene para los niveles de los factores los EFECTOS PRINCIPALES respecto a la gran media (promedio de todos los niveles de los factores) ## los coeficientes de los niveles de los factores miden los cambios respecto a la gran media. ## contr.sum asegura que todos los contrastes sumen cero, de manera que el "intercepto" sea la gran media. ## Los efectos se resumen con coeficientes que representan el número de niveles de factor (k) menos 1, eliminando el último nivel. ## Los coeficientes en contr.sum representan lna diferencia media respecto a la "gran media" para los primeros niveles del factor (del 1 al k???1). ## ved la diferencia para los contrastes en un factor con 5 niveles ## el baseline es el último nivel, y las columnas suman CERO (ortogonalidad) contr.sum(5) mat.desv1 <- matrix(c(1,0,0,0,-1, 0,1,0,0,-1, 0,0,1,0,-1, 0,0,0,1,-1), ncol=4) ## igual que el anterior mat.desv2 <- matrix(c(4,-1,-1,-1,-1, -1,4,-1,-1,-1, -1,-1,4,-1,-1, -1,-1,-1,4,-1), ncol=4) ## equivalente al anterior pero con diferentes coeficiente ## el baseline es el primer nivel, y las columnas NO suman CERO (ausencia de ortogonalidad) contr.treatment(5) ## para ver la secuencia ordenada de los niveles de un factor levels(datos$especie) ## por si hay problemas de estimas y convergencia ## introducid en lmer(..., control=control.lmer) lme(..., control=control.lme) control.lme <- lmeControl(maxIter=200, msMaxIter=200, msVerbose=TRUE, opt="optim") ## el optimizador opt="nlminb" hay veces que da problemas control.lmer <- lmerControl(check.conv.grad=.makeCC(action ="ignore", tol=1e-6, relTol=NULL), optimizer="bobyqa", optCtrl=list(maxfun=100000)) ## control.lmer <- lmerControl(check.conv.grad=.makeCC(action ="ignore", tol=1e-6, relTol=NULL), optCtrl=list(maxfun=100000)) ## control.lmer <- lmerControl(check.conv.grad=.makeCC(action ="ignore", tol=1e-6, relTol=NULL), optimizer="Nelder_Mead", optCtrl=list(maxfun=100000)) ## MODELOS ## desarrollados con: lme4::lmer (modelo.) ## El método ML produce estimas (negativamente) sesgadas de las componentes de la varianza. ## ML produce menores valores MSE (mean-squared errors) que REML ## El método REML es preferido para bajos tamaños muestrales donde el problema de ML se ve incrementado ## REML se usa para comparar modelos que comparten exactamente la misma parte fija ## REML estima mejores valores medios de los coeficientes, pero con mayores variabilidades en sus estimas ## en modelos lmer(...) frecuentemente aparece esta alarma: boundary (singular) fit ## en este caso se obtiene un ajuste "singular", indicativo de que el modelo está sobreajustado ## o que los efectos aleatorios son muy pequeños (<1.0e-3) ## (i.e., la estructura de efectos aleatorios es demasiado compleja para ser apoyada por los datos) ## ésto sugiere eliminar la parte más compleja de la estructura de los efectos aleatorios (generalmente "random slopes"), ## que conduce a un modelo más parsimonioso, no sobreajustado. ## ¡¡¡ PONED AQUI VUESTROS DATOS !!! (las variables respuesta y predictoras a utilizar) ## Intraclass correlation: the "unconditional means model" (or "null model"); Random intercept with fixed mean ## diferencias entre los niveles del factor aleatorio (individuo) lmer.0 <- lmer(a_por_b ~ 1 + (1|id), data=datos, control=control.lmer, REML=TRUE) ## lme.0 <- lme(a_por_b ~ 1, random=list(id=~1), data=datos, control=control.lme, method="REML") ## modelo de RANDOM INTERCEPT & FIXED SLOPES lmer.1 <- lmer(a_por_b ~ t_sonda+fuente*especie + (1|id), data=datos, control=control.lmer, REML=TRUE) ## lme.1 <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=~1), data=datos, control=control.lme, method="REML") ## modelo de FIXED INTERCEPT & RANDOM SLOPES, with UNCORRELATED random slopes for different covariates lmer.2 <- lmer(a_por_b ~ t_sonda+fuente*especie + (0+t_sonda|id)+(0+fuente|id), data=datos, control=control.lmer, REML=TRUE) ## lme.2 <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=~0+t_sonda, id=~0+fuente), data=datos, control=control.lme, method="REML") lme.2b <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=pdDiag(~0+t_sonda+fuente)), data=datos, control=control.lme, method="REML") ## modelo de UNCORRELATED RANDOM INTERCEPT & RANDOM SLOPES, with UNCORRELATED random slopes for different covariates; elimina la covarianza intercepto-pendientes lmer.3 <- lmer(a_por_b ~ t_sonda+fuente*especie + (1|id)+(0+t_sonda|id)+(0+fuente|id), data=datos, control=control.lmer, REML=TRUE) ## lme.3 <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=~1, id=~0+t_sonda, id=~0+fuente), data=datos, control=control.lme, method="REML") lme.3b <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=pdDiag(~1+t_sonda+fuente)), data=datos, control=control.lme, method="REML") ## modelo de CORRELATED RANDOM INTERCEPT & RANDOM SLOPES, with CORRELATED random slopes for different covariates lmer.4 <- lmer(a_por_b ~ t_sonda+fuente*especie + (1+t_sonda+fuente|id), data=datos, control=control.lmer, REML=TRUE) ## lme.4 <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=pdSymm(~1+t_sonda+fuente)), data=datos, control=control.lme, method="REML") ## modelo de CORRELATED RANDOM INTERCEPT & RANDOM SLOPES, with UNCORRELATED random slopes for different covariates lmer.5 <- lmer(a_por_b ~ t_sonda+fuente*especie + (1+t_sonda|id)+(1+fuente|id), data=datos, control=control.lmer, REML=TRUE) ## lme.5 <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=pdSymm(~1+t_sonda), id=pdSymm(~1+fuente)), data=datos, control=control.lme, method="REML") lme.5b <- lme(a_por_b ~ t_sonda+fuente*especie, random=list(id=~1+t_sonda, id=~1+fuente), data=datos, control=control.lme, method="REML") ## modelo de EFECTOS ALEATORIOS ANIDADOS, RAMDOM INTERCEPT (lme no puede gestionar, de manera directa, dos efectos aleatorios cruzados) lmer.anid.1 <- lmer(a_por_b ~ t_sonda+fuente + (1|especie/id), data=datos, control=control.lmer, REML=TRUE) ## los dos modelos lme de abajo son equivalentes (importa el orden de los efectos aleatorios: especie/id ---> especie, id; especie/id xxx> id, especie) lme.anid.1 <- lme(a_por_b ~ t_sonda+fuente, random=~1|especie/id, data=datos, control=control.lme, method="REML") lme.anid.1b <- lme(a_por_b ~ t_sonda+fuente, random=list(especie=~1, id=~1), data=datos, control=control.lme, method="REML") ## modelo de EFECTOS ALEATORIOS ANIDADOS, CORRELATED RAMDOM INTERCEPT & SLOPES lmer.anid.2 <- lmer(a_por_b ~ t_sonda+fuente + (1+t_sonda+fuente|especie/id), data=datos, control=control.lmer, REML=TRUE) ## en los modelos lme importa el orden de los efectos aleatorios: especie/id ---> especie, id; especie/id xxx> id, especie lme.anid.2 <- lme(a_por_b ~ t_sonda+fuente, random=list(especie=pdSymm(~1+t_sonda+fuente), id=pdSymm(~1+t_sonda+fuente)), data=datos, control=control.lme, method="REML") ## modelo de EFECTOS ALEATORIOS ANIDADOS, CORRELATED RANDOM SLOPES AND INTERCEPTS FOR TWO UNCORRELATED FIXED EFFECTS lmer.anid.3 <- lmer(a_por_b ~ t_sonda+fuente + (1+t_sonda|especie/id)+(1+fuente|especie/id), data=datos, control=control.lmer, REML=TRUE) ## en los modelos lme importa el orden de los efectos aleatorios: especie/id ---> especie, id; especie/id xxx> id, especie lme.anid.3 <- lme(a_por_b ~ t_sonda+fuente, random=list(especie=pdSymm(~1+t_sonda), especie=pdSymm(~1+fuente), id=pdSymm(~1+t_sonda), id=pdSymm(~1+fuente)), data=datos, control=control.lme, method="REML") ## modelo de EFECTOS ALEATORIOS ANIDADOS, UNCORRELATED RANDOM SLOPES AND INTERCEPTS FOR TWO UNCORRELATED FIXED EFFECTS lmer.anid.4 <- lmer(a_por_b ~ t_sonda+fuente + (1|especie/id)+(0+t_sonda|especie/id)+(0+fuente|especie/id), data=datos, control=control.lmer, REML=TRUE) ## en los modelos lme importa el orden de los efectos aleatorios: especie/id ---> especie, id; especie/id xxx> id, especie lme.anid.4 <- lme(a_por_b ~ t_sonda+fuente, random=list(especie=~1, especie=~0+t_sonda, especie=~0+fuente, id=~1, id=~0+t_sonda, id=~0+fuente), data=datos, control=control.lme, method="REML") ## comparación de modelos con IDENTICA PARTE FIJA y diferente aleatoria ## A la hora de comparar diferentes modelos anidados, el método REML lo podremos utilizar si, y sólo si, ## ¡¡¡se mantienen constantes los efectos fijos!!! (sean factores o covariantes) y SOLO cambia la estructura de los efectos aleatorios. ## NUNCA deberemos utilizar el método REML si en los modelos anidados que se comparan cambian los efectos fijos. ## Podemos actualizar los modelos usando ML o REML de la siguiente manera: ## update(modelo.lmer, REML=TRUE) o update(modlme, method="REML") ## AICc modelos con un factor aleatorio AICc(lmer.1, lmer.2, lmer.3, lmer.4, lmer.5) AICc(lme.1, lme.2, lme.3, lme.4, lme.5) ## AICc modelos con dos factores aleatorios anidados AICc(lmer.anid.1, lmer.anid.2, lmer.anid.3, lmer.anid.4) AICc(lme.anid.1, lme.anid.2, lme.anid.3, lme.anid.4) ## pesos de los modelos con un factor aleatorio round(Weights(AICc(lmer.1, lmer.2, lmer.3, lmer.4, lmer.5)), 5) exp(-0.5 * (AICc(lmer.1) - AICc(lmer.3))) ## ¿cuántas veces es mejor un modelo que otro? (primero el de menor valor de AICc) ## pesos de los modelos con dos factores aleatorios anidados round(Weights(AICc(lmer.anid.1, lmer.anid.2, lmer.anid.3, lmer.anid.4)), 5) exp(-0.5 * (AICc(lmer.anid.1) - AICc(lmer.anid.3))) ## ¿cuántas veces es mejor un modelo que otro? (primero el de menor valor de AICc) ## comparamos modelos anidados obtenidos mediante ML (no-REML) utilizando anova(..., ...) anova(update(lmer.1, REML=F), update(lmer.3, REML=F)) anova(update(lme.1, method="ML"), update(lme.3, method="ML")) anova(update(lmer.3, REML=F), update(lmer.5, REML=F)) anova(update(lme.3, method="ML"), update(lme.5, method="ML")) ## equivalente a utilizar likelihood ratio tests (usando estimas ML y no-REML) lrtest(update(lmer.1, REML=F), update(lmer.3, REML=F)) lrtest(update(lme.1, method="ML"), update(lme.3, method="ML")) ## elegimos un modelo modelo <- lmer.1 modlme <- lme.1 ## para sacar la fórmula de lo que hemos hecho: modelo@call; formula(modelo) modlme$call; formula(modlme) ## los datos con los que hemos trabajado están ahora en (respuesta y luego las columnas de las predictoras): modelo@frame modlme$data ## con modelos lme aparecen todos los datos del dataframe usado ## RESIDUOS DEL MODELO ## lecturas: https://arxiv.org/pdf/1502.06988.pdf ## https://besjournals.onlinelibrary.wiley.com/doi/pdfdirect/10.1111/2041-210X.13434 ## https://peerj.com/articles/9522.pdf ## https://link.springer.com/article/10.3758/s13428-021-01587-5 ## Exploración visual de la normalidad de los residuos del modelo ## modelo lmer hist(residuals(modelo), density=5, freq=FALSE, main="residuos del modelo", ylim=c(0,0.3)) curve(dnorm(x, mean=mean(residuals(modelo)), sd=sd(residuals(modelo))), col="red", lwd=2, add=TRUE, yaxt="n") ## modelo lme hist(residuals(modlme), density=5, freq=FALSE, main="residuos del modelo", ylim=c(0,0.3)) curve(dnorm(x, mean=mean(residuals(modlme)), sd=sd(residuals(modlme))), col="red", lwd=2, add=TRUE, yaxt="n") ## ## modelo lmer qqnorm(residuals(modelo), main="residuos del modelo") qqline(residuals(modelo), col="red", lwd=2) ## modelo lme qqnorm(residuals(modlme), main="residuos del modelo") qqline(residuals(modlme), col="red", lwd=2) ## ## test de normalidad ## https://en.wikipedia.org/wiki/Shapiro%E2%80%93Wilk_test ## En vez del test de Kolmogorov-Smirnov (que asume que muestra y población son coincidentes) ## (en la mayoría de los casos tenemos una muestra que no coincide con la población) ## https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test shapiro.test(residuals(modelo)) ## modelo lmer shapiro.test(residuals(modlme)) ## modelo lme ## cuidado con este test: http://www.statisticalmisses.nl/index.php/frequently-asked-questions/77-what-is-wrong-with-tests-of-normality ## ## parametrización de la normalidad de los residuos de devianza ## kurtosis: puntiagudez de la distribución (Ho= 3 ó 0) ## https://en.wikipedia.org/wiki/Kurtosis ## mayor efecto sobre las P's en la violación de la normalidad ## si K>0 (leptokurtosis) hay un mayor error tipo II (aceptar la Ho [nula] cuando de hecho es falsa) ## si K<0 (platikurtosis) hay un mayor error tipo I (rechazar la Ho [nula] cuando es cierta) kurtosis(residuals(modelo)) ## kurtosis de Pearson, Ho = 3 anscombe.test(residuals(modelo)) ## ## sesgo: simetría de la distribución (Ho= 0) ## https://en.wikipedia.org/wiki/Bias ## poco efecto sobre las P's en la violación de la normalidad ## leve aumento del error de tipo I (rechazar la Ho [nula] cuando es cierta) skewness(residuals(modelo)) ## sesgo Ho = 0 agostino.test(residuals(modelo)) ## para plots identificando valores con residuos extremos, ponemos en id la proporción de residuos extremos qqmath(modelo, id=0.05) ## para modelos lmer plot(modlme, id=0.05) ## para modelos lme ## ## "outliers"; solo para modelos lmer outliers.cook <- check_outliers(modelo, method="cook") plot(outliers.cook) ## otra opción más elaborada: distancia de Cook y leverage plot(hatvalues(modelo), cooks.distance(modelo)) abline(h=4/length(modelo@frame[,1]), col="red") abline(v=2*length(fixef(modelo))/length(residuals(modelo)), col="red") identify(hatvalues(modelo), cooks.distance(modelo)) ## si no se ven dos líneas rojas es que no se supera un umbral de preocupación ## otra opción visual con el paquete {car} ## nos fijamos en las observaciones con residuos studentizados (aprox. Student-t) >2 y valores altos de leverage y distancia de Cook influencePlot(modelo, main="distancia de Cook: tamaño del círculo") ## con un listado de las observaciones potencialmente problemáticas ## heterrocedasticidad de los residuos ## con lmer plot(fitted(modelo), residuals(modelo), xlab="valores predichos", ylab="residuos", main="¿HAY HETEROCEDASTICIDAD?") abline(h=0, lty=2, lwd=2, col="red") lines(smooth.spline(fitted(modelo), residuals(modelo), spar=0.95), col="blue") ## ## otro modo más directo plot(modelo) ## con lmer plot(modlme) ## con lme ## podemos abordar la heterocedasticidad de los residuos usando modelos mixtos lme(...) ¡¡¡ no con lmer(...) !!! ## modelo Heteroskedasticity-Corrected (hay otras opciones: mirad en Help varPower, varExp, varIdent, varFunc) ## https://quantdev.ssri.psu.edu/sites/qdev/files/ILD_Ch06_2017_MLMwithHeterogeneousVariance.html plot(modlme, main="modelo original") modlme.HC0 <- update(modlme, weights=varExp()) plot(modlme.HC0, main="modelo HC0") modlme.HC1 <- update(modlme, weights=varExp(form=~t_sonda)) plot(modlme.HC1, main="modelo HC1") modlme.HC2 <- update(modlme, weights=varExp(form=~t_sonda|id)) plot(modlme.HC2, main="modelo HC2") modlme.HC3 <- update(modlme, weights=varFixed(~t_sonda)) plot(modlme.HC3, main="modelo HC3") modlme.HC4 <- update(modlme, weights=varIdent(~1|id)) plot(modlme.HC4, main="modelo HC4") modlme.HC5 <- update(modlme, weights=varComb(varIdent(form=~1|id), varExp(form=~t_sonda))) plot(modlme.HC5, main="modelo HC5") ## comparamos los modelos que intentan corregir la heterocedasticidad residual AICc(modlme.HC0, modlme.HC1, modlme.HC2, modlme.HC3, modlme.HC4, modlme.HC5) ## comparamos con el modelo lme original frente a ## los modelos que en los plots ya no muestran heterocedasticidad residual AICc(modlme, modlme.HC0, modlme.HC2, modlme.HC5) ## ANALISIS GLOBAL DE RESIDUOS de manera sintética para la revisión de los modelos lmer: mcp.fnc(modelo) ## con {LMERConvenienceFunctions} check_model(modelo) ## con {performance} ## con modelos lmer es posible efectuar una diagnosis utilizando residuos DHARMa standarizados en 0-1 ## normalidad, heterocedasticidad de los residuos y existencia de outliers testUniformity(modelo) testQuantiles(modelo) testOutliers(modelo) ## sintéticamente residuos.lmer <- simulateResiduals(modelo) plot(residuos.lmer) ## valoramos la correspondencia entre residuos clásicos y los de DHARMA plot(residuos.lmer$scaledResiduals, residuals(modelo)) plot(rank(fitted(modelo))/dim(modelo@frame)[1], residuos.lmer$scaledResiduals, ylab="residuos DHARMa", xlab="predicciones rango-transformadas", main="¿HAY HETEROCEDASTICIDAD?") ## HOMOGENEIDAD DE VARIANZAS ENTRE LOS GRUPOS DE LOS PREDICTORES NOMINALES (factores) ## con el paquete {DHARMa} para los residuos simulados de modelos lmer testCategorical(modelo, datos$fuente) ### ¡¡ ojo !! poned aquí VUESTROS DATOS testCategorical(modelo, datos$especie) ### ¡¡ ojo !! poned aquí VUESTROS DATOS testCategorical(modelo, datos$fuente:datos$especie) ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## varianzas de las diferentes celdas de interacción entre factores ### ¡¡ ojo !! poned aquí VUESTROS DATOS miformula <- as.formula(a_por_b ~ fuente*especie) aggregate(miformula, data=datos, FUN=var) ## varianza aggregate(miformula, data=datos, FUN=mean) ## media aggregate(miformula, data=datos, FUN=length) ## tamaño muestral ## 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(a_por_b ~ especie, data=datos, center=median) ### ¡¡ ojo !! poned aquí VUESTROS DATOS leveneTest(a_por_b ~ especie, data=datos, center=mean) ### ¡¡ ojo !! poned aquí VUESTROS DATOS bartlett.test(a_por_b ~ especie, data=datos) ### ¡¡ ojo !! poned aquí VUESTROS DATOS fligner.test(a_por_b ~ especie, data=datos) ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## para dos factores o más leveneTest(a_por_b ~ especie:fuente, data=datos, center=median) ### ¡¡ ojo !! poned aquí VUESTROS DATOS leveneTest(a_por_b ~ especie:fuente, data=datos, center=mean) ### ¡¡ ojo !! poned aquí VUESTROS DATOS bartlett.test(a_por_b ~ interaction(especie, fuente), data=datos) ### ¡¡ ojo !! poned aquí VUESTROS DATOS fligner.test(a_por_b ~ interaction(especie, fuente), 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 miformula <- as.formula(a_por_b ~ fuente*especie) num.col <- dim(aggregate(miformula, data=datos, FUN=mean))[2] medias.originales <- aggregate(miformula, data=datos, FUN=mean)[num.col] varianzas.originales <- aggregate(miformula, data=datos, FUN=var)[num.col] 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=paste("r entre medias y varianzas =",cor.media.varianza)) abline(lm((varianzas.originales[,1]) ~ medias.originales[,1]), col="red", lwd=2) par(mfcol=c(1,1)) ## MULTICOLINEARIDAD ENTRE PREDICTORES ## para modelos lmer check_collinearity(modelo) plot(check_collinearity(modelo)) ## para modelos lme solo para random intercept models check_collinearity(modlme) plot(check_collinearity(modlme)) ## OMNIBUS TESTS de la "significación" global del modelo ## Utilizamos la aproximación de los "likelihood ratio tests" (lrtest) utilizando Maximum Likelihood (ML) ## https://en.wikipedia.org/wiki/Likelihood-ratio_test ## Se comparan modelos anidados, como son el modelo de interés y el MODELO NULO ## el MODELO NULO es el que tiene la misma parte aleatoria, pero en la fija no tiene ningún efecto (sólo el intercepto) ## También podríamos comparar dos modelos de interés que tuviesen diferentes subconjuntos de efectos fijos, y misma parte aleatoria ## modelos lmer modelo.ML <- update(modelo, REML=FALSE) modelo.nulo.ML <- update(modelo.ML, .~1+(1|id), REML=FALSE) ## ¡¡¡PONED AQUI VUESTROS DATOS!!! con la misma parte aleatoria ## modelos lme modlme.ML <- update(modlme, method="ML") modlme.nulo.ML <- update(modlme.ML, .~1, random=list(id=~1), method="ML") ## ¡¡¡PONED AQUI VUESTROS DATOS!!! con la misma parte aleatoria ## ponemos antes el modelo más sencillo (i.e., menos grados de libertad, Df); también podemos hacerlo con anova(..., ...) lrtest(modelo.nulo.ML, modelo.ML) lrtest(modlme.nulo.ML, modlme.ML) ## ## Aproximación basada en métodos de remuestreo para comparar modelos anidados. Solo para modelos lmer(...). ## ponemos antes el modelo más complejo y luego el más sencillo ## con nsim definimos el número de procesos de bootstrapping ## con seed podemos cambiar los procesos de aleatorización-remuestreo PBmodcomp(modelo.ML, modelo.nulo.ML, nsim=1000, seed=123, details=0) ## SIGNIFICACION DE EFECTOS DEL MODELO ## lo siguiente es lo mismo que summary(modelo, ddf="Satterthwaite") ## sobre todo nos fijaremos en los coeficientes de los "Fixed effects" ## también es muy interesante la matriz "Correlation of Fixed Effects" ## cuyos valores de correlación entre los efectos fijos deberían ser próximos a cero summary(modelo) summary(modlme) ## solo con la tabla de efectos fijos ## para modelos lmer con aproximación "Satterthwaite" (por defecto) y "Kenward-Roger" ## para los modelos lme son los grados de libertad marginalñes round(summary(modelo, ddf="Kenward-Roger")$coefficients, 5) round(summary(modelo)$coefficients, 5) round(summary(modlme)$tTable, 5) ## COEFCIENTES DE REGRESIÓN ESTANDARIZADOS PARA LAS COVARIANTES (predictores continuos) ## zeta-standarizamos (a MEDIA=0 Y SD=1) ## para obtener los coeficientes estandarizados del modelo (hace comparables los coeficientes para valorar su importancia relativa EN VALOR ABSOLUTO) ## modelos lmer std.coef(modelo, partial.sd=T) ## coeficientes de regresión para la parte fija ESTANDARIZADOS fixef(modelo) ## coeficientes de regresión para la parte fija SIN ESTANDARIZAR ## modelos lme ## no funciona en los modelos con más de un término random en la lista de random=list(tu.efecto.random=~...) std.coef(modlme, partial.sd=T) ## coeficientes de regresión para la parte fija ESTANDARIZADOS fixef(modlme) ## coeficientes de regresión para la parte fija SIN ESTANDARIZAR ## TEST DE EFECTOS ALEATORIOS: valores de Chi2 y sus p-values para los efectos aleatorios ## modelos lmer rand(modelo) ## modelos lme anova(gls(formula(modlme), data=datos, method="REML"), modlme) ## VARIANZAS de los términos aleatorios summary(modelo)$varcor ## modelo lmer VarCorr(modlme) ## modelo lme ## VALORES de los efectos aleatorios qqmath(ranef(modelo)) ## solo para modelos con un único término aletorio plot(ranef(modlme)) ## TEST DE EFECTOS FIJOS ## opción ideal cuando tengamos como predictores a factores con tres niveles o más ## significaciones de los Fixed Effects usando la aproximación Kenward-Roger ## proporciona la Sum Sq de los efectos fijos ## los DenDf son los mismos que los obtenidos en df de summary ## los valores de la F son idénticos a los cuadrados de los t value de las salidas summary ## los valores de significaciones (p) son idénticos a los obtenidos en summary ## no proporciona los interceptos ## Kenward-Roger modification of the F-statistic for some linear mixed models fitted with lmer ## http://web.warwick.ac.uk/statsdept/useR-2011/abstracts/290311-halekohulrich.pdf ## utilizamos una versión de "anova" (anova.merModLmerTest) aplicada a objetos merMod, no equivalente a la usada en objetos lm o glm anova(modelo, ddf="Kenward-Roger", type=3) anova(modelo, ddf="Satterthwaite", type=3) anova.lme(modlme, type="marginal")[-1,] ## ## con el comando Anova de {car} se obtienen directamente los resultados de significación usando la aproximación Kenward-Roger ## proporciona las significaciones para el intercepto ## los resultados de F, p y Df son idénticos a los de anova(..., ddf="Kenward-Roger", type=3) Anova(modelo, type=3, test="F") ## no lo aplicamos a modelos lmer o lme con el test="Chisq"; produce resultados erróneos ## INTERVALOS DE CONFIANZA DE LOS COEFICIENTES ## Podemos efectuar parametrizaciones de nuestro modelo mixto lmer utilizando ## procedimientos de remuestreo (bootstrapping), que también serán de utilidad ## para obtener, de otro modo, estimas de significación (si el intervalo no incluye el valor cero). ## Para cuando nuestro modelo no proporciona significaciones. ## cambiamos level=0.99 para p=0.01. Más rápido con la estima asintótica de "profile": confint.merMod(modelo, parm="beta_", level=0.95, method="profile") confint.fixed <- confint.merMod(modelo, parm="beta_", level=0.95, method="boot", boot.type="perc", nsim=500) print(confint.fixed, digits=3) ## para modelos lme en el caso de que haya un solo término en random=list(...); es igual que: intervals(modlme, levels=0.95)$fixed ## BOOTSTRAPPING DEL MODELO ## Model-based (Semi-)Parametric Bootstrap for Mixed Models; para modelos lmer mcmc.fixed <- bootMer(modelo, FUN=fixef, nsim=500, type="parametric", use.u=FALSE, seed=1) coef.mcmc.fixed <- as.data.frame(mcmc.fixed$t) ## resumen con coeficientes del modelo de interés (original), la media de los coeficientes del bootstrap (bootMed), ## su error estándard (bootSE) y el bias (bootBias) summary(mcmc.fixed) ## tabla de resultados con cuantiles para alfa 0.05 y 0.01 tabla.mixedmodel.simulada <- describe(coef.mcmc.fixed, quant=c(0.025, 0.975, 0.005, 0.995))[,c(1:4, 14:17)] colnames(tabla.mixedmodel.simulada)[4] <- "std_error" tabla.mixedmodel.simulada$coeficiente.modelo <- round(fixef(modelo), 5) tabla.mixedmodel.simulada$std_error.modelo <- round(summary(modelo, ddf="Kenward-Roger")$coefficients[,2], 5) tabla.mixedmodel.simulada$P.modelo <- round(summary(modelo, ddf="Kenward-Roger")$coefficients[,5], 5) ## comparar los coeficientes obtenidos por bootstrapping (en mean) con los originales del modelo (coeficiente.modelo) ## si los valores de kurtosis y sesgo no son cercanos a cero eso es indicativo de la existencia de valores extremos (influyentes, perdidos, ...) ## utilizaremos los intervalos de confianza derivados de cuantiles (95%: Q0.025 y Q0.975; 99%: Q0.005 y Q0.995) ## que no incluyan el valor "cero" print(tabla.mixedmodel.simulada, digits=5) ## Valoración del grado de correlación entre los coeficientes del modelo; los valores de r obtenidos se deberían aproximar a "cero" si los efectos fijos fuesen independientes entre si cor(coef.mcmc.fixed) ## ¿CÚANTO EXPLICA EL MODELO? ## VALORES OBSERVADOS Y PREDICHOS ## con modelos lmer predicho <- fitted(modelo) plot(modelo@frame[,1] ~ predicho, ylab="RESPUESTA", xlab="PREDICCIONES") abline(lm(modelo@frame[,1] ~ predicho), col="green", lwd=2) ## con modelos lme (no efectuamos esta estima con modelos lme corregidos por heterocedasticidad) predicho <- fitted(modlme) plot(get.response(modlme) ~ predicho, ylab="RESPUESTA", xlab="PREDICCIONES") abline(lm(get.response(modlme) ~ predicho), col="green", lwd=2) ## ## PARTICIÓN DE LA VARIANZA ## R2m: marginal R2 (the proportion of variance explained by the fixed factor(s) alone) ## R2c: conditional R2 (the proportion of variance explained by both the fixed and random factors; i.e. the entire model) ## consultad: https://jonlefcheck.net/2013/03/13/r2-for-linear-mixed-effects-models/comment-page-1/ ## http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210x.2012.00261.x/epdf ## http://onlinelibrary.wiley.com/doi/10.1111/2041-210X.12225/epdf r.squaredGLMM(modelo) ## modelo lmer performance(modelo) r.squaredGLMM(modlme) ## modelo lme (no efectuamos esta estima con modelos lme corregidos por heterocedasticidad) ## PARTICION DE LA VARIACION EN LA RESPUESTA ENTRE EFECTOS FIJOS con lmer ## tabla sintética de los resultados con la partición de la variabilidad de la respuesta (y magnitudes de efectos parciales) ## utilizando la estima de grados de libertad de Kenward-Roger ## eta2 es el tanto por uno de la variación en la respuesta explicado exclusivamente por ese efecto fijo SStotal <- sum((modelo@frame[,1]-mean(modelo@frame[,1]))^2) ## incluir aquí manualmente el nombre de la variable respuesta SSerror <- SStotal*(1-as.numeric(r.squaredGLMM(modelo)[2])) tabla.SS <- as.data.frame(anova(modelo, ddf="Kenward-Roger", type=3)) eta2 <- tabla.SS[,1]/SStotal tabla.SS <- data.frame(eta2, tabla.SS) names(tabla.SS)[7] <- "Pr.F" print(tabla.SS, digits=4) var.compartida <- (as.numeric(r.squaredGLMM(modelo)[1])-sum(tabla.SS$eta2)) print(c("proporción de la varianza atribuible a los efectos fijos =", round(as.numeric(r.squaredGLMM(modelo)[1]), 4)), quote=FALSE) print(c("proporción de la varianza compartida por los efectos fijos =", round(var.compartida, 4)), quote=FALSE) ## TESTS POST-HOC library(emmeans) ## para tests post hoc library(phia) ## para medias de interacciones de factores ## otros métodos de ajuste son: "scheffe", "fdr", "holm" ## ¡¡¡ PONED AQUI VUESTROS DATOS !!! (nombre del factor en c("...")) ## modelo lmer pairs(emmeans(modelo, c("especie")), adjust="tukey") pairs(emmeans(modelo, c("fuente", "especie")), adjust="tukey") ## modelo lme pairs(emmeans(modlme, c("especie")), adjust="tukey") pairs(emmeans(modlme, c("fuente", "especie")), adjust="tukey") ## MEDIAS MARGINALES DE LOS NIVELES DE LOS FACTORES ## (controlando por el efecto de las covariantes) interactionMeans(modelo) ## modelo lmer interactionMeans(modlme) ## modelo lme ## ## otra posibilidad ## ¡¡¡ PONED AQUI VUESTROS DATOS !!! emmeans(modelo, ~fuente*especie) ## modelo lmer emmeans(modlme, ~fuente*especie) ## modelo lme ## VISUALIZACION DE EFECTOS ## valores medios para los niveles de los factores ajustados por las covariantes ## ¡¡¡ PONED AQUI VUESTROS DATOS !!! ## modelo lmer plot(emmeans(modelo, ~fuente*especie)) ## modelo lme plot(emmeans(modlme, ~fuente*especie)) detach("package:emmeans") detach("package:phia") ## Partial effects plots con {effects} library(effects) ## elegimos un modelo lmer de interés que contenga a los predictores cuyos efectos parciales sobre la respuesta deseamos estimar plot(allEffects(modelo), rotx=45, rug=F) ## intervalos de confianza del 95% ## para plots de efectos particulares de factores (si queremos ver los residuos parciales: partial.residuals=T) plot(Effect("especie", partial.residuals=F, modelo), lwd=3, lty=1, rotx=45, partial.residuals=list(col="blue", pch=1, cex=1.25), rug=F) plot(Effect(c("fuente", "especie"), partial.residuals=F, modelo), lwd=3, lty=1, rotx=45, partial.residuals=list(col="blue", pch=1, cex=1.25), rug=F) ## para covariantes plot(Effect("t_sonda", partial.residuals=F, modelo), lwd=3, lty=1, rotx=0, partial.residuals=list(col="blue", pch=1, cex=1.25), rug=F) ## ## elegimos un modelo lme(...) plot(allEffects(modlme), rotx=45, rug=F) ## intervalos de confianza del 95% ## para plots de efectos particulares plot(Effect("especie", partial.residuals=F, modlme), lwd=3, lty=1, rotx=45, partial.residuals=list(col="blue", pch=1, cex=1.25), rug=F) plot(Effect(c("fuente", "especie"), partial.residuals=F, modlme), lwd=3, lty=1, rotx=45, partial.residuals=list(col="blue", pch=1, cex=1.25), rug=F) plot(Effect("t_sonda", partial.residuals=F, modlme), lwd=3, lty=1, rotx=0, partial.residuals=list(col="blue", pch=1, cex=1.25), rug=F) ## mirad la ayuda de "plot.effects" en Detailed Argument Descriptions ## mediante el paquete {sjPlot} ## en modelos lme solo cuando son random intercepts ## type="int" para representar solo las interacciones en modelos factoriales ## para mostrar ca. un error estándard en vez del intervalo al 95%: ci.lvl=0.68 ## "std2" muestra los efectos estandarizados tanto para covariantes como para los predictores binarios de los factores ## plot de valores parciales plot_model(modelo, show.data=FALSE, type="eff", ci.lvl=0.95, terms=NULL) ## plot de coeficientes estandarizados (comparables entre si en su magnitud) plot_model(modelo, show.data=FALSE, type="std2", ci.lvl=0.95, terms=NULL) plot_model(modelo, show.data=FALSE, type="std2", ci.lvl=0.95, terms=NULL, show.values=TRUE, value.offset=0.4) ## ejemplos por variables concretas plot_model(modelo, show.data=FALSE, type="eff", ci.lvl=0.95, terms="t_sonda") plot_model(modelo, show.data=FALSE, type="eff", ci.lvl=0.95, terms="especie") plot_model(modelo, show.data=FALSE, type="eff", ci.lvl=0.95, terms=c("t_sonda", "fuente")) plot_model(modelo, show.data=FALSE, type="eff", ci.lvl=0.95, terms=c("fuente", "especie")) plot_model(modelo, show.data=FALSE, type="int", ci.lvl=0.95, terms=NULL) plot_model(modelo, show.data=FALSE, type="est", ci.lvl=0.95, terms=NULL) ## para diagnosis del modelo plot_model(modelo, show.data=TRUE, type="resid", ci.lvl=0.95, terms=NULL) plot_model(modelo, show.data=TRUE, type="diag", ci.lvl=0.95, terms=NULL) ## detach(package:sjPlot) ## HIPOTESIS NULAS para los coeficientes con las que contrastar nuestro modelo mixto de interés lmer(...) ## Randomization (Permutation) test under the null hypothesis that coefficients are cero. ## Then we compare the observed coefficients in the model of interest with the 1000 permuted coefficient values. ## The justification of the randomization test derives from the fact that under the null hypothesis of no treatment effect, ## the random assignment procedure produces a random shuffle of the responses. ## The justification of the permutation test derives from the fact that under the null hypothesis of identical distributions, ## all permutations of the responses are equally likely. ## Randomization tests consider every possible permutation of the labels, permutation tests take a random sample of permutations of the labels ## para ello aleatorizamos los valores de la variable respuesta por remuestreo-sin-reemplazo: sample(respuesta, replace=FALSE) ## conveniente cuando ha habido desvíos de los supuestos canónicos de los residuos del modelo ## https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2687965/ ## https://link.springer.com/content/pdf/10.1007%2F978-1-4614-1365-3_10.pdf ## https://lmcarrascal.eu/cursos/randtests.pdf ## https://www.sciencedirect.com/topics/mathematics/randomization-test/pdf ## https://en.wikipedia.org/wiki/Permutation_test ## https://www.researchgate.net/publication/321385778_Randomization_Tests_or_Permutation_Tests_A_Historical_and_Terminological_Clarification ## https://stats.stackexchange.com/questions/55742/difference-between-randomization-test-and-permutation-test ## https://www.r-bloggers.com/2021/03/randomization-tests-make-fewer-assumptions-and-seem-pretty-intuitive/ ## https://uoftcoders.github.io/rcourse/lec12-randomization-tests.html ## https://towardsdatascience.com/how-to-use-permutation-tests-bacc79f45749 ## https://measuringu.com/randomization-test/ ## ¡¡¡ MUY LENTO !!! library(bayestestR) ## MODELOS lmer N_boots <- 1000 N_coefs <- length(summary(modelo)$coefficients[,1]) N_Fs <- dim(anova(modelo, ddf="Satterthwaite", type=3))[1] null_coeficientes <- as.data.frame(matrix(9999, nrow=N_boots, ncol=N_coefs)) null_Fs <- as.data.frame(matrix(9999, nrow=N_boots, ncol=N_Fs)) null_R2 <- as.data.frame(matrix(9999, nrow=N_boots, ncol=2)) respuesta <- get.response(modelo) for (i in 1:N_boots) { print(i) datos$respuesta.null <- sample(respuesta, replace=FALSE) lmer_null <- update(modelo, respuesta.null~., REML=TRUE) null_coeficientes[i,] <- summary(lmer_null)$coefficients[,1] null_Fs[i,] <- anova(lmer_null, ddf="Satterthwaite", type=3)[,5] null_R2[i,] <- r.squaredGLMM(lmer_null) } ## comprobad que en la matriz de null_coeficientes no hay valores 9999 indicativos de errores por problemas de estima (singularidad) colnames(null_coeficientes) <- rownames(summary(lmer_null)$coefficients) colnames(null_Fs) <- rownames(anova(modelo, ddf="Satterthwaite", type=3)) colnames(null_R2) <- c("R2m", "R2c") tabla.null <- describe(null_coeficientes, quant=c(0.025, 0.975, 0.005, 0.995))[, c(-1, -5:-10, -13)] tabla.null$Estimate <- summary(modelo)$coefficients[,1] colnames(tabla.null)[3] <- "std_error" F.null <- describe(null_Fs, quant=c(0.90, 0.95, 0.99, 0.995))[,c(2, 14:17)] F.null$Estimate <- anova(modelo, ddf="Satterthwaite", type=3)[,5] R2.null <- describe(null_R2, quant=c(0.90, 0.95, 0.99, 0.995))[,c(2, 14:17)] R2.null$Estimate <- r.squaredGLMM(modelo)[1,] ## que el Estimate del modelo mixto de interés quede fuera de los intervalos al 95% (Q0.025 - Q0.975) ó 99% (Q0.005 - Q0.995) print(tabla.null, digits=3) plot(ecdf(null_coeficientes$t_sonda)); abline(h=c(0.025, 0.975), col="red") ## ¡¡¡ AQUÍ VUESTROS DATOS !!! ## que las F y R2 queden a la derecha de los cuantiles correspondientes a alpha = 1 - quantile print(F.null, digits=3) 1-ecdf(null_Fs$especie)(2.290) ## ¡¡¡ AQUÍ VUESTROS DATOS !!!; para sacar una p para una F observada en el modelo mixto de interés print(R2.null, digits=3) 1-ecdf(null_R2$R2m)(0.475) ## ¡¡¡ AQUÍ VUESTROS DATOS !!!; para sacar una p para una R2 observada en el modelo mixto de interés plot(density(null_R2$R2c), lwd=2, col="red", xlab="R2 total del análisis permulacional", main="", xlim=c(0, 1)); abline(v=R2.null[2,6], col="blue") plot(density(null_R2$R2m), lwd=2, col="red", xlab="R2 efectos fijos del análisis permulacional", main="", xlim=c(0, 1)); abline(v=R2.null[1,6], col="blue") ## ## INTERVALOS con los métodos ETI (Equal-Tailed Interval; por defecto) y HDI (Highest Density Interval) ## usando el método ETI plot(eti(null_R2, ci=0.95)) print(eti(null_coeficientes, ci=0.95), digits=3) ## igual que la salida previa plot(eti(null_coeficientes, ci=0.95)) ## usando el método HDI plot(hdi(null_R2, ci=0.95)) print(hdi(null_coeficientes, ci=0.95), digits=3) plot(hdi(null_coeficientes, ci=0.95)) ## MODELOS lme N_boots <- 1000 N_coefs <- length(summary(modlme)$coefficients$fixed) N_Fs <- dim(anova.lme(modlme, type="marginal"))[1] null_coeficientes_lme <- as.data.frame(matrix(9999, nrow=N_boots, ncol=N_coefs)) null_Fs_lme <- as.data.frame(matrix(9999, nrow=N_boots, ncol=N_Fs)) null_R2_lme <- as.data.frame(matrix(9999, nrow=N_boots, ncol=2)) respuesta <- get.response(modlme) for (i in 1:N_boots) { print(i) datos$respuesta.null <- sample(respuesta, replace=FALSE) lme_null <- update(modlme, respuesta.null~., method="REML") null_coeficientes_lme[i,] <- summary(lme_null)$coefficients$fixed null_Fs_lme[i,] <- anova.lme(lme_null, type="marginal")[,3] null_R2_lme[i,] <- r.squaredGLMM(lme_null) } ## comprobad que en la matriz de null_coeficientes no hay valores 9999 indicativos de errores por problemas de estima (singularidad) colnames(null_coeficientes_lme) <- names(summary(lme_null)$coefficients$fixed) colnames(null_Fs_lme) <- rownames(anova.lme(modlme, type="marginal")) colnames(null_R2_lme) <- c("R2m", "R2c") tabla.null.lme <- describe(null_coeficientes_lme, quant=c(0.025, 0.975, 0.005, 0.995))[, c(-1, -5:-10, -13)] tabla.null.lme$Estimate <- summary(modlme)$coefficients$fixed colnames(tabla.null.lme)[3] <- "std_error" F.null_lme <- describe(null_Fs_lme, quant=c(0.90, 0.95, 0.99, 0.995))[,c(2, 14:17)] F.null_lme$Estimate <- anova.lme(modlme, type="marginal")[,3] R2.null_lme <- describe(null_R2_lme, quant=c(0.90, 0.95, 0.99, 0.995))[,c(2, 14:17)] R2.null_lme$Estimate <- r.squaredGLMM(modlme)[1,] ## que el Estimate del modelo mixto de interés quede fuera de los intervalos al 95% (Q0.025 - Q0.975) ó 99% (Q0.005 - Q0.995) print(tabla.null.lme, digits=3) plot(ecdf(null_coeficientes_lme$t_sonda)); abline(h=c(0.025, 0.975), col="red") ## ¡¡¡ AQUÍ VUESTROS DATOS !!! ## que las F y R2 queden a la derecha de los cuantiles correspondientes a alpha = 1 - quantile print(F.null_lme[-1,], digits=3) 1-ecdf(null_Fs_lme$especie)(2.290) ## ¡¡¡ AQUÍ VUESTROS DATOS !!!; para sacar una p para una F observada en el modelo mixto de interés print(R2.null_lme, digits=3) 1-ecdf(null_R2_lme$R2m)(0.475) ## ¡¡¡ AQUÍ VUESTROS DATOS !!!; para sacar una p para una R2 observada en el modelo mixto de interés plot(density(null_R2_lme$R2c), lwd=2, col="red", xlab="R2 total del análisis permulacional", main="", xlim=c(0, 1)); abline(v=R2.null_lme[2,6], col="blue") plot(density(null_R2_lme$R2m), lwd=2, col="red", xlab="R2 efectos fijos del análisis permulacional", main="", xlim=c(0, 1)); abline(v=R2.null_lme[1,6], col="blue") ## ## ## INTERVALOS con los métodos ETI (Equal-Tailed Interval; por defecto) y HDI (Highest Density Interval) ## usando el método ETI tabla.eti <- as.data.frame(eti(null_coeficientes_lme, ci=0.95)) tabla.eti$Estimate <- tabla.null.lme$Estimate tabla.eti plot(eti(null_coeficientes_lme, ci=0.95)) ## usando el método HDI tabla.hdi <- as.data.frame(hdi(null_coeficientes_lme, ci=0.95)) tabla.hdi$Estimate <- tabla.null.lme$Estimate tabla.hdi plot(hdi(null_coeficientes_lme, ci=0.95)) #################################### ## Estima ROBUSTA del modelo lmer ## #################################### library(robustlmm) ## https://cran.r-project.org/web/packages/robustlmm/vignettes/rlmer.pdf ## métodos: ## DAStau: (default) For this method, the consistency factors are computed using numerical ## quadrature. This is slower but yields more accurate results. This is the direct analogue to ## the DAS-estimate in robust linear regression. ## DASvar: This method computes the consistency factors using a direct approximation ## which is faster but less accurate. For complex models with correlated random effects ## with more than one correlation term, this is the only method available. eqt <- modelo@call ## modelos rlmer ## para modelos mixtos complejos con términos aleatorios correlacionados el método DASvar es el único posible mm.rvar2 <- rlmer(eqt, data=datos, rel.tol=1e-08, max.iter=1000, method="DASvar") mm.rtau2 <- rlmer(eqt, data=datos, rel.tol=1e-08, max.iter=1000, method="DAStau") mm.rvar <- update(mm.rvar2, rho.sigma.e=psi2propII(smoothPsi, k=2.28), rho.sigma.b=psi2propII(smoothPsi, k=2.28)) mm.rtau <- update(mm.rtau2, rho.sigma.e=psi2propII(smoothPsi, k=2.28), rho.sigma.b=psi2propII(smoothPsi, k=2.28)) ## los PESOS DE LAS UNIDADES MUESTRALES en los modelos se obtienen para las observaciones ("w_e") y para los efectos aleatorios ("w_b") getME(mm.rvar, "w_e") getME(mm.rvar, "w_b") par(mfcol=c(1,1)) ## fija un sólo panel gráfico plot(residuals(modelo, type="response"), getME(mm.rvar, "w_e")) plot(scale(residuals(modelo, type="response")), getME(mm.rvar, "w_e")) plot(cooks.distance(modelo), getME(mm.rvar, "w_e")) plot(hatvalues(modelo), getME(mm.rvar, "w_e")) ## RESIDUOS del modelo mcp.fnc(mm.rvar) mcp.fnc(mm.rtau) mcp.fnc(modelo) ## plots de los residuos y pesos del modelo ## demanda en la consola: Hit to see next plot plot(mm.rvar, which=c(1:3)) plot(mm.rtau, which=c(1:3)) ## RESULTADOS del modelo ## varianza de los términos aleatorios summary(mm.rvar)$varcor summary(mm.rtau)$varcor ## ## los modelos rlmer no proporcionan significaciones summary(mm.rvar)$coefficients summary(mm.rtau)$coefficients summary(modelo, ddf="lme4")$coefficients ## SIGNIFICACIONES de los coeficientes del modelo ## como rlmer no proporciona significaciones, ## podríamos utilizar los grados de libertad proporcionados por el modelo lmer de interés ## ## modelo robusto rlmer DASvar con significaciones asociadas a los grados de libertad tabla.mm.rvar <- as.data.frame(summary(mm.rvar)$coefficients) tabla.mm.rvar$dfKR <- round(summary(modelo, ddf="Kenward-Roger")$coefficients[,3], 1) tabla.mm.rvar$P_dfKR <- (1-pt(abs(tabla.mm.rvar[,3]), df=tabla.mm.rvar[,4])) round(tabla.mm.rvar, 4) ## ## modelo robusto rlmer DAStau con significaciones asociadas a los grados de libertad tabla.mm.rtau <- as.data.frame(summary(mm.rtau)$coefficients) tabla.mm.rtau$dfKR <- round(summary(modelo, ddf="Kenward-Roger")$coefficients[,3], 1) tabla.mm.rtau$P_dfKR <- (1-pt(abs(tabla.mm.rtau[,3]), df=tabla.mm.rtau[,4])) round(tabla.mm.rtau, 4) ## ## modelo original lmer sin pesar las observaciones: round(summary(modelo, ddf="Kenward-Roger")$coefficients, 4) #################################### ## OTRO JUEGO DE DATOS Y EJEMPLOS ## #################################### ## datos de internet meconecto <- url("http://www.lmcarrascal.eu/cursos/nwayMIXED_ANOVA_1.RData") ## abro una conexión con internet load(meconecto) ## cargo el archivo de la conexión close(meconecto); rm(meconecto) ## cierro la conexión options(contrasts=c(factor="contr.sum", ordered="contr.poly")) control.lme <- lmeControl(maxIter=200, msMaxIter=200, msVerbose=TRUE, opt="optim") ## el optimizador opt="nlminb" hay veces que da problemas control.lmer <- lmerControl(check.conv.grad=.makeCC(action ="ignore", tol=1e-6, relTol=NULL), optimizer="bobyqa", optCtrl=list(maxfun=100000)) ## Intraclass correlation: the "unconditional means model" (or "null model"); Random intercept with fixed mean ## diferencias entre los niveles del factor aleatorio (individuo) lmer.0 <- lmer(ln_uso10h ~ 1 + (1|individuo), data=datos, control=control.lmer, REML=FALSE) ## lme.0 <- lme(ln_uso10h ~ 1, random=list(individuo=~1), data=datos, control=control.lme, method="ML") ## modelo de RANDOM INTERCEPT & FIXED SLOPES lmer.1 <- lmer(ln_uso10h ~ temperatura+spp*distancia + (1|individuo), data=datos, control=control.lmer, REML=FALSE) ## lme.1 <- lme(ln_uso10h ~ temperatura+spp*distancia, random=list(individuo=~1), data=datos, control=control.lme, method="ML") ## modelo de FIXED INTERCEPT & RANDOM SLOPES, with UNCORRELATED random slopes for different covariates lmer.2 <- lmer(ln_uso10h~temperatura+spp*distancia + (0+temperatura|individuo)+(0+distancia|individuo), data=datos, control=control.lmer, REML=FALSE) ## lme.2 <- lme(ln_uso10h ~ temperatura+spp*distancia, random=list(individuo=~0+temperatura, individuo=~0+distancia), data=datos, control=control.lme, method="ML") lme.2b <- lme(ln_uso10h ~ temperatura+spp*distancia, random=list(individuo=pdDiag(~0+temperatura+distancia)), data=datos, control=control.lme, method="ML") ## modelo de UNCORRELATED RANDOM INTERCEPT & RANDOM SLOPES, with UNCORRELATED random slopes for different covariates; elimina la covarianza intercepto-pendientes lmer.3 <- lmer(ln_uso10h ~ temperatura+spp*distancia + (1|individuo)+(0+temperatura|individuo)+(0+distancia|individuo), data=datos, control=control.lmer, REML=FALSE) ## lme.3 <- lme(ln_uso10h ~ temperatura+spp*distancia, random=list(individuo=~1, individuo=~0+temperatura, individuo=~0+distancia), data=datos, control=control.lme, method="ML") lme.3b <- lme(ln_uso10h ~ temperatura+spp*distancia, random=list(individuo=pdDiag(~1+temperatura+distancia)), data=datos, control=control.lme, method="ML") ## modelo de CORRELATED RANDOM INTERCEPT & RANDOM SLOPES, with CORRELATED random slopes for different covariates lmer.4 <- lmer(ln_uso10h ~ temperatura+spp*distancia + (1+temperatura+distancia|individuo), data=datos, control=control.lmer, REML=FALSE) ## lme.4 <- lme(ln_uso10h ~ temperatura+spp*distancia, random=list(individuo=pdSymm(~1+temperatura+distancia)), data=datos, control=control.lme, method="ML") ## modelo de CORRELATED RANDOM INTERCEPT & RANDOM SLOPES, with UNCORRELATED random slopes for different covariates lmer.5 <- lmer(ln_uso10h ~ temperatura+spp*distancia + (1+temperatura|individuo)+(1+distancia|individuo), data=datos, control=control.lmer, REML=FALSE) ## lme.5 <- lme(ln_uso10h ~ temperatura+spp*distancia, random=list(individuo=pdSymm(~1+temperatura), individuo=pdSymm(~1+distancia)), data=datos, control=control.lme, method="ML") ###################################################### ## TRANSFORMACION DE BOX-COX ## ## búsqueda automatizada del valor óptomo de lambda ## ###################################################### ## TRANSFORMACIÓN BOX-COX ## función boxCox del paquete {car} ## 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)) ## corremos un modelo general lineal frecuentista (lm) formula(modelo) ## eliminamos los términos aleatorios ## ¡¡¡ PONED AQUÍ VUESTROS DATOS !!! eqt <- as.formula(I(uso10h+1) ~ temperatura + spp * distancia) mod_bc <- lm(eqt, data=datos) ## comenzamos el proceso automatizado de búsqueda de la lamba para la transformación óptima ## si no observamos un máximo, modificamos el rango de lambdas en lambda=seq(-2,2, 1/100) print(boxCox(mod_bc, lambda=seq(-2,2, 1/100))) lambda.bc <- with(boxCox(mod_bc, 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 incluida en datos datos$respuesta.bc <- ((mod_bc$model[,1]^lambda.bc)-1)/(lambda.bc) ## para explorar la relación entre la respuesta original y su transformación Box-Cox plot(datos$respuesta.bc, mod_bc$model[,1], xlab="variable respuesta transformada", ylab="variable respuesta original") ## renombramos y guardamos el modelo stan_lmer previo modelo.previo <- modelo ## ahora rehacemos el modelo stan_lmer ## sustituyendo la variable respuesta del modelo por datos$respuesta.bc modelo <- update(modelo, respuesta.bc~., na.action=na.omit) ## EXPLICACIONES, RESUMEN Y SUGERENCIAS ## Los modelos mixtos lme(...) tiene algunas limitaciones a la hora de especificar términos aleatorios complejos (efectos aleatorios cruzados). ## Los modelos mixtos lmer(...) gestionan de manera sencilla esos efectos aleatorios cruzados (e.g.: 1|factor_aleatorio_1 + 1|factor_aleatorio_2) ## Los modelos mixtos lme(...) utilizan una aproximación clásica de grados de libertad de diseños balanceados ANOVA n-way ## Los modelos mixtos lmer(...) utilizan dos aproximaciones diferentes: Satterthwaite y Kenward-Roger. ## La de Kenward-Roger es la más restrictiva respecto a los grados de libertad del denominador de los efectos fijos al usar "random slopes" ## Los modelos mixtos lmer(...) no gestionan la heterocedasticidad de los residuos. ## Los modelos mixtos lme(...) gestionan la heterocedasticidad de los residuos a través del argumento weights=varExp(), varFixed(), varIdent(), ... ## Los modelos mixtos lme(...) también gestionan situaciones de autocorrelación (e.g., espacial y temporal) mediante el argumento cor=corAR1(), corARMA(), corEXp, corSpher, ... ## Debemos evitar situaciones en las cuales la varianza asociada a los términos aleatorios sea muy pequeña, tendente a cero. ## Este asunto conlleva estimas singulares, al ser la estructura de efectos aleatorios es demasiado compleja para ser apoyada por los datos. ## Eliminaremos la parte más compleja de la estructura de los efectos aleatorios (generalmente "random slopes"), evitando sobreparametrizarlos-ajustarlos ## Para ello miraremos el resultado de varianzas de la parte aleatoria: ## Std.Dev. de summary(lmer(...)) o StdDev de summary(lme(...)) en la parte de Random effects: ## Si no hemos tenido un buen ajuste a los supuestos canónicos de los modelos mixtos (normalidad, homocedasticidad de los residuos, puntos outliers) ## contruiremos hipótesis nulas para examinar la magnitud de los coeficientes de los efectos fijos (cuyos valores deben estar fuera de los intervalos de confianza nulos) ## obtendremos los parámetros de los modelos mixtos mediante bootstraps (cuyos intervalos de confianza no deben incluir el valor cero) ## Anova{car}: (Type III Wald F tests with Kenward-Roger df) ## http://ir.library.oregonstate.edu/xmlui/bitstream/handle/1957/5262/mydissertation.pdf ## Instead of approximating the Wald-type statistic itself by an F distribution as in the Satterthwaite approximation, ## Kenward and Roger considered a scaled form of the statistic that makes the approximation more flexible since ## they approximate the denominator degrees of freedom and the scale as well. Second, Kenward and Roger ## modified the procedure in such a way that they get the right values for both the denominator degrees of freedom ## and the scale for the two special cases considered in chapter 3. Third, the adjusted estimator of the ## variance-covariance matrix of the fixed effects estimator, which is less biased than the conventional estimator, ## as we saw in chapter 2, is used in constructing the Wald-type statistic. ## The K-R, the Satterthwaite, and the proposed methods produce the same estimate of the denominator degrees of ## freedom, and the scale estimate is one. Even though the estimates are the same, the Satterthaite and the K-R tests ## are not necessarily identical, and this is true because in the K-R statistic, we use the adjusted estimator of the ## variance-covariance matrix of the fixed effects estimator, whereas the conventional estimator is used in the Satterthwaite test. ## We conducted a simulation study for three kinds of block designs; partially balanced block designs, balanced incomplete block ## designs, and complete block designs with missing data. The K-R, the proposed, the Satterthwaite, and the Containment methods ## were compared in the study. Three factors were considered in the study: the variance components ratio, the efficiency factor, ## and the sample size. In most cases, we found the K-R and the proposed methods performed as well as or better than the other ## methods. For designs with small values for all the three factors, the Satterthwaite method performed poorly. ## Kenward-Roger modification of the F-statistic for some linear mixed models fitted with lmer ## http://web.warwick.ac.uk/statsdept/useR-2011/abstracts/290311-halekohulrich.pdf