modelo lineal » G. R. Serrano

GEBR 5: Regresión con series temporales (IV)

Docencia, Econometría Grados No Responses »

may 162012

5.4 Autocorrelación y matriz de covarianzas consistente

Aunque la varianza de las estimaciones de los coeficientes del modelo se suele presentar como $\mathrm{V}(\hat{\boldsymbol{\beta}})=\sigma^2(X'X)^{-1}$ , se trata de un caso particular. En general, la varianza de los estimadores por MCO es:

$\mathrm{V}(\hat{\boldsymbol{\beta}})=(X'X)^{-1} X'\Omega X (X'X)^{-1},$

de donde se obtiene el caso particular si $\mathrm{V}(\boldsymbol{\varepsilon})=\Omega=\sigma^2 I_n$ . Sin embargo, existen situaciones en las que el supuesto ideal de perturbaciones esféricas no se cumple.

El primer caso, llamado heterocedasticidad, ocurre cuando las perturbaciones no tienen idéntica varianza, es decir V(ε_i)=σ²_i ≠ σ². En esta situación Ω es una matriz diagonal, pero los elementos no nulos son diferentes. El segundo caso, que denominamos autocorrelación, consiste en que las perturbaciones no son independientes (Cov(ε_t,ε _s) ≠ 0) y por lo tanto, los elementos de fuera de la diagonal de la matriz Ω son distintos de cero. En cualquiera de los casos mencionados, la estimación de $\mathrm{V}(\hat{\boldsymbol{\beta}})$ por mínimos cuadrados produce un estimador inconsistente, que invalida los contrastes realizados sobre los coeficientes estimados.

Se ha propuesto una familia general de estimadores consistentes de la matriz de covarianzas, de manera que distintas elecciones para estimar X’ΩX producen estimaciones consistentes ante la presencia de heterocedasticidad y/o autocorrelación, siendo pioneros el de White para heterocedasticidad y el de Newey-West para el caso de heterocedasticidad y autocorrelación.

La presencia de autocorrelación, además, si se incluyen retardos de la variable explicada como explicativa, puede dar lugar a estimaciones inconsistentes del vector de parámetros. Para detectar la presencia de autocorrelación en el modelo lineal general, además de las funciones de autocorrelación simple y parcial de los residuos del modelo, se puede utilizar el contraste de Breusch-Godfrey, donde hay que especificar el orden de la autocorrelación de las perturbaciones que se desea contrastar o el más genérico contraste de Ljung-Box, cuya hipótesis nula es que la serie (de residuos en este caso) sigue un proceso de ruido blanco. Volviendo a los modelos de consumo de tabaco:

> bgtest(tabm1, 1)

	Breusch-Godfrey test for serial correlation of order up to 1

data:  tabm1
LM test = 8.8763, df = 1, p-value = 0.002889

> Box.test(residuals(tabm1), 6,  type="Ljung-Box")

	Box-Ljung test

data:  residuals(tabm1)
X-squared = 14.3463, df = 6, p-value = 0.026

> bgtest(tabm2, 1)

	Breusch-Godfrey test for serial correlation of order up to 1

data:  tabm2
LM test = 2.4486, df = 1, p-value = 0.1176

> Box.test(residuals(tabm2), 6, type="Ljung-Box")

	Box-Ljung test

data:  residuals(tabm2)
X-squared = 4.9471, df = 6, p-value = 0.5506

En el modelo tabm1 existe evidencia de autocorrelación de primer orden de acuerdo con el contraste de Breusch-Godfrey y se rechaza que los residuos sean ruido blanco según el contraste de Ljung-Box, por lo que los estadísticos t de significación individual (o cualquier otro contraste) no son válidos. Sin embargo, no aparecen problemas de autocorrelación en el modelo tabm2.

Para estimar la matriz de varianzas covarianzas consistente se puede usar la función vcovHAC() (sería vcovHC() si el problema fuese sólo de heterocedasticidad). Como existen distintos estimadores consistentes, la función vcovHAC admite un argumento type para seleccionar un estimador concreto o alternativamente, utilizar la función NeweyWest() que devuelve el estimador de Newey-West de $\mathrm{V}(\hat{\boldsymbol{\beta}})$ . Por su parte, la función coeftest() realiza los contrastes t habituales admitiendo diferentes estimaciones de la matriz de covarianzas de los coeficientes estimados:

> library(sandwich)
> coeftest(tabm1) ## Estimador de la varianza MCO

t test of coefficients:

                     Estimate Std. Error t value Pr(>|t|)
(Intercept)         0.0267053  0.0074777  3.5713 0.001031 **
d(log(publicidad)) -0.2627841  0.0872122 -3.0132 0.004712 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> coeftest(tabm1, vcov.=vcovHAC(tabm1))

t test of coefficients:

                    Estimate Std. Error t value Pr(>|t|)
(Intercept)         0.026705   0.012354  2.1617  0.03737 *
d(log(publicidad)) -0.262784   0.126461 -2.0780  0.04490 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> coeftest(tabm1, vcov.=NeweyWest(tabm1, lag=2))

t test of coefficients:

                    Estimate Std. Error t value Pr(>|t|)
(Intercept)         0.026705   0.011400  2.3426  0.02479 *
d(log(publicidad)) -0.262784   0.107444 -2.4458  0.01947 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

La inconsistencia de la matriz de covarianzas no sólo invalida los contrastes t, también cualquier contraste de hipótesis lineales presentados en la Sección 4.2 utilizando la forma de Wald, no así la forma basada en sumas residuales, que no se ven alteradas. En su momento vimos la función linearHypothesis(), que también admite un argumento vcov.= para especificar una función que calcule la matriz de covarianzas consistente o bien la propia matriz de covarianzas.

GEBR 5: Regresión con series temporales (III)

Docencia, Econometría Grados No Responses »

may 082012

5.3 Modelo lineal con series temporales

En las secciones anteriores se ha visto como cargar y crear un objeto de clase zooreg para tratar datos de series temporales. A continuación, se presentan algunas funciones específicas para la estimación del modelo lineal general en incluyendo elementos habituales para esta clase de datos como retardos y diferencias. Tomando el conjunto de datos sobre consumo de tabaco y publicidad, cuya descripción se puede encontrar al final, se cargan los datos y se crea un objeto zooreg excluyendo la primera columna cargada (que contiene el año):

> library(zoo)
> tabd  <- read.csv("http://grserrano.es/datos/datos_tabaco.csv")
> tabdz <- zooreg(tabd[, -1], start=1930, frequency=1)

La falta de estacionariedad (media no constante) de las variables del ejemplo es evidente a la vista del gráfico en logaritmos que se presentó en la entrada sobre gráficos de series temporales, por lo que especificaremos los modelos en tasas de variación logarítmica. Para justificar formalmente tomar diferencias, se pueden utilizar los contrastes de Dickey-Fuller (adf.test(), que se ilustra un poco más adelante) o el de Phillips-Perron (pp.test()) aplicados a las series en logaritmos.

Para la estimación, el paquete dynlm estima modelos lineales por mínimos cuadrados, pero facilita la especificación de regresiones con retardos y diferencias en las variables. No insistiré aquí en la importancia de no estimar regresiones con variables que comparten tendencias, puesto que casi siempre serán regresiones espurias.

En el siguiente fragmento se estima un modelo en tasas de variación (diferencias de logartimos, usando la función d()) explicando el consumo de tabaco en función del gasto en publicidad utilizando datos desde 1940 (argumento data con la función window()) que se muestran en el Gráfico 1.

Gráico 1: Evolución del consumo de tabaco y el gasto en publicidad en tasas logarítmicas de variación interanual.

> library(dynlm)
> tabdz1 <- diff(log(window(tabdz[, c("consumo", "publicidad")], start=1940)))
> tmp  <- data.frame(date=as.Date(paste(1941:1978, "-6-30", sep="")),
+                    consumo=tabdz1[, "consumo"],
+                    publicidad=tabdz1[, "publicidad"])
> tmpl <- melt(tmp, id.vars="date")
> qplot(date, value*100, data=tmpl, geom=c("point", "line"),
+       ylab="T.V. Interanual (%)", xlab="Año", colour=variable)
> tabm1 <- dynlm(d(log(consumo))~d(log(publicidad)),
+                data=window(tabdz, start=1940))
> summary(tabm1)

Time series regression with "zooreg" data:
Start = 1941, End = 1978

Call:
dynlm(formula = d(log(consumo)) ~ d(log(publicidad)), data = window(tabdz,
    start = 1940))

Residuals:
      Min        1Q    Median        3Q       Max
-0.077950 -0.032970 -0.004458  0.012886  0.098907

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)
(Intercept)         0.026705   0.007478   3.571  0.00103 **
d(log(publicidad)) -0.262784   0.087212  -3.013  0.00471 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04285 on 36 degrees of freedom
Multiple R-squared: 0.2014,	Adjusted R-squared: 0.1792
F-statistic: 9.079 on 1 and 36 DF,  p-value: 0.004712

Como se aprecia en la estimación, el consumo está inversamente relacionado con el gasto publicitario, resultado anómalo sólo explicable por la progresiva concienciación de los perjuicios del tabaco para la
salud.

Para introducir un retardo de la tasa de variación del consumo de tabaco como variable explicativa se usa la función L() (dentro de la fórmula de dynlm()):

> tabm2 <- dynlm(d(log(consumo))~L(d(log(consumo))) + d(log(publicidad)),
+                data=window(tabdz, start=1940))
> summary(tabm2)

Time series regression with "zooreg" data:
Start = 1942, End = 1978

Call:
dynlm(formula = d(log(consumo)) ~ L(d(log(consumo))) + d(log(publicidad)),
    data = window(tabdz, start = 1940))

Residuals:
     Min       1Q   Median       3Q      Max
-0.07860 -0.02460 -0.00015  0.02212  0.08261

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)
(Intercept)         0.011739   0.007496   1.566  0.12659
L(d(log(consumo)))  0.430304   0.138768   3.101  0.00386 **
d(log(publicidad)) -0.136513   0.081429  -1.676  0.10282
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03638 on 34 degrees of freedom
Multiple R-squared: 0.3658,	Adjusted R-squared: 0.3285
F-statistic: 9.806 on 2 and 34 DF,  p-value: 0.0004341

Ahora es la tasa de variación del consumo en el año anterior lo que mejor explica la variación en el presente y el gasto publicitario ha perdido significación. No obstante, estos resultados deben tomarse con precaución hasta realizar los diagnósticos del modelo, ya que es frecuente la presencia de autocorrelación, lo que puede hacer que la estimación de las varianzas de los coeficientes resulte inconsistente (Sección 5.4).

El gráfico de residuos (Gráfico 2(a)), como siempre, es una primera herramienta para el diagnóstico de un modelo estimado. Primero hay que preparar un data.frame con las fechas y los residuos tipificados para usar ggplot2:

> tmp <- data.frame(date=as.Date(paste(1942:1978, "-6-30", sep="")),
+                   resid=scale(residuals(tabm2)),
+                   precio=subset(tabd, subset=(year>=1942), select="precio"))
> g6 <- qplot(date, resid, data=tmp, color=precio, geom=c("point", "line"),
+             ylab="Residuos tipificados", xlab="Año")

A la vista de las funciones de autocorrelación simple y parcial del Gráfico 2(b) (calculadas con acf() y pacf()) no se sospecha la presencia de autocorrelación en las perturbaciones del modelo, aunque sí se observa cierta variabilidad en la dispersión en el gráfico de residuos. Veremos contrastes formales en la Sección 5.4. Sobre el código del ejemplo, en la práctica no hay que programar el gráfico de las funciones de autocorrelación cada vez, sino que habría que empaquetarlo en una función que estaría siempre disponible (véase el Apéndice C).

Gráfico de residuos de tabm2

Funciones de autocorrelación simple y parcial

Gráfico 2: Residuos del modelo tabm2.

Para especificar más retardos de una variable se puede usar el segundo argumento de la función de retardo L() dentro de dynlm(). A continuación se estima un modelo autorregresivo de orden 2 para el consumo de tabaco:

> tabm2.ar <- dynlm(d(log(consumo))~L(d(log(consumo)), 1:2),
+                data=window(tabdz, start=1940))
> summary(tabm2.ar)

Time series regression with "zooreg" data:
Start = 1943, End = 1978

Call:
dynlm(formula = d(log(consumo)) ~ L(d(log(consumo)), 1:2), data = window(tabdz,
    start = 1940))

Residuals:
      Min        1Q    Median        3Q       Max
-0.067064 -0.018912  0.002138  0.019484  0.089672

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)
(Intercept)              0.003018   0.006310   0.478   0.6356
L(d(log(consumo)), 1:2)1 0.275740   0.156556   1.761   0.0875 .
L(d(log(consumo)), 1:2)2 0.213038   0.145727   1.462   0.1532
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03445 on 33 degrees of freedom
Multiple R-squared: 0.2615,	Adjusted R-squared: 0.2167
F-statistic: 5.843 on 2 and 33 DF,  p-value: 0.006726

