Introduccion a los modelos lineales generales (LMs)

Instituto de Investigaciones en Ecosistemas y Sustentabilidad, Universidad Nacional Autonoma de Mexico. Morelia, Michoacan.
fmora@iies.unam.mx

Como se mencion? en la secci?n anterior, la tabla de datos de emisiones de CO₂ fue ensamblada con el fin de establecer cu?l es la importancia relativa del tama?o de la poblaci?n y de la econom?a en la generaci?n de emisiones de CO₂. Esta pregunta parte de las premisas de que ambas variables constituyen una causalidad de las emisiones: se espera que la cantidad de CO₂ emitida por un pa?s se incremente a medida que crece su poblaci?n (posiblemente bajo el argumento de que personas adicionales implican el uso de recursos adicionales) y a medida que crece el tama?o de la actividad econ?mica (a mayor actividad de producci?n, mayor cantidad de procesos que generan emisiones). La Figura 1 sugiere una asociaci?n entre las emisiones de CO₂ y la poblaci?n, en concordancia con una de las premisas del ejercicio. ?Podemos sintetizar esta relaci?n empleando un modelo estad?stico? Por su puesto que s?. De hecho, existen diferentes tipos de modelos estad?sticos para resumir esta relaci?n. En la presente sesi?n haremos uso de una estrategia (herramienta) de modelaci?n muy com?n: la regresi?n lineal simple, la cual constituye un caso particular de los modelos lineales generales.

Figura 1. Relaci?n entre las emisiones de CO₂ y la poblaci?n. La l?nea corresponde al ajuste de un modelo de regresi?n lineal.

1 Recordando la l?nea recta

Cualquier l?nea recta puede describirse a trav?s de la siguiente funci?n matem?tica:

\[ \begin{aligned} y = f(x) = \beta_{0} + \beta_{1}x \end{aligned} \]

En esta f?rmula, la variable \(y\) est? en funci?n de la variable \(x\) a trav?s de dos par?metros: el intercepto (\(\beta_{0}\)) y la pendiente (\(\beta_{1}\)). El intercepto establece cu?l es el valor de la variable \(y\) cuando \(x = 0\). Por su parte, la pendiente establece en cu?ntas unidades se incrementa el valor de Y por cada unidad de incremento en la variable X. Dependiendo del valor espec?fico que tomen estos par?metros, la l?nea recta se trazar? en una posici?n y/o con una inclinaci?n diferente, como se ilustra en la Figura 2.

Figura 2. Ejemplos de l?neas rectas con diferentes interceptos y misma pendiente (a), o diferentes pendientes con el mismo intercepto (b). La l?nea vertical punteada indica el valor de x = 0 y las horizontales los valores del intercepto.

2 La l?nea recta aplicada a “respuestas ruidosas”: Regresi?n lineal

La l?nea recta constituye un modelo completamente determin?stico: cada valor de \(x\) predice un ?nico valor posible de \(y\) (Figura 2). Sin embargo, la mayor?a de los datos reales no se comportan de esta forma: un mismo valor de \(x\) puede estar asociado a uno o m?s valores de \(y\), la mayor?a de las veces claramente diferentes al valor de \(y\) establecido por la l?nea recta. Tal es el caso de las emisiones de CO₂: si bien las emisiones siguen un patr?n claro respecto a la poblaci?n, ?stas pueden tomar valores contrastantes y/o muy alejados de los valores esperados a partir del valor de poblaci?n correspondiente (Figura 1). Esto se debe a que las emisiones responden no solo al valor de la poblaci?n, sino de otros factores no considerados hasta el momento y cuya acci?n parece ser completamente azarosa (aleatoria).

La regresi?n lineal simple es un modelo estad?stico que permite conocer el valor espec?fico que toma una variable \(y\) al incorporar tanto el componente determin?stico (relaci?n entre \(x\) y \(y\), ambas cont?nuas), como la variaci?n aleatoria en torno a esta relaci?n. El componente determin?stico describe cu?l ser?a el valor esperado de la variable \(y\) para cada valor de la variable \(x\). En otras palabras, este componente “idealiza” el comportamiento de \(y\) en relaci?n a \(x\). Dicha idealizaci?n o valor esperado de la variable \(y\) para cada valor de la variable \(x\) se denota como \(\hat{y}\), por lo que estrictamente el modelo de la l?nea recta debe re-escribirse como:

\[ \begin{aligned} \hat{y} = f(x) = \beta_{0} + \beta_{1}x \end{aligned} \]

