##################################################################### ## MODELOS MIXTOS GENERALIZADOS MULTINOMIALES LINEALES (ordinales) ## ##################################################################### ## Luis M. Carrascal ## https://lmcarrascal.eu ## revisión 20/02/2022 ## https://cran.r-project.org/web/packages/ordinal/vignettes/clm_article.pdf ## https://cran.r-project.org/web/packages/ordinal/vignettes/clmm2_tutorial.pdf ## datos de ejemplo 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 ## CARGAMOS PAQUETES library(ordinal) ## para modelos mixtos multinomiales library(MASS) ## para usar dropterm y polr library(MuMIn) ## para calcular AICc library(lmtest) ## para obtener la p del modelo, en sustitución de anova(nulo, model, test="Chisq") library(sandwich) ## para estima robusta, estimando sandwich standard errors library(car) ## para Anova.clm ## PREPARACIÓN ## cargamos tablas de contrastes compartidas con las de otros programas estadísticos (como SPSS, STATISTICA) options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## ## a modo de ejemplo recodificamos una variable continua para convertirla en una multinomial ordinal ## creamos ahora una multinomial ordinal, poniendo niveles de corte a la variable fuente en <=0 (de -1 a 0), 1-60, 61-120 y 121-1000 datos$uso10h_mnord <- as.ordered(cut(datos$uso10h, c(-1,0,60,120,1000), labels=c("nada","poco","algo","mucho"))) ## ESTABLECEMOS EL MODELO ## ESTABLECEMOS EL MODELO ## creamos nuestras ecuaciones ## modelo nulo ("unconditional means model") eqt.nulo <- as.formula(uso10h_mnord ~ 1 + (1|individuo)) ## random intercept & fixed slopes eqt.randint <- as.formula(uso10h_mnord ~ temperatura+distancia*spp + (1|individuo)) ## random intercept & random (uncorrelated) slopes eqt.randintslo <- as.formula(uso10h_mnord ~ temperatura+distancia*spp + (temperatura|individuo)+(ln_luz|individuo)+(distancia|individuo)) ## random intercept & random (correlated) slopes eqt.randintslo2 <- as.formula(uso10h_mnord ~ temperatura+distancia*spp + (temperatura+ln_luz+distancia|individuo)) ## ## creamos MODELOS GENERALIZADOS MIXTOS LINEALES MULTINOMIALES ORDINALES ## link puede ser "logit" y "cloglog" ## nAGQ=0 no da problemas, aunque utiliza una forma más rápida pero menos exacta de estimación ## de parámetros para los GLMM (usando penalized iteratively reweighted least squares) ## al optimizar los efectos aleatorios y los coeficientes de efectos fijos en la estima penalizada IrWLS, ## y no es válida para modelos con más de un término aleatorio ## los modelos random intercept plus random slopes son muy lentos al ejecutarse modelo.nulo <- clmm(eqt.nulo, data=datos, link="logit", threshold="flexible", nAGQ=1) modelo.randint <- clmm(eqt.randint, data=datos, link="logit", threshold="flexible", nAGQ=1) modelo.randintslo <- clmm(eqt.randintslo, data=datos, link="logit", threshold="flexible", nAGQ=1) modelo.randintslo2 <- clmm(eqt.randintslo2, data=datos, link="logit", threshold="flexible", nAGQ=1) ## COMPARACION DE MODELOS ## Aproximación basada en "información" usando AICc (válida porque con estamos operando con ML) ## seleccionamos el modelo más adecuado utilizando el criterio de información de Akaike valores_AICc <- na.omit(AICc(modelo.nulo, modelo.randint, modelo.randintslo, modelo.randintslo2)) valores_AICc Weights(valores_AICc) ## OMNIBUS TEST comparando el modelo de interés con el modelo nulo ## Utilizamos la aproximación de los "likelihood ratio tests" (lrtest) habiendo usado ML ## Se comparan modelos anidados, poniendo antes el modelo más sencillo (i.e., menos grados de libertad, df) ## test para la relevancia de los efectos fijos (usando como modelo base el de random intercepts) ## OMNIBUS TEST comparando el modelo de interés con el modelo nulo lrtest(modelo.nulo, modelo.randint) anova(modelo.nulo, modelo.randint) ## seleccionamos el modelo de interés modelo <- modelo.randint ## RESIDUOS ## leed al autor del paquete: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2013q4/021281.html ## RESULTADOS summary(modelo) ## solo la tabla de regresión de los efectos fijos asumiendo Inf grados de libertad summary(modelo)$coefficients coeftest(modelo) ## intervalo de confianza de los coeficientes round(confint(modelo, level=0.95),3) ## To convert the coefficients into odds ratios, we just exponentiate the estimates and confidence intervals: ## mirad https://en.wikipedia.org/wiki/Odds_ratio ## descartad por dificultad de interpretación las primeras líneas de los interceptos (xxx|yyy) round(exp(confint(modelo, level=0.95)), 3) ## Likelihood ratio test con devianzas para efectos (indicada para cuando existan factores) ## No estima los efectos principales usando el equivalente a SS de tipo III (efectos parciales totales) si hay interacción. ## lo hacemos directamente con Anova de {car} usando el equivalente a SS type-II Anova(modelo, type=2, test="Chisq") ## Si no se ejecuta podemos utilizar la aproximación de likelihood ratio tests drop1(modelo, test="Chisq") ## de estos resultados solo nos quedamos con la estima del efcto principal ## como no produce resultados para los efectos principales de factores que interaccionan ## construimos un modelo de efectos principales drop1(update(modelo, .~.-distancia:spp), test="Chisq") ## VARIACION EXPLICADA por el modelo (usando devianzas) d2 <- round(100*((-2*modelo.nulo$logLik)-(-2*modelo$logLik))/(-2*modelo.nulo$logLik), 2) print(c("D2 de MCFadden (%)=", d2), quote=FALSE) ## ANÁLISIS DE VALORES OBSERVADOS Y PREDICHOS POR EL MODELO plot(modelo$model[,1], 1-modelo$fitted.values, xlab="NIVELES DE LA RESPUESTA", ylab="PROBABILIDAD PREDICHA") ## VISUALIZACION DE LAS MEDIAS DE LOS NIVELES DE LOS FACTORES EN LA ESCALA LOGIT library(emmeans) ## ¡¡¡ poned aquí VUESTROS DATOS !!! tabla_medias <- as.data.frame(emmeans(modelo, ~spp:distancia)) ## probabilidades de ocurrencia y de los intervalos de confianza asintóticos tabla_medias$probabilidad <- round(exp(tabla_medias$emmean), 3) tabla_medias$limINF95 <- round(exp(tabla_medias$asymp.LCL), 3) tabla_medias$limSUP95 <- round(ifelse(exp(tabla_medias$asymp.UCL)>1, 1, exp(tabla_medias$asymp.UCL)), 3) tabla_medias[,c(1:2,8:10)]