Datos de consumo de tabaco y publicidad

El archivo datos_tabaco.csv contiene datos anuales desde 1930 a 1978 para EE.UU. de las siguientes variables:

consumo de cigarrillos per cápita,
renta per cápita,
precio medio de la cajetilla de tabaco en términos reales (es decir, deflactado por un índice de precios) en centavos,
un indicador del gasto en publicidad del sector tabaquero (una variable stock que se construye a partir del gasto corriente en publicidad de años anteriores, aplicando una cierta tasa de depreciación).

Además, se conoce la siguiente información:

En 1953, la American Cancer Society y el British Medical Research Council emitieron sendos informes que mostraban cómo las tasas de mortalidad eran más elevadas entre fumadores.
En 1964, el US General Surgeon Report concluía que fumar cigarrillos estaba ligado al cáncer de laringe.
En 1970 se aprobó la legislación que obligaba a las compañías a mostrar en las cajetillas mensajes de advertencia sobre los peligros del tabaco.

GEBR 6.3: Error de especificación y no linealidades

Docencia, Econometría Grados No Responses »

may 082012

Uno de los supuestos más difíciles de verificar en la práctica es que el modelo lineal está correctamente especificado y que tanto la forma funcional es correcta como que no falta ninguna variable relevante. Se ha visto que si disponemos de datos, es posible contrastar la relevancia de cada variable explicativa en el modelo, sin embargo, el problema real se produce cuando no disponemos de información (no observamos) una variable posiblebemente relevante o cuando la especificación del modelo debería incluir transformaciones no lineales de las variables.

El contraste RESET (Regression Equation Specification Error Test, contraste de especificación del error del modelo de regresión) trata de arrojar alguna luz sobre los tipos de error de especificación mencionados. El supuesto es que si faltan términos no lineales o variables no observadas, eso quedaría reflejado en que ciertas no linealidades de los valores ajustados del modelo resultan significativas en la explicación de la variable dependiente.

La versión de valores ajustados del contraste RESET requiere estimar por MCO el modelo auxiliar:

$y_i = \beta_1 + \beta_2 x_{i2} + \ldots + \beta_k x_{ik} + \delta_2 \hat{y}_i^2 + \ldots + \delta_p \hat{y}_i^p + \nu_i,$
donde $\hat{y}_i$ son los valores ajustados del modelo en y el orden p es decisión del analista, y suele ser un valor pequeño como 3. La hipótesis nula del contraste RESET es H₀:δ₂=… = δ_p=0 y para contrastarla se puede utilizar un estadístico F y puede que sea necesario utilizar la estimación de la matriz de covarianzas consistente a la heterocedasticidad. Existe también la versión de regresores del contraste RESET, donde en la regresión auxiliar se incluyen potencias y productos cruzados de las variables explicativas en lugar de los valores ajustados del modelo de partida.

El contraste RESET ha sido criticado por su escasa potencia, es decir, aún cuando la nula es falsa y son necesarios términos no lineales u otras variables en el modelo, el contraste tiende a no recharzala.

Utilizando datos de 88 viviendas unifamiliares, se intenta explicar el precio (price) en función del tamaño de la parcela (lotsize), los pies cuadrados de superficie (sqrft) y el número de dormitorios (bdrms). A continuación se estima el modelo por MCO, se realiza el contraste RESET con p=3 (argumento power=2:3 donde deben especificarse todas las potencias hasta p) y versión de valores ajustados (type="fitted"):