?C?mo saber cu?l l?nea es la que mejor representa la tendencia general de la relaci?n? En el caso particular de nuestro ejemplo, ?cu?les deben de ser los valores de los par?metros para obtener la mejor idealizaci?n posible de la relaci?n entre emisiones y poblaci?n? Existen diferentes m?todos matem?ticos para establecer cu?les deben ser los valores de los par?metros para obtener el modelo m?s adecuado, de los cuales los m?s comunes son el m?todo de m?nimos cuadrados y el de m?xima verosimilitud (ver Cap?tulo 6 en Bolker [-@Bolker2008]). Por ahora, basta con saber que las estimaciones de los par?metros de la l?nea recta en la Figura 1 son \(\beta_{0} = 0.4101\) (intercepto) y \(\beta_{1} = 0.8815\) (pendiente) y que el componente determin?stico del modelo de regresi?n lineal puede ser escrito entonces como:

\[ \begin{aligned} \log_{10}(\hat{CO_{2}}) = 0.4101 + 0.8815* \log_{10}(pob) \end{aligned} \]

Para convertir este modelo determin?stico en uno probabil?stico (y por ende estad?stico), se requiere la inclusi?n de un t?rmino que represente la variaci?n de las emisiones reales en torno a esta tendencia central o valor esperado. A ese t?rmino se le denomina error o residuo y no es m?s que la diferencia entre el valor real de \(CO_{2}\) para cada observaci?n y el valor esperado o idealizado de las emisiones (\(\hat{CO_{2}}\)) dado el valor de \(x\):

\[ \begin{aligned} error_{i} = \log_{10}(CO_{2_{i}}) - \log_{10}(\hat{CO_{2_{i}}}) \ \ \ \ \ \ \ \ i: 1 \dots n \end{aligned} \]

El sub?ndice \(i\) indica la observaci?n particular a la cual se hace referencia y por ende var?a entre 1 (la primera observaci?n) y \(n\) (la ?ltima observaci?n). Remplazando el valor de \(\hat{CO_{2_{i}}}\) de la Eq.1 en la Eq. 2 y despejando para \(CO_{2_{i}}\) se obtiene :

\[ \begin{aligned} \log_{10}(CO_{2_{i}}) &= 0.4101 + 0.8815*\log_{10}(pob_{i})\ + \ error_{i} \end{aligned} \]

?sta es la expresi?n del modelo de regresi?n lineal aplicado a la relaci?n entre las emisiones de CO₂ y la poblaci?n. Esta expresi?n permite obtener con precisi?n el valor de emisiones de \(CO_{2_{i}}\) de un pa?s cualquiera en la tabla de datos a partir de su valor de poblaci?n (ambas en escala logar?tmica de base 10) y de su error particular. El error corresponde a la porci?n aleatoria, es decir, aquella porci?n del valor de \(CO_{2}\) que no es explicable o predecible por el tama?o de la poblaci?n. Podr?an existir pa?ses que tengan exactamente la misma poblaci?n pero no las mismas emisiones de CO₂. La raz?n para dicha variaci?n es hasta este momento desconocida, y por ende, se entiende que su valor es resultado ?nicamente del azar.

Veamos la representaci?n gr?fica de esta expresi?n (Figura 3). En este caso, se resalta el dato correspondiente a los Estados Unidos de Am?rica, cuyo n?mero de registro en la tabla de datos es el 201. El punto rojo con relleno azul claro indica el valor observado de emisiones (\(CO_{2_{201}} = 3.72\), ver la tabla de datos emi). El punto rojo relleno corresponde al valor “idealizado”, es decir, el valor predicho al remplazar el valor de poblaci?n en la Eq.1. (\(\hat{CO_{2_{201}}} = 2.62\)).

Figura 3. Representaci?n gr?fica del modelo de la ecuaci?n 3. Relaci?n entre las emisiones de CO₂ y la poblaci?n. La l?nea corresponde al ajuste de un modelo de regresi?n lineal.

La poblaci?n de los Estados Unidos de Am?rica (en escala log10) es de 2.50. Empleando la Eq. 2 podemos calcular el valor predicho o esperado de emisiones de \(CO_{2}\) para un pa?s con dicha poblaci?n:

\[ \begin{aligned} \log_{10}(\hat{CO_{2_{201}}}) &= 0.4101 + 0.8815 * 2.50 \\ \log_{10}(\hat{CO_{2_{201}}}) &= 2.62 \end{aligned} \]

Las emisiones reales de Estados Unidos de Am?rica son superiores al valor esperado a partir del tama?o de su poblaci?n. Las razones por las cuales esto ocurre son hasta ahora desconocidas. La diferencia entre estos dos valores est? indicada por la l?nea roja en la Figura 3, y su valor (error) est? dado por:

