########################################## ## Modelos Generalizados Aditivos ## ## con diferentes distribuciones ## ## modelos GAM con mgcv ## ########################################## ## Luis M. Carrascal ## https://lmcarrascal.eu ## revisión 02/03/2020 ## cargamos las siguientes librerías library(nlme) ## es necesario para mgcv library(mgcv) ## este es el paquete para efectuar modelo GAM library(MuMIn) library(lmtest) library(sandwich) ## para correcciones por heterocedasticidad en los residuos del modelo ## ## ## ## IMPORTACIÓN DE DATOS ## para importar / exportar los datos fuente ## datos <- read.table("clipboard", header=T, sep="\t", dec=".") ## importar desde portapapeles ## write.table(datos, "kk.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/regravanz/datosregravanz.RData") ## abro una conexión con internet load(meconecto) ## cargo el archivo de la conexión close(meconecto); rm(meconecto) ## cierro la conexión ## names(datos) ## ## ## ## Podemos trabajar con diferentes tipos de distribuciones de la variable respuesta ## family= poisson(link="log"); binomiales negativas: nb(); gaussian(link="identity"); binomial(link="logit"); quasipoisson(link="log"); quasibinomial(link="logit") ## ## Los términos no lineales se pueden definir de diferentes modos ## consultad el capítulo 5 de http://reseau-mexico.fr/sites/reseau-mexico.fr/files/igam.pdf (paginas 224, 226) ## thin plate splines (http://en.wikipedia.org/wiki/Thin_plate_spline); k define la complejidad máxima alcanzable ## por defecto, al utilizar s(...) la base del spline (bs) es thin plate (tp) ## ejemplos: s(altmed, bs="tp", k=10) s(altmed, bs="tp", k=5); k por defecto es 10 (si no se escribe) ## cubic regression splines, bs="cr": ## ejemplos: s(altmed, bs="cr", k=10) s(altmed, bs="cr", k=5); k por defecto es 10 (si no se escribe) ## ti(...) para "interaction splines" cuando las dos predictoras se miden en la misma escala; suelen requerir un valor de k más alto ## ejemplos: ti(variable_1, variable_2, bs="tp", k=15) podemos modificar bs="tp" por bs="cr" ## te(...) para "tension product spline" cuando las dos predictoras no se miden en la misma escala ¡¡¡ mayoritariamente utilizaremos éste !!! ## ejemplos: te(variable_1, variable_2, bs="tp", k=15) ## ## ## Creamos la fórmula eqt <- as.formula(turmer~te(longitud, latitud, bs="tp", k=15) + s(altmed, bs="tp", k=10) + s(rangoalt, bs="tp", k=10) + s(shannon, bs="tp", k=10) + s(tempmin, bs="tp", k=10) + s(precip, bs="tp", k=10) + sizepa + distmar) ## ## ¡IMPORTANTE CORRER ESTA LINEA! en el caso de que existan factores options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## ## ## Creamos el modelo ## ## Método GCV.Cp: método "general cross validation for unknown scale parameter" ## añadiendo scale=-1 corrige por sobredispersión en variables respuestas BINOMIALES o POISSON ¡¡¡ hacedlo sólo con ellas !!! ## entonces no es necesario crear un modelo quasbinomial o quasipoisson ## el argumento gamma: "The argument gamma=1.4, forces each model effective degree of freedom to count as ## 1.4 degrees of freedom in the GCV score, which forces models to be a little smoother than they might otherwise be, ## and is an ad hoc way of avoiding overfitting" modelo <- gam(eqt, data=datos, family=poisson(link="log"), method="GCV.Cp", gamma=1.4, scale=0) ## ## modelos ML (Maximum Likelihood) o REML ("REML estimation, including the unknown scale") ## por defecto REML en BINOMIALES NEGATIVAS (family=nb(link ="log")) y en modelos mixtos modelo.ml <- gam(eqt, data=datos, family=poisson(link ="log"), method="ML", optimizer=c("outer","newton")) modelo.reml <- gam(eqt, data=datos, family=poisson(link ="log"), method="REML", optimizer=c("outer","newton")) ## ## ## ## seleccionar el modelo que queramos modelo <- modelo ## ## ## ## SUPUESTOS CANÓNICOS (muy importante en distribuciones poisson y binomiales) ## ## estimamos el coeficiente de sobredispersión texto <- paste("coef dispersión modelo", modelo$call$method, "=") print(c(texto, sum(residuals(modelo, type="pearson")^2)/modelo$df.residual), quote=FALSE) ## corregir por sobredispersión alterará la significación y la complejidad de los términos spline: o con family=quasi... o con el argumento scale=-1 en el modelo ## ## primeros resultados gráficos (supuestos canónicos con los residuos -3 paneles- y relación observado vs. predicho) ## también una salida de texto con indicación de si necesitamos elevar el parámetro de complejidad k de los splines ## k' es el máximo k obtenible menos uno (k'= k-1); edf son los grados de libertad efectivos; en te o ti k' = k*k-1 ## lo ideal es que edf sea menor que k'; si no es así ... volvemos a construir la ecuación eqt elevando k (sobre todo si p<0.01) par(mfcol=c(1,1)) gam.check(modelo, k.rep=1000) ## ## más detalle del análisis de la normalidad de los residuos par(mfcol=c(1,1)) qqnorm(residuals(modelo,type="deviance")) qqline(residuals(modelo,type="deviance"), col="red") shapiro.test(residuals(modelo,type="deviance")) ## ## más detalle para la homocedasticidad de los residuos plot(modelo$linear.predictors, residuals(modelo,type="deviance"), main="¿Hay heterocedasticidad?") abline(h=0, col="red") ## ## concurvidad (equivalente a colinealidad en los modelos lineales) concurvity(modelo) ## para se refiere a los términos paramétricos no curvilíneos ## ## ## ## Estima de la significación global del modelo; comparación de nuestro modelo con el nulo (respuesta ~ 1) lrtest(modelo) ## ## ## ## El resumen clásico de los resultados ## proporciona la variación observada en la variable respuesta (dev expl %) y la estima de R2 ajustada por la complejidad del modelo (R-sq.(adj)) ## el parámetro UBRE es una estima de parsimonia equivalente a AIC pero cuando trabajamos con el método "GCV.Cp" ## mejor cunato menor sea su valor al comparar varios modelos "anidados" (también lo haremos con AICc) summary(modelo) ## Resumen con equivalente a SS-III. ## Ideal para examinar efectos globales de factores, en vez de a través de sus columnas de contraste. anova.gam(modelo) ## ## R2 sin ajustar (el ajustado lo proporciona summary(...) en R-sq.(adj) print(c("R2 (%) =", round(100*cor(modelo$fitted, modelo$y)^2,2)), quote=FALSE) ## ## de nuevo, cómo sale la relación final entre observado y predicho en la escala original de medida de la variable respuesta par(mfcol=c(1,1)) plot(modelo$fitted, modelo$y, main="OBSERVADO vs PREDICHO") abline(lm(modelo$y~modelo$fitted), col="red", lwd=2) ## ## ## ## REPRESENTACIÓN GRÁFICA DE EFECTOS plot(modelo, scale=0, pages=1, shade=T, shade.col="gray70", se=2, all.terms=T) plot(modelo, pages=1, residuals=TRUE, shade=T, shade.col="gray70", se=2, all.terms=TRUE) plot(modelo, pages=1, shade=T, shade.col="gray70", se=2, all.terms=T, scheme=2) plot(modelo, scale=0, pages=1, shade=T, shade.col="gray70", se=2, all.terms=T, scheme=1) plot(modelo, shade=T, shade.col="gray70", all.terms=F, residuals=T, pch=19, cex=0.1, scheme=1, col="red", pages=1, rug=F) ## otro modo de verlo con los datos individuales plot(modelo, shade=T, shade.col="gray70", se=2, scheme=2, rug=T, too.far=0.05, select=5) plot(modelo, shade=T, shade.col="gray70", se=2, scheme=2, residuals=T, pch=19, col="red", rug=F, too.far=0.05, select=5) plot(modelo, shade=T, shade.col="gray70", se=2, scheme=2, residuals=T, pch=19, col="red", cex=0.1, rug=F, too.far=0.05, select=5) plot(modelo, shade=T, shade.col="gray70", se=2, all.terms=T, scheme=2, select=1) ## ## Representaciones de posibles interacciones ## seleccionar variables para vis.gam y escribidlas en view (preservad las comillas) ## modificamos en Zlim los límites de representación de la variación de la variable respuesta vis.gam(modelo, view=c("altmed", "shannon"), type = "response", zlim=c(0,80), plot.type = "contour", n.grid=100) ## lo mismo de lo anterior pero tridimensionalmente; con phi y theta modificamos los ángulos de giro para una mejor visualización vis.gam(modelo, view=c("altmed", "shannon"), type = "response", zlim=c(0,80), plot.type = "persp", phi=30, theta=30, n.grid=100, border=NA) par(mfcol=c(1,1)) ## ## ## ## generación de un nuevo modelo revisado a partir de otro previo ## quitamos (con - delante) o añadimos (con + delante) términos modelo.modificado <- update(modelo, .~. -s(shannon, bs="tp", k=10) +s(shannon, bs="tp", k=5) -distmar +s(distmar, bs="tp", k=5) -sizepa) summary(modelo.modificado) anova.gam(modelo.modificado) ## COMPARACIÓN DE MODELOS AICc(modelo, modelo.modificado) ## ## Muchas de estas cosas ya las hemos visto con otros tipos de modelos