####################################### ## DISEÑOS FACTORIALES GENERALIZADOS ## ## n-way AN(c)OVA ## ## CON DISTRIBUCIONES POISSON ## ## Y BINOMIALES NEGATIVAS ## ## 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 ponerlo en uso en "environment" meconecto <- url("http://www.lmcarrascal.eu/cursos/nwayGENERALIZADO.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 MUCHAS OTRAS COSAS OS REMITO AL script "GLM factorial designs - n-way AN(c)OVA.R" ## CARGAMOS PAQUETES library(MuMIn) ## para calcular AICc library(car) ## para obtener VIF's y usar Anova, bootCase library(MASS) ## para usar dropterm library(lmtest) ## para obtener la p del modelo, en sustitución de anova(nulo, model, test="Chisq") library(moments) ## para obtener el sesgo, kurtosis y sus significaciones: kurtosis, anscombe.test, skewness, agostino.test library(sandwich) ## para estima robusta, estimando sandwich standard errors library(robustbase) ## para estimas robustas usando Mallows-Huber quasilikelihhod estimation library(robust) ## para estimas robustas library(psych) ## para uso del comando describe library(fit.models) ## para usar el comando leverage library(phia) ## para plots de interacciones library(emmeans) ## para tests a posteriori library(effects) ## para visualizar efectos parciales ## PARA TRABAJAR CON VUESTROS DATOS ... cargar vuestro data frame y después poner su nombre abajo en XXXXXX datos <- XXXXXXX names(datos) ## cear fórmula de nuestro modelo; ejemplo: eqt <- as.formula(abundancia ~ covariante + insolacion*tratamiento) ## establecemos los tipos de contrastes para las variables predictoras nominales (factores) ## 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")) ## crear MODELO POISSON (y su nulo) modelo.p <- glm(eqt, data=datos, family=poisson(link="log")) ## otra opción sería [si no cargamos la línea de options(contrasts=...)] modelo.p <- glm(eqt, data=datos, family=poisson(link="log"), contrasts=c(factor="contr.sum", ordered="contr.poly")) modelo.p.nulo <- update(modelo.p, .~1) ## crear MODELO asumiendo una BINOMIAL NEGATIVA (y su nulo -con diferente theta-) modelo.nb <- glm.nb(eqt, data=datos, link=log) modelo.nb.nulo <- update(modelo.nb, .~1) ## para correr (la mayoría) de las líneas siguientes hacemos (seleccionamos la instrucción y la corremos): ## modelo <- modelo.nb ; modelo.nulo <- modelo.nb.nulo ## seguimos trabajando con el modelo Poisson modelo <- modelo.p ; modelo.nulo <- modelo.p.nulo ## variabilidad explicada por el modelo (usando devianzas) d2 <- round(100*(summary(modelo)$null.deviance-modelo$deviance)/summary(modelo)$null.deviance,2) print(c("D2 de MCFadden (%) =",d2), quote=FALSE) ## representación de los valores observados y predichos en la respuesta ## usando la escala original de medida (habiendo destransformado los datos desde la transformación de la link function) plot(modelo$y~fitted(modelo), ylab="respuesta original") abline(lm(modelo$y~fitted(modelo)), col="red", lwd=2) ## esta R2 es diferente a D2 de MCFadden porque opera en la escala original de la respuesta ## sin aplicarle la transformación asociada a la link function (log en el caso de la Poisson y binomial negativa) print(c("R2 (%) =",round(100*cor(modelo$model[,1],fitted(modelo))^2,2)), quote=FALSE) ## estima de significación global del modelo comparándolo con el modelo nulo ## no hemos creado el modelo nulo, pero glm lo estima automáticamente (modelo$null.deviance) lrtest(modelo) ## likelihood ratio test waldtest(modelo) ## test de Wald; usa test="F"; otra opción es test="Chisq" (la Chisq = F*Df) ## es preferible el LR test; consultad: https://en.wikipedia.org/wiki/Wald_test#Alternatives_to_the_Wald_test ## sólo procederemos valorando los resultados de este modelo, SÍ Y SÓLO SÍ, ## este "omnibus test" ha sido significativo ## valores AICc del modelo de interés y el modelo nulo AICc(modelo, modelo.nulo) veces.mejor <- exp(-0.5*(AICc(modelo)-AICc(modelo.nulo))) print(c("Veces que MI MODELO es mejor que el modelo NULO =", veces.mejor), quote=FALSE) ## exploración de los residuos del modelo ## unos gráficos generales que ya nos dan muchas pistas: ## normalidad de los residuos de devianza (¡¡no en la escala original de la respuesta!!) ## homocedasticidad de los residuos a través de las predicciones del modelo (habiendo aplicado la transformación de la link function) ## existencia de datos influyentes y perdidos con la distancia de Cook y el Leverage ## 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(modelo, c(1:2,4,5)) ## saca los cuatro paneles; en otras ocasiones puede ser interesante: par(mfcol=c(1,1)); par(mfcol=c(1,2)); plot(modelo, 1:2); par(mfcol=c(1,1)) par(mfcol=c(1,1)) ## volvemos al modo gráfico de un solo panel ## o algunos plots individuales hist(residuals(modelo, type="deviance"), main="residuos de devianza") qqnorm(residuals(modelo, type="deviance"), main="residuos de devianza") qqline(residuals(modelo, type="deviance"), main="residuos de devianza", col="red") plot(modelo$linear.predictors, residuals(modelo, type="deviance"), ylab="residuos de devianza", xlab="predicciones del modelo con link function", main="¿hay heterocedasticidad en los residuos?") abline(h=0, col="red") ## test de Shapiro-Wilk de la normalidad de los residuos de devianza shapiro.test(residuals(modelo, type="deviance")) ## puntos influeyentes y perdidos ## leverage con sus límites "críticos" ## si no vemos la línea roja es que no hay problemas con puntos con gran leverage leverage.modelo <- hatvalues(modelo) plot(leverage.modelo) ## niveles críticos 2*(g.l. del modelo)/(Número de datos) abline(h=2*(length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals), col="red") ## dffits con sus límites "críticos" plot(dffits(modelo)) ## niveles críticos 2*raiz((g.l. del modelo)/(Número de datos)) 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") ## distancia de Cook con sus límites "críticos" plot(cooks.distance(modelo)) ## menor que 4/(número de datos) abline(h=4/length(modelo$residuals), col="red") identify(cooks.distance(modelo)) ## terminad dando clik en Finish del panel Plots (esquina sup. derecha) ## plot(modelo, 6) plot(leverage.modelo, cooks.distance(modelo)) abline(v=2*(length(modelo$residuals)-modelo$df.residual-1)/length(modelo$residuals), col="blue") abline(h=4/length(modelo$residuals), col="red") identify(cooks.distance(modelo)) ## terminad dando clik en Finish del panel Plots (esquina sup. derecha) ## INDEPENDENCIA ENTRE LAS PREDICTORAS - 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)) ## Vemos los valores medios de la respuesta (en su escala original) ## ¡¡ojo!! los valores son los ajustados si tuviésemos covariantes ## se dan los errores estándard; ## modificad con legend.margin la posición de las etiquetas y el ancho de la ventana Plots plot(interactionMeans(modelo), legend.margin=0.3) # TESTS POST HOC ## 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("insolacion", "tratamiento")), adjust="tukey") ### ¡¡ ojo !! poned aquí VUESTROS DATOS pairs(emmeans(modelo, c("insolacion", "tratamiento")), adjust="holm") ### ¡¡ ojo !! poned aquí VUESTROS DATOS pairs(emmeans(modelo, c("insolacion", "tratamiento")), adjust="scheffe") ### ¡¡ ojo !! poned aquí VUESTROS DATOS pairs(emmeans(modelo, c("insolacion", "tratamiento")), adjust="fdr") ### ¡¡ ojo !! poned aquí VUESTROS DATOS ## resultados del modelo y SIGNIFICACIONES de efectos parciales (de las predictoras) ## tabla exhaustiva del modelo summary(modelo) ## TEST DE WALD de coeficientes del modelo versus Ho=0 ## https://en.wikipedia.org/wiki/Wald_test ## si sólo queremos la tabla de efectos mirando los coeficientes summary(modelo)$coefficients ## asumiendo Inf grados de libertad coeftest(modelo) ## asumiendo Inf grados de libertad Anova(modelo, type=3, test="Wald") ## da las mismas P's asumiendo Inf grados de libertad coeftest(modelo, df=modelo$df.residual) ## asumiendo los grados de libertad REALES del modelo ## INTERVALOS DE CONFIANZA de los coeficientes asumiendo normalidad asintótica de los mismos: coefci(modelo, level=0.95) ## con Inf grados de libertad coefci(modelo, df=modelo$df.residual, level=0.95) ## con los grados de libertad REALES del modelo ## COEFCIENTES DE REGRESIÓN ESTANDARIZADOS ## zeta-standarizamos (a MEDIA=0 Y SD=1) std.coef(modelo, partial.sd=T) ## INTERPRETACIÓN DE LOS COEFICIENTES ## si hemos utilizado como la función de vínculo el logaritmo (link="log") ## El coeficiente b en antilogaritmo, exp(b), mide el cambio multiplicativo ## en la variable respuesta "Y" cuando esa variable predictora cambia en una unidad. ## O dicho de otro modo, el coeficiente b es el cambio esperado en el log(Y) ## cuando la variable predictora aumenta una unidad. exp(coefficients(modelo)) exp(coefci(modelo, df=modelo$df.residual, level=0.95)) ## LIKELIHOOD RATIO TEST ## de especial relevancia para diseños con factores ## https://en.wikipedia.org/wiki/Likelihood-ratio_test ## dos funciones diferentes con el mismo resultado drop1(modelo, ~., test="Chisq", sorted=FALSE) Anova(modelo, type=3, test="LR") ## ¡¡¡NUNCA utilizaremos anova(...)!!! ## magnitudes de efectos usando las devianzas tabla.devianza <- drop1(modelo, ~., test="LRT")[,c(1:3)] tabla.devianza[1,1] <- sum(tabla.devianza[-1,1]) tabla.devianza$Dev.retenida <- tabla.devianza$Deviance-tabla.devianza[1,2] concomitancias <- sum(tabla.devianza[,4])/modelo$null.deviance tabla.devianza[1,4] <- modelo$null.deviance-modelo$deviance rownames(tabla.devianza)[1] <- "modelo" tabla.devianza$proporcion.explicada <- tabla.devianza$Dev.retenida/modelo$null.deviance ## la diferencia entre la "proporcion.explicada" para el modelo y la suma de esos valores ## para los efectos mide las concomitancias entre los términos predictores round(tabla.devianza, 3) print(c("Concomitancias (proporción) =",round(tabla.devianza[1,5]-concomitancias, 3)), quote=FALSE) ## VISUALIZACION DE EFECTOS PARCIALES, incluyendo covariante(s) ## 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). ## ## con el paquete effects ## proporciona intervalos de confianza al 95% plot(allEffects(modelo)) ## VALIDACIÓN CRUZADA DE MODELOS GLM Poisson ## elegimos un modelo glm(...) modelo <- modelo.p ## ## 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 <- glm(respuesta~.-id, family=modelo$family, 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] ## estima de la SOBREDISPERSION del modelo ## https://data.princeton.edu/wws509/notes/c4a.pdf ## lo hacemos SOLO PARA la Poisson y Binomial ## nunca lo haremos para la Binomial-Negativa, Gamma, y Gaussiana ## este valor canónico debería de ser igual a la unidad phi <- sum((residuals(modelo, type="pearson"))^2)/modelo$df.residual print(c("Pearson overdispersion =", round(phi, 5)), quote=FALSE) ## si phi está muy desviado de 1 recalculamos las estimas de significación ## teniendo en cuenta la sobredispersión. ## no se recomienda corregir la sobredispersión si phi < 1 ## F = (diferencias en Devianza) / (d.f. x phi) ## directamente con: ## estima asintótica (con la distribución Z que asume Inf df) summary(modelo, dispersion=phi) ## estima basada en el test de la F que incluye la sobredispersión Anova(modelo, type=3, test="F") ## la F previa es igual que la división de LR Chisq / (Df * phi) de la siguiente salida Anova(modelo, type=3, test="LR") ## otra manera basada en las familias QUASI ## da lo mismo que el modelo basado en quasipoisson; e.g.: family=quasipoisson(link="log") modelo.sobredisp <- update(modelo.p, family=quasipoisson(link="log")) summary(modelo.sobredisp) ## estima con la distribución de la t de Student que asume los df reales coeftest(modelo.sobredisp, df=modelo$df.residual) ## asumiendo los grados de libertad REALES del modelo ## 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 ## https://www.researchgate.net/publication/225715335_A_New_Heteroskedasticity-Consistent_Covariance_Matrix_Estimator_for_the_Linear_Regression_Model ## nueva significación de coeficientes (no cambian los coeficientes; solo los error stándard y las P) coeftest(modelo, df=modelo$df.residual) ## sin corrección coeftest(modelo, df=modelo$df.residual, vcovHC(modelo, type="const")) ## asumiendo ausencia de heterocedasticidad coeftest(modelo, df=modelo$df.residual, vcovHC(modelo, type="HC5")) ## con corrección HC5 de heterocedasticidad ## no funciona Anova(modelo, type=3, test="Wald", white.adjust="hc4") ## intervalos de confianza coefci(modelo, df=modelo$df.residual, level=0.95) ## sin corrección coefci(modelo, df=modelo$df.residual, level=0.95, vcov=vcovHC(modelo, type="HC5")) ## RECÁLCULO DEL MODELO QUITANDO ALGUNOS PUNTOS modelo.sin_outliers <- update(modelo, .~., data=datos[c(-23, -28, -46, -54, -55, -56),]) summary(modelo.sin_outliers) ## ESTIMACIONES ROBUSTAS (no funciona con modelos glm.nb para Binomiales Negativas) ## volvemos al modelo Poisson modelo <- modelo.p ## ESTIMAS ROBUSTAS CON MODELOS GENERALIZADOS LINEALES ## (usando IWLS: iteratively (re)weighted least squares) ## sin quitar datos del modelo, aproximación robusta que estima nuevos coeficientes y errores standard ## Mqle es el método Mallows-Hubber quasi-liquelihood. ## además podemos trabajar con weights.on.x en función del "hat" (i.e., los potenciales "outliers" se penalizan de acuerdo con el leverage) ## los métodos method="WBY" y method="BY" están disponibles solo para regresión logistic (family=binomial) ## tcc es la constante de ajuste c para Huber's psi-function: ## a mayor valor de tcc menos peso a valores extremos considerados ## https://cran.r-project.org/web/packages/robustbase/vignettes/psi_functions.pdf ## tenemos en cuenta que nuestro modelo utiliza una ecuación-función eqt y un dataframe llamado datos ## en familia definimos la distribución de la respuesta y la función de vínculo: ## family=gaussian(link ="identity") ## family=poisson(link="log") ## family=binomial(link="logit") ó family=binomial(link ="cloglog") ## family=Gamma(link="log") ## hacemos el modelo robusto sin infrapesar por el leverage robusto <- glmrob(eqt, data=datos, family=poisson(link="log"), method="Mqle", control=glmrobMqle.control(tcc=1.2, maxit=100)) ## hacemos el modelo robusto infrapesabdo por el leverage, de manera inversamente proporcional a su valor robusto.hat <- glmrob(eqt, data=datos, family=poisson(link="log"), weights.on.x="hat", method="Mqle", control=glmrobMqle.control(tcc=1.2, maxit=100)) ## seleccionamos el modelo robusto modelo.rob <- robusto.hat ## valoramos los RESULTADOS exhaustivos summary(modelo.rob) ## tabla de Anova para coeficientes, con estima asintótica de significación (asumiendo Inf grados de libertad con la distribución Z) round(summary(modelo)$coefficients, 5) ## modelo original de interés round(summary(modelo.rob)$coefficients, 5) ## tabla de Anova para coeficientes, con estima de significación con los grados de libertad correctos del modelo de interés round(coeftest(modelo, df=modelo$df.residual), 5) ## modelo original de interés round(coeftest(modelo.rob, df=modelo$df.residual), 5) ## modelo robusto ## tabla de Anova para efectos (útil con factores), con estima asintótica de significación Anova(modelo, type=3, test="Wald") ## modelo original de interés Anova(modelo.rob, type=3, test="Wald") ## aproximación a la proporción de la variabilidad explicada (pseudo-R2) devianza.robusto <- sum(residuals(modelo.rob, type="deviance")^2) devianza.robusto.nula <- sum(residuals(update(modelo.rob, .~1), type="deviance")^2) d2.robusto <- (devianza.robusto.nula - devianza.robusto) / devianza.robusto.nula print(c("D2 de McFadden (%) modelo robusto=", round(100*d2.robusto, 2)), quote=FALSE) ## robustez de los datos individuales plot(modelo.rob$w.r, modelo.rob$w.x) robustez <- modelo.rob$w.r*modelo.rob$w.x plot(robustez, ylab="robustez de las observaciones") identify(modelo.rob$w.r) ## finalizar dando clik a "Finish" (esquina superior derecha de Plots) ## ## relación de los pesos del modelo robusto con la distancia de Cook del modelo original plot(robustez ~ cooks.distance(modelo), ylab="robustez de las observaciones") ## relación de los pesos del modelo robusto con los valores absolutos de dffit del modelo original plot(robustez ~ abs(dffits(modelo)), ylab="robustez de las observaciones") ## relación de los pesos del modelo robusto con los valores de leverage del modelo original plot(robustez ~ hatvalues(modelo), ylab="robustez de las observaciones") ## aproximación package "robust" (método conditionally unbiased bounded influence estimator) modelo.rob2 <- glmRob(modelo$formula, data=datos, family=poisson(), method="cubif") ## valoramos los RESULTADOS exhaustivos summary.glmRob(modelo.rob2, correlation=FALSE) ## omnibus test df.rob2 <- summary.glmRob(modelo.rob2, correlation=FALSE)$df[1] chi2.rob2 <- summary.glmRob(modelo.rob2, correlation=FALSE)$null.deviance - summary.glmRob(modelo.rob2, correlation=FALSE)$deviance P.rob2 <- pchisq(chi2.rob2, df=df.rob2, lower.tail=FALSE) print(paste("Chi2 =", round(chi2.rob2,2), " df =", df.rob2, " P =", round(P.rob2, 5)), quote=F) ## aproximación a la proporción de la variabilidad explicada (pseudo-R2) print(c("D2 de MCFaddeen (%) =",round(100*(modelo.rob2$null.deviance-modelo.rob2$deviance)/modelo.rob2$null.deviance,2)), quote=FALSE) ## algunos plots individuales hist(residuals(modelo.rob2, type="deviance"), main="residuos de devianza") qqnorm(residuals(modelo.rob2, type="deviance"), main="residuos de devianza") qqline(residuals(modelo.rob2, type="deviance"), main="residuos de devianza", col="red") plot(modelo.rob2$linear.predictors, residuals(modelo.rob2, type="deviance"), ylab="residuos de devianza", xlab="predicciones del modelo con link function", main="¿hay heterocedasticidad en los residuos?") abline(h=0, col="red") ## plot(modelo.rob2$weights, ylab="robustez de las observaciones") ## robustez de los datos individuales identify(modelo.rob2$weights) ## terminad dando clik en Finish del panel Plots (esquina sup. derecha) plot(modelo.rob2$weights, cooks.distance(modelo), xlab="robustez de las observaciones") identify(cooks.distance(modelo)~modelo.rob2$weights) ## terminad dando clik en Finish del panel Plots (esquina sup. derecha) plot(modelo.rob2$weights, hatvalues(modelo), xlab="robustez de las observaciones") ## frente a leverage ## PRIORIZACION DE MODELOS SEGUN CRITERIOS DE AIC AICc(modelo.p, modelo.nb) ## Multimodel Inference options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## volvemos a escribir la ecuación de nuestro modelo dentro de glm(...) ## para POISSON ### ¡¡¡ PONED AQUI VUESTROS DATOS !!! modelo <- glm(abundancia ~ covariante + insolacion*tratamiento, data=datos, family=poisson(link="log")) ## para BINOMIAL NEGATIVO ### ¡¡¡ PONED AQUI VUESTROS DATOS !!! modelo <- glm.nb(abundancia ~ covariante + insolacion*tratamiento, data=datos, link=log) ## options(na.action = "na.fail") dredge.modelo <- dredge(modelo, beta="sd", rank="AICc") ## 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)