\[ \begin{aligned} error_{201} &= 3.72 - 2.62\\ error_{201} &= 1.10 \end{aligned} \]

La Eq. 4. expresada para los Estados Unidos de Am?rica corresponde a:

\[ \begin{aligned} 3.72 &= 0.4101 + 0.8815 * 2.50 + 1.10 \end{aligned} \]

3 El modelo lineal general

El modelo de la Eq.4 puede ser generalizado para el caso de cualquier variable respuesta aleatoria Y y cualquier predictor x:

\[ \begin{aligned} y_{i} = \beta_{0} + \beta_{1}x_{i} + \ \epsilon_{i} \end{aligned} \]

A este modelo se le conoce como modelo lineal general. El modelo es lineal no porque la relaci?n se describa usando una l?nea recta, sino porque la variable \(y\) es una funci?n lineal de los par?metros del modelo \(\beta_{0},\beta_{1}\). Dicho de forma m?s simple, las par?metros aparecen como m?ltiplos en cada t?rmino. Los modelos lineales pueden incluir relaciones no lineales, como en el caso de un polinomio de segundo orden.

\[ \begin{aligned} y_{i} = \beta_{0} + \beta_{1}x^{2}_{i} + \ \epsilon_{i} \end{aligned} \]

En este caso la variable \(y\) se incrementa de manera no lineal respecto a \(x\). Sin embargo, el modelo es lineal en la medida en que sus par?metros \(\beta_{0},\beta_{1}\) aparecen como m?ltiplos en sus t?rminos correspondientes. Por el contrario, el modelo

\[ \begin{aligned} y_{i} = \beta_{1}x^{\beta_{2}}_{i} + \ \epsilon_{i} \end{aligned} \]

no es un modelo lineal, pues \(\beta_{2}\) aparece como potencia de \(x\) y no como m?ltiplo de ?sta.

Los modelos lineales generales tienen una serie de supuestos fundamentales:

1. Los valores \(y_{i}\) son independientes entre s?.
2. \(\epsilon\) es una variable aleatoria con distribuci?n normal con media cero.
3. La varianza de \(\epsilon\) es constante a trav?s de los valores predichos \(\hat{y}\).
4. la variable \(x\) se mide sin error.

Todos estos supuestos surgen del m?todo empleado para la estimaci?n de los par?metros del modelo y por ende, para que dichas estimaciones sean v?lidas (y en particular su error est?ndar asociado), se requiere que los supuestos se cumplan.

4 Ajuste del modelo lineal en R

Se le denomina ajuste al proceso de estimaci?n de los par?metros de un modelo estad?stico cuando ?ste se aplica a un conjunto de datos espec?fico. El proceso parte de que la formulaci?n del modelo describe de manera adecuada la relaci?n que se quiere modelar. En el caso de los datos de emisiones, vimos previamente que la asociaci?n entre emisiones y poblaci?n (ambos en escala logar?tmica de base 10) puede modelarse adecuadamente empleando una l?nea recta (Figura 1). Ajustaremos entonces el modelo lineal general definido en la Eq. 4 para modelar la relaci?n entre estas dos variables. La funci?n m?s ampliamente utilizada para ajustar modelos lineales generales en R es lm.

lm(log10(CO2.emissions.Tg) ~ log10(Population.millions), data = emi)

## 
## Call:
## lm(formula = log10(CO2.emissions.Tg) ~ log10(Population.millions), 
##     data = emi)
## 
## Coefficients:
##                (Intercept)  log10(Population.millions)  
##                     0.4101                      0.8815

La funci?n se emplea usando el operador ~, que se lee “en funci?n de”. A la izquierda del operador se coloca la variable respuesta y a la derecha la variable predictora (o variables predictoras, como veremos m?s adelante). El argumento data sirve para especificar donde se encuentran almacenadas las variables. El llamado de la funci?n genera la estimaci?n de los par?metros del modelo, llamados Coefficients. El primer par?metro, llamado Intercept, corresponde a la estimaci?n de \(\beta_{0}\), es decir, el intercepto de la regresi?n lineal. El segundo par?metro, llamado log10(Population.millions), corresponde a la estimaci?n de \(\beta_{1}\), que es la pendiente de la relaci?n entre emisiones y poblaci?n. El nombre se debe a que ese par?metro define el efecto de la poblaci?n sobre la variable respuesta. Con estos valores se defini? la l?nea de regresi?n de la Figura 1.

Existen una serie de funciones adicionales que permiten obtener m?s informaci?n sobre el modelo ajustado. Para ello es conveniente guardar el ajuste del modelo en un objeto.