> library(lmtest)
> hprice   <- read.csv("hprice1.csv")
> hprice.1 <- lm(price ~ lotsize + sqrft + bdrms, data=hprice)
> summary(hprice.1)

Call:
lm(formula = price ~ lotsize + sqrft + bdrms, data = hprice)

Residuals:
     Min       1Q   Median       3Q      Max
-120.026  -38.530   -6.555   32.323  209.376

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.177e+01  2.948e+01  -0.739  0.46221
lotsize      2.068e-03  6.421e-04   3.220  0.00182 **
sqrft        1.228e-01  1.324e-02   9.275 1.66e-14 ***
bdrms        1.385e+01  9.010e+00   1.537  0.12795
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 59.83 on 84 degrees of freedom
Multiple R-squared: 0.6724,	Adjusted R-squared: 0.6607
F-statistic: 57.46 on 3 and 84 DF,  p-value: < 2.2e-16

> resettest(hprice.1, power=2:3, type="fitted")

	RESET test

data:  hprice.1
RESET = 4.6682, df1 = 2, df2 = 82, p-value = 0.01202

Los gráficos de diagnósitico de residuos, cuya creación con plot.lm() se presentó en la sección sobre atípicos e influyentes, se muestran a continuación:

Gráfico de residuos

Matriz sombrero, residuos y distancia de Cook

Gráfico 1: Gráficos de diagnósitico del modelo hprice.1.

Como se puede apreciar, aparece una observación influyente, la vivienda 77 que podría estar alterando los resultados. Así, primero se elimina la observación problemática, se reestima el modelo y se realiza de nuevo el contraste:

> hprice   <- hprice[-77,]
> hprice.1 <- lm(price~lotsize+sqrft+bdrms, data=hprice)
> summary(hprice.1)

Call:
lm(formula = price ~ lotsize + sqrft + bdrms, data = hprice)

Residuals:
     Min       1Q   Median       3Q      Max
-115.072  -30.340   -0.739   31.913  211.138

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -28.301170  25.149383  -1.125  0.26370
lotsize       0.009195   0.001363   6.748 1.88e-09 ***
sqrft         0.086495   0.012949   6.680 2.55e-09 ***
bdrms        20.481920   7.767117   2.637  0.00998 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 51 on 83 degrees of freedom
Multiple R-squared: 0.7646,	Adjusted R-squared: 0.7561
F-statistic: 89.89 on 3 and 83 DF,  p-value: < 2.2e-16

> resettest(hprice.1, power=2:3, type="fitted")

	RESET test

data:  hprice.1
RESET = 4.2617, df1 = 2, df2 = 81, p-value = 0.01739

La estimación sin la vivienda 77 es claramente distinta de la obtenida con todas las observaciones. Sin embargo, el contraste RESET sigue rechazando la hipótesis nula de especificación correcta con un 5% de significación. Como alternativa, se propone una especificación en logaritmos (excepto el número de dormitorios), se estima y se realiza el contraste RESET sobre los residuos del nuevo modelo:

> hprice.2 <- lm(log(price) ~ log(lotsize) + log(sqrft) + bdrms, data = hprice)
> summary(hprice.2)

Call:
lm(formula = log(price) ~ log(lotsize) + log(sqrft) + bdrms,
    data = hprice)

Residuals:
     Min       1Q   Median       3Q      Max
-0.67967 -0.08873 -0.00435  0.10679  0.66830

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  -1.33721    0.64787  -2.064   0.0421 *
log(lotsize)  0.20376    0.04553   4.475 2.41e-05 ***
log(sqrft)    0.66183    0.09611   6.886 1.02e-09 ***
bdrms         0.04140    0.02754   1.503   0.1365
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1835 on 83 degrees of freedom
Multiple R-squared: 0.6508,	Adjusted R-squared: 0.6382
F-statistic: 51.57 on 3 and 83 DF,  p-value: < 2.2e-16

> resettest(hprice.2, power=2:3, type="fitted")

	RESET test

data:  hprice.2
RESET = 2.3629, df1 = 2, df2 = 81, p-value = 0.1006

Con la transformación logaritmica, el número de dormitorios deja de ser significativo. Según el contraste RESET, no se puede rechazar la hipótesis nula de especificación correcta con un 5% de significación.

GEBR 6.2: Observaciones atípicas e influyentes

Docencia, Econometría Grados No Responses »

may 072012

Un problema que se presenta en los modelos de regresión, pero no ligado al incumplimiento de un supuesto concreto, es la presencia de observaciones extremas e influyentes. Se dice que una observación es atípica si el residuo asociado es grande. Por su parte, una observación es extrema (o potencialmente influyente o apalancada) si se encuentra apreciablemente alejada del resto de observacines de la muestra. Por su parte, se habla de observación influyente si la presencia de dicha observación en la muestra altera significativamente algún aspecto de la estimación del modelo.

Debe quedar claro que tanto la presencia de observaciones atípicas como extremas no implica, necesariamente, que influyan en la estimación, aunque deberían ser analizadas las causas de la presencia de dichas observaciones en la muestra. De no pertenecer a la población objeto de estudio, podrían eliminarse sin más. Más problemáticas resultan las observaciones influyentes, puesto que su peso en la estimación es considerable, pero si no son extremas, no se podría justificar su exclusión.

Si bien el gráfico de residuos es una primera aproximación y suficiente para la detección de atípicos, los residuos no son la herramienta adecuada para la detección de observaciones extremas e influyentes, ya que por la naturaleza de estas, podrían alterar la estimación de forma que el residuo asociado resultase pequeño.

Para detectar la presencia de observaciones potencialmente influyentes (también se les llama extremas o apalancadas) se utilizan los elementos de la diagonal principal de la matriz sombrero. Partiendo de la notación del modelo lineal general, se define H=X(X′X)⁻¹X′. El curioso nombre de la matriz H se debe a que cumple $\hat{\mathbf{y}}=H\mathbf{y}$ , es decir, le “pone el sombrero” al vector de valores de la variable dependiente y. Los elementos de la diagonal de matriz H, que se denotan h_ii, son una distancia de cada observación a la nube de puntos y por lo tanto, valores elevados (superiores a 2k/n) señalan observaciones relativamente alejadas del conjunto y podrían influir sustancialmente en la estimación.

Para detectar observaciones influyentes se han definido estadísticos que miden el efecto de cada observación en diferentes elementos de la estimación. El más utilizado es la distancia de Cook (o estadístico de Cook), que mide el cambio en el vector de coeficientes estimados $\hat{\boldsymbol{\beta}}$ al eliminar sucesivamente cada observación de la muestra. Si se denota por $\hat{\boldsymbol{\beta}}_{(i)}$ la estimación por mínimos cuadrados sin la i-ésima observación, la distancia de Cook para la observación i es:

$D_{i} = \frac{1}{k}(\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}}_{(i)})' \left [ \mathrm{V}(\hat{\boldsymbol{\beta}}) \right ]^{-1} (\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}}_{(i)}),$

donde k es la dimensión del vector $\hat{\boldsymbol{\beta}}$ . Aunque el estadístico de Cook no sigue una distribución concreta, para determinar aquellos valores grandes se suele utilizar 1 como indicador de problemas y de forma más refinada el valor crítico de una F_k,n−k.

En el contexto de las medidas de influencia, también se suelen considerar los residuos estudentizados, que se obtienen dividiendo el residuo MCO entre la desviación típica residual obtenida al eliminar la i-ésima observación de la estimación, de manera que dicha observación no tenga impacto en la estimación de σ². Los residuos estudentizados se definen:

$r_i = \frac{\hat{\varepsilon}_i}{\hat{\sigma}\sqrt{1-h_{ii}}},$

donde $\hat{\varepsilon}_i$ y $\hat{\sigma}^2$ se obtienen de la estimación del modelo por MCO y h_ii es el i-ésimo elemento de la diagonal de la matriz sombrero.

Las medidas presentadas están relacionadas con la distancia de Cook a través de las siguientes expresiones:

$D_i = \frac{r_i^2}{k} \frac{h_{ii}}{1-h_{ii}} = \frac{1}{k} \left ( \frac{\hat{\varepsilon}_i}{\hat{\sigma}\sqrt{1-h_{ii}}} \right )^2 \left ( \frac{h_{ii}}{1-h_{ii}} \right ).$

