#################################### ## GENERALIZED LEAST SQUARES ## ## MODELOS LINEALES CON CONTROL ## ## DE LA VARIACION RESIDUAL ## #################################### ## Luis M. Carrascal ## https://lmcarrascal.eu ## revisión 03/03/2022 ## referencias relevantes: ## https://en.wikipedia.org/wiki/Generalized_least_squares ## http://www.biostat.jhsph.edu/~iruczins/teaching/jf/ch5.pdf ## https://socialsciences.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Timeseries-Regression.pdf ## commented examples of gls use: https://goo.gl/ETBda4 library(nlme) ## para gls y Variogram library(MuMIn) ## para calcular AICc library(car) ## para obtener VIF library(lmtest) ## para obtener la p del modelo, y comparar varios modelos library(effects) ## para representación de efectos parciales library(corrgram) ## Operaremos con modelos lineales y una variable respuesta gausiana ## 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 ###################################### ## CONTROL DE LA HETEROCEDASTICIDAD ## ###################################### ## https://fukamilab.github.io/BIO202/03-C-heterogeneity.html ## ## IMPORTAR EXPORTAR DATOS RData ## 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) ## CONSTRUCCIÓN DE LA ECUACIÓN DEL MODELO GLS ## fórmula ("sizepa" es un factor) eqt <- as.formula(nspp ~ altmed+rangoalt+shannon+tempmin+precip+sizepa) ## si tenemos un factor, correremos la siguiente línea para establecer unos contrastes correctos options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## nos fijamos en la gran heterocedasticidad que existe en las relaciones respuesta - predictores ## representación gráfica de las relaciones entre variables pairs(eqt, data=datos) ## los pie charts denotan la proporción de varianza compartida (relaciones: roja-negativa, verde-positiva) ## los numeritos en la diagonal denotan los valores mínimos y máximos de cada variable corrgram(lm(eqt, data=datos)$model, order=FALSE, upper.panel=panel.pts, lower.panel=panel.pie, text.panel=panel.txt, diag.panel=panel.minmax, cex.labels = 2, col.regions = colorRampPalette(c("red", "salmon", "white", "mediumseagreen", "darkgreen")), main="relación entre variables") ## EL MODELO ## para poder comparar modelos usando AICc usaremos method="ML" ## también se pueden comparar modelos con method="REML" si los efectos fijos no cambian ## para la estima del modelo final utilizaremos method="REML" más precisa modelo.nulo <- gls(nspp ~ 1, data=datos, method="ML") ## el siguiente modelo converge con lm(eqt, data=datos) modelo <- gls(eqt, data=datos, method="ML") modelo.lm <- lm(eqt, data=datos) ## valoramos la heterocedasticidad de los residuos plot(modelo) plot(modelo.lm, 1) ## heterocedasticidad en relación con predictores ## poned en la X los diferentes predictores del modelo datos$... eqt ## comparad las asociaciones diferentes de altmed y rangoalt respecto a shannon, tempmin, sizepa plot(datos$altmed, residuals(modelo)) plot(datos$tempmin, residuals(modelo)) ## para corregir la heterocedasticidad en los residuos podemos usar como pesos (weights) ## varFixed, varExp o varPower para covariantes (si las predictoras no incluyen el valor cero) ## varIdent para factores modelo.hdc.exp <- gls(eqt, data=datos, method="ML", weights=varExp()) modelo.hdc.pow <- gls(eqt, data=datos, method="ML", weights=varPower()) modelo.hdc.fix <- gls(eqt, data=datos, method="ML", weights=varFixed(~altmed+rangoalt)) modelo.hdc.fac <- gls(eqt, data=datos, method="ML", weights=varIdent(form = ~1|sizepa)) ## para corregir la heterocedasticidad a través de los niveles de un factor varIdent ## y combinándolo con otra correción varComb(...) modelo.hdc.expfac <- gls(eqt, data=datos, method="ML", weights=varComb(varIdent(form = ~1|sizepa), varExp())) ## comparación de modelos habiendo utilizado la estima ML (Maximum_Likelihood) AICc(modelo.nulo, modelo.lm, modelo, modelo.hdc.fix, modelo.hdc.exp, modelo.hdc.pow, modelo.hdc.fac, modelo.hdc.expfac) ## comparación de dos modelos entre si anova(modelo.hdc.exp, modelo) anova(modelo.hdc.fix, modelo.hdc.exp) anova(modelo.hdc.exp, modelo.hdc.expfac) ## lrtest(modelo, modelo.hdc.exp) lrtest(modelo.hdc.fix, modelo.hdc.exp) lrtest(modelo.hdc.exp, modelo.hdc.expfac) ## elegimos el modelo con menor AICc ## y lo corremos usando method="REML" mimodelo <- update(modelo.hdc.expfac, .~., method="REML") ## valoración de los RESIDUOS del modelo ## normalidad par(mfcol=c(1,2)) hist(residuals(modelo), main="¡¡modelo sin correcciones!!") hist(residuals(mimodelo), density=5, freq=FALSE, main="¡¡modelo con correcciones!!") curve(dnorm(x, mean=mean(residuals(mimodelo)), sd=sd(residuals(mimodelo))), col="red", lwd=2, add=TRUE, yaxt="n") par(mfcol=c(1,1)) ## par(mfcol=c(1,2)) qqnorm(residuals(modelo), main="¡¡modelo sin correcciones!!") qqline(residuals(modelo), col="red", lwd=2) qqnorm(residuals(mimodelo), main="¡¡modelo con correcciones!!") qqline(residuals(mimodelo), col="red", lwd=2) par(mfcol=c(1,1)) ## ## heterocedasticidad de los residuos plot(modelo, main="¡¡modelo sin correcciones!!") plot(mimodelo, main="¡¡modelo con correcciones!!") ## INDEPENDENCIA ENTRE VARIABLES PREDICTORAS (colinealidad de los predictores) ## VIF = 1 / (1 - R2), de cada predictora por las restantes ## VIF specifically indicates the magnitude of the inflation in the standard errors ## associated with a particular beta weight that is due to multicollinearity pairs(eqt, data=datos) ## representación gráfica de las relaciones entre variables vif(mimodelo) sqrt(vif(mimodelo)) sqrt(4.070960) ## esas veces se inflan los errores standard de esa predictora ## OMNIBUS TEST: significación global del modelo lrtest(mimodelo) ## VARIACION EXPLICADA: pseudo R2 pseudoR2.gls <- cor(getResponse(mimodelo), predict(mimodelo))^2 print(c("Pseudo R2 =", round(pseudoR2.gls, 3)), quote=F) ## RESULTADOS summary(mimodelo) summary(mimodelo)$tTable coeftest(mimodelo) ## estima asintótica asumiendo infinitos grados de libertad ## estima con los grados de libertad correctos idéntica a summary(mimodelo)$tTable df.residual <- mimodelo$dims$N-mimodelo$dims$p coeftest(mimodelo, df=df.residual) ############################################ ## CONTROL DE LA AUTOCORRELACION ESPACIAL ## ############################################ plot(datos$longitud, datos$latitud) ## fórmula ("sizepa" es un factor) eqt <- as.formula(nspp ~ altmed+rangoalt+shannon+tempmin+precip+sizepa) ## MODELOS ## http://rfunctions.blogspot.com/2017/06/how-to-identify-and-remove-spatial.html ## corrección de la autocorrelación espacial en los residuos ## (i.e., diferente distancia entre unidades muestrales) ## tenemos las clases corExp, corGaus, ... ## help(corClasses) ## lleva mucho tiempo modelo <- gls(eqt, data=datos, method="ML") modelo.esp.lin <- gls(eqt, data=datos, method="ML", correlation=corLin(form=~longitud+latitud, nugget=T)) modelo.esp.exp <- gls(eqt, data=datos, method="ML", correlation=corExp(form=~longitud+latitud, nugget=T)) modelo.esp.gau <- gls(eqt, data=datos, method="ML", correlation=corGaus(form=~longitud+latitud, nugget=T)) modelo.esp.sph <- gls(eqt, data=datos, method="ML", correlation=corSpher(form=~longitud+latitud, nugget=T)) ## semivariogramas ## https://en.wikipedia.org/wiki/Variogram ## http://gisgeography.com/semi-variogram-nugget-range-sill/ ## http://pro.arcgis.com/es/pro-app/help/analysis/geostatistical-analyst/understanding-a-semivariogram-the-range-sill-and-nugget.htm plot(Variogram(modelo, form=~longitud+latitud, resType="normalized")) plot(Variogram(modelo.esp.lin, form=~longitud+latitud, resType="normalized")) plot(Variogram(modelo.esp.exp, form=~longitud+latitud, resType="normalized")) plot(Variogram(modelo.esp.gau, form=~longitud+latitud, resType="normalized")) plot(Variogram(modelo.esp.sph, form=~longitud+latitud, resType="normalized")) ## comparación de modelos AICc(modelo, modelo.esp.lin, modelo.esp.exp, modelo.esp.gau, modelo.esp.sph) model.sel(modelo, modelo.esp.lin, modelo.esp.exp, modelo.esp.gau, modelo.esp.sph) ## ## comparación de dos modelos entre si anova(modelo.esp.exp, modelo.esp.lin) anova(modelo.esp.exp, modelo) lrtest(modelo.esp.exp, modelo.esp.lin) lrtest(modelo.esp.exp, modelo) ## elegimos el modelo con menor AICc ## y lo corremos usando method="REML" mimodelo <- update(modelo.esp.exp, .~., method="REML") ## valoración de los RESIDUOS del modelo ## normalidad par(mfcol=c(1,2)) hist(residuals(modelo), main="¡¡modelo sin correcciones!!") hist(residuals(mimodelo), density=5, freq=FALSE, main="¡¡modelo con correcciones!!") curve(dnorm(x, mean=mean(residuals(mimodelo)), sd=sd(residuals(mimodelo))), col="red", lwd=2, add=TRUE, yaxt="n") par(mfcol=c(1,1)) ## par(mfcol=c(1,2)) qqnorm(residuals(modelo), main="¡¡modelo sin correcciones!!") qqline(residuals(modelo), col="red", lwd=2) qqnorm(residuals(mimodelo), main="¡¡modelo con correcciones!!") qqline(residuals(mimodelo), col="red", lwd=2) par(mfcol=c(1,1)) ## ## heterocedasticidad de los residuos plot(modelo, main="¡¡modelo sin correcciones!!") plot(mimodelo, main="¡¡modelo con correcciones!!") ## podemos actualizar el modelo de interés añadiendo otra estructura de correlación de residuos ## y lo corremos usando method="REML" ## muy lento mimodelo <- gls(eqt, data=datos, method="REML", correlation=corExp(form=~longitud+latitud, nugget=T), weights=varExp()) ## plot(mimodelo, main="¡¡modelo con correcciones!!") ## plot(Variogram(mimodelo, form=~longitud+latitud)) ## INDEPENDENCIA ENTRE VARIABLES PREDICTORAS (colinealidad de los predictores) ## VIF = 1 / (1 - R2), de cada predictora por las restantes ## VIF specifically indicates the magnitude of the inflation in the standard errors ## associated with a particular beta weight that is due to multicollinearity pairs(eqt, data=datos) ## representación gráfica de las relaciones entre variables vif(mimodelo) sqrt(vif(mimodelo)) sqrt(3.098340) ## esas veces se inflan los errores standard de esa predictora ## VARIACION EXPLICADA: pseudo R2 ## con predicciones del modelo incluyendo la autocorrelación espacial pseudoR2.gls <- cor(getResponse(mimodelo), predict(mimodelo))^2 print(c("Pseudo R2 =", round(pseudoR2.gls, 3)), quote=F) ## RESULTADOS summary(mimodelo) summary(mimodelo)$tTable df.residual <- mimodelo$dims$N-mimodelo$dims$p coeftest(mimodelo, df=df.residual) coeftest(mimodelo) ## estima asintótica ## tabla de ANOVA Anova(mimodelo, type=3, test="Chisq") ## ## efectos parciales (con un std.err) plot(allEffects(mimodelo)) ## para efectos principales; poner en "..." lo que sea plot(Effect(mod=mimodelo, focal.predictors="...")) ############################################## ## CONTROL DE LA AUTOCORRELACION TEMPORAL ## ############################################## ## IMPORTAMOS DATOS ## de internet (web en lmcarrascal.eu) meconecto <- url("http://www.lmcarrascal.eu/cursos/regravanz/anfibios.RData") load(meconecto) close(meconecto); rm(meconecto) ## names(datos) ## CONSTRUCCIÓN DE LA ECUACIÓN DEL MODELO GLS eqt <- as.formula(I(log(trimar)) ~ TMINmed.6 + TMEDmayago.6) ## ## ¡IMPORTANTE CORRER ESTA LINEA! en el caso de que existan factores options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## plots plot(datos$year, datos$TMINmed.6, main="temperetura media de las mínimas durante 6 años"); lines(datos$year, datos$TMINmed.6, type="l") plot(datos$year, datos$TMEDmayago.6, main="temperetura media mayo-agosto durante 6 años"); lines(datos$year, datos$TMEDmayago.6, type="l") plot(datos$year, datos$trimar, main="larvas por año"); lines(datos$year, datos$trimar, type="l") ## 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 años previos ## help(corClasses) modelo.lin <- gls(eqt, data=datos, method="ML") modelo.car1 <- gls(eqt, data=datos, method="ML", corr=corCAR1(form=~year)) ## equivalente al siguiente modelo.arma.10 <- gls(eqt, data=datos, method="ML", corr=corARMA(form=~year, p=1, q=0)) modelo.arma.01 <- gls(eqt, data=datos, method="ML", corr=corARMA(form=~year, p=0, q=1)) modelo.arma.11 <- gls(eqt, data=datos, method="ML", corr=corARMA(form=~year, p=1, q=1)) modelo.arma.20 <- gls(eqt, data=datos, method="ML", corr=corARMA(form=~year, p=2, q=0)) modelo.arma.02 <- gls(eqt, data=datos, method="ML", corr=corARMA(form=~year, p=0, q=2)) modelo.arma.22 <- gls(eqt, data=datos, method="ML", corr=corARMA(form=~year, p=2, q=2)) ## comparación de modelos AICc(modelo.lin, modelo.car1, modelo.arma.10, modelo.arma.01, modelo.arma.11, modelo.arma.20, modelo.arma.02, modelo.arma.22) ## elegimos un modelo mimodelo <- modelo.arma.10 ## análisis de los RESIDUOS plot(mimodelo) ## qqnorm(residuals(mimodelo)) ## RESULTADOS summary(mimodelo)$tTable ## VARIACION EXPLICADA: pseudo R2 pseudoR2.gls <- cor(getResponse(mimodelo), predict(mimodelo))^2 print(c("Pseudo R2 =", round(pseudoR2.gls, 3)), quote=F) ## rehacemos el modelo corrigiendo la heterocedasticidad mimodelo <- update(mimodelo, .~., weights=varExp()) summary(mimodelo)$tTable ####################################################### ## MODELOS MIXTOS CON ## ## CONTROL DE LA AUTOCORRELACION TEMPORAL ## ## PARA DATOS TEMPORALES NO UNIFORMEMENTE ESPACIADOS ## ####################################################### ## IMPORTAMOS DATOS ## de internet (web en lmcarrascal.eu) meconecto <- url("http://www.lmcarrascal.eu/cursos/regravanz/datosgamm.RData") load(meconecto) close(meconecto); rm(meconecto) ## convertimos una variable carácter en factor class(datos$factor) datos$factor <- as.factor(datos$factor) ## CONSTRUCCIÓN DE LA ECUACIÓN DEL MODELO GLS ## definimos la posibilidad de establecer relaciones no lineales meniante el uso de polinomios ortogonales ## haciendo uso del comando poly help(poly) eqt <- as.formula(humedad ~ poly(precip_24h, degree=2) + poly(ndays_dry, degree=2) + poly(temperatura, degree=2) + seno_fecha + coseno_fecha + factor) ## ¡IMPORTANTE CORRER ESTA LINEA! en el caso de que existan factores options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## estructura de los contrastes utilizada levels(datos$factor) contr.sum(2) contr.sum(levels(datos$factor)) ## tenemos el día calendario comenzando por 1 - 1-Enero plot(datos$fecha, datos$humedad); lines(datos$fecha, datos$humedad, type="l") ## corrección de la autocorrelación temporal en los residuos ## tenemos la clase corCAR1 aplicado a una fecha con datos irregularmente espaciados ## help(corClasses) ## además contamos con un factor aleatorio denominado LOCATION que controla la pseudorreplicación muestral modelo.car1.lme <- lme(eqt, data=datos, method="REML", corr=corCAR1(form=~fecha), random=~1|LOCATION) modelo.car1.gls <- gls(eqt, data=datos, method="REML", corr=corCAR1(form=~fecha|LOCATION)) ## ANALISIS DE LOS RESIDUOS plot(modelo.car1.lme) plot(modelo.car1.gls) ## qqnorm(residuals(modelo.car1.lme), main="Modelo lme") qqnorm(residuals(modelo.car1.gls), main="Modelo gls") ## RESULTADOS ## modelo mixto de intercepta aleatoria summary(modelo.car1.lme)$tTable ## modelo que no tiene en cuenta el factor aleatorio summary(modelo.car1.gls)$tTable ## omnibus test: significación global del modelo lrtest(modelo.car1.lme) lrtest(modelo.car1.gls) ## pseudo R2 pseudoR2.lme <- cor(getResponse(modelo.car1.lme), predict(modelo.car1.lme))^2 print(c("Pseudo R2 =", round(pseudoR2.lme, 3)), quote=F) ## pseudoR2.gls <- cor(getResponse(modelo.car1.gls), predict(modelo.car1.gls))^2 print(c("Pseudo R2 =", round(pseudoR2.gls, 3)), quote=F) ## rehacemos el modelo corrigiendo la heterocedasticidad modelo.car1.lme.hc <- update(modelo.car1.lme, .~., weights=varComb(varIdent(form=~1|factor), varExp(form=~temperatura+ndays_dry))) plot(modelo.car1.lme.hc) ## modelo.car1.gls.hc <- update(modelo.car1.gls, .~., weights=varComb(varIdent(form=~1|factor), varExp(form=~temperatura+ndays_dry))) plot(modelo.car1.gls.hc)