mod1 <- lm(log10(CO2.emissions.Tg) ~ log10(Population.millions), data = emi)
summary(mod1)

## 
## Call:
## lm(formula = log10(CO2.emissions.Tg) ~ log10(Population.millions), 
##     data = emi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6440 -0.4687  0.1179  0.4675  1.2917 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.41011    0.05379   7.624 9.41e-13 ***
## log10(Population.millions)  0.88152    0.04422  19.936  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6433 on 202 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.663,  Adjusted R-squared:  0.6613 
## F-statistic: 397.4 on 1 and 202 DF,  p-value: < 2.2e-16

La funci?n summary sirve para obtener un resumen ampliado del ajuste del modelo: 1) genera la f?rmula empleada en el ajuste; 2) un resumen de los errores (Residuals) del modelo ; 3) un cuadro con estad?sticas relacionadas a los par?metros del modelo (Coefficients), donde se observa de nuevo su estimaci?n en la columna Estimate, el error est?ndar asociado a cada estimaci?n (Std.Error) y valores del estad?stico t (t value) y su probabilidad asociada (Pr(>|t|)); y 4) una serie de estad?sticas adicionales del modelo, tales como el coeficiente de determinaci?n R² y el estad?stico F y valor de p asociado a la prueba de hip?tesis del modelo completo.

Es posible extraer los valores estimados de los par?metros, los valores de emisi?n predichos por el modelo para cada pa?s con base en su poblaci?n (\(\hat{CO_{2_{i}}}\)), o el error asociado a cada pa?s (\(\epsilon_{i}\)).

# Para extraer la estimaci?n de los par?metros
coefficients(mod1)
# Para extraer las predicciones del modelo
predict(mod1)
# Para extraer los residuos del modelo
residuals(mod1)

5 Validaci?n del ajuste del modelo lineal

Como se mencion? previamente, todo modelo lineal tiene cuatro supuestos b?sicos: 1) independencia de datos, 2) normalidad de errores, 3) homogeneidad de varianza de los errores (homocedasticidad), y 4) variable independiente medida sin error. La funci?n plot aplicada al modelo lineal genera cuatro gr?ficos que, adem?s ayudar a comprobar el cumplimiento de los supuestos 2 y 3, permite diagnosticar otro tipo de situaciones relacionadas con el ajuste del modelo, tales como la idoneidad del modelo empleado (que la forma que describe es adecuada para el patr?n de relaci?n entre las variables) o la presencia de datos muy influyentes en el ajuste del modelo.

par(mfrow=c(2,2))
plot(mod1, col="blue", pch=21)

Figura 4. Gr?ficos diagn?sticos del modelo lineal obtenidos usando la funci?n plot.

En la Figura 4, el gr?fico en la esquina superior izquierda se denomina com?nmente gr?fico de residuos o gr?fico de errores y sirve para diagnosticar la homogeneidad de varianza, as? como la ideneidad del modelo ajustado para describir la forma de la relaci?n entre las variables. En ?l los valores predichos por el modelo (\(\hat{CO_{2}}\)) aparecen en el eje horizontal y los errores correspondientes (\(\epsilon\)) en el vertical. El gr?fico a su derecha se denomina gr?fico Q-Q y sirve principalmente para diagnosticar si los errores siguen la distribuci?n de probabilidades esperada (en este caso, normal). En ?ste se emplean los residuos estandarizados, es decir, divididos por la desviaci?n est?ndar. El gr?fico inferior izquierdo es un scale-location plot, que muestra la raiz cuadrada del valor absoluto de los residuos estandarizados en el eje vertical y los valores predichos en el horizontal. Este gr?fico tambi?n permite diagnosticar la homogeneidad de varianza. Por ?ltimo, el gr?fico inferior derecho es un gr?fico de residuos contra leverage. El leverage (algunas veces traducido como apalancamiento) se refiere a la distancia de un dato respecto al resto de datos en el eje de la variable predictora. Un dato en el eje predictor que se encuentra muy alejado de los dem?s puede llegar a influir de forma importante el resultado del ajuste. Por ello, esta gr?fica ayuda a detectar datos altamente influyentes.

?Qu? se deber?a ver en estos gr?ficos? Para darse una idea, en la Figura 5 se presentan los mismos gr?ficos diagn?sticos pero para diferentes situaciones particulares. En el escenario ideal la relaci?n entre un par de variables \(x\) y \(y\) es adecuadamente descrita por la l?nea recta asociada al modelo lineal empleado hasta el momento (gr?fica superior de la columna Ideal); . adicionalmente, los errores o residuos del modelo siguen una distribuci?n normal con varianza homog?nea. En este caso, la distribuci?n de los residuos respecto a los valores predichos por el modelo (segunda l?nea de la primera columna) es homog?nea a lo largo del rango de valores predichos, el gr?fico-qq muestra los puntos dispuestos a lo largo de una l?nea recta punteada, y la gr?fica de residuos contra leverage no se?ala ninguna tendencia en los residuos ni alguna observaci?n cuyo valor de distancia de Cook sea cercano a 0.5. este es, como su nombre lo dice, el escenario ideal.

