Manejo de datos en el tidyverse

Los objetos fundamentales para la modelización son los objetos de datos. En R hay varias clases asociadas: data.frame, tibble. Ambos son similares. Vamos a ver algunos de los objetos que emplearemos en nuestros ejemplos:

⊕Un objeto de datos tiene siempre la misma estructura: Las columnas son variables, las filas son observaciones y cada celda es un dato. Para cada estructura de datos diferente, tenemos un objeto diferente.

⊕Cada columna puede ser un objeto de tipo diferente. En este caso tenemos columnas numéricas (con números en doble precisión y enteros) y columnas con factores (variables cualitativas)

⊕La tiabla (tibble) tiene la ventaja de que al imprimirse sólo se muestran la información básica de la estructura y las primeras filas

library(tidyverse)
library(gamlss.data)
data(rent99)
rent99 = rent99 %>% as_tibble
rent99

## # A tibble: 3,082 x 9
##     rent rentsqm  area yearc location bath  kitchen cheating district
##    <dbl>   <dbl> <int> <dbl> <fct>    <fct> <fct>   <fct>       <int>
##  1  110.    4.23    26  1918 2        0     0       0             916
##  2  243.    8.69    28  1918 2        0     0       1             813
##  3  262.    8.72    30  1918 1        0     0       1             611
##  4  106.    3.55    30  1918 2        0     0       0            2025
##  5  133.    4.45    30  1918 2        0     0       1             561
##  6  339.   11.3     30  1918 2        0     0       1             541
##  7  215.    6.94    31  1918 1        0     0       0             822
##  8  323.   10.4     31  1918 1        0     1       1            1713
##  9  216.    6.76    32  1918 1        0     0       0            1812
## 10  245.    7.43    33  1918 2        0     0       0             152
## # ... with 3,072 more rows

⊕Datos sobre taquillaje de películas en los años 90. Incluyen el número de pantallas en que se mostró la primera semana y los ingresos de la primera semana, los ingresos tras la primera semana (nuestra variable de interés) y el tipo de distribuidora.

data(film90)
film90= film90 %>% as_tibble
film90

## # A tibble: 4,031 x 4
##    lnosc lboopen lborev1 dist       
##    <dbl>   <dbl>   <dbl> <fct>      
##  1  6.93   14.1    13.9  Independent
##  2  1.79    9.72    9.16 Independent
##  3  7.29   14.2    14.3  Independent
##  4  1.39    8.22    8.44 Independent
##  5  0       9.24   12.1  Independent
##  6  0       7.92   10.7  Independent
##  7  1.39    7.56    9.38 Independent
##  8  1.95    9.03   10.5  Independent
##  9  0       6.64    7.12 Independent
## 10  0       8.87   11.8  Independent
## # ... with 4,021 more rows

La siguiente tabla resume los paquetes principales para cargar datos en R:

gather

spread

separate

Para manejar los datos una vez están en R, las herramientas claves son los paquetes tidyr y dplyr que forman parte del tidyverse.

tidyr lo empleamos cuando queremos cambiar la estructura de los datos para llevarlos a la forma ordenada: que las columnas se correspondan con las variables de nuestro análisis. Las tres funciones principales de tidyr son las siguientes:

Estas funciones son muy útiles también para el manejo de los resultados o la pre-preparación de gráficos.

El tidyverse incluye numerosas funciones para manipular objetos de datos. Recordamos las principales funciones del paquete dplyr. La filosofia del tidyverse para manipular datos es aplicar de forma sucesiva las funciones elementales que sean necesarias separadas con el operador después %>%.

En la selección de observaciones con filter tenemos que emplear condiciones lógicas como las siguientes:

group_by y summarize según Ismay and Kim (⊕2018Ismay, Chester, and Albert Y. Kim. 2018. ModernDive: An Introduction to Statistical and Data Sciences via R. moderndive.com. https://moderndive.com/.)

Los tipos de join en Wickham and Grolemund (⊕2017Wickham, Hadley, and Garrett Grolemund. 2017. R for Data Science: Import, Tidy, Transform, Visualize, and Model Data. 1st ed. O’Reilly Media, Inc. https://r4ds.had.co.nz/.)

La potencia de estas funciones tan básicas se deriva de la combinación con la declaración de grupos mediante group_by. group_by declara que la tabla está compuesta por distintas subtablas definidas a partir de la variable o variables de grupo. Después de declarados los grupos, tanto mutate como arrange o summarize operan dentro de cada subtabla. En particular, summarize cambia la estructura de los datos calculando una fila para cada combinación de las variables de grupo.

Por último, una última operación de gran utilidad la constituyen los joins. Mediante los joins, y en particular el left_join: estos permiten pegar a nuestro objeto de datos variables que procedan de otro objeto relacionado. Los requisitos para que funcione bien es que las variables relacionadas tengan el mismo nombre (se pueden renombrar antes del join) y que contengan información del mismo tipo (ejemplo: numérica, factores, …)

Podéis encontrar más información sobre manejo de datos en el tidyverse, con ejemplos, en este enlace de José Antonio Ortega o, mucho más detallado, en Ginolhac, Koncina, and Krause (⊕2017Ginolhac, Aurélien, Eric Koncina, and Roland Krause. 2017. Data Processing with R Tidyverse. University of Luxembourg. https://lsru.github.io/tv_course/.). Para una introducción rápida podéis consultar esta simpática presentación de Omaima Said. El texto de Kaplan and Shaw (⊕2016Kaplan, Daniel, and Frank Shaw. 2016. Statistical Modeling: Computational Technique. mosaic-web.org. http://mosaic-web.org/go/SM2-technique/.) proporciona una muy buena introducción a R orientada hacia la modelización utilizando el paquete mosaic

Gráficos ggplot2

La herramienta básica para visualización con el tidyverse es el paquete ggplot2. Los gráficos ggplot2 están compuestos por distintas capas separadas por +. La información necesaria mínima para un gráfico es:

• Los datos. Normalmente estarán en forma larga. data, primer argumento.
• Los atributos estéticos, aes, del gráfico: Qué representan los ejes, para qué utilizo el color, …
• Los objetos geométricos, geom, que añado al gráfico. Separados por +

Ingresos de las películas tras la primera semana en función de los ingresos de la primera semana (en logaritmos). Aparece superpuesta la estimación de un modelo de regresión no paramétrico gam

film90 %>% ggplot(aes(y=lborev1,x=lboopen)) + 
  geom_point() +
  geom_smooth()

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Del gráfico básico anterior puedo pasar a uno más completo cambiando lo que desée del aspecto del gráfico. En particular, es posible establecer para cada capa geométrica parámetros como el tamaño, el color, el tipo de línea o la transparencia que pueden asociarse a variables dentro de aes() o ser establecidos de manera directa fuera de aes(). Las declaraciones de ggplot() son globales, las de cada geom_xxx() son locales.

film90 %>% ggplot(aes(y=lborev1,x=lboopen)) + 
  geom_point(aes(color=dist,size=lnosc),alpha=0.3) + 
  geom_smooth() +
  labs(color="Productor",
       y="log de ingresos resto",
       x="log de ingresos en primera semana",
       size="log de número de pantallas") +
  theme_classic(20) +
  theme(legend.position="bottom",legend.box="vertical") 

Ingresos de las películas de los 90 utilizando las posibilidades de ggplot: Área proporcional al número de pantallas (log), y color asociado al tipo de productor.

plotly::ggplotly()

Como referencia podemos consultar la chuleta de RStudio sobre ggplot, o las secciones correspondientes de los libros online de Wickham and Grolemund (⊕2017Wickham, Hadley, and Garrett Grolemund. 2017. R for Data Science: Import, Tidy, Transform, Visualize, and Model Data. 1st ed. O’Reilly Media, Inc. https://r4ds.had.co.nz/.), Ismay and Kim (⊕2018Ismay, Chester, and Albert Y. Kim. 2018. ModernDive: An Introduction to Statistical and Data Sciences via R. moderndive.com. https://moderndive.com/.) o Anderson and Severson (⊕2018Anderson, Brooke, and Rachel Severson. 2018. R Programming for Research. Colorado State University. https://geanders.github.io/RProgrammingForResearch/.). Para un análisis más detallado se puede consultar Healy (⊕2018Healy, Kieran. 2018. Data Visualization: A Practical Introduction. Forthcoming, Princeton University Press. https://socviz.co/.) o esta presentación muy completa de Ginolhac, Koncina, and Krause (2017).

Teoría de la estimación

Para introducirnos en la modelización vamos a centrarnos ahora en modelos de la relación entre dos únicas variables. La teoría se puede encontrar más desarrollada en James et al. (⊕2014James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2014. An Introduction to Statistical Learning: With Applications in R. Springer Publishing Company, Incorporated. http://www-bcf.usc.edu/~gareth/ISL/.) y Shalizi (⊕2018Shalizi, Cosma Rohilla. 2018. Advanced Data Analysis from an Elementary Point of View. Forthcoming, Cambridge University Press. http://www.stat.cmu.edu/~cshalizi/ADAfaEPoV/.). Es extensible al caso en que \(x\) sea un vector sin más que cambiar \(x\) por \(X\).

La función \(f(x)\) capta la parte sistemática de \(y\), mientras que el término de error o componente no observado, \(\varepsilon\), capta la parte no explicada por el modelo.

\[y= f(x) + \varepsilon\]

Ajuste de un modelo dado por \(\hat{f}\) y residuos como distancia vertical a la función. La suma de cuadrados de los residuos es proporcional a las áreas de los puntos, cuyo radio es proporcional al residuo.

• El supuesto básico para estimar los modelos es el que garantiza que la función \(f\) represente la esperanza condicionada de la variable \(y\). El motivo de hacer este supuesto es que garantiza que cometamos el error cuadrático mínimo al predecir utilizando \(f\):

Supuesto básico: Esperanza condicionada nula \(E(\varepsilon/x)=0\)

Desde el punto de vista de la estimación, en todos los casos tendremos una función de R que escoge de entre qué familia de funciones escogemos en el ajuste.

• Ajuste: Se basa en optimizar una función objetivo. La más común es la suma de cuadrados de los residuos que da lugar a los estimadores mínimo cuadráticos.

• Valor ajustado y predicción: La función resultante \(\hat{f}\) me define, cuando la utilizo dentro de la muestra, el valor ajustado, \(\hat{y}\), y cuando la utilizo para observaciones fuera de la muestra, la predicción puntual.

\[\hat{y}=\hat{f}(x)\]

• Error reducible e irreducible: Otra forma de ver el término de error es como componente no reducible de error en la predicción. En la medida en que mi método de estimación sea adecuado y se adapte a la naturaleza del problema puedo aspirar a que mi función \(\hat{f}\) se aproxime a la verdadera función \(f\). Llamamos error reducible a la relación \(f-\hat{f}\) que determina en qué medida mi estimación está sesgada. Si consigo reducir el sesgo conseguiré predecir mejor, pero no acertaré debido a la variabilidad que tiene la \(y\) para un valor dado de la \(x\). Este término de error no es reducible utilizando sólo la información de la variable \(x\), por ello lo llamamos error no reducible. Podemos conseguir reducir el error si existe información relevante fuera del modelo que se puede utilizar. Esto implica cambiar de modelo.

• Residuos: El residuo es el equivalente muestral del término de error, \(e=y-\hat{f}(x)\). Puede contener dos partes: una sistemática asociada a los sesgos en que hemos incurrido al estimar, \(f-\hat{f}\), y otra aleatoria asociada al término de error, \(\varepsilon\).

• Error cuadrático medio: Supongamos que utilizamos una función \(\hat{f}\) fija para estimar, que no es la verdadera \(f\). El valor esperado de los errores de estimación al cuadrado se puede descomponer como

\[ECM_{\hat{f}}(x)=E_{\hat{f}}(y - \hat{f}(x))^2 = [f(x) - \hat{f}(x)]^2 + Var_x(\varepsilon)= \text{Sesgo}^2+ Varianza\]

Conflicto sesgo-varianza en la estimación según Fortmann-Roe(2012).

• Conflicto sesgo-varianza: La descomposición anterior condicionaba a una \(\hat{f}\) concreta. Pero la \(\hat{f}\) estimada dependerá de los datos: unas veces me equivocaré más al estimar, otras veces menos. Si consideramos las distintas posibilidades de estimación, podemos descomponer el error de estimación en dos componentes: el error sistemático de estimación dado por el sesgo esperado \(f(x)-E{\hat{f}(x)}\) y la varianza de la estimación \(Var(\hat{f}(x))\). Se demuestra que el error cuadrático medio es igual a:

\[ECM(x)=\left(f(x)-E[\hat{f}(x)]\right)^2 + Var\left(\hat{f}(x)\right) + Var_x(\varepsilon)\] \[(\text{Sesgo de la predicción})^2+ \text{Varianza de la predicción} + \text{Error irreducible}\]

Complejidad del modelo, subajuste y sobreajuste en un ejemplo Fuente.

De este modo tenemos que lo bien que puede funcionar nuestro modelo se topa con tres limitaciones:

• La variabilidad intrinseca asociada al error irreducible: Hay fenómenos más difíciles de predecir que otros. Dónde caen los dardos una vez en concreto.
• Nuestro estimador puede estar equivocándose de manera sistemática (estimadores sesgados).
• Nuestro estimador puede ser más o menos dependiente de los datos concretos (Varianza de la predicción).

¿Cómo tiran los dardos el que sobreajusta y el que subajusta? Fuente.

• Complejidad del modelo, y conflicto sesgo-varianza: Cada modelo tiene una forma de ver la realidad. Los estimadores más complejos tienden a adaptarse mejor a las posibles situaciones, pero al coste de variar mucho en función de los datos:

Con modelos más complejos podemos reducir el sesgo, pero el coste es, en general, aumentar la varianza del estimador.

Para cada problema habrá un grado de complejidad óptimo. Si la relación entre las variables es aproximadamente lineal, los modelos lineales serán aproximaciones difíciles de mejorar, pero si la relación no es lineal, habrá otros estimadores más complejos que serán mejores, porque el modelo lineal será sesgado.

Sobreajuste y complejidad del modelo según Fortmann-Roe(2012).

• Sobreajuste en la muestra y validación predictiva: Si sólo trabajamos con una muestra, el riesgo es que nuestro modelo capte exactamente la relación entre las variables en la muestra, pero que ésta no sea generalizable a otras muestras posibles. Esto no es lo que queremos. Queremos modelos que predigan bien fuera de la muestra. Las técnicas de validación predictiva nos permiten protegernos frente al sobreajuste. Consisten en, para elegir entre modelos alternativos, contemplar cómo funcionan no en la muestra que se utilizado para estimar (muestra de trabajo), sino en otra/s muestra/s no empleadas en la estimación: muestra de test. Cuando un modelo ajusta bien en la muestra de trabajo y mal en la muestra de test, indica un problema de sobreajuste. El sobreajuste tiende a ser mayor cuanto más complejo sea el modelo. Para cada problema de predicción habrá un grado de complejidad que minimice el error de predicción fuera de la muestra, mientras que con modelos más complejos podré seguir reduciendo el error en la muestra de trabajo (training error)

Modelos de regresión con una variable explicativa en R

Vamos a centrarnos en un modelo para los datos de ingresos de las películas tras la primera semana utilizando como única variable los ingresos de la primera semana. En este caso los datos ya están preparados para poder estimar los modelos. Por razones de conveniencia y claridad, únicamente renombraremos las variables: lboopen, el logaritmo de la taquilla de la primera semana, como taq0, y lborev1, el logaritmo del taquillaje en el resto de semanas como taq1.

film90 = film90 %>% rename(taq0=lboopen,taq1=lborev1)

Desde el punto de vista del comportamiento predictivo tenemos dos pasos que hacer con cada modelo: estimación y predicción. Para comparar inicialmente cómo funcionan los modelos emplearemos la función predict. La función predict es una función genérica de R que generalmente toma como primer argumento el modelo y como segundo, opcional, un objeto de datos en el que predecir. A menudo también es posible obtener intervalos de predicción analíticos o errores estándar de predicción analiticos especificando los argumentos adecuados.

⊕Las funciones genéricas de R examinan la clase del objeto del primer argumento. La clase se puede obtener como class(modelo). Por ejemplo, para modelos lineales es lm. Esto implica que la función que realmente se ejecuta es predict.lm. Para cada tipo de modelos existirá la función predict.estimador correspondiente. Una de las partes desafortunadas es que los argumentos y los resultados de la función pueden ser distintos dependiendo del modelo.

• Estimación (o ajuste): mxxx = estimador(lborev1~lboopen, data=film90)
• Resumen del modelo: summary(mxxx)
• Predicción: predict(mxxx,tibble()) o con la función de mosaic, mod_fun(modelo)

Estimadores paramétricos: lm

⊕poly(var,k) añade un polinomio de orden k para la variable var. Por defecto los polinomios se parametrizan utilizando una base ortogonal. Si se desea utilizar el polinomio con la fórmula habitual especificamos poly(var,k,raw=TRUE). ns(var,df), del paquete splines, añade splines cúbicos naturales. bs(var,df,degree=k) permite añadir splines de distinto orden. Es posible controlar la posición de los nodos con knots y Boundary.knots

## Estimación de modelos
# Paramétricos lineales
ml=lm(taq1~taq0,film90) # lineal
mp2=lm(taq1~poly(taq0,2),film90) #cuadrático
mp3=lm(taq1~poly(taq0,3),film90) #cúbico
mp4=lm(taq1~poly(taq0,4),film90) #pol. grado 4º
mp6=lm(taq1~poly(taq0,6),film90) #pol. grado 6º
mp8=lm(taq1~poly(taq0,8),film90) #pol. grado 8º
library(splines)
msl=lm(taq1~bs(taq0,df=4,degree=1),film90)
msc5=lm(taq1~ns(taq0,5),film90) #Splines cúbicos
msc6=lm(taq1~ns(taq0,6),film90) # con más grados de libertad
msc7=lm(taq1~ns(taq0,7),film90) 
msc8=lm(taq1~ns(taq0,8),film90) 

Los splines definen funciones polinomiales locales en tramos definidos por los nodos. Desde el punto de vista de la estimación son modelos lineales donde hay que estimar el peso de cada elemento de la base. En el diagrama la base de splines del modelo msc6

Todos los modelos estimados son paramétricos y lineales en los parámetros. Por ello se han estimado todos con lm. Que los modelos son lineales en los parámetros quiere decir que la función de predicción consiste en una combinación lineal de funciones que se pueden calcular con los datos, no que el modelo sea lineal en la variable. Únicamente sería lineal el ml según ese criterio. En lo que difieren los modelos estimados es en su complejidad. La complejidad aumenta con el orden del polinomio o con los grados de libertad de los splines.

Podemos “ver” las estimaciones de uno de estos modelos. Los gráficos de diagnóstico en dos formatos, y el ajuste del modelo msc6 comparado con uno de complejidad similar, el mp6, y con el modelo lineal. Para esto último utilizamos una función del paquete mosaic que es más cómoda que predict en este contexto: mod_fun(modelo) da como resultado una función que evalua la predicción del modelo para los valores de las x que se especifiquen.

library(mosaic); library(mosaicModel)
film90 %>% ggplot(aes(y=taq1,x=taq0))+ geom_point(size=0.2) +
  geom_smooth(aes(color="ml"),method="lm",se=FALSE,size=1.2) +
  geom_line(aes(y=model_output,color="mp6"),size=1.2,data=mod_fun(mp6)(taq0=seq(4,18,by=.1))) +
  geom_line(aes(y=model_output,color="msc6"),size=1.2,data=mod_fun(msc6)(taq0=seq(4,18,by=.1))) +
  labs(color="Modelo")+
  theme_bw(18)

En este caso podemos ver como el modelo lineal es un caso típico de subajuste: no es capaz de adaptarse a la naturaleza de los datos y será un estimador sesgado. El modelo polinómico, por su parte, presenta síntomas de sobreajuste e inestabilidad, como puede ser la caída fuerte para valores muy altos. Esto ocurre porque al tratarse de un ajuste global, y de haber pocos puntos en la zona alta, al modelo “no le importa” que pasa con esos puntos si puede mejorar el ajuste en otra región. El modelo de splines parece adaptarse a los datos sin que se observe sobreajuste.

Gráficos de diagnóstico por defecto, con R base.

Gráficos de diagnóstico por defecto, con R base.

Gráficos de diagnóstico por defecto, con R base.

Gráficos de diagnóstico por defecto, con R base.

plot(msc6)

summary(msc6)

## 
## Call:
## lm(formula = taq1 ~ ns(taq0, 6), data = film90)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4842 -0.7057 -0.0786  0.6799  7.8536 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.5235     0.3315   13.64   <2e-16 ***
## ns(taq0, 6)1   8.0214     0.3094   25.92   <2e-16 ***
## ns(taq0, 6)2  10.3332     0.3856   26.80   <2e-16 ***
## ns(taq0, 6)3   8.9515     0.3458   25.89   <2e-16 ***
## ns(taq0, 6)4  11.6965     0.2199   53.20   <2e-16 ***
## ns(taq0, 6)5  16.6452     0.7248   22.97   <2e-16 ***
## ns(taq0, 6)6  12.8291     0.3184   40.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.397 on 4024 degrees of freedom
## Multiple R-squared:  0.8161, Adjusted R-squared:  0.8159 
## F-statistic:  2977 on 6 and 4024 DF,  p-value: < 2.2e-16

Vemos que el ajuste dentro de la muestra es muy bueno, de un 81.6135444% de la variabilidad de la y. Por otro lado, la alta significatividad de los parámetros estimados es también síntoma de que no estamos sobreajustando. En todo caso podríamos utilizar modelos más complejos.

library(ggfortify)
autoplot(msc6,size=.2) + theme_bw(14)

Gráficos de diagnóstico ggplot con ggfortify::autoplot

Las gráficas de diagnóstico nos muestran que capta bastante bien la forma funcional sin que se detecten no linealidades en el gráfico de residuos frente a valores ajustados. Esta es la hipótesis más fundamental del modelo. Sin embargo, se detectan problemas en los otros tres gráficos:

• El gráfico qqplot (normal Q-Q) indica que la distribución del término de error dista bastante de la normal. En particular tiene colas anchas tanto en la parte alta como baja de la distribución. Esto afecta también a los intervalos de predicción que no tendrían el nivel nominal si se han construido para una normal.

• El gráfico scale-location sirve para detectar heterocedasticidad. En este caso se ve qu hay un nivel (en torno a 14) en que cambia la varianza. Las implicaciones no son en principio graves: no afecta a la consistencia o insesgadez de los estimadores, pero exige que si queremos realizar inferencia sobre los parámetros deberíamos utilizar un estimador robusto frente a heterocedasticidad. La fórmula por defecto no es válida. También implica que las bandas de predicción puedan calcularse mal cuando se utilice predict. Por otro lado, si combinamos la información del qqplot con ésta y el gráfico del modelo, podríamos decir que el problema es que hay dos tipos de películas de entre las que hicieron una taquilla media la primera semana: algunas que consiguen mejorar y hacer taquillajes estupendos, y otras que se quedan más abajo. Esto no ocurre cuando taq0 es muy alta o baja. El problema, visto de esta manera, es que la distribución de la variable es distinta para distintos niveles de taq0. lm, aunque pueda funcionar cuando la distribución no es normal o no es homogénea, deja de ser el mejor estimador posible, y es posible que existan estimadores mejores que se adapten a esta distribución complicada.

• El gráfico del apalancamiento frente a los residuos. nos muestra que hay ciertas observaciones influyentes, y otras con residuos estandarizados muy elevados, por encima de 6.

Los paquetes mosaic y mosaicModel tienen herramientas útiles para visualizar deprisa el ajuste de un modelo. Hay alternativas más sofisticadas que veremos como visreg, LinRegInteractive o ggeffects:

• modelPlot(modelo): Gráfico de los puntos y del ajuste. Sólo para modelos lineales y glm.
• plot_mod(): Gráfico del ajuste del modelo. Se le pueden añadir capas ggplot con +.

Gráficos de modelo con mosaic y mosaicModel

plotModel(msc8)

Gráficos de modelo con mosaic y mosaicModel

mod_plot(msc8)

Gráficos de modelo con mosaic y mosaicModel

mod_plot(msc8) + 
  geom_point(aes(y=taq1),data=film90,color="orange",alpha=0.4,size=0.4) +
  labs(title="Modelo de splines cúbicos `msc8`") + theme_dark(20)

mloess1=loess(taq1~taq0,data=film90,span=0.1)
mloess2=loess(taq1~taq0,data=film90,span=0.5)
tibble(taq0=seq(4.4,18,by=0.2)) %>% 
  mutate(mloess1=predict(mloess1,.),
         mloess2=predict(mloess2,.)) %>% 
  ggplot(aes(x=taq0))+
  geom_point(aes(y=taq1),color="grey80",size=0.4,alpha=0.4,data=film90)+
  geom_line(aes(y=mloess1,color="mloess1"))+
  geom_line(aes(y=mloess2,color="mloess2"))+
  labs(title="Modelos loess con `span` 0.1 y 0.50",color="Modelo")+
  theme_bw(20)

library(np)
mkernel = npregbw(taq1~taq0,data=film90 %>% sample_n(1000))

## 
Multistart 1 of 1 |
Multistart 1 of 1 |
Multistart 1 of 1 |
Multistart 1 of 1 /
Multistart 1 of 1 |
Multistart 1 of 1 |
                   

plot(mkernel,plot.errors.method="asymptotic",
     plot.errors.style="band")

library(rpart); library(randomForest)
mrpart=rpart(taq1~taq0, data=film90)
mod_plot(mrpart) + geom_point(aes(y=taq1),data=film90,size=.4,alpha=.4,color="orange") + 
  labs(title="Modelo `rpart`") +   theme_bw(20)

mrf=randomForest(taq1~taq0, data=film90)
mod_plot(mrf) + geom_point(aes(y=taq1),data=film90,size=.4,alpha=.4,color="orange") + 
  labs(title="Modelo `randomForest`") +   theme_bw(20)

mod_plot(mp3,bootstrap=10)+labs(title="mp3")

mod_plot(msc5,bootstrap=10)+labs(title="msc5")

mod_plot(msc8,bootstrap=10)+labs(title="msc8")

mod_plot(mrpart,bootstrap=10)+labs(title="mrpart")

Regresión con más de una variable

mod0=lm(taq1~taq0+lnosc+dist, data=film90)
summary(mod0)

## 
## Call:
## lm(formula = taq1 ~ taq0 + lnosc + dist, data = film90)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7037 -0.6093 -0.0370  0.5604  6.6343 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.31064    0.15243  -34.84   <2e-16 ***
## taq0         1.88574    0.01837  102.67   <2e-16 ***
## lnosc       -1.06824    0.01715  -62.30   <2e-16 ***
## distMajor    0.60189    0.04691   12.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.101 on 4027 degrees of freedom
## Multiple R-squared:  0.8856, Adjusted R-squared:  0.8855 
## F-statistic: 1.039e+04 on 3 and 4027 DF,  p-value: < 2.2e-16

autoplot(mod0,size=.3,alpha=0.4)

mod_plot(mod0)

mod1=lm(taq1~(ns(taq0, 5)+poly(lnosc,2)+dist)^2, data=film90)
summary(mod1)

## 
## Call:
## lm(formula = taq1 ~ (ns(taq0, 5) + poly(lnosc, 2) + dist)^2, 
##     data = film90)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6952 -0.5563 -0.0497  0.4834  6.7567 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                     22.8016    15.4458   1.476   0.1400  
## ns(taq0, 5)1                   -11.2774    15.3548  -0.734   0.4627  
## ns(taq0, 5)2                    -6.8699    15.5001  -0.443   0.6576  
## ns(taq0, 5)3                     1.6300    10.0970   0.161   0.8718  
## ns(taq0, 5)4                   -18.3129    33.6021  -0.545   0.5858  
## ns(taq0, 5)5                     4.6894    22.3991   0.209   0.8342  
## poly(lnosc, 2)1               1177.3814  1056.3393   1.115   0.2651  
## poly(lnosc, 2)2                242.3497   266.9729   0.908   0.3641  
## distMajor                        0.6146     0.9832   0.625   0.5319  
## ns(taq0, 5)1:poly(lnosc, 2)1 -1423.5809  1050.7854  -1.355   0.1756  
## ns(taq0, 5)2:poly(lnosc, 2)1 -1434.3901  1057.7521  -1.356   0.1752  
## ns(taq0, 5)3:poly(lnosc, 2)1  -900.7836   709.3884  -1.270   0.2042  
## ns(taq0, 5)4:poly(lnosc, 2)1 -2257.1053  2077.4965  -1.086   0.2773  
## ns(taq0, 5)5:poly(lnosc, 2)1  -481.6520  1266.7160  -0.380   0.7038  
## ns(taq0, 5)1:poly(lnosc, 2)2  -239.4557   266.6321  -0.898   0.3692  
## ns(taq0, 5)2:poly(lnosc, 2)2  -225.6316   265.8145  -0.849   0.3960  
## ns(taq0, 5)3:poly(lnosc, 2)2  -170.5628   178.4988  -0.956   0.3394  
## ns(taq0, 5)4:poly(lnosc, 2)2  -494.5523   551.4341  -0.897   0.3699  
## ns(taq0, 5)5:poly(lnosc, 2)2  -239.4729   298.6366  -0.802   0.4227  
## ns(taq0, 5)1:distMajor          -0.5207     0.9114  -0.571   0.5678  
## ns(taq0, 5)2:distMajor           0.3292     1.0916   0.302   0.7630  
## ns(taq0, 5)3:distMajor          -1.3705     0.8277  -1.656   0.0979 .
## ns(taq0, 5)4:distMajor           1.1083     2.4349   0.455   0.6490  
## ns(taq0, 5)5:distMajor           1.3627     2.1451   0.635   0.5253  
## poly(lnosc, 2)1:distMajor        0.9000    11.3919   0.079   0.9370  
## poly(lnosc, 2)2:distMajor        9.3479     4.8776   1.916   0.0554 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.039 on 4005 degrees of freedom
## Multiple R-squared:  0.8987, Adjusted R-squared:  0.8981 
## F-statistic:  1421 on 25 and 4005 DF,  p-value: < 2.2e-16

autoplot(mod1,size=.3,alpha=0.4)

mod_plot(mod1)

anova(mod0,mod1)

## Analysis of Variance Table
## 
## Model 1: taq1 ~ taq0 + lnosc + dist
## Model 2: taq1 ~ (ns(taq0, 5) + poly(lnosc, 2) + dist)^2
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   4027 4883.9                                  
## 2   4005 4323.4 22    560.55 23.603 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod_plot(mod1,bootstrap=10)

library(car)
Anova(mod1)

## Anova Table (Type II tests)
## 
## Response: taq1
##                             Sum Sq   Df   F value    Pr(>F)    
## ns(taq0, 5)                12390.0    5 2295.5192 < 2.2e-16 ***
## poly(lnosc, 2)              3120.7    2 1445.4339 < 2.2e-16 ***
## dist                          78.0    1   72.2793 < 2.2e-16 ***
## ns(taq0, 5):poly(lnosc, 2)   100.5   10    9.3135 1.995e-15 ***
## ns(taq0, 5):dist              12.3    5    2.2876   0.04356 *  
## poly(lnosc, 2):dist            5.8    2    2.6947   0.06769 .  
## Residuals                   4323.4 4005                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(Rcmdr)

library(mgcv)
modgam1=gam(taq1~s(taq0,by=dist)+s(lnosc,dist,bs="fs"),data=film90)
summary(modgam1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## taq1 ~ s(taq0, by = dist) + s(lnosc, dist, bs = "fs")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.1832     0.4771   27.63   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df     F p-value    
## s(taq0):distIndependent  5.514  6.560 942.2  <2e-16 ***
## s(taq0):distMajor        6.054  7.229 661.2  <2e-16 ***
## s(lnosc,dist)           11.893 19.000 158.5  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.898   Deviance explained = 89.8%
## GCV = 1.0892  Scale est. = 1.0826    n = 4031

plot(modgam1)

library(mgcViz)

## Loading required package: qgam

## Loading required package: rgl

## 
## Attaching package: 'mgcViz'

## The following object is masked from 'package:lattice':
## 
##     qq

## The following objects are masked from 'package:stats':
## 
##     qqline, qqnorm, qqplot

pgam1=getViz(modgam1)
plot(pgam1)

## Hit <Return> to see next plot:

## Hit <Return> to see next plot:

## Hit <Return> to see next plot:

check(pgam1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 24 iterations.
## The RMS GCV score gradient at convergence was 1.123817e-07 .
## The Hessian was positive definite.
## Model rank =  39 / 39 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                            k'   edf k-index p-value    
## s(taq0):distIndependent  9.00  5.51    0.98    0.12    
## s(taq0):distMajor        9.00  6.05    0.98    0.08 .  
## s(lnosc,dist)           20.00 11.89    0.87  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

gamlss: Generalized additive models for location, scale and shape

Los gamlss generalizan el modelo GAM en tres direcciones:

• Permiten especificar funciones lineales para todos los momentos de la distribución. En particular los cuatro primeros argumentos son cuatro posibles fórmulas para la esperanza condicionada, la varianza, la asimetría y la curtosis.

library(gamlss)
args(gamlss)

## function (formula = formula(data), sigma.formula = ~1, nu.formula = ~1, 
##     tau.formula = ~1, family = NO(), data = sys.parent(), weights = NULL, 
##     contrasts = NULL, method = RS(), start.from = NULL, mu.start = NULL, 
##     sigma.start = NULL, nu.start = NULL, tau.start = NULL, mu.fix = FALSE, 
##     sigma.fix = FALSE, nu.fix = FALSE, tau.fix = FALSE, control = gamlss.control(...), 
##     i.control = glim.control(...), ...) 
## NULL

• Es mucho más rico en la variedad de familias que permite, de modo que se puede utilizar una distribución que se adapta a las características de nuestros datos. Incluye distribuciones mixtas, truncadas, …

• Es también más rico en la variedad de términos suaves que se pueden incluir para captar las relaciones entre las variables.

La ventaja de tener un modelo que se adapta a las cararcterísticas de nuestros datos es que facilita la predicción, al proporcionar intervalos de predicción creibles, por ejemplos. Estos se pueden extraer a partir de los centiles de la distribución ajustada. En esta exposición seguiremos el siguiente texto de introducción a gamlss

Veamos el ejemplo primero con una variable utilizando smoothing splines penalizadas, pb y una distribución normal:

library(broom)
m_1 = gamlss(taq1~pb(taq0), data=film90, family=NO)

## GAMLSS-RS iteration 1: Global Deviance = 14109.58 
## GAMLSS-RS iteration 2: Global Deviance = 14109.58

augment(m_1) %>% ggplot(aes(y=taq1,x=pb.taq0.))+geom_point(size=0.4,alpha=0.4,color="grey80")+geom_line(aes(y=.fitted),color="red")+ theme_bw()

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = taq1 ~ pb(taq0), family = NO, data = film90) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.347147   0.087053   26.96   <2e-16 ***
## pb(taq0)    0.928889   0.007149  129.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33120    0.01114   29.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  4031 
## Degrees of Freedom for the fit:  12.73672
##       Residual Deg. of Freedom:  4018.263 
##                       at cycle:  2 
##  
## Global Deviance:     14109.58 
##             AIC:     14135.05 
##             SBC:     14215.32 
## ******************************************************************

## Don't know how to automatically pick scale for object of type smooth. Defaulting to continuous.

Vemos como el resumen nos proporciona información de los distintos momentos. También nos proporciona criterios estadísticos que nos pueden ayudar en la selección de modelos, en particular el AIC y el SBC o BIC. El AIC tiende a seleccionar el modelo que mejor comportamiento predictivo tiene. El BIC lleva a escoger modelos más parsimoniosos al penalizar más cuando tenemos más observaciones.

Podemos estimar y comparar con el modelo que también modeliza la desviación típica, sigma.

m_2 = gamlss(taq1~pb(taq0), ~pb(taq0),data=film90, family=NO)

## GAMLSS-RS iteration 1: Global Deviance = 12263.21 
## GAMLSS-RS iteration 2: Global Deviance = 12263.55 
## GAMLSS-RS iteration 3: Global Deviance = 12263.54 
## GAMLSS-RS iteration 4: Global Deviance = 12263.54

augment(m_2) %>% ggplot(aes(y=taq1,x=pb.taq0.))+geom_point(size=0.4,alpha=0.4,color="grey80")+geom_line(aes(y=.fitted),color="red")+ theme_bw()

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = taq1 ~ pb(taq0), sigma.formula = ~pb(taq0),  
##     family = NO, data = film90) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.343498   0.081359    28.8   <2e-16 ***
## pb(taq0)    0.929188   0.005519   168.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.905753   0.044196   43.12   <2e-16 ***
## pb(taq0)    -0.153053   0.003629  -42.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  4031 
## Degrees of Freedom for the fit:  22.82189
##       Residual Deg. of Freedom:  4008.178 
##                       at cycle:  4 
##  
## Global Deviance:     12263.54 
##             AIC:     12309.19 
##             SBC:     12453 
## ******************************************************************

## Don't know how to automatically pick scale for object of type smooth. Defaulting to continuous.

Podemos ver los gráficos de diagnóstico con plot

plot(m_2)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  0.0006979142 
##                        variance   =  1.000248 
##                coef. of skewness  =  0.5907226 
##                coef. of kurtosis  =  3.940587 
## Filliben correlation coefficient  =  0.9909749 
## ******************************************************************

Y una herramienta muy útil para poner la lupa de aumento en la distribución del término de error: los gráficos de gusano, **worm plots*, que consisten en una versión de qqplot quitando la tendencia lineal.

wp(m_2, ylim.all=3)

Vemos que la distribución dista mucho de estar bien capturada debido a la presencia de mucha probabilidad en las colas. Vamos a utilizar distribuciones más generales que permiten modelizar la asimetría o la kurtosis. La distribución Box-Cox-Cole-Green (BCCG) permite especificar tambíén la asimetría

m_3 <-gamlss(taq1~pb(taq0), ~pb(taq0), ~pb(taq0), data=film90, family=BCCG)

## GAMLSS-RS iteration 1: Global Deviance = 11888.56 
## GAMLSS-RS iteration 2: Global Deviance = 11818.61 
## GAMLSS-RS iteration 3: Global Deviance = 11811.08 
## GAMLSS-RS iteration 4: Global Deviance = 11809.64 
## GAMLSS-RS iteration 5: Global Deviance = 11809.64

wp(m_3, ylim.all=0.5)

Ajusta bastante mejor sobre todo por la cola de la izquierda. La familia Box-Cox power exponential (BCPE) permite modelizar los cuatro momentos:

m_4 <-gamlss(taq1~pb(taq0), ~pb(taq0), ~pb(taq0), ~pb(taq0),data=film90, family=BCPE,start.from=m_3)

## GAMLSS-RS iteration 1: Global Deviance = 11738.54 
## GAMLSS-RS iteration 2: Global Deviance = 11735.35 
## GAMLSS-RS iteration 3: Global Deviance = 11734.89 
## GAMLSS-RS iteration 4: Global Deviance = 11734.67 
## GAMLSS-RS iteration 5: Global Deviance = 11734.41 
## GAMLSS-RS iteration 6: Global Deviance = 11734.08 
## GAMLSS-RS iteration 7: Global Deviance = 11733.81 
## GAMLSS-RS iteration 8: Global Deviance = 11733.66 
## GAMLSS-RS iteration 9: Global Deviance = 11733.6 
## GAMLSS-RS iteration 10: Global Deviance = 11733.59 
## GAMLSS-RS iteration 11: Global Deviance = 11733.59 
## GAMLSS-RS iteration 12: Global Deviance = 11733.6 
## GAMLSS-RS iteration 13: Global Deviance = 11733.6 
## GAMLSS-RS iteration 14: Global Deviance = 11733.61 
## GAMLSS-RS iteration 15: Global Deviance = 11733.62 
## GAMLSS-RS iteration 16: Global Deviance = 11733.62 
## GAMLSS-RS iteration 17: Global Deviance = 11733.62 
## GAMLSS-RS iteration 18: Global Deviance = 11733.63 
## GAMLSS-RS iteration 19: Global Deviance = 11733.63 
## GAMLSS-RS iteration 20: Global Deviance = 11733.63

## Warning in RS(): Algorithm RS has not yet converged

# Hay problemas de convergencia. Le decimos que continue donde lo dejó. Y en una iteración más termina.
m_4 <-gamlss(taq1~pb(taq0), ~pb(taq0), ~pb(taq0), ~pb(taq0),data=film90, family=BCPE,start.from=m_4)

## GAMLSS-RS iteration 1: Global Deviance = 11733.63

wp(m_4, ylim.all=0.5)

Que nos avisa de que el algoritmo no ha conseguido converger en 20 iteraciones, por lo que reestimamos hasta conseguir convergencia.

La selección entre estos modelos en base al AIC sería:

AIC(m_1,m_2,m_3,m_4)

##           df      AIC
## m_4 44.97853 11823.59
## m_3 36.06436 11881.77
## m_2 22.82189 12309.19
## m_1 12.73672 14135.05

Sugiere que el mejor modelo es el m_4, que además, capta la distribución del término de error.

Podemos comparar los modelos m_3 y m_4 con:

fittedPlot(m_3, m_4, x=film90$taq0, line.type = TRUE)

Funciones estimadas para cada uno de los momentos en los modelos 3, con distribución BCCG, y 4, con distribución BCPE, en verde.

Que nos muestra las funciones suaves asociadas a los cuatro momentos.

En cuanto a la distribución de la variable, podemos captarla con los centiles correspondientes.

centiles.fan(m_4, xvar=film90$taq0, cent=c(3,10,25,50,75,90,97),
colors="terrain",ylab="taq1", xlab="taq0")

Percentiles 3%, 10%, 25%, 50% y simétricos del modelo gamlss preferido

Podemos también estimar por separado las distribuciones condicionales en distintos valores de la variable. Para ello se utiliza PlotSimpleGamlss del paquete gamlss.util.

library(gamlss.util)
library(colorspace)
plotSimpleGamlss(taq1,taq0, model=m_4, data=film90,
x.val=seq(6,16,2), val=5, N=1000, ylim=c(0,25),
cols=heat_hcl(100))

Distribución condicionada del ingreso en el resto de semanas en función del ingreso en la primera semana

## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3)

Podemos explorar las distribuciones que nos ofrece gamlss con esta demostración

library(gamlss.demo)
gamlss.demo()

Asimismo, podemos estimar un modelo completo que incluya a las otras variables con

m_f <- gamlss(taq1 ~ pb(taq0)+pb(lnosc) + dist,
sigma.fo = ~ pb(taq0)+pb(lnosc) + dist,
nu.fo = ~ pb(taq0)+pb(lnosc) + dist,
tau.fo = ~ pb(taq0)+pb(lnosc) + dist,
family = BCPE, data = film90, trace=FALSE)

# Vemos las características de las funciones ajustadas con `tidy` del paquete broom
tidy(m_f)

## # A tibble: 16 x 6
##    parameter term        estimate std.error statistic  p.value
##    <chr>     <chr>          <dbl>     <dbl>     <dbl>    <dbl>
##  1 mu        (Intercept)  -5.66      0.106    -53.3   0.      
##  2 mu        pb(taq0)      1.92      0.0126   152.    0.      
##  3 mu        pb(lnosc)    -1.04      0.0126   -82.2   0.      
##  4 mu        distMajor     0.211     0.0316     6.66  3.04e-11
##  5 sigma     (Intercept)  -0.0872    0.0978    -0.891 3.73e- 1
##  6 sigma     pb(taq0)     -0.225     0.0119   -18.8   5.21e-76
##  7 sigma     pb(lnosc)    -0.0125    0.0113    -1.11  2.68e- 1
##  8 sigma     distMajor     0.0617    0.0314     1.97  4.91e- 2
##  9 nu        (Intercept)  -5.56      0.826     -6.73  1.90e-11
## 10 nu        pb(taq0)      0.742     0.108      6.90  5.99e-12
## 11 nu        pb(lnosc)    -0.905     0.126     -7.18  8.50e-13
## 12 nu        distMajor    -1.06      0.329     -3.24  1.21e- 3
## 13 tau       (Intercept)   1.27      0.306      4.16  3.19e- 5
## 14 tau       pb(taq0)     -0.0710    0.0366    -1.94  5.23e- 2
## 15 tau       pb(lnosc)     0.0224    0.0339     0.661 5.09e- 1
## 16 tau       distMajor     0.142     0.0918     1.54  1.22e- 1

wp(m_f, xvar=~taq0+lnosc, ylim.worm=1) # Worm plot por tramos

## number of missing points from plot= 1  out of  835

## Warning in panel(x[id], y[id], col = col[id], pch = pch[id], ...): Some points are missed out 
## increase the y limits using ylim.worm

## number of missing points from plot= 0  out of  382

## number of missing points from plot= 0  out of  12

## number of missing points from plot= 0  out of  318

## number of missing points from plot= 0  out of  689

## number of missing points from plot= 0  out of  143

## number of missing points from plot= 2  out of  5

## Warning in panel(x[id], y[id], col = col[id], pch = pch[id], ...): Some points are missed out 
## increase the y limits using ylim.worm

## number of missing points from plot= 0  out of  146

## number of missing points from plot= 0  out of  688

## number of missing points from plot= 0  out of  174

## number of missing points from plot= 0  out of  174

## number of missing points from plot= 0  out of  834

Ajuste de la distribución en el modelo final con todas las variables por tramos.

# Gráficos de las funciones suaves
term.plot(m_f,pages=1,main="mu")

Gráficos de las funciones suaves para los momentos en el gamlss final.

term.plot(m_f,what="sigma",pages=1,main="sigma")

Gráficos de las funciones suaves para los momentos en el gamlss final.

term.plot(m_f,what="nu",pages=1,main="nu")

Gráficos de las funciones suaves para los momentos en el gamlss final.

term.plot(m_f,what="tau",pages=1,main="tau")

Gráficos de las funciones suaves para los momentos en el gamlss final.

Una vez estimado el modelo general, podemos comprobar que mejora sustancialmente respecto al AIC al mejor de los modelos que habíamos estimado hasta ahora (el gam). Es también posible realizar validación cruzada con gamlss para comprobar qué modelos funcionan mejor en la predicción.

Para un modelo tan complejo como éste, es dificil la interpretación, incluso a partir de los gráficos anteriores. Existe un paquete, distreg.vis, que permite trabajar de forma interactiva a través de una aplicación shiny. Lo podeís iniciar con

library(distreg.vis)
vis()

Sin embargo, todavía no incluye las distribuciones que hemos empleado BCCG y BCCP. Podéis seleccionar m_2. Siguiendo las instrucciones podéis crear escenarios (condicionar a valores de las variables), y observar cómo cambia la distribución condicionada. Además tienes la opción de obtener el código para replicar el análisis desde R.

Validación predictiva de modelos

La idea de la validación consiste en escoger los modelos que mejor predicen fuera de la muestra. Para saberlo, el procedimiento general consiste en dividir nuestra muestra en dos partes: la muestra de trabajo en la que estimamos, y la muestra de test en la que predecimos. De este modo podemos obtener una idea de cómo funcionan los modelos en la práctica y evitar los problemas ligados al sobreajuste.

La parte mala de validar de este modo es que perdemos una parte de los datos para la validación,y también que nada nos garantiza que la muestra de test que estemos utilizando no sea “rara”. Un método que evita ambos problemas es la validación cruzada k-veces (k-fold cross-validation). La validación cruzada se basa en dividir la muestra en k-partes, y estimar el modelo k veces dejando como muestra de test sucesivamente cada una de las partes. Podemos entonces comparar el ajuste del modelo con toda la muestra con las predicciones obtenidas por la validación cruzada. En general será peor el ajuste. La diferencia entre ambas mide el sobreajuste que, como decíamos, es en general mayor cuanto más complejo es el modelo.

Para muchas funciones de R existen equivalentes que escogen los parámetros por validación cruzada. Un ejemplo era el paquete np que comentamos. Otro es

Lógica de la validación cruzada Fuente wikipedia

Lógica de la validación cruzada Fuente wikipedia

Hay muchas funciones en R que realizan automáticamente la validación cruzada. Por ejemplo, para loess, tenemos la función loess.wrapper que evalua loess con distintos valores de span para ver cuál se comporta mejor en validación cruzada.

library(bisoreg)
mloesscv=loess.wrapper(film90$taq0, film90$taq1, span.vals = seq(0.05, 0.5, by = 0.05), folds = 10)
summary(mloesscv)

## Call:
## loess(formula = y ~ x, span = span)
## 
## Number of Observations: 4031 
## Equivalent Number of Parameters: 30.66 
## Residual Standard Error: 1.39 
## Trace of smoother matrix: 33.9  (exact)
## 
## Control settings:
##   span     :  0.1 
##   degree   :  2 
##   family   :  gaussian
##   surface  :  interpolate      cell = 0.2
##   normalize:  TRUE
##  parametric:  FALSE
## drop.square:  FALSE

En nuestro caso tenemos que, dada la gran cantidad de puntos disponibles, y lo complejo de la función subyacente, el span óptimo está en torno a 0.1. Si quisieramos podríamos estimar con más precisión utilizando span.vals en un entorno de 0.1.

Esto que hemos hecho a pequeña escala, es lo que hacen automáticamente gam y gamlss cuando estiman el grado de complejidad de las funciones suaves.

Existen distintos paquetes para ayudar en la validación cruzada. Un instrumento muy sencillo es mod_error y mod_cv del paquete mosaic. Para las familias que lo tienen definido, que incluyen los modelos lineales, podemos obtener el ECM (rmse) dentro de la muestra (o prediciendo si se especifica newdata) con mod_error y la predicción fuera de la muestra por validación cruzada k-veces (por defecto 10) con mod_cv(). Además, para evitar problemas de varianza al estimar, lo hace por defecto 5 veces (se puede controlar con ntrials)

mod_error(mp8)

##      mse 
## 1.947198

mod_cv(mp8)

##        mse model
## 1 1.954641   mp8
## 2 1.957472   mp8
## 3 1.961952   mp8
## 4 1.958341   mp8
## 5 1.957114   mp8

mod_cv(ml,mp2,mp3,mp4,mp6,mp8,mrf,mrpart,msl,msc5,msc6,msc7,msc8,mod0,mod1,ntrials=1,k=5) %>% 
  arrange(mse)

##         mse  model
## 1  1.093432   mod1
## 2  1.214031   mod0
## 3  1.951899   msc8
## 4  1.953624    mp8
## 5  1.955262   msc6
## 6  1.959387   msc7
## 7  1.965561    mp6
## 8  1.968247   msc5
## 9  2.015312    msl
## 10 2.048654    mp4
## 11 2.150606    mp3
## 12 2.252271    mp2
## 13 2.257199 mrpart
## 14 2.468782     ml
## 15 2.601622    mrf

Manejo de resultados con broom

El paquete broom, que ya hemos utilizado, tiene tres funciones genéricas:

• augment: Medidas asociadas a las observaciones.
• tidy: Medidas asociadas a los coeficientes en los modelos paramétricos. Con conf.int=TRUE proporciona intervalo de confianza al nivel de confianza conf.level (por defecto 95%)
• glance: medidas asociadas al modelo.

# Augment para distintos tipos de modelos
augment(mod0) %>% head() #lm

## # A tibble: 6 x 11
##    taq1  taq0 lnosc dist  .fitted .se.fit  .resid    .hat .sigma .cooksd
##   <dbl> <dbl> <dbl> <fct>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>   <dbl>
## 1 13.9  14.1   6.93 Inde~   13.8   0.0462  0.130  1.76e-3   1.10 6.12e-6
## 2  9.16  9.72  1.79 Inde~   11.1   0.0263 -1.95   5.70e-4   1.10 4.47e-4
## 3 14.3  14.2   7.29 Inde~   13.7   0.0484  0.569  1.93e-3   1.10 1.29e-4
## 4  8.44  8.22  1.39 Inde~    8.70  0.0333 -0.262  9.15e-4   1.10 1.29e-5
## 5 12.1   9.24  0    Inde~   12.1   0.0347 -0.0387 9.92e-4   1.10 3.06e-7
## 6 10.7   7.92  0    Inde~    9.62  0.0287  1.12   6.81e-4   1.10 1.77e-4
## # ... with 1 more variable: .std.resid <dbl>

augment(mloesscv) %>% head() #loess

## # A tibble: 6 x 5
##       y     x .fitted .se.fit .resid
##   <dbl> <dbl>   <dbl>   <dbl>  <dbl>
## 1 13.9  14.1    14.6    0.117 -0.710
## 2  9.16  9.72   12.2    0.118 -3.03 
## 3 14.3  14.2    14.8    0.113 -0.533
## 4  8.44  8.22    9.16   0.110 -0.719
## 5 12.1   9.24   11.5    0.118  0.614
## 6 10.7   7.92    8.70   0.121  2.04

augment(modgam1) %>% head() #gam

## # A tibble: 6 x 10
##    taq1  taq0 dist    lnosc .fitted .se.fit  .resid    .hat .sigma .cooksd
##   <dbl> <dbl> <fct>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <lgl>    <dbl>
## 1 13.9  14.1  Indepe~  6.93   14.3   0.121  -0.345  0.0134  NA     6.21e-5
## 2  9.16  9.72 Indepe~  1.79   11.5   0.0594 -2.38   0.00326 NA     7.04e-4
## 3 14.3  14.2  Indepe~  7.29   14.3   0.163  -0.0349 0.0245  NA     1.19e-6
## 4  8.44  8.22 Indepe~  1.39    8.66  0.0722 -0.222  0.00481 NA     9.07e-6
## 5 12.1   9.24 Indepe~  0      11.8   0.0547  0.242  0.00277 NA     6.13e-6
## 6 10.7   7.92 Indepe~  0       9.01  0.0481  1.73   0.00213 NA     2.41e-4

augment(m_4) %>% head() #gamlss

## ******************************************************************
## Family:  c("BCPE", "Box-Cox Power Exponential") 
## 
## Call:  gamlss(formula = taq1 ~ pb(taq0), sigma.formula = ~pb(taq0),  
##     nu.formula = ~pb(taq0), tau.formula = ~pb(taq0),  
##     family = BCPE, data = film90, start.from = m_4) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.40437    0.07743   18.14   <2e-16 ***
## pb(taq0)     0.97076    0.00528  183.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.145623   0.041874   3.478 0.000511 ***
## pb(taq0)    -0.220777   0.003568 -61.871  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.54502    0.43420  -1.255    0.209
## pb(taq0)     0.04146    0.04428   0.936    0.349
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.38463    0.14240   9.724  < 2e-16 ***
## pb(taq0)    -0.04542    0.01144  -3.969 7.34e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  4031 
## Degrees of Freedom for the fit:  44.97853
##       Residual Deg. of Freedom:  3986.021 
##                       at cycle:  1 
##  
## Global Deviance:     11733.63 
##             AIC:     11823.59 
##             SBC:     12107.03 
## ******************************************************************

## # A tibble: 6 x 5
##    taq1 pb.taq0. .fitted .se.fit .resid
##   <dbl>    <dbl>   <dbl>   <dbl>  <dbl>
## 1 13.9     14.1    14.5   0.0345 -1.23 
## 2  9.16     9.72   12.2   0.0656 -1.73 
## 3 14.3     14.2    14.7   0.0288 -0.933
## 4  8.44     8.22    9.14  0.0603 -0.485
## 5 12.1      9.24   11.3   0.0620  0.446
## 6 10.7      7.92    8.56  0.0587  1.17

# Tidy para modelos paramétricos o semiparamétricos
broom::tidy(mp4)

## # A tibble: 5 x 5
##   term           estimate std.error statistic  p.value
##   <chr>             <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)        13.3    0.0225     590.  0.      
## 2 poly(taq0, 4)1    181.     1.43       127.  0.      
## 3 poly(taq0, 4)2    -29.6    1.43       -20.7 8.08e-91
## 4 poly(taq0, 4)3     20.4    1.43        14.2 6.87e-45
## 5 poly(taq0, 4)4     20.5    1.43        14.4 1.35e-45

tidy(modgam1)

## # A tibble: 3 x 5
##   term                      edf ref.df statistic p.value
##   <chr>                   <dbl>  <dbl>     <dbl>   <dbl>
## 1 s(taq0):distIndependent  5.51   6.56      942.       0
## 2 s(taq0):distMajor        6.05   7.23      661.       0
## 3 s(lnosc,dist)           11.9   19         159.       0

tidy(m_4)

## # A tibble: 8 x 6
##   parameter term        estimate std.error statistic  p.value
##   <chr>     <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 mu        (Intercept)   1.40     0.0685     20.5   7.48e-89
## 2 mu        pb(taq0)      0.971    0.00471   206.    0.      
## 3 sigma     (Intercept)   0.146    0.0405      3.60  3.26e- 4
## 4 sigma     pb(taq0)     -0.221    0.00344   -64.1   0.      
## 5 nu        (Intercept)  -0.545    0.432      -1.26  2.07e- 1
## 6 nu        pb(taq0)      0.0415   0.0442      0.939 3.48e- 1
## 7 tau       (Intercept)   1.38     0.146       9.49  3.87e-21
## 8 tau       pb(taq0)     -0.0454   0.0117     -3.87  1.12e- 4

# Glance para distintos tipos de modelo
glance(mp4)

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df logLik    AIC
## *     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl>  <dbl>
## 1     0.807         0.807  1.43     4213.       0     5 -7159. 14329.
## # ... with 3 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>

glance(modgam1)

## # A tibble: 1 x 6
##      df logLik    AIC    BIC deviance df.residual
##   <dbl>  <dbl>  <dbl>  <dbl>    <dbl>       <dbl>
## 1  24.5 -5867. 11786. 11946.    4338.       4007.

glance(m_4)

## # A tibble: 1 x 6
##      df logLik    AIC    BIC deviance df.residual
##   <int>  <dbl>  <dbl>  <dbl>    <dbl>       <dbl>
## 1     0 -5867. 11824. 12107.   11734.       3986.

# Comparación de los modelos gam y gamlss
list(gam1=modgam1,gamlss1=m_1,gamlss2=m_2,gamlss3=m_3,gamlss4=m_4) %>% 
  map_df(glance,.id="Model") %>% arrange(AIC)

## # A tibble: 5 x 7
##   Model      df logLik    AIC    BIC deviance df.residual
##   <chr>   <dbl>  <dbl>  <dbl>  <dbl>    <dbl>       <dbl>
## 1 gam1     24.5 -5867. 11786. 11946.    4338.       4007.
## 2 gamlss4   0   -5867. 11824. 12107.   11734.       3986.
## 3 gamlss3   0   -5905. 11882. 12109.   11810.       3995.
## 4 gamlss2   0   -6132. 12309. 12453.   12264.       4008.
## 5 gamlss1   0   -7055. 14135. 14215.   14110.       4018.

Es lógico que sea mejor el modelo gam puesto que incluía información del resto de variables. Sería interesante comparar incorporando el resto de variables al gamlss.

Tablas comparadas con texreg

texreg es un paquete que incluye funciones para mostrar de manera sencilla tablas de modelos que se pueden incluir en artículos de investigación. Funciona casi con todos los modelos. Un ejemplo:

library(texreg)
list(gam1=modgam1,gamlss1=m_1,gamlss2=m_2,gamlss3=m_3,gamlss4=m_4) %>% htmlreg()

Al meterlo en un .Rmd hay que especificar la opción results=\asis`en el encabezamiento del *chunk*. Si queremos tablas interactivas empleamosscreenreg`

screenreg(list(ml,mp2,mod0))

## 
## =====================================================
##                 Model 1      Model 2      Model 3    
## -----------------------------------------------------
## (Intercept)        2.35 ***    13.29 ***    -5.31 ***
##                   (0.10)       (0.02)       (0.15)   
## taq0               0.93 ***                  1.89 ***
##                   (0.01)                    (0.02)   
## poly(taq0, 2)1                180.94 ***             
##                                (1.50)                
## poly(taq0, 2)2                -29.65 ***             
##                                (1.50)                
## lnosc                                       -1.07 ***
##                                             (0.02)   
## distMajor                                    0.60 ***
##                                             (0.05)   
## -----------------------------------------------------
## R^2                0.77         0.79         0.89    
## Adj. R^2           0.77         0.79         0.89    
## Num. obs.       4031         4031         4031       
## RMSE               1.57         1.50         1.10    
## =====================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

Variabilidad bootstrap con mosaic y bcaboot

mosaictiene una forma muy sencilla para realizar inferencia bootstrap. Por ejemplo, si nos interesa obtener los intervalos de confianza de los parámetros del mod0 podemos hacer

mod0bs=do(500)*lm(formula = taq1 ~ taq0 + lnosc + dist, data = resample(film90))

## Using parallel package.
##   * Set seed with set.rseed().
##   * Disable this message with options(`mosaic:parallelMessage` = FALSE)

mod0bs %>% head

##   Intercept     taq0     lnosc distMajor    sigma r.squared         F
## 1 -5.480574 1.906551 -1.083093 0.5510458 1.092290 0.8831479 10145.119
## 2 -5.279149 1.874485 -1.053040 0.6312589 1.093452 0.8855268 10383.843
## 3 -5.286818 1.890631 -1.077345 0.5756066 1.086774 0.8862589 10459.320
## 4 -5.261512 1.884019 -1.084497 0.6800503 1.094319 0.8859579 10428.175
## 5 -5.307785 1.881260 -1.065819 0.6603023 1.123564 0.8801964  9862.118
## 6 -5.346134 1.888138 -1.075851 0.6523747 1.121774 0.8802147  9863.831
##   numdf dendf .row .index
## 1     3  4027    1      1
## 2     3  4027    1      2
## 3     3  4027    1      3
## 4     3  4027    1      4
## 5     3  4027    1      5
## 6     3  4027    1      6

confint(mod0bs)

##        name        lower         upper level     method      estimate
## 1 Intercept   -5.6989615    -4.9449678  0.95 percentile    -5.3106371
## 2      taq0    1.8417811     1.9300239  0.95 percentile     1.8857356
## 3     lnosc   -1.1062116    -1.0274760  0.95 percentile    -1.0682367
## 4 distMajor    0.5082691     0.7031641  0.95 percentile     0.6018866
## 5     sigma    1.0697557     1.1319340  0.95 percentile     1.1012705
## 6 r.squared    0.8796699     0.8927855  0.95 percentile     0.8855833
## 7         F 9813.0887111 11177.7414357  0.95 percentile 10389.6342338

confint(mod0bs,method="se")

## Warning: confint: Using df=Inf.

##        name        lower         upper level method      estimate
## 1 Intercept   -5.6889488    -4.9579960  0.95 stderr    -5.3106371
## 2      taq0    1.8439679     1.9302900  0.95 stderr     1.8857356
## 3     lnosc   -1.1078467    -1.0313439  0.95 stderr    -1.0682367
## 4 distMajor    0.5000552     0.7046069  0.95 stderr     0.6018866
## 5     sigma    1.0663044     1.1321089  0.95 stderr     1.1012705
## 6 r.squared    0.8791340     0.8930040  0.95 stderr     0.8855833
## 7         F 9732.8228461 11169.1766966  0.95 stderr 10389.6342338
##   margin.of.error
## 1    3.654764e-01
## 2    4.316107e-02
## 3    3.825139e-02
## 4    1.022758e-01
## 5    3.290228e-02
## 6    6.934971e-03
## 7    7.181769e+02

var(mod0bs %>% dplyr::select(1:4)) #Matriz de varianzas-covarianzas estimada por el método bootstrap.

##              Intercept          taq0         lnosc     distMajor
## Intercept  0.034771420 -0.0040327474  0.0032107708  0.0030271412
## taq0      -0.004032747  0.0004849402 -0.0004004957 -0.0004249267
## lnosc      0.003210771 -0.0004004957  0.0003808889  0.0001268334
## distMajor  0.003027141 -0.0004249267  0.0001268334  0.0027230146

library(sandwich)
vcovBS(mod0) # Se puede obtener directamente a partir del paquete sandwich, pero sólo hace 100 réplicas (`R`) por defecto.

##              (Intercept)          taq0         lnosc     distMajor
## (Intercept)  0.032400390 -0.0036809123  2.809613e-03  2.984384e-03
## taq0        -0.003680912  0.0004346287 -3.479184e-04 -3.934319e-04
## lnosc        0.002809613 -0.0003479184  3.325880e-04  6.946027e-05
## distMajor    0.002984384 -0.0003934319  6.946027e-05  2.760402e-03

La ventaja de los intervalos bootstrap es que replican la distribución de los estimadores al remuestrar las observaciones.

Los métodos de bootstrap de mosaic son facilmente trazables, pero no incorporan todos los avances de la técnica bootstrap. Esto es lo que hace el paquete bcaboot que estima intervalos bootstrap con corrección de sesgo. Veamos un ejemplo prediciendo para taq=8 con el modelo msc5.

library(bcaboot)

## 
## Attaching package: 'bcaboot'

## The following object is masked from 'package:bootstrap':
## 
##     diabetes

Xy <- film90 %>% dplyr::select(x=taq0,y=taq1)
rfun <- function(Xy) {
    lm(y ~ ns(x,5),data=Xy) %>% predict(data.frame(x=8))
}
set.seed(1234)
result <- bcajack(x = Xy, B = 200, func = rfun, verbose = FALSE)


kable(result$lims, digits = 3) # Percentiles de la distribución bootstrap y su error estándar.

kable(result$stats, digits = 3) # Estadísticos del ajuste. El z0 mide el sesgo. La segunda fila, los errores estándar de los estimadores

Percentiles bootstrap de acuerdo al estimador naive (verdes), y corregidos de sesgo (negros)

bcaplot(result)