Para ilustrar el uso los estadísticos vamos a usar un conjunto de datos de un estudio de 1976 sobre la concentración media de nitrógeno en mg/litro (nitrogen) en 20 ríos del estado de Nueva York en EE.UU. Las variables explicativas son el porcentaje de tierra de la cuenca de cada río destinados a agricultura (agriculture), bosque (forests), residencial (residential) y comercial o industrial (comm.indust).

> nyrivers <- read.csv("http://grserrano.es/datos/nyrivers.csv",
+                      row.names="river" )
> m.nit <- lm(nitrogen~agriculture+forest+residential+com.indust, data=nyrivers)
> summary(m.nit)

Call:
lm(formula = nitrogen ~ agriculture + forest + residential +
    com.indust, data = nyrivers)

Residuals:
     Min       1Q   Median       3Q      Max
-0.49404 -0.13180  0.01951  0.08287  0.70480

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.722214   1.234082   1.396   0.1832
agriculture  0.005809   0.015034   0.386   0.7046
forest      -0.012968   0.013931  -0.931   0.3667
residential -0.007227   0.033830  -0.214   0.8337
com.indust   0.305028   0.163817   1.862   0.0823 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2649 on 15 degrees of freedom
Multiple R-squared: 0.7094,	Adjusted R-squared: 0.6319
F-statistic: 9.154 on 4 and 15 DF,  p-value: 0.0005963

Como se puede apreciar, ninguna de las variables es significativa al 5%, si bien el porcentaje de la tierra destinado a usos comerciales o industriales lo es al 10%. La función plot.lm(), a la se se le pasa un modelo estimado y una selección de gráfico mediante el argumento which= nos ofrece la información necesaria.

> pdf(file="infl%02d.pdf", onefile=FALSE)
> plot.lm(m.nit, which=1)
> plot.lm(m.nit, which=4)
> plot.lm(m.nit, which=5)
> dev.off()

Valores ajustados y residuos.

Distancia de Cook.

Matriz sombrero, residuos y distancia de Cook.

Gráfico 1: Gráficos de diagnósitico de atípicos e influyentes.

Como se puede apreciar en los gráficos, los ríos Fishkill y Oswegatchie presentar residuos elvados, sin embargo, no son influyentes. A la vista de las distancias de Cook y del gráfico combinado, son el Hackensack y el Neversink los influyentes, es decir, su presencia en la estimación puede estar introduciendo alteraciones importantes en los resultados. El origen de la influencia es la gran distancia que los separa del resto de observaciones, con valores de h_ii próximos a la unidad.

Los valores numéricos de las medidas de influencia presentadas más arriba (junto con otras) se pueden obtener mediante la función influence.measures(). A continuación se utiliza dicha función, que devuelve una estructura algo compleja, y se crea una variable que indica si la observación es influyente o potencialmente influyente usando la distancia de Cook y la matriz sombrero. Después se listan las filas del data.frame original que podrían resultar problemáticas. Puede mostrar las medidas relevantes para todas las observaciones listando infl.m.nit.

> infl.m.nit <- influence.measures(m.nit)
> esinf      <- apply(infl.m.nit$is.inf[, c("cook.d", "hat")], 1, any)
> cbind(nyrivers[esinf,], infl.m.nit$infmat[esinf, c("cook.d", "hat")])

           agriculture forest residential com.indust nitrogen   cook.d
Neversink            2     84         1.9       1.98     1.00 13.21955
Hackensack           3     27        29.4       3.11     1.99 65.42632
                 hat
Neversink  0.8967834
Hackensack 0.9720820

A la vista de los resultados, en ambos casos la matriz sombrero indica posibles problememas, puesto que los valores de h_ii superan con mucho la frontera de 2k/n que en este caso es 0.5. También la distancia de Cook es enorme comparada con el valor de referencia, el valor que deja a la izquierda un 95% de probabilidad en una F_{5, 20}, que es 2.71.

> nyrivers.noinf <- nyrivers[!esinf, ]
> m.nit.2 <- lm(nitrogen~agriculture+forest+residential+com.indust,
+               data=nyrivers.noinf)
> summary(m.nit.2)

Call:
lm(formula = nitrogen ~ agriculture + forest + residential +
    com.indust, data = nyrivers.noinf)

Residuals:
     Min       1Q   Median       3Q      Max
-0.36990 -0.05218 -0.00979  0.09231  0.24847

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.379765   0.774262   1.782   0.0981 .
agriculture  0.005197   0.009431   0.551   0.5909
forest      -0.010462   0.008668  -1.207   0.2489
residential  0.078339   0.084329   0.929   0.3698
com.indust   0.674612   0.437107   1.543   0.1467
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1619 on 13 degrees of freedom
Multiple R-squared: 0.8816,	Adjusted R-squared: 0.8452
F-statistic: 24.21 on 4 and 13 DF,  p-value: 6.365e-06

A pesar de haber eliminado las observaciones influyentes, el modelo todavía presenta algún problema, puesto que las variables no resultan individualmente significativas pero sí en conjunto, se puede sospechar la presencia de multicolinealidad, que se comprueba mediante los factores de inflación de varianza (usando la función vif()). Dicha sospecha está bien justificada en la naturaleza de las variables, que son el porcentaje de la cuenca de cada río destinado a cierto uso, por lo que cabe suponer que sumaran aproximadamente 100. Si eliminamos las variables residential y forest, vemos que las restantes variables sí resultan claramente significativas y con los signos esperados:

> vif(m.nit.2)

agriculture      forest residential  com.indust
  11.832121   12.335352    7.533501    7.638111

> m.nit.3 <- lm(nitrogen~agriculture+com.indust, data=nyrivers.noinf)
> summary(m.nit.3)

Call:
lm(formula = nitrogen ~ agriculture + com.indust, data = nyrivers.noinf)

Residuals:
      Min        1Q    Median        3Q       Max
-0.270138 -0.101492 -0.002511  0.075934  0.275115

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.428827   0.083290   5.149 0.000119 ***
agriculture 0.016706   0.002782   6.005 2.41e-05 ***
com.indust  1.157635   0.160485   7.213 3.01e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1641 on 15 degrees of freedom
Multiple R-squared: 0.8596,	Adjusted R-squared: 0.8409
F-statistic: 45.93 on 2 and 15 DF,  p-value: 4.022e-07

GEBR 6.1: Multicolinealidad

Docencia, Econometría Grados 2 Responses »

may 042012

La multicolinealidad se produce cuando algunas variables explicativas de un modelo lineal son combinación lineal de otras. La consecuencia es que el determinante de X′X se encuentra próximo a cero y su inversa, por lo tanto, puede ser inestable (pequeños cambios en X producen grandes cambios en la inversa). En un modelo lineal, a menudo, la primera señal de multicolinealidad es que las varianzas estimadas de los coeficientes resultan grandes, llevando a aceptar la nula de no significatividad, aunque la regresión en su conjunto resulte significativa.
Continue reading »

GEBR A.1: Cálculos matriciales y tablas estadísticas

Docencia, Econometría Grados No Responses »

may 032012

Aunque en el trabajo aplicado se dispone de datos originales, en la enseñanza de econometría es frecuente proponer problemas donde es necesario realizar ciertos cálculos a mano y posteriormente consultar tablas estadísticas. Vamos a ver cómo realizar esas tareas con R mediante un ejemplo.
Continue reading »

GEBR 5: Regresión con series temporales (II)

Docencia, Econometría Grados, Introducción al Análisis de Datos No Responses »

abr 162012

5.2 Análisis simple de series temporales

Un supuesto clásico es que una serie temporal (Y_t) se puede descomponer en ciclo-tendencia (C_t), componente estacional (S_t) y componente irregular o ruido (v_t):

$Y_t = C_t + S_t + v_t, \qquad t=1,\ldots,T.$

Si la descomposición anterior se aplica a una serie en logaritmos, estaríamos suponiendo que los componentes de la serie original son multiplicativos (Y_t = C_t ·S_t ·e^v_t) en lugar de aditivos, lo que resulta conveniente si la varianza no es constante.