Figura 5. Gr?ficos diagn?sticos para datos con caracter?sticas particulares. La columna Ideal corresponde a una relaci?n entre variables x-y que cumple los supuestos de normalidad y homocedasticidad de los errores, y para la que el modelo lineal describe adecuadamente la relaci?n. La columna Heterocedasticidad presenta una situaci?n con errores heteroced?sticos, cuya varianza se incrementa con el valor de x. En No linealidad se presenta el caso de una relaci?n en la que las variables x y y presentan una asociaci?n no lineal. Por ?ltimo, en la columna de Datos extremos se presenta un caso en el que dos de los datos incluidos presentan valores extremos de la variable y.

Ahora veamos que pasa cuando la situaci?n no es la ideal. En la columna Heterocedasticdad de la Figura 5 se hace referencia a una situaci?n en la cual los errores no presentan varianza homog?nea, sino que por el contrario la varianza de los residuos se incrementa a medida que los valores predichos se incrementan (gr?fico en la segunda l?nea). Asociado a ello, el gr?fico–qq muestra una desviaci?n importante respecto a la l?nea recta punteada. Este patr?n es bastante com?n al ajustar modelos lineales generales e implica una violaci?n a uno de los supuestos fundamentales de ?stos, la homocedaticidad. En estos casos, se requiere hacer alguna modificaci?n durante la modelaci?n, ya sea atrav?s de la transformaci?n de la variable respuesta, o mediante la modelaci?n de la varianza (veremos algunos casos de ello m?s adelante).

Otra situaci?n no ideal ocurre cuando el modelo lineal empleado no representa adecuadamente la relaci?n entre las variables, como se muestra en el primer gr?fico de la columna Bo linealidad (Figura 5). En este caso, la relaci?n es claramente no lineal, tal vez de tipo polinomial de segundo orden (que puede ser modelada empleando un modelo lineal), de tipo exponencial, en cuyo caso existen modelos m?s adecuados (como veremos en los modelos lineales generalizados, o con modelos no lineales). El ?ltimo caso de situaci?n no-ideal presentado ac? es el de la presencia de datos extremos, como se muestra en la ?ltima columna de la Figura 5. En este ejemplo hay dos datos cuyo valor para la variable \(y\) es muy alto o muy bajo dado su valor de \(x\). En este caso, tanto el gr?fico de residuos contra predicho como el gr?fico qq muestran claramente que estos dos datos se salen de la masa central de datos. El gr?fico de Leverage presenta claramente como la distancia de Cook para estos dos datos tiene valores cercanos o incluso superiores a 0.5, indicando que son datos altamente influyentes. En estos casos, es adecuado revisar la tabla de datos para ver si se trata de un error en la transcripci?n de la informaci?n. Si no lo es, hay que considerar si es adecuado o no remover dichos datos del an?lisis; aunque esta opci?n no es la ideal, el ajuste del modelo y las conclusiones que se sacan del mismo pueden depender dr?sticamente de su inclusi?n. Si la decisi?n es excluirlos, se debe de reportar tal situaci?n durante la presentaci?n de los resultados del ajuste del modelo.

6 Inferencia estad?stica respecto al modelo

Si el diagn?stico gr?fico de los modelos es positivo, de tal forma que se considera que el modelo es adecuado y que se cumplen los supuestos del mismo, es adecuado proceder entonces a la inferencia estad?stica. La inferencia estad?stica se refiere a la toma de decisi?n respecto al efecto hipotetizado de la(s) variable(s) predictoras sobre la respuesta. En t?rminos estad?sticos, la inferencia se refiere a llegar a conclusiones respecto a propiedades de una poblaci?n (estad?stica) a partir de una muestra. Existen al menos tres aproximaciones diferentes a la inferencia estad?stica en la actualidad: 1) las pruebas de hip?tesis basadas en estad?stica frecuentista, 2) la selecci?n de modelos empleando criterios de informaci?n, y 3) la estad?stica bayesiana. En este curso abordaremos las dos primeras.

6.1 Pruebas de hip?tesis

