###################################### ## MODELOS GENERALIZADOS ADITIVOS ## ## (con control de autocorrelación) ## ###################################### ## Luis M. Carrascal ## https://lmcarrascal.eu ## revisión 20/02/2022 ## PAQUETES DE ANALISIS library(car) library(mgcv) library(MuMIn) library(nlme) library(lmtest) library(MASS) library(lattice) ## 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 ## de internet y ponderlo en uso en "environment" meconecto <- url("http://www.lmcarrascal.eu/cursos/regravanz/datosgamm.RData") load(meconecto) close(meconecto); rm(meconecto) ## convertimos los niveles del factor aleatorio en una variable que sea FACTOR datos$LOCATION <- as.factor(datos$LOCATION) ## cuidado con variables "character" o "factor" class(datos$factor) datos$factor <- as.factor(datos$factor) ## ¡IMPORTANTE CORRED ESTA LINEA! en el caso de que existan factores ## establecemos la tabla de contrastes ortogonales para los factores options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## FORMULACION DE LOS MODELOS gam Y gamm ## Podemos trabajar con diferentes tipos de distribuciones de la variable respuesta ## mirad: ?family.mgcv ## family=gaussian(link="identity"); family=Gamma(link="identity"); family=poisson(link="log"); binomial(link="logit"); ## family=quasipoisson(link="log"); family=quasibinomial(link="logit") ## para binomiales negativas con estima de theta: family=nb(link="log") con gam, no con gamm ## gamm produce resultados "pobres" con distribuciones binomiales ## mirad WARNINGS al final de ?gamm ## 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 222, 224) ## thin plate splines (http://en.wikipedia.org/wiki/Thin_plate_spline); k define la complejidad máxima alcanzable ## consultad (tecleando en la consola): help(smooth.terms) ## 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) ## s(x,z) para "interaction splines" cuando las dos predictoras se miden en la misma escala; se asume un mismo valor de k para x y z ## ejemplos: ti(variable_1, variable_2, bs="tp", k=15) podemos modificar bs="tp" por bs="cr" ## te(x,z) para "tension product spline" cuando las dos predictoras no se miden en la misma escala; ## se puede estimar un diferente valor de k para x y z ## ejemplos: te(variable_1, variable_2, bs="tp", k=15) ## te(x,z) y s(x,z) controlan por los efectos principales x y z ## ti(x,z) solo estima el efecto interacción, pero sin controlar por los efectos principales s(x) y s(z) ## para corregir la HETEROCEDASTICIDAD EN LOS RESIDUOS podemos usar como pesos (weights) ## para covariantes: varFixed, varExp o varPower (si las predictoras no incluyen el valor cero) ## para factores: varIdent ## con varComb podemos combinar estructuras de corrección considerando factores y covariantes ## corrección de la AUTOCORRELACIÓN TEMPORAL en los residuos ## tenemos la clase corCAR1 aplicado a una fecha con datos irregularmente espaciados ## (equivalente a ARMA p=1, q=0; si los datos están regularmente espaciados) ## o la más compleja ARMA de AutoRegressive and MovingAverage, con los parámetros p y q respectivamente ## p y q hacen referencia al número de unidades temporales previas uniformemente espaciadas ## help(corClasses) ## fórmula del modelo ## humedad ~ s(temperatura, k=4) + s(precip_24h, k=4) + factor ## usando LOCATION como efecto aleatorio): random=list(LOCATION=~1) ## otra manera de definir términos "Random_Effect" ("re") es con s(factor_aletorio, bs="re") ## con autocorrelación temporal por separado para los niveles del factor aleatorio: correlation=corCAR1(form=~fecha|factor) ## (usando la fecha y asumiendo efectos diferentes dentro de cada LOCATION) ## controlando la heterocedasticidad en los residuos: weights=varFixed(~temperatura+precip_24h), weights=varExp(), ... ## construimos los modelos con Maximum Likelihood (ML) para compatrarlos con Akaike AICc ## para modelos Gaussianos podemos utilizar REML; para los demás generalizados (poisson, binomial, ...) usaremos siempre ML ## EMPLEANDO LA FUNCION gam QUE NO ACEPTA ESTRUCTURAS DE CONTROL DE HETEROCEDASTICIDAD Y AUTOCORRELACION RESIDUAL ## haciendo uso de un spline complejo para controlar el efecto del tiempo (i.e., fecha, autocorrelación temporal) ## en vez de con estructuras correlation=corCAR1(...): s(fecha, k=12) ## modelo RANDOM INTERCEPT & FIXED SLOPES con una gausiana usando gam gam.1 <- gam(humedad ~ s(fecha, k=12) + s(temperatura, k=4) + s(precip_24h, k=4) + factor + s(LOCATION, bs="re"), data=datos, family=gaussian(link="identity"), method="ML") ## modelo RANDOM INTERCEPT & FIXED SLOPES con una gausiana usando gam ## con temperatura como efecto lineal y no usando thin plate smooth gam.11 <- gam(humedad ~ s(fecha, k=12) + temperatura + s(precip_24h, k=4) + factor + s(LOCATION, bs="re"), data=datos, family=gaussian(link="identity"), method="ML") ## modelo RANDOM INTERCEPT para todos los predictores fijos con una gausiana usando gam ## añadiendo UNCORRELATED RANDOM INTERCEPT & RANDOM SLOPE PARA LA TEMPERATURA como efecto lineal (no usando thin plate smooth) gam.21 <- gam(humedad ~ s(fecha, k=12) + temperatura + s(precip_24h, k=4) + factor + s(LOCATION, bs="re") + s(temperatura, LOCATION, bs="re"), data=datos, family=gaussian(link="identity"), method="ML") ## modelo RANDOM INTERCEPT para todos los predictores fijos con una gausiana usando gam ## añadiendo un RANDOM SMOOTH PARA LA TEMPERATURA ## asume un parámetro de suavizado común para la temperatura, aunque puede tener diferente forma en los niveles de random gam.3 <- gam(humedad ~ s(fecha, k=12) + s(precip_24h, k=4) + factor + s(LOCATION, bs="re") + s(temperatura, LOCATION, bs="fs"), data=datos, family=gaussian(link="identity"), method="ML") ## modelo RANDOM SMOOTH PARA LA TEMPERATURA con una gausiana usando gam ## asume un parámetro de suavizado común para la temperatura, aunque puede tener diferente forma en los niveles de random gam.31 <- gam(humedad ~ s(fecha, k=12) + s(precip_24h, k=4) + factor + s(temperatura, LOCATION, bs="fs"), data=datos, family=gaussian(link="identity"), method="ML") ## modelo RANDOM SMOOTH AND POPULATION-LEVEL EFFECT PARA LA TEMPERATURA con una gausiana usando gam ## asume un parámetro de suavizado común, aunque puede tener diferente forma en los niveles de random ## y añade un suavizado global (population level curve s(X)) al Random smooth gam.41 <- gam(humedad ~ s(fecha, k=12) + s(temperatura, k=4) + s(precip_24h, k=4) + factor + s(temperatura, LOCATION, bs="fs"), data=datos, family=gaussian(link="identity"), method="ML") ## EMPLEANDO LA FUNCION gamm QUE ACEPTA ESTRUCTURAS DE CONTROL DE HETEROCEDASTICIDAD Y AUTOCORRELACION RESIDUAL ## modelo RANDOM INTERCEPT & FIXED SLOPES con una gausiana usando gamm ## controlando la heterocedasticidad residual (weights=...) y la autocorrelación temporal (correlation=...) gamm.1 <- gamm(humedad ~ s(temperatura, k=4) + s(precip_24h, k=4) + factor + s(LOCATION, bs="re"), data=datos, family=gaussian(link="identity"), method="ML", correlation=corCAR1(form=~fecha|LOCATION), weights=varFixed(~temperatura+precip_24h)) ## los otros modelos serían similares a los presentados para gam.11, gam.21, gam.3, gam.31, gam.41 ## controlando la heterocedasticidad residual y la autocorrelación temporal ## como consecuencia de ellos son ¡¡muy lentos!! ## COMPARACION DE MODELOS AICc(gam.1, gam.11, gam.21, gam.3, gam.31, gam.41) ## elegimos el modelo con menor valor de AICc modelo <- gam.1 ## construimos el MODELO NULO (SIN EFECTOS FIJOS) para el mejor modelo mixto modelo.nulo <- update(modelo, .~ 1 + s(LOCATION, bs="re")) ## OMNIBUS TEST ## comparamos los dos modelos ML que han utilizado la estructura aleatoria: random=list(LOCATION=~1) AICc(modelo.nulo, modelo) ## comparamos los dos modelos anidados ML mediante lrtest lrtest(modelo.nulo, modelo) ## VOLVEMOS A GENERAR EL MODELO CON LA ESTIMA REML (solamente para los modelos Gaussianos) modelo <- update(modelo, method="REML") ## no es posible con update(...) para los modelos gamm ## hay que volverlos a escribir desde el principio añadiendo method="REML" ## PARA LOS MODELOS OBTENIDOS HACIENDO USO DE gamm ## la parte condicional asociada a los efectos fijos está en modelo.gamm$gam ## si tenemos un modelo gamm llamado "modelo" hay que sustituir "modelo" de las líneas siguientes por "modelo$gam" ## PRUEBAS CANÓNICAS DEL MODELO ## NORMALIDAD de los residuos hist(residuals(modelo, type="deviance")) ## HETEROCEDASTICIDAD de los residuos plot(modelo$linear.predictors, residuals(modelo, type="deviance")) abline(h=0, col="red", lwd=2) ## de manera global y valorando (en la salida de la consola) la NECESIDAD DE ESTABLECER UN MAYOR VALOR DE K para los efectos curvos s(...) ## el modelo ha sido construido intencionadamente suavizado con bajos valores de K par(mfcol=c(1,1)) par(mfcol=c(2,2)) gam.check(modelo) par(mfcol=c(1,1)) ## concurvidad (equivalente a colinealidad en los modelos lineales) ## "para" se refiere a los términos paramétricos no curvilíneos concurvity(modelo) ## RESULTADOS Y SIGNIFICACIÓN DE EFECTOS ## significación del tipo SS-III de la parte fija summary(modelo) ## tabla de ANOVA (ideal para el caso de factores, con tres niveles o más) anova.gam(modelo) ## VARIACION EXPLICADA ## pseudo-R2 ajustado por los grados de libertad de los efectos summary(modelo)$r.sq ## pseudo-R2 explicado sin ajustar cor(fitted(modelo), modelo$y)^2 plot(fitted(modelo), modelo$y) abline(a=0, b=1, col="red", lwd=2) ## VISUALIZACION DE RESULTADOS ## la etiqueta para el eje Y son los RESIDUOS PARCIALES DE LA RESPUESTA para esa predictora, controlando por las demás en el modelo plot(modelo, shade=T, shade.col="gray70", all.terms=T, scheme=2, pages=1, se=2) ## sin predictoras no sometidas a thin plate spline plot(modelo, shade=T, shade.col="gray70", all.terms=F, scheme=2, pages=1, se=2) ## mostrando los valores de las unidades muestrales (en azul) plot(modelo, shade=T, shade.col="gray70", all.terms=F, residuals=TRUE, pch=19, se=2, cex=0.2, col="blue", pages=1, rug=TRUE) ## representación de un solo plot (seleccionándolo de la secuencia de paneles gráficos) plot(modelo, shade=T, shade.col="gray70", se=2, scheme=2, rug=T, too.far=0.05, select=1) ## Otro ejemplo de EFECTOS ALEATORIOS sin incluir un RANDOM FACTOR ## pero con control de la AUTOCORRELACION ESPACIAL que definiría la parte aleatoria (psudo-replicación por autocorrelación) ## para una distribución no gausiana ## IMPORTAMOS DATOS DE TRABAJO meconecto <- url("http://www.lmcarrascal.eu/cursos/regravanz/gamm_pinnas.RData") load(meconecto) close(meconecto); rm(meconecto) ## levels(datos$lugar) levels(datos$tipo) ## ¡IMPORTANTE CORRED ESTA LINEA! en el caso de que existan factores options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## MODELO gamm CON CONTROL DE LA AUTOCORRELACION ESPACIAL ## gamm no puede abordar distribuciones binomiales negativas; lo corremos asumiendo una distribución (quasi)poisson ## sin existir un factor random, al menos hay que definir un término s(...) ## intencionadamente se definen splines poco complejas ## ¡¡¡ muy lento !!! mod.gamm <- gamm(n_pinnas ~ lugar + tipo + s(altitud, k=5) + s(alt_pinos, k=5) + s(cob_pinos, k=5) + s(radiac_jul_ago, k=5), data=datos, family=quasipoisson(link="log"), niterPQL=200, correlation=corExp(form=~longitud+latitud)) ## MODELO gam MODELIZANDO LA AUTOCORRELACION ESPACIAL ## en gam sí podemos usar distribuciones binomiales negativas ## utilizamos un plano de interacción de longitud x latitud: s(longitud, latitud, k=10) ## ¡¡¡ muy rápido !!! mod.gam <- gam(n_pinnas ~ s(longitud, latitud, k=10) + lugar + tipo + s(altitud, k=5) + s(alt_pinos, k=5) + s(cob_pinos, k=5) + s(radiac_jul_ago, k=5), data=datos, family=nb(link="log"), method="ML") ## los comparamos con AICc AICc(mod.gamm, mod.gam) ## bondad del modelo según su ajuste a supuestos canónicos par(mfcol=c(1,1)) par(mfcol=c(2,2)) gam.check(mod.gamm$gam) par(mfcol=c(1,1)) ## par(mfcol=c(1,1)) par(mfcol=c(2,2)) gam.check(mod.gam) par(mfcol=c(1,1)) ## concurvidad concurvity(mod.gamm$gam) concurvity(mod.gam) ## variación (devianza) explicada summary(mod.gamm$gam)$r.sq summary(mod.gam)$r.sq ## tabla de coeficientes paramétricos (i.e., términos fijos no-spline) round(summary(mod.gamm$gam)$p.table, 3) round(summary(mod.gam)$p.table, 3) ## tabla de coeficientes para los términos spline round(summary(mod.gamm$gam)$s.table, 3) round(summary(mod.gam)$s.table, 3) ## tabla de ANOVA anova.gam(mod.gamm$gam) anova.gam(mod.gam) ## ilustración de efectos plot(mod.gamm$gam, shade=T, shade.col="gray70", all.terms=T, scheme=2, pages=1, se=2) plot(mod.gam, shade=T, shade.col="gray70", all.terms=T, scheme=2, pages=1, se=2) ## plot(mod.gamm$gam, shade=T, shade.col="gray70", se=2, scheme=2, rug=T, too.far=0.05, select=1) plot(mod.gam, shade=T, shade.col="gray70", se=2, scheme=2, rug=T, too.far=0.05, select=2)