Un supuesto simple (aunque casi siempre inadecuado) es que el componente tendencial es determinista y se puede ajustar a una recta (o polinomio) en función del tiempo: C_t = α+ βt + w_t. Una vez estimado el modelo anterior por mínimos cuadrados, el componente estacional puede extraerse promediando, para cada mes o trimestre, las diferencias entre el valor observado y el componente tendencial ajustado por el modelo (equivalente a añadir variables ficticias estacionales).

> ipc09  <- window(ipcz, end=c(2009, 12))       ## Datos hasta diciembre 2009
> tiempo <- time(ipc09)                         ## Variable de tiempo
> mlipc  <- lm(log(ipc09)~tiempo)               ## Tendencia lineal
> mlipc2 <- lm(log(ipc09)~tiempo+I(tiempo^2))   ## Tendencia cuadrática
> ## Componente estacional
> factor.mes.lm <- aggregate(log(ipc09)-fitted(mlipc), cycle(ipc09), mean)
> factor.mes.lm

            1             2             3             4             5
-0.0054713293 -0.0066538299 -0.0028545010  0.0053092476  0.0056314714
            6             7             8             9            10
 0.0055681477 -0.0014711077 -0.0012383893 -0.0021442988  0.0017199817
           11            12
 0.0013317933  0.0002728143

> ## Preparación para el gráfico con ggplot2, cada dato se asocia al día 15
> ipc.df      <- data.frame(log.IPC=log(ipc09), Tend.Lin.=fitted(mlipc),
+                           Tend.Cuad.=fitted(mlipc2))
> ipc.df$date <- as.Date(paste(floor(time(ipc09)), cycle(ipc09), 15, sep="-"))
> ipc.dfl     <- melt(ipc.df, id="date")
> qplot(date, value, data=ipc.dfl, geom="line", color=variable,
+       ylab="logaritmos", xlab="Año")
> ggsave("desc-ts-1.pdf", last_plot())

Como se puede apreciar (Gráfico 8), las tendencias deterministas tanto lineal como cuadrática son incapaces de ajustarse al cambio que experimenta la evolución del IPC en 2008, con un crecimiento por encima del esperado, y al contrario en 2009. Dado que estos cambios se producen de manera habitual en datos económicos, el uso de tendencias deterministas sólo es razonable en contextos muy controlados.

Gráfico 8: Extracción de tendencia determinista lineal y cuadrática.

El cálculo de previsiones, bajo el supuesto de tendencia determinista, tan sólo requiere la extrapolación de la tendencia a la que habría que añadir el componente estacional. A continuación se ilustra el cálculo de previsiones para 2010 con la tendencia lineal, y se observa el error sistemático que comete el modelo de tendencia determinista:

> tend.prev <- predict(mlipc, data.frame(tiempo=2010 + 0:11/12))
> ipc.prev  <- zooreg(as.numeric(exp(tend.prev+factor.mes.lm)),
+                     start=c(2010,1), frequency=12)
> window(merge(ipcz, ipc.prev), start=c(2010, 1), end=c(2010, 12))

            ipcz ipc.prev
2010(1)  106.680 110.0170
2010(2)  106.480 110.1650
2010(3)  107.270 110.8641
2010(4)  108.420 112.0556
2010(5)  108.660 112.3753
2010(6)  108.850 112.6525
2010(7)  108.360 112.1453
2010(8)  108.640 112.4552
2010(9)  108.710 112.6376
2010(10) 109.705 113.3598
2010(11) 110.300 113.6025
2010(12) 110.979 113.7694

La función stl() devuelve una descomposición en ciclo-tendencia, componente estacional e irregular mucho más refinada, utilizando regresiones locales para la estimación del componente tendencial. El paquete mFilter contiene distintas funciones para la extracción de tendencias, entre ellas el conocido filtro de Hodrick-Prescott. Finalmente, la función StructTS() también realiza esta clase de descomposiciones permitiendo especificar qué componentes se consideran incluidos en la variable utilizando modelos estructurales de series temporales.

A continuación se ilustra el uso de stl() para la descomposición del IPC. Como se puede apreciar en el Gráfico 9 la tendencia es adaptativa, por lo que recoge el cambio en la evolución del IPC a partir de 2008.

> lipc.stl <- stl(log(ipc09), s.window="periodic")
> lipc.stl <- as.zooreg(lipc.stl$time.series) ## Series resultantes
> tend.stl <- lipc.stl[, "trend"]
> factor.mes.stl <- aggregate(lipc.stl[, "seasonal"],
+                             cycle(lipc.stl[, "seasonal"]), mean)
> factor.mes.stl

            1             2             3             4             5
-0.0073789123 -0.0085894080 -0.0048180723  0.0036906824  0.0043579095
            6             7             8             9            10
 0.0050215278 -0.0012907873 -0.0003226822 -0.0004932039  0.0036591368
           11            12
 0.0035590064  0.0026048053

> lipc.stl.df <- data.frame(log.IPC=log(ipc09), Tend.loc.=tend.stl)
> lipc.stl.df$date <- as.Date(paste(floor(time(ipc09)),
+                                   cycle(ipc09), 15, sep="-"))
> ipc.stl.dfl      <- melt(lipc.stl.df, id="date")
> qplot(date, value, data=ipc.stl.dfl, geom="line", color=variable,
+       ylab="logaritmos", xlab="Año")
> ggsave("desc-ts-2.pdf", last_plot())

Gráfico 9: Extracción de tendencia mediante regresiones locales.

Un método más simple consiste en utilizar medias móviles para extraer el componente ciclo-tendencia de forma adaptativa. Sin embargo, la media móvil centrada tiene el inconveniente de que no se puede calcular para las últimas observaciones, por eso es recomendable usar stl(), que básicamente calcula una media móvil ponderada robusta y se puede obtener la tendencia desde la primera a la última observación. Para calcular la media móvil centrada de orden 12 de los datos de IPC se haría:

> tend.mm <- rollapply(log(ipc09), 12, mean)
> ## Comparación tendencia stl y media móvil centrada en 2009
> window(merge(ipc09, exp(tend.stl), exp(tend.mm)), start=c(2009,1))

          ipc09 exp(tend.stl) exp(tend.mm)
2009(1)  105.59      106.8233     106.7913
2009(2)  105.60      106.7330     106.7191
2009(3)  105.78      106.6429     106.6277
2009(4)  106.81      106.6611     106.5691
2009(5)  106.77      106.6794     106.5963
2009(6)  107.24      106.7155     106.6667
2009(7)  106.33      106.7516           NA
2009(8)  106.70      106.7965           NA
2009(9)  106.45      106.8415           NA
2009(10) 107.21      106.8961           NA
2009(11) 107.79      106.9508           NA
2009(12) 107.76      107.0139           NA

En cualquiera de los casos anteriores, es necesario definir alguna estrategia de extrapolación de la tendencia si se pretenden calcular previsiones. Una opción simple, aunque no óptima, sería ajustar splines cúbicos a la tendencia, por lo que la extrapolación se basaría en los últimos datos.

GEBR 4. Modelo lineal general III

Docencia, Econometría Grados No Responses »

mar 052012

4.3 Predicción

Una vez estimado el modelo (1) se puede utilizar para prever introduciendo nuevos valores de las variables explicativas en la estimación. Formalmente $\hat{y}_{N+1} = \mathbf{x}_{N+1}'\hat{\boldsymbol{\beta}}$ , que es idéntico al cálculo de los valores ajustados por la regresión para las observaciones en la muestra de estimación.

Los valores estimados por el modelo para la variable dependiente se pueden obtener con fitted() o con la función predict() pasándole como único argumento el modelo estimado. La función predict() admite, además, un data.frame como segundo argumento con nuevos valores de las variables explicativas para los que calcula previsiones. En el fragmento que sigue se crea un data.frame llamado new con nuevos valores de las variables cigs y etnic y se utiliza el modelo m1.a para obtener previsiones. Como se aprecia en el ejemplo, en new sólo son necesarias las variables que entran en el modelo:

> new <- data.frame(cigs=c(10, 15, 20, 20),
+                 etnic=factor(c("Blanco", "Blanco", "Blanco", "No blanco")))
> new

  cigs     etnic
1   10    Blanco
2   15    Blanco
3   20    Blanco
4   20 No blanco

> predict(m1.a, new)

       1        2        3        4 