Las pruebas de hip?tesis consisten, en t?rminos generales, en calcular la probabilidad asociada a observar un evento espec?fico si se toma por cierta una hip?tesis particular, que generalmente es la de no existencia de efecto o relaci?n (hip?tesis nula). Si dicha probabilidad calculada es baja (< 0.05), entonces se presume que es debido la hip?tesis de partida no es cierta en realidad, por lo cual se rechaza dicha hip?tesis nula y se toma como cierta la hip?tesis alterna. En el caso de los modelos lineales existen dos tipos generales de pruebas de hip?tesis. El primero corresponde a las pruebas sobre los par?metros espec?ficos del modelo, las cuales se obtienen en la secci?n de coefficientes del resumen del modelo.

coef(summary(mod1))

##                             Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)                0.4101107 0.05379027  7.624253 9.405818e-13
## log10(Population.millions) 0.8815223 0.04421830 19.935687 1.335593e-49

En estas pruebas de hip?tesis se prueba si cada par?metro estimado es diferente de cero (hip?tesis nula) empleando para ello el estad?stico t de Student (t value) y su probabilidad asociada. (Pr(>|t|)). A estas pruebas se les conoce como pruebas t de Wald. La conclusi?n en este caso es que hay evidencia suficiente para rechaz<ar la hip?tesis nula de que tanto el intercepto como la pensiente son diferentes de cero.

El segundo tipo de prueba de hip?tesis corresponde a la prueba del modelo en su conjunto: se prueba si la capacidad explicativa del modelo ajustado es superior a la de un modelo nulo que no posee variables predictoras y en la que por ende la variaci?n en la respuesta es completamente debida al azar. El estad?stico F y su probabilidad p asociada que se muestran en el resumen del modelo presenta corresponden a dicha prueba de hip?tesis.

# Extracci?n del estad?stico F asociado al modelo completo y su valor de *p*
f <- summary(mod1)$fstatistic
f

##    value    numdf    dendf 
## 397.4316   1.0000 202.0000

pf(f[1],f[2],f[3], lower.tail = FALSE)

##        value 
## 1.335593e-49

La conclusic?n es este caso es que el modelo tiene mayor capacidad explicativa de las emisiones de CO₂ que un modelo nulo que no tiene predictores. Es posible generar el mismo resultado si se ajusta el modelo nulo (sin predictores) y se comparan el modelo original y el nulo empelando la funci?n anova

#Ajuste del modelo nulo sin predictores
mod0 <- lm(log10(CO2.emissions.Tg) ~ 1, data = emi)

# Comparaci?n de los modelos 
anova(mod1,mod0)

## Analysis of Variance Table
## 
## Model 1: log10(CO2.emissions.Tg) ~ log10(Population.millions)
## Model 2: log10(CO2.emissions.Tg) ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1    202  83.593                                  
## 2    203 248.060 -1   -164.47 397.43 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

N?tese que el valor de F obtenido y su probabilidad asociada es el mismo en ambos casos. A este tipo de prueba se le conoce como prueba de raz?n de verosimilitud y puede ser usada para contrastar dos o m?s modelos anidados. Un modelo est? anidado en otro cuando su conjunto de variables predictores es subconjunto de los predictores incluidos en el modelo m?s grande. En este caso, el modelo nulo est? anidado en el modelo original, pues sus predictores (conjunto vac?o), son un subconjunto de los predictores en el modelo original (poblaci?n). Usaremos este tipo de prueba de manera recurrente m?s adelante para modelos m?s complejos.

6.2 Selecci?n de modelos basada en criterios de informaci?n

Cuando tenemos un conjunto de explicaciones alternativas (hip?tesis), representadas en diferentes modelos estad?sticos, es necesario poder comparar el desempe?o de dichos modelos en la recreaci?n de la realidad (datos). Los criterios de informaci?n proveen estimaciones de la evidencia emp?rica relativa en favor de cada modelo, es decir, cu?l de los modelos est?, en t?rminos relativos, mejor sustentado por los datos. El criterio de informaci?n m?s ampliamente utilizado es el de Akaike (AIC), el cual est? constuido a partir de dos t?rminos, uno que mide qu? tan bien el modelo reproduce la realidad y otro que mide la complejidad del modelo:

\[ \begin{equation} AIC = -2log(\mathcal{L}(\hat{\theta})|datos) + 2K \end{equation} \]

