############################### ## Regresiones por Quantiles ## ############################### ## Luis M. Carrascal ## https://lmcarrascal.eu ## revisión 22/11/2021 ## corrido bajo R-3.6.3 ## REFERENCIAS INTERESANTES ## http://dx.doi.org/10.1890/1540-9295(2003)001[0412:AGITQR]2.0.CO;2 ## https://www.r-bloggers.com/2019/01/quantile-regression-in-r-2/ ## https://cran.r-project.org/web/packages/quantreg/vignettes/rq.pdf ## https://data.library.virginia.edu/getting-started-with-quantile-regression/ ## http://dx.doi.org/10.1080/01621459.1999.10473882 ## cargamos librermas library(quantreg) ## para regresiones de cuantiles library(MuMIn) ## para AICc library(psych) ## para estadmsticos basicos de las variables usando "describe" library(visreg) ## para visualizar efectos parciales ## 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 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 ## definimos nuestro modelo de interes (eqt) y el modelo nulo (eqt.n) sin predictores ## ¡¡¡ PONED AQUÍ VUESTROS DATOS !!! eqt <- as.formula(stuuni ~ rangoalt+tempmin+precip) eqt.n <- as.formula(stuuni ~1) eqt; eqt.n ## si hay muchas observaciones en la variable respuesta con el mismo valor podrmamos querer desempatarlos al azar ## ¡¡¡ PONED AQUÍ VUESTROS DATOS !!! ... tras datos$ datos$respuesta.t <- dither(datos$stuuni, type="right", value=1) plot(datos$respuesta.t, datos$stuuni) ## si este fuese el caso, re-definimos las ecuaciones del modelo ## ¡¡¡ PONED AQUÍ VUESTROS DATOS !!! ... en eqt eqt <- as.formula(respuesta.t ~ rangoalt+tempmin+precip) eqt.n <- as.formula(respuesta.t ~ 1) eqt; eqt.n ## RELACIONES VISUALES ## creamos antes un modelo OLS lineal que utilizaremos y contiene los datos de analisis relaciones <- lm(eqt, data=datos) pairs.panels(relaciones$model, smooth=T, scale=F, density=T, ellipses=T, digits=2, method="pearson", pch=20, rug=T, cex.cor=2) ## MODELO de quantiles lineales ## existen diferentes mitodos (ved la ayuda de rq) help(rq) ## definimos las taus referentes a los percentiles para los que se estiman las pendientes taus <- c(0.1, 0.25, 0.5, 0.75, 0.90, 0.95) ## otra forma: tau=seq(0.1, 0.9, by=0.1) (primer cuantil, ultimo cuantil, incremento con by=) ## modelo de interes definido por la ecuacion eqt modelo <- rq(eqt, data=datos, method="br", tau=taus) ## modelo nulo para pseudo-R2 modelo.n <- rq(eqt.n, data=datos, method="br", tau=taus) ## VALORAMOS LOS RESULTADOS ## los valores al cuantil 50% (0.50) se denominan regresisn a la mediana ## AIC (no es aplicable AICc de MuMIn) AIC.mod <- AIC(modelo); AIC.nul <- AIC(modelo.n) tabla.AIC <- data.frame(taus, AIC.mod, AIC.nul) tabla.AIC$veces.mejor.mod <- exp(-0.5*(tabla.AIC$AIC.mod-tabla.AIC$AIC.nul)) ## cuantas veces es mejor el modelo que el nulo ## veces que es mejor el modelo a ese cuantil que el modelo nulo round(tabla.AIC, 2) ## PSEUDO-R2; como aproximacion no exactamente equivalente a R2 y D2 de modelos Generales y Generalizados ## "rho" is the weighted distance used to calculate the objective function in the minimisation algorithm ## as tau * pmax(resid, 0) + (1 - tau) * pmin(resid, 0) ## "rho" es una medida analoga a la devianza de los modelos glm con la que hacemos D2 <- 1 - modelo$deviance / modelo$null.deviance rho.nulo <- as.vector(modelo.n$rho) rho.modelo <- as.vector(modelo$rho) pseudoR2 <- 1-rho.modelo/rho.nulo rho.mod.quant <- data.frame(taus, pseudoR2, rho.nulo, rho.modelo) round(rho.mod.quant,4) ## ## sintesis de las dos tablas previas tabla.rho.AIC <- data.frame(round(rho.mod.quant,4), round(tabla.AIC[,c(2:4)], 2)) tabla.rho.AIC ## resultados de los COEFICIENTES tabla.coefs <- as.data.frame(t(modelo$coefficients)) ## tabla sintetica de coeficientes a los cuantiles definidos; t para trasponer round(tabla.coefs, 5) ## vemos lo que sale con precision de cinco decimales ## ## resultados de errores standard, limites y significaciones; mirad "summary.rq" y "rq.fit.br" en la ayuda summary(modelo, se="rank", iid=FALSE, alpha=0.05) ## produces confidence intervals for the estimated parameters by inverting a rank test as described in Koenker (1994) summary(modelo, se="ker") ## poco conservador; uses a kernel estimate of the sandwich as proposed by Powell(1990) summary(modelo, se="nid") ## computes a Huber sandwich estimate using a local estimate of the sparsity ## ## Otra opcion muy lenta pero robusta utilizando bootstrapping de los datos; ## to construct standard errors, confidence intervals and tests of hypotheses using bootstrapping with R replicates ## como es muy lento, para no correrlo cada vez que lo necesitamos, guardamos sus resultados ## en R definimos el numero de replicas bootstrap modelo.boot <- summary(modelo, se="boot", R=2000) ## los resultados de std.error, t y Pr cambian (muy poco) de vez en vez porque son remuestreos modelo.boot ## compaeremos los diferentes resultados para la secuencia de taus taus ## la secuencia de las tau la especificamos en [[...]]; [[6]] para la sexta tau ## las tau las hemos definido previamente con: taus <- c(0.1, 0.25, 0.5, 0.75, 0.90, 0.95) modelo.boot[[6]] summary(modelo, se="nid")[[6]] summary(modelo, se="ker")[[6]] ## RESULTADOS GRAFICOS ## necesitamos realizar quantile regression si las bandas grises del modelo cuantil (IC al 95%) ## quedan fuera de las lineas rojas (IC al 95%) indicativas de los coeficientes de un modelo OLS ## los efectos son significativos si no incluyen la linea gis horizontal (coeficiente = 0) ## las lineas rojas son los coeficientes del modelo lineal OLS (relaciones) coefficients(relaciones) ## plot(summary(modelo, se="rank", iid=FALSE, alpha=0.05), ols=TRUE, lcol="red") plot(summary(modelo, se="ker"), lcol="red", ols=TRUE) plot(summary(modelo, se="nid"), lcol="red", ols=TRUE) plot(modelo.boot, lcol="red", ols=TRUE) ## opciones ## sin lineas rojas de OLS plot(modelo.boot, ols=FALSE) ## sacando solo el panel 3 (param=3) plot(summary(modelo, se="nid"), parm=3, lcol="red") ## otra manera es mediante el uso del paquete {visreg} ## se muestra un panel grafico ara cada variable predictora ## segun la secuencia de percentiles establecida en "taus" ## los puntos grises se refieren a cada unidad muestral ## las lineas de regresion definen efectos parciales ## controlando cada predictora pos las otras predictoras dentro de cada cuantil ## dad para sacar el panel de cada coeficiente visreg(modelo) ## PREDECIMOS USANDO UN NUEVO CONJUNTO DE DATOS ## la matriz que contiene esos datos (con los mismos nobres de las variables usadas) ## la ponemos en newdata=... ## ¡¡ OJo !!! aqui denominada "predecir" predichos <- as.data.frame(predict(modelo, newdata=predecir)) ## asignamos a los valores predichos la variable respuesta original en la matriz predecir ## ¡¡¡ PONED AQUI VUESTROS DATOS !!! e.g.: predecir$stuuni predichos$respuesta <- predecir$stuuni ## plot de predichos y observados en un juego de dato externo ## para el cuantil 0.95 plot(predichos$`tau= 0.95`, predichos$respuesta) ##################################### ### PARA REGRESIONES POLINOMIALES ### ##################################### ## TRANSFORMACIONES POLINOMIALES (de orden 3 a cada predictor) ## pero antes un ejemplo del uso de la transformacion polinomial correcta z <- c(1:10) z.pol3 <- poly(z, degree=3) ## podemos reducirla y poner solo: poly(z, 3) ## contiene una matriz de datos con columnas 1 (lineal), 2 (cuadratico) y 3 (cubico) z.pol3 ## los vectores de la transformación polinomial tienen dos propiedades: ## tienen media "cero" mean(z.pol3[,1]); mean(z.pol3[,2]); mean(z.pol3[,2]) ## mantienen entre ellos correlaciones "cero" (i.e., ortogonales) plot(z.pol3[,1], z.pol3[,2]) plot(z.pol3[,1], z.pol3[,3]) cor(z.pol3[,1], z.pol3[,2]) cor(z.pol3[,1], z.pol3[,3]) ## DEFINIMOS LA FORMULA DE NUESTRO MODELO DE REGRESION POR CUANTILES eqt.pol3 <- as.formula(stuuni ~ poly(rangoalt, 3) + poly(tempmin, 3) + poly(precip, 3)) eqt.n <- as.formula(stuuni ~ 1) ## definimos las TAUS referentes a los percentiles para los que se estiman las pendientes taus <- c(0.1, 0.25, 0.5, 0.75, 0.90, 0.95) ## ESTABLECEMOS EL MODELO DE REGRESION POR CUANTILES ## modelo polinomial de interes modelo.pol3 <- rq(eqt.pol3, data=datos, method="br", tau=taus) ## modelo polinomial nulo polinomial para pseudo-R2 modelo.n <- rq(eqt.n, data=datos, method="br", tau=taus) ## modelo polinomial de relaciones OLS relaciones <- lm(eqt.pol3, data=datos) ## miramos como hemos construido terminos polinomiales que no inflan el VIF (Variance Inflation Factor) ## https://en.wikipedia.org/wiki/Variance_inflation_factor library(car) vif(relaciones) ## asignamos esos modelos a los nombres del apartado anterior ## para poder correr las mismas lineas de codigo modelo <- modelo.pol3 ## poned aqui el modelo deseado "no nulo" ## ## el modelo "nulo" ya viene definido por modelo.n ## VOLVEMOS AL INICIO DE LOS ANALISIS PREVIOS EN: ## linea 81 ## VALORAMOS LOS RESULTADOS ############################################### ### QUANTILE NON-PARAMETRIC ADDITIVE MODELS ### ############################################### library(qgam) library(mgcViz) ## referencias ## https://cran.r-project.org/web/packages/qgam/vignettes/qgam.html uso general ## https://mfasiolo.github.io/mgcViz/articles/qgam_mgcViz.html para efectos aleatorios con s(..., bs="re") ## DEFINIMOS LA FORMULA DE NUESTRO MODELO DE REGRESION ADITIVO POR CUANTILES ## sequimos el mismo esquema utilizado por {mgcv} ## definimos nuestro modelo de interes (eqt) ## los basis splines pueden ser bs="tp" (thin plate), bs="cr" (cubic regression), bs="re" (para factores aleatorios; mixed quantile models) ## ¡¡¡ PONED AQUÍ VUESTROS DATOS !!! eqt <- as.formula(stuuni ~ s(rangoalt, k=10, bs="cr")+s(tempmin, k=10, bs="cr")+s(precip, k=10, bs="cr")) ## MODELOS qgam (para un solo valor de tau definido en qu=), o mqgam (para varios valores de tau) modelo.qg <- qgam(eqt, data=datos, qu=0.9, err=0.05) tau.g <- c(0.5, 0.75, 0.9) modelo.mqg <- mqgam(eqt, data=datos, qu=tau.g, err=0.05) ## MODEL CHECKING ## valoramos si el criterio de calibración usando el "learning rate" ha tenido éxito evaluando la "loss function" ## con esta otra función podemos ver si existe un único mínimo check.learnFast(modelo.qg$calibr) check.learnFast(modelo.mqg$calibr) ## valoración de los residuos negativos para modelos qgam de múltiples taus (produce dos paneles gráficos) ## los circulos deberían ser negros, estar dentro de las líneas horizontales y estar cercano al quantil qu ## también es muy útil para saber si es necesario incrementar la curvilinealidad del spline definida por k check.qgam(modelo.qg) ## para modelos qgam qdo(modelo.mqg, qu=0.9, check) ## para modelos mqgam definiendo con qu cada cuantil ## también podemos hacer una exploración a priori calibracion <- tuneLearn(eqt, data=datos, qu=0.9, lsig=seq(1, 4, length.out=40), control=list("progress"="none", loss="cal", sam="boot")) check(calibracion) ## cantidad de observaciones bajo los cuantiles para cada predictor ## los círculos deberían ser negros y estar dentro de las líneas horizontales y cercanos a la línea ---- del cuantil cqcheck(modelo.qg, v="rangoalt", X=datos, y=modelo.qg$y) cqcheck(modelo.qg, v="tempmin", X=datos, y=modelo.qg$y) ## visualizado para la combinación de dos variables ## los valores anaranjados indican situaciones no bien representadas, tanto más cuanto los valores se desvían del cuantil ## los asteriscos (*) tras los números denotan diferencias significativas respecto a los cuantiles esperados cqcheck(modelo.qg, v=c("rangoalt", "tempmin"), X=datos, y=modelo.qg$y) ## igual que antes pero de manera más visual combinado con {shiny} library(shiny) cqcheckI(modelo.qg, v="rangoalt", X=datos, y=modelo.qg$y) ## RESULTADOS DEL MODELO ## los resultados son solo aproximados, ya que derivan de aplicar la aproximación paramétrica de mgcv, y no de cuantiles ## resultados NUMERICOS ## para modelos qg con un solo cuantil summary(modelo.qg) ## para modelo mqg qdo(modelo.mqg, qu=tau.g, summary) ## resultados GRAFICOS para modelos de un solo cuantil ## la más sencilla de una-en-una plot(modelo.qg) ## más elaboradas (definiendo con se el múltiplo del std.err; 1 ca. 68%, 2 ca. 95%, 5.58 ca. 99%, 3.3 ca. 99.9%) ## todos los predictores en diferente escala en el eje Y plot(modelo.qg, pages=1, shade=T, shade.col="gray70", se=2, all.terms=T, scale=0) ## todos los predictores con la misma escala en el eje Y para hacer comparables sus efectos plot(modelo.qg, pages=1, shade=T, shade.col="gray70", se=2, all.terms=T) ## un solo plot de los N paneles resultantes (uno por predictor), eligiendo uno de ellos con select=... plot(modelo.qg, shade=T, shade.col="gray70", se=2, rug=T, select=3) ## usando el paquete {mgcViz} ## volvemos a correr los modelos qgam o mqgam haciendo uso de las funciones mgcViz::mqgamV o mgcViz::qgamV modelo.mqg.v <- mqgamV(eqt, data=datos, qu=tau.g, err=0.05) modelo.qg.v <- qgamV(eqt, data=datos, qu=0.9, err=0.05) ## presentamos los resultados gráficos sintéticos para modelos mqgamV de múltiples cuantiles print(plot(modelo.mqg.v, allTerms=TRUE), pages = 1) ## con sus intervalos de confianza para cada uno de los cuantiles del modelo [[1]]: primer cuantil, [[2]]: segundo cuantil ... print(plot(modelo.mqg.v[[1]], allTerms=TRUE), pages = 1) print(plot(modelo.mqg.v[[2]], allTerms=TRUE), pages = 1) print(plot(modelo.mqg.v[[3]], allTerms=TRUE), pages = 1) ## y para uno solo de los predictores print(plot(modelo.mqg.v[[3]], allTerms=TRUE, select=3), pages = 1) ## para el check de los residuos negativos a ese cuantil check1D(modelo.qg.v, datos$precip) check1D(modelo.mqg.v[[3]], datos$precip) ## para modelos qgam con PLANOS DE INTERACCION ## ¡¡¡ PONED AQUÍ VUESTROS DATOS !!! eqt2 <- as.formula(stuuni ~ s(rangoalt, k=10, bs="cr") + s(tempmin, precip, k=10, bs="tp")) eqt3 <- as.formula(stuuni ~ s(rangoalt, k=10, bs="cr") + te(tempmin, precip, k=10, bs="tp")) eqt4 <- as.formula(stuuni ~ s(rangoalt, k=10, bs="cr") + ti(tempmin, precip, k=10, bs="tp")) ## modelos usando mgcViz::mqgamV modelo.mqg.v2 <- mqgamV(eqt2, data=datos, qu=tau.g) ## modelo s(x,z) asume isotropía (controla por s(x)+s(z)) modelo.mqg.v3 <- mqgamV(eqt3, data=datos, qu=tau.g) ## modelo te(x,z) no asume isotropía (controla por s(x)+s(z)) modelo.mqg.v4 <- mqgamV(eqt4, data=datos, qu=tau.g) ## modelo ti(x,z) no asume isotropía (no controla por s(x)+s(z)) ## modelo de interés modelo.int <- modelo.mqg.v3 ## resultados numéricos para el cuantil predefinido en las qu=tau.g según su secuencia (primer cuantil [[1]], segundo cuantil [[2]] ...) summary(modelo.int[[3]]) ## resultados gráficos sintéticos para modelos mqgamV de múltiples cuantiles print(plot(modelo.int, allTerms=TRUE), pages = 1) ## para cada uno de los cuantiles del modelo [[1]]: primer cuantil, [[2]]: segundo cuantil ... print(plot(modelo.int[[3]], allTerms=TRUE), pages = 1) ## y para uno solo de los predictores print(plot(modelo.int[[3]], allTerms=TRUE, select=1), pages = 1) print(plot(modelo.int[[3]], allTerms=TRUE, select=2), pages = 1) ## para el check de los residuos negativos a ese cuantil check1D(modelo.int[[3]], datos$precip) ## PREDICCIONES DE LOS MODELOS ## para modelos de un solo cuantil; qgam::qgam y mgcViz::qgamV ## para trabajar con los mismos datos usados para modelizar, con estima del std.error predicciones <- predict(modelo.qg, se=TRUE) predicciones <- predict(modelo.qg.v, se=TRUE) ## plot de valores reales y líneas al valor del cuantil qu predefinido ## CI95%: ca. 1.97; CI99%: ca. 2.60; CI99.5%: ca. 2.84; CI99.9%: ca. 3.34; predecir <- data.frame(matrix(ncol=4, nrow=dim(datos)[1])) colnames(predecir) <- c("y", "x", "fit", "se.fit") predecir[,1] <- modelo.qg$y predecir[,3] <- predicciones$fit predecir[,4] <- predicciones$se.fit ## ¡¡¡ PONED AQUI EN x VUESTROS DATOS!!! predecir[,2] <- datos$tempmin ## no debe fallar esta entrada por error ## ## el plot queda muy bien con modelos qgam de un solo predictor predecir <- predecir[order(predecir$x),] plot(predecir$x, predecir$y, col="grey", ylab="RESPUESTA", xlab="PREDICTOR SELECCIONADO") lines(predecir$x, predecir$fit, col="blue") lines(predecir$x, predecir$fit+1.97*predecir$se.fit, col="red") lines(predecir$x, predecir$fit-1.97*predecir$se.fit, col="red") ## ## para modelos de varios cuantiles; sólo para mgcViz::mqgamV ## seleccionando uno de ellos de la serie puesta en qu=...; [[1]] para el primer cuantil, [[2]] para el segundo cuantil, ... predicciones <- predict(modelo.mqg.v[[3]], se=TRUE) ## y volvemos a la línea 361