3.289126 3.217409 3.145691 2.971390

Algunos argumentos opcionales controlan los resultados que devuelve la función predict(). En particular, si a la llamada se le añade se.fit=TRUE, los resultados incluyen las desviaciones típicas de las previsiones.

Además, se pueden obtener intervalos de confianza utilizando el argumento interval, que puede tomar los valores "prediction" y "confidence" en función de que se opte por calcular los intervalos de confianza de las previsiones puntuales o de la esperanza de las mismas, respectivamente. Ligado a interval esta el argumento level, que sirve para especificar el nivel de confianza de los intervalos (por defecto level=0.95). Compare:

> predict(m1.a, new, interval="prediction")

       fit      lwr      upr
1 3.289126 2.176970 4.401283
2 3.217409 2.104078 4.330739
3 3.145691 2.030629 4.260754
4 2.971390 1.855093 4.087688

> predict(m1.a, new, interval="confidence")

       fit      lwr      upr
1 3.289126 3.236987 3.341266
2 3.217409 3.144394 3.290424
3 3.145691 3.049819 3.241563
4 2.971390 2.862088 3.080692

4.4 Regresión local

El modelo lineal general representa la mejor aproximación lineal a todo el conjunto de datos. Sin embargo, cabe plantearse si dicha aproximación es adecuada para todo el rango de valores de cada variable explicativa o sería mejor utilizar aproximaciones lineales algo diferentes para distintos tramos. Los algoritmos de LOWESS (locally weighted exponential smoothing) estiman regresiones ponderadas utilizando un subconjunto de h observaciones contiguas denominado ancho de ventana (que se suele especificar como un porcentaje del total de observaciones). Para cada observación de la muestra, el algoritmo LOWESS devuelve el valor ajustado mediante la regresión local y a la secuencia gráfica de dichos valores ajustados se le denomina regresograma (aunque no es una definición única).

Si bien dichas estimaciones no proporcionan una forma funcional como la del modelo lineal de (1), sirven para analizar la presencia de no linealidades en la relación entre la variable dependiente y cada una de las explicativas en el modelo. La estimación LOWESS se suele presentar de forma gráfica, y una vez obtenida es posible calcular interpolaciones para nuevos valores de la variable explicativa.

En R, las funciones loess() y lowess() devuelven esta clase de ajuste, que es sencillo añadir a un gráfico y se denominan regresogramas. Además, el paquete ggplot2 incluye opciones para incluir un ajuste mediante regresión local a una nube de puntos. A continuación se añade un ajuste lineal (en azul) y otro local (en rojo) al gráfico XY del peso de los recién nacidos y cigarrillos fumados por la embarazada (Gráfico 6):

> bwd  <- subset(bwd, bwght < 5)  ## quitamos un dato atípico
> g.sm  <- qplot(cigs, bwght, data=bwd, geom="point") +
+                geom_smooth(color="red", se=FALSE) +
+                geom_smooth(method="lm", se=FALSE)
> ggsave('smooth-plot.pdf', g.sm)

Figure 6: Gráfico XY con ajuste lineal y local.

Como se puede apreciar, el ajuste lineal indica un descenso sistemático del peso de los recién nacidos, sin embargo, el ajuste local parece indicar que la media de peso deja de descender a partir de 20 cigarrillos (si bien hay pocas observaciones).

Si queremos los resultados numéricos, la función loess() es similar a lm(). Se le pasa una fórmula, un data.frame con los datos y admite argumentos adicionales para controlar el algoritmo de suavizado. En particular el porcentaje de observaciones que se usan en un ajuste local se especifica con span. En el siguiente ejemplo hace falta un porcentaje muestral alto para que el algoritmo funcione:

> m.loess <- loess(bwght~cigs, data=bwd, span=0.9)
> predict(m.loess, data.frame(cigs=seq(0, 40, by=5)))

       1        2        3        4        5        6        7        8 
3.398268 3.061770 2.846092 2.730101 2.692666 2.712654 2.768933 2.840370 
       9 
2.905834

Las funciones predict() y fitted() también operan sobre modelos de ajuste local. En el ejemplo se ha creado directamente un data.frame en el segundo argumento para obtener los valores ajustados para cierto número de cigarrillos. Si el modelo se estima con el objetivo de extrapolar (utilizar valores de la variable explicativa fuera del rango muestral), es necesario especificar el argumento surface="direct" en la llamada a loess().

En la Sección 6.4 veremos cómo contrastar más formalmente la posible presencia de no linealidades mediante el test RESET.

GEBR 4. Modelo lineal general II

Docencia, Econometría Grados No Responses »

feb 282012

4.2 Contraste de hipótesis lineales

Para fijar la notación, se quiere contrastar un conjunto de m hipótesis lineales sobre el modelo que se formulan: H₀:Aβ = c, donde A (m ×k) y c (m ×1) están formados por constantes conocidas. La hipótesis alternativa es genérica: H₁: Aβ ≠ c. Existen dos formas equivalentes de realizar el contraste, aunque con resultados idénticos:

$F = \frac{1}{m} (A\hat{\boldsymbol{\beta}}-\mathbf{c})' \left [ A V(\hat{\boldsymbol{\beta}}) A \right ]^{-1} (A\hat{\boldsymbol{\beta}}-\mathbf{c}) = \frac{n-k}{m} \frac{SCR_R - SCR_{NR}}{SCR_{NR}} \sim F_{m, n-k}.$

El uso de una u otra forma dependerá de la información disponible. La primera versión sólo requiere estimar el modelo sin restricciones (contraste de Wald), mientras que para usar la segunda hacen falta las sumas de cuadrados residuales del modelo que incorpora las restricciones (SCR_R) y del que no las incorpora (SCR_NR), que es la formulación de razón de verosimilitudes.

Consideremos el modelo de peso de los bebés donde, al apreciar un comportamiento no lineal de los cigarrillos fumados por la madre sobre el peso (Gráfico 3(c) y modelo m2.b), se decide introducir dicha variable por tramos basados en los cuartiles (para tener suficientes observaciones en cada tramo):

> quantile(bwd$cigs[bwd$cigs>0])

  0%  25%  50%  75% 100%
   1    8   10   20   50

> bwd$cigsint <- cut(bwd$cigs, c(0, 1, 5, 10, 20, 60), right=FALSE)
> table(bwd$cigsint)

  [0,1)   [1,5)  [5,10) [10,20) [20,60)
   1176      23      35      79      75

> m2.int <- lm(bwght~log(faminc)+sexo+etnic + cigsint, data=bwd)
> summary(m2.int)

Call:
lm(formula = bwght ~ log(faminc) + sexo + etnic + cigsint, data = bwd)

Residuals:
    Min      1Q  Median      3Q     Max
-2.6662 -0.3350  0.0135  0.3696  4.2728

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)     3.39114    0.06326  53.604  < 2e-16 ***
log(faminc)     0.02733    0.01799   1.519  0.12889
sexoMujer      -0.08366    0.03045  -2.748  0.00608 **
etnicNo blanco -0.15798    0.03951  -3.999 6.71e-05 ***
cigsint[1,5)   -0.15307    0.11928  -1.283  0.19960
cigsint[5,10)  -0.16568    0.09706  -1.707  0.08805 .
cigsint[10,20) -0.24400    0.06606  -3.694  0.00023 ***
cigsint[20,60) -0.28276    0.06788  -4.165 3.30e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5644 on 1380 degrees of freedom
Multiple R-squared: 0.04816,	Adjusted R-squared: 0.04333
F-statistic: 9.975 on 7 and 1380 DF,  p-value: 3.351e-12

A la vista de la estimación, cabe preguntarse si el efecto sobre el peso de fumar entre 1 y 4 cigarrillos es igual al de fumar entre 5 y 9, y también si el efecto de fumar entre 10 y 19 es igual al de fumar 20 o más. Para hacer el contraste en R podemos utilizar la función linearHypothesis() del paquete car. Además de pasarle el modelo sin restricciones, las hipótesis se pueden formular utilizando las etiquetas de las variables tal y como aparecen al hacer el summary() del modelo:

> library(car)
> linearHypothesis(m2.int, c("cigsint[1,5)=cigsint[5,10)",
+                            "cigsint[10,20)=cigsint[20,60)"),
+                  test="F")

Linear hypothesis test

Hypothesis:
cigsint[5) - cigsint[5,10) = 0
cigsint[10,20) - cigsint[20,60) = 0

Model 1: restricted model
Model 2: bwght ~ log(faminc) + sexo + etnic + cigsint

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1382 439.63
2   1380 439.57  2  0.059924 0.0941 0.9102

No se puede rechazar la hipótesis nula, por lo que reformulamos el modelo con nuevos intervalos:

> bwd$cigsint <- cut(bwd$cigs, c(0, 1, 10, 60), right=FALSE)
> m3.int <- lm(bwght~log(faminc)+sexo+etnic + cigsint, data=bwd)
> #mtable(m2.int, m3.int)  # Omito la salida

El contraste anterior también se puede realizar especificando la matriz A y el vector c:

> A <- matrix(0, nrow=2, ncol=8)
> A[1, 5] <- 1
> A[1, 6] <- -1
> A[2, 7] <- 1
> A[2, 8] <- -1
> linearHypothesis(m2.int, hypothesis.matrix=A, rhs=c(0,0))

En ocasiones no es sencillo especificar la hipótesis nula como un conjunto de restricciones lineales en la forma Aβ=c, en casos así se puede recurrir a la estimación del modelo restringido y realizar el contraste mediante sumas residuales. Por ejemplo, a la vista de las estimaciones de m2.int y m3.int, cabría contrastar si los coeficientes asociados a los distintos niveles de consumo de cigarrillos son todos iguales. Es decir, si es el hecho de fumar y no la cantidad de cigarrillos fumados lo que tiene un efecto significativo sobre el peso del bebé. En este caso, queremos contrastar H₀: β₅ = β₆ = β₇ = β₈, pero no se puede especificar en la forma Aβ=c, excepto que contrastemos la igualdad de pares de coeficientes. En un caso como ese, una alternativa consiste en estimar el modelo imponiendo la restricción y realizar el contraste utilizando las sumas de cuadrados residuales del modelo restringido y no restringido. Este contraste se puede realizar en R mediante la función anova() a la que se le pasan ambos modelos:

> m4.int <- lm(bwght~log(faminc)+sexo+etnic + mothsmoke, data=bwd)
> anova(m4.int, m2.int)

Analysis of Variance Table

Model 1: bwght ~ log(faminc) + sexo + etnic + mothsmoke
Model 2: bwght ~ log(faminc) + sexo + etnic + cigsint
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1383 440.07
2   1380 439.57  3   0.49767 0.5208  0.668

Como se puede apreciar, no podemos rechazar la hipótesis nula de que m2.int y m4.int son estadísticamente iguales: apenas existe ganancia al utilizar los cuatro niveles de consumo de cigarrillos frente a considerar, sencillamente, que la madre es fumadora. La comparación de los tres modelos que se han estimado en este apartado aparece a continuación.

> mtable(m2.int, m3.int, m4.int)

Calls:
m2.int: lm(formula = bwght ~ log(faminc) + sexo + etnic + cigsint, data = bwd)
m3.int: lm(formula = bwght ~ log(faminc) + sexo + etnic + cigsint, data = bwd)
m4.int: lm(formula = bwght ~ log(faminc) + sexo + etnic + mothsmoke,
    data = bwd)

======================================================
                          m2.int    m3.int    m4.int
------------------------------------------------------
(Intercept)               3.391***  3.391***  3.153***
                         (0.063)   (0.063)   (0.067)
log(faminc)               0.027     0.027     0.028
                         (0.018)   (0.018)   (0.018)
sexo: Mujer/Hombre       -0.084**  -0.083**  -0.081**
                         (0.030)   (0.030)   (0.030)
etnic: No blanco/Blanco  -0.158*** -0.158*** -0.157***
                         (0.040)   (0.039)   (0.039)
cigsint: [1,5)/[0,1)     -0.153
                         (0.119)
cigsint: [5,10)/[0,1)    -0.166
                         (0.097)
cigsint: [10,20)/[0,1)   -0.244***
                         (0.066)
cigsint: [20,60)/[0,1)   -0.283***
                         (0.068)
cigsint: [1,10)/[0,1)              -0.161*
                                   (0.076)
cigsint: [10,60)/[0,1)             -0.263***
                                   (0.049)
mothsmoke: No fuma/Fuma                       0.235***
                                             (0.043)
------------------------------------------------------
R-squared                    0.048     0.048     0.047
adj. R-squared               0.043     0.045     0.044
sigma                        0.564     0.564     0.564
F                            9.975    13.946    17.083
p                            0.000     0.000     0.000
Log-likelihood           -1171.513 -1171.607 -1172.298
Deviance                   439.572   439.631   440.069
AIC                       2361.025  2357.214  2356.596
BIC                       2408.146  2393.864  2388.009
N                         1388      1388      1388
======================================================

GEBR 4. Modelo lineal general I

Docencia, Econometría Grados 5 Responses »

feb 242012

La especificación del modelo lineal general que usaremos es:

$y_i = \beta_1 + \beta_2 x_{i2} + \ldots + \beta_k x_{ik} + \varepsilon_i = \mathbf{x}_i ' \boldsymbol{\beta} + \varepsilon_i, \qquad i=1, \ldots, n,$

que en notación matricial se expresa:

$\mathbf{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}.$

Necesitamos algunos supuestos sobre la especificación anterior, que resumidos son:

La relación entre la variable dependiente y_i y las explicativas x_i1,…, x_ik:
1. Es una relación lineal (véase Sección 6.4).
2. Los parámetros desconocidos β_i son constantes (véase Sección 6.3).
3. Entre las variables explicativas no falta ninguna relevante y no se incluyen variables irrelevantes (véase Sección 7.3).
Sobre las perturbaciones ε_i se supone que (véase Sección 5.5):
1. La esperanza es nula: E(ε _i)=0 ∀i.
2. La varianza es constante: E(ε_i)=σ² ∀i.
3. Son independientes entre sí, por lo que Cov(ε_i, ε_j)=0 ∀j ≠ i.
Las hipótesis anteriores, en notación matricial, se expresan: E(ε) = 0_n y V(ε) = σ² I_n. Aunque no es estrictamente necesario, a menudo se supone normalidad: ε ∼ N(0_n, σ²I_n ).
Sobre las variables explicativas se supone que:
1. Se obtendrían los mismos datos en muestreos repetidos (son deterministas) o, que al menos, son independientes de las perturbaciones.
2. Son linealmente independientes, rango(X)= rango(X′X) = k (véase Sección 6.1).
3. Se dispone de muchas más observaciones que parámetros a estimar: n >> k.

4.1 Estimación por mínimos cuadrados ordinarios

Los estimadores de de los parámetros desconocidos por mínimos cuadrados ordinarios (MCO) son:

$\hat{\boldsymbol{\beta}} = (X'X)^{-1} X' \mathbf{y} \qquad$ y $\qquad \hat{\sigma}^2 = \frac{\hat{\boldsymbol{\varepsilon}} ' \hat{\boldsymbol{\varepsilon}}}{n-k},$

donde $\hat{\boldsymbol{\varepsilon}}= \mathbf{y} - X \hat{\boldsymbol{\beta}}$ . Si las hipótesis formuladas sobre X y ε se cumplen, entonces se puede probar que $E(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta}$ y que $\widehat{V(\hat{\boldsymbol{\beta}})} = \hat{\sigma}^2 (X'X)^{-1}$ , y bajo condiciones bastante generales que $\hat{\boldsymbol{\beta}} \rightarrow N_k [ \boldsymbol{\beta}, \hat{\sigma}^2 (X'X)^{-1} ]$ .

La función para estimar un modelo lineal como el por MCO en R es lm(), que requiere dos argumentos: una fórmula de R que especifica el modelo y un data.frame que contiene las variables. Por ejemplo, un modelo simple que explica el peso en función de los cigarrillos que fuma la madre durante el embarazo se estima con:

> m1 <- lm(bwght~cigs, data=bwd)
> summary(m1)

Call:
lm(formula = bwght ~ cigs, data = bwd)

Residuals:
    Min      1Q  Median      3Q     Max
-2.7434 -0.3337  0.0084  0.3750  4.2872

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.395476   0.016226 209.267  < 2e-16 ***
cigs        -0.014565   0.002565  -5.678 1.66e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5706 on 1386 degrees of freedom
Multiple R-squared: 0.02273,	Adjusted R-squared: 0.02202
F-statistic: 32.24 on 1 and 1386 DF,  p-value: 1.662e-08

En el resumen de la estimación que se obtiene con la función summary(), los coeficientes de cada variable explicativa aparecen etiquetados con el nombre de la variable. La etiqueta (Intercept) corresponde al término constante que R incluye automáticamente. Además de los coeficientes estimados, en la tabla aparecen sus desviaciones típicas, contrastes t para la nula de que el coeficiente es 0 y el p-valor del contraste para la alternativa de dos colas. En el ejemplo, partiendo de un peso medio de 3.395 kg. para un recién nacido de madre no fumadora, el peso medio se reduce en 14.6 gramos por cada cigarrillo diario que fuma la recién parida durante el embarazo.

Se pueden añadir más variables explicativas y las transformaciones se pueden especificar directamente en la fórmula:

> m2 <- lm(bwght~cigs+log(faminc), data=bwd)
> summary(m2)

Call:
lm(formula = bwght ~ cigs + log(faminc), data = bwd)

Residuals:
    Min      1Q  Median      3Q     Max
-2.7429 -0.3322  0.0199  0.3744  4.2544

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.231582   0.055109  58.640  < 2e-16 ***
cigs        -0.013246   0.002592  -5.110 3.68e-07 ***
log(faminc)  0.052467   0.016865   3.111   0.0019 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5689 on 1385 degrees of freedom
Multiple R-squared: 0.02951,	Adjusted R-squared: 0.02811
F-statistic: 21.06 on 2 and 1385 DF,  p-value: 9.796e-10

Algunas transformaciones deben definirse detro de la función I() porque tienen una interpretación especial dentro de una fórmula. En particular el cuadrado, que genera interacciones entre variables:

> m2.b <- lm(bwght~cigs+I(cigs^2)+log(faminc), data=bwd)
> summary(m2.b)

Call:
lm(formula = bwght ~ cigs + I(cigs^2) + log(faminc), data = bwd)

Residuals:
    Min      1Q  Median      3Q     Max
-2.7477 -0.3359  0.0314  0.3787  4.2503

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.2393400  0.0552460  58.635  < 2e-16 ***
cigs        -0.0225151  0.0058922  -3.821 0.000139 ***
I(cigs^2)    0.0003635  0.0002076   1.751 0.080083 .
log(faminc)  0.0515053  0.0168613   3.055 0.002296 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5684 on 1384 degrees of freedom
Multiple R-squared: 0.03166,	Adjusted R-squared: 0.02956
F-statistic: 15.08 on 3 and 1384 DF,  p-value: 1.16e-09

> m2.c <- lm(bwght~(cigs+log(faminc))^2, data=bwd)
> summary(m2.c)

Call:
lm(formula = bwght ~ (cigs + log(faminc))^2, data = bwd)

Residuals:
    Min      1Q  Median      3Q     Max
-2.7436 -0.3366  0.0246  0.3757  4.2482

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)       3.205332   0.057641  55.609  < 2e-16 ***
cigs             -0.003480   0.006829  -0.510 0.610389
log(faminc)       0.061113   0.017761   3.441 0.000597 ***
cigs:log(faminc) -0.003752   0.002428  -1.545 0.122463
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5686 on 1384 degrees of freedom
Multiple R-squared: 0.03118,	Adjusted R-squared: 0.02908
F-statistic: 14.85 on 3 and 1384 DF,  p-value: 1.616e-09

Las variables explicativas cualitativas (factores en la notación de R) se introducen en el modelo igual que las continuas. R se encarga de generar las variables ficticias necesarias y hacer el tratamiento adecuado de las mismas (excluyendo una de ellas por defecto). Si en el modelo queremos introducir, por ejemplo, el grupo étnico de los bebés como explicativo del peso además de los cigarrilos, las instrucciones para estimar son:

> m1.a <- lm(bwght~cigs+etnic, data=bwd)
> summary(m1.a)

En cambio, si queremos introducir el grupo étnico más su interacción con los cigarrillos fumados por la madre (ficticia multiplicativa) se haría:

> m1.m <- lm(bwght~cigs*etnic, data=bwd)
> summary(m1.m)

Call:
lm(formula = bwght ~ cigs * etnic, data = bwd)

Residuals:
    Min      1Q  Median      3Q     Max
-2.5885 -0.3489  0.0197  0.3672  4.2438

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)          3.438962   0.018206 188.888  < 2e-16 ***
cigs                -0.017497   0.003020  -5.795 8.47e-09 ***
etnicNo blanco      -0.198418   0.038983  -5.090 4.08e-07 ***
cigs:etnicNo blanco  0.010864   0.005604   1.939   0.0528 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5658 on 1384 degrees of freedom
Multiple R-squared: 0.04076,	Adjusted R-squared: 0.03868
F-statistic:  19.6 on 3 and 1384 DF,  p-value: 1.898e-12

Como se puede apreciar, lm() incluye automáticamente una constante en los modelos, pero si se desea excluirla se puede añadir +0 a la fórmula. La función lm() acepta otros argumentos, como subset que permite seleccionar un subconjunto de observaciones para la estimación. Por ejemplo, para estimar el modelo de peso en función de cigarrillos sólo para las madres fumadoras:

> m1.mfuma <- lm(bwght~cigs, data=bwd, subset=(cigs>0))
> summary(m1.mfuma)

Además del comando summary() que resume la estimación de un modelo, el paquete memisc facilita la comparación de las estimaciones con la función mtable(), que calcula algunos estadísticos adicionales. Por ejemplo:

> library(memisc)
> mtable(m1.a, m1.m)

Calls:
m1.a: lm(formula = bwght ~ cigs + etnic, data = bwd)
m1.m: lm(formula = bwght ~ cigs * etnic, data = bwd)

===================================================
                                  m1.a      m1.m
---------------------------------------------------
(Intercept)                      3.433***  3.439***
                                (0.018)   (0.018)
cigs                            -0.014*** -0.017***
                                (0.003)   (0.003)
etnic: No blanco/Blanco         -0.174*** -0.198***
                                (0.037)   (0.039)
cigs x etnic: No blanco/Blanco             0.011
                                          (0.006)
---------------------------------------------------
R-squared                           0.038     0.041
adj. R-squared                      0.037     0.039
sigma                               0.566     0.566
F                                  27.472    19.604
p                                   0.000     0.000
Log-likelihood                  -1178.769 -1176.887
Deviance                          444.192   442.989
AIC                              2365.537  2363.774
BIC                              2386.480  2389.952
N                                1388      1388
===================================================

Los residuos y valores ajustados de un modelo estimado los podemos extraer con las funciones residuals() y fitted(), respectivamente. Vamos a añadir los residuos tipificados (usando la función scale() para tipificar) y los valores ajustados de la variable dependiente al data.frame que contiene los datos para hacer un gráfico (coloreando los puntos según los cigarrillos que fuma la madre):

> bwd$e    <- scale(residuals(m1))
> bwd$yhat <- fitted(m1)
> resplot  <- qplot(yhat, e, data=bwd, geom="point", colour=cigs,
+                 xlab=expression(hat(y)[i]), ylab=expression(hat(epsilon)[i]))

En el Gráfico 1 se aprecia la mayor dispersión de los residuos en los pesos más elevados comparada con la dispersión en los pesos más bajos. Sin embargo la comparación visual es delicada porque hay, comparativamente, pocas observaciones con pesos bajos.

Gráfico 1: Residuos del modelo m1

Además de residuos y valores ajustados, es posible acceder a otros resultados de la estimación con coef(), vcov() y los distintos componentes del objeto summary(). Por ejemplo:

> coef(m1)

(Intercept)        cigs
 3.39547626 -0.01456519

> vcov(m1)

              (Intercept)          cigs
(Intercept)  2.632696e-04 -1.373601e-05
cigs        -1.373601e-05  6.581147e-06

> summary(m1)$sigma

[1] 0.5706356

> summary(m1)$adj.r.squared

[1] 0.02202402

> summary(m1)$fstatistic

     value      numdf      dendf
  32.23524    1.00000 1386.00000