El t?rmino \(-2log(\mathcal{L}(\hat{\theta})|datos)\) es sencillamente el logaritmo del valor de la funci?n de verosimilitud. Como lo vimos previamente, la verosimilitud es una medida relativa de la posibilidad de obtener un modelo con los par?metros especificados, dados los datos que tenemos a la mano. Un modelo que describa mejor un conjunto de datos tiene un valor m?s grande para la funci?n de verosimilitud en relaci?n a otros modelos, de tal forma que el t?rmino \(-2log(\mathcal{L}(\hat{\theta})|datos)\) se hace m?s peque?o (m?s negativo) a medida que incrementamos la calidad del modelo. Sin embargo, es conocido que entre m?s parametros tiene un modelo mejor se ajusta a los datos y por ende mayor es su valor de verosimilud. El t?rmino \(2K\) adiciona 2 unidades por cada par?metro adicional en el modelo, penalizando modelos cada vez m?s complejos. El valor de AIC es entonces un indicador del ajuste relativo de los modelos al conjunto de datos, de tal forma que modelos con AIC m?s bajos (mayor valor de verosimilitud y menor n?mero de par?metros) son considerados mejores en t?rminos relativos [@Anderson2008].

El criterio de Akaike puede ser calculado empleando la funci?n AIC. Sin embargo, en la pr?ctica es muy com?n que el conjunto de datos analizado sea peque?o (por ejemplo, menos de 50 observaciones diferentes). Cuando la muestra es peque?a el valor de AIC tiene un sesgo que favorece a modelos con m?s par?metros, por lo que com?nmente se utiliza una formulaci?n del AIC corregida para muestras peque?as, AICc:

\[ AIC = -2log(\mathcal{L}(\hat{\theta})|datos) + 2K\Big(\frac{n}{n-k-1}\Big) \]

El AICc puede ser calculado empleando varias funciones en R. En adelante emplearemos la funci?n AICc de la librer?a MuMIn. Veamos entonces cu?les son los valores de AIC para los modelos ajustados a los datos de emisiones de CO₂:

library(MuMIn)
AICc(mod0,mod1)

##      df     AICc
## mod0  2 622.8789
## mod1  3 403.0455

Claramente el valor de AICc para el modelo que incluye como predictor a la poblaci?n es inferior al modelo nulo que no incluye predictores. En este caso, existe mucha m?s evidencia relativa en favor del modelo alternativo que del modelo nulo. Otra funci?n ?til para comparar dos o m?s modelos empleando el AIcc es la funci?n model.sel:

model.sel(mod0,mod1)

## Model selection table 
##       (Int) l10(Ppl.mll) df   logLik  AICc  delta weight
## mod1 0.4101       0.8815  3 -198.463 403.0   0.00      1
## mod0 0.9964               2 -309.410 622.9 219.83      0
## Models ranked by AICc(x)

Esta funci?n despliega una tabla que presenta cada modelo incluido en una fila diferente y en cada columna informaci?n sobre los mismos. Se presentan los valores de los par?metros estimados (columnas (Int) y l10(Ppl.mll)), los grados de libertad asociados a cada modelo, su valor de verosimilitud (logLik), el valor del criterio de Akaike corregido para muestras peque?as (AICc), la diferencia en el valor de (AICc) entre cada modelo y el modelo con el AICc m?s bajo (delta) y el peso de evidencia relativo asociado a cada modelo (weight). Estas dos ?ltimas columnas son importantes para hacer la inferencia respecto a los modelos.

El \(\Delta AICc\) indica qu? tanto apoyo tiene un modelo en relaci?n al modelo con mayor evidencia a su favor. Entre m?s cercano sea este valor a cero, quiere decir que dicho modelo se parece mucho al mejor modelo y que existe poca evidencia que favorezca a uno de los dos en particular. Si bien se ha usado como regla pr?ctica que cualquier modelo con \(\Delta AIC > 2\) puede ser descartado como un modelo interesante o informativo, no existe en realidad un punto de corte para descartar modelos. La posibilidad de descartar un modelo debe ser evaluada a la luz de otros argumentos [@Anderson2008]. Por su parte, la columna de pesos de evidencia \(w\) representa una medida de la probabilidad relativa de los modelos. El modelo con menor valor de AIC obtiene el valor m?s alto de \(w\). Los valores de \(w\) suman en total 1, por lo que proveen evidencia de la probabilidad relativa de cada modelo de representar el proceso que gener? los datos.

7 Visualizaci?n del ajuste del modelo

Una vez se ha tomado una decisi?n respecto a la pertinencia de un modelo para representar la relaci?n entre las variables es ideal visualizar el ajuste del modelo, ya sea como parte del reporte de los resultados, o simplemente para confirmar que el modelo constituye una representaci?n adecuada de los datos y as? evaluar si es suficiente o si se requiere emplear otras aproximaciones para la modelaci?n. En el caso de la regresi?n lineal simple, la visualizaci?n del ajuste es relativamente sencilla. En la secci?n anterior vimos como aplicando un geom_smooth en una gr?fica de ggplot podemos obtener r?pidamente una representaci?n del modelo de regresi?n:

ggplot(emi, aes(log10(Population.millions), log10(CO2.emissions.Tg))) + 
  geom_point() +
  geom_smooth(method = "lm")

La opci?n “lm” en el argumento method permite especificar que el modelo que se quiere representar es un modelo lineal, en este caso una regresi?n lineal simple. La l?nea azul representa la estimaci?n derivada de la porci?n determin?stica del modelo, mientras que el ?rea gris a su alrededor el intervalo de confianza del 95% para dicha estimaci?n. Si bien esta opci?n es muy pr?ctica, tiene restricciones en cuanto al tipo y detalles de los modelos que se pueden desplegar (ver la ayuda de la funci?n, ?geom_smooth). Una forma m?s general de graficar el ajuste de un modelo es mediante la predicci?n de nuevos valores a partir del modelo empleando la funci?n predict.

# Generaci?n de una tabla de datos con nuevos valores de poblaci?n y PIB
nuevos <- data.frame(Population.millions = seq(from = min(emi$Population.millions, na.rm = TRUE), to = max(emi$Population.millions, na.rm = TRUE), length.out = 1000),
                     GDP.million.US.dollar = seq(from = min(emi$GDP.million.US.dollar, na.rm = TRUE), to = max(emi$GDP.million.US.dollar, na.rm = TRUE), length.out = 1000))

# Generaci?n de las estimaciones y su intervalo de confianza
predichos <- predict(mod1, newdata = nuevos, interval = c("confidence"))

# Integraci?n de las variables predictoras y los valores estimados en una tabla
emi.nue <- data.frame(nuevos, predichos)

# Generaci?n de la gr?fica
ggplot(emi, aes(log10(Population.millions), log10(CO2.emissions.Tg))) + 
  geom_point() +
  geom_line(data=emi.nue, aes(log10(Population.millions), fit), col="blue3", size=1) +
  geom_line(data=emi.nue, aes(log10(Population.millions), lwr), col="blue1", size=1, linetype=2) +
  geom_line(data=emi.nue, aes(log10(Population.millions), upr), col="blue1", size=1, linetype=2)

Una ?ltima opci?n muy pr?ctica, particularmente para modelos m?s complejos es la funci?n visregen la biblioteca del mismo nombre:

library(visreg)
visreg(mod1, xtrans = log)

El uso del argumento xtrans es necesario en este caso porque la funci?n visreg grafica por descarte el predictor en su escala original, no transformada.

8 ?Qu? reportar?

En el caso de un modelo de regresi?n lineal simple, posiblemente lo m?s sencillo es presentar una gr?fica que muestre el ajuste del modelo y que incluya informaci?n adicional, como la ecuaci?n del modelo y el valor del coeficiente de determinaci?n \(R^2\). Esto permite adem?s evaluar la calidad del ajuste, al ver el comportamiento del modelo en relaci?n ala nube de puntos. Si el objetivo es tan solo mencionar que existe el efecto (por ejemplo, cuando el an?lisis no es central en los resultados), se puede incluir directamente en el texto un resumen del procedimiento de inferencia (valores del estad?stico de prueba y probabilidad asociados, o la probabilidad \(w\) asociada al modelo si se emple? selecci?n de modelos por criterios de informaci?n). En cualquier caso, siempre es altamente recomendable reportar si el modelo cumpli? con los supuestos (en primera instancia s?, o de lo contrario ser?a un error reportarlo), o mejor a?n, presentar como mqaterial suplementario dicha evidencia. Esto con el fin de mostrar que el ajuste del modelo es adecuado y que las inferencias que de ?l se derivan son correctas [@Zuur2010]

8.1 Ejercicios

1. ?Existe relaci?n entre las emisiones de CO₂ y el PIB? Ajuste un modelo de regresi?n lineal para la evaluar dicha relaci?n siguiendo la serie de pasos desarrollados en esta secci?n.

2. La instalaci?n b?sica de R trae integradas algunas tablas de datos, entre las que se encuentra airquality. Dicha tabla contiene informaci?n de cinco variables ambientales medidas en la ciudad de Nueva York. El objetivo original era probar como var?a la concentraci?n de ozono en el aire (gas contaminante cuando se presenta a nivel del suelo) en relaci?n a las condiciones de radiaci?n solar, velocidad del viento y la temperatura. Ajuste un modelo de regresi?n lineal simple a la relaci?n entre ozono y velocidad del viento. ?Es este un modelo adecuado? ?C?mo se podr?a mejorar el modelo?

9 Bibliograf?a