4 Comparación de Modelos

Objetivos del capítulo

Al finalizar este capítulo, el lector estará en capacidad de:

Explicar en sus propias palabras la intuición detrás de las métricas \(R^2\) ajustado, el \(AIC\) y el \(BIC\) para comparar modelos.
Emplear el \(R^2\) ajustado, el \(AIC\) y el \(BIC\) para seleccionar el mejor modelo.
Explicar en sus propias palabras la intuición detrás de las pruebas estadísticas que permitan seleccionar entre modelos anidados y no anidados empleando R.
Realizar pruebas estadísticas que permitan seleccionar entre modelos anidados empleando R.
Realizar pruebas estadísticas que permitan seleccionar entre modelos no anidados empleando R.

4.1 Introducción

Con frecuencia, el científico de datos se enfrenta con el problema de tener que comparar diferentes modelos y escoger uno de ellos. En este capítulo nos concentraremos en dos tipos de aproximaciones, que no son excluyentes, para seleccionar modelos. La primera aproximación emplea medidas de bondad de ajuste y la segunda pruebas estadísticas. En lo que resta de este capítulo supondremos que queremos comparar modelos cuya variable dependiente es la misma y que emplean la misma muestra.

Pero antes de entrar en detalle, es importante anotar que al momento de comparar modelos que expliquen por ejemplo la variable \(y_i\) podemos enfrentar dos situaciones. Por ejemplo, consideremos los siguientes modelos para explicar las unidades vendidas de un SKU de una bebida carbonizada en el mes \(t\) (\(V_t\)):

\[\begin{equation} V_t = \beta_1 + \beta_2 p_t + \beta_3 p{c_t} + \beta_4 tem{p_t} + \beta_5 PubliCom{p_t} + \varepsilon_t \tag{4.1} \end{equation}\]

\[\begin{equation} V_t = \beta_1 + \beta_2 p_t + \beta_3 p{c_t} + \varepsilon_t \tag{4.2} \end{equation}\]

\[\begin{equation} V_t = \beta_1 + \beta_6 PubliT{V_t} + \beta_7 PubliRadi{o_t} + \beta_8 PubliInt{r_t}+ \varepsilon_t \tag{4.3} \end{equation}\] donde \(p_t\) y \(p{c_t}\) denotan el precio por mililitro en el mes \(t\) de la bebida bajo estudio y del competidor más cercano. \(tem{p_t}\) representa la temperatura promedio en el mes \(t\) en grados centígrados. \(PubliCom{p_t}\), \(PubliT{V_t}\), \(PubliRadi{o_t}\) y \(PubliInt{r_t}\) corresponden a la inversión en el mes \(t\) en millones de pesos en publicidad total del competidor más cercano, la inversión propia en televisión, radio e Internet, respectivamente.

Noten que el modelo (4.2) es una versión restringida del modelo (4.1); pues si \({\beta_4}= {\beta_5} =0\) en el modelo (4.1) entonces se obtiene el modelo restringido (4.2). Al modelo (4.2) se le denomina modelo anidado (en el modelo (4.1)) . Por otro lado, el modelo (4.3) no se encuentra anidado en los modelos (4.2) o (4.1). Es decir, no hay forma de restringir los modelos (4.1) o (4.2) para obtener el modelo (4.3).

4.2 Comparación de modelos empleando medidas de bondad de ajuste

En el capítulo anterior (ver sección 3.3) discutimos las limitaciones del \(R^2\) al comparar modelos. Concluimos que el \(R^2\) solo permite comparar modelos si estos tienen la misma variable dependiente y el mismo número de variables explicativas. Mostramos además cómo el \(R^2\) se infla con la inclusión de variables.

Para tener en cuenta la relación directa entre el \(SSR\) y el número de regresores (\(k-1\)), se ha diseñado el \(R^2\) ajustado (\({\bar R^2}\)). El \({\bar R^2}\) penaliza la inclusión de nuevas variables explicativas al modelo, esa penalización se da de la siguiente forma: \[\begin{equation*} {\bar R^2} = 1 - \left( {1 - {R^2}} \right)\frac{{n - 1}}{{n - k}} \end{equation*}\] \(\frac{{n - 1}}{{n - k}}\) es el factor que penaliza el \(R^2\) al incluir más variables explicativas, \(\frac{{n - 1}}{{n - k}}\) decrece con el aumento de \(k\) (inclusión de variables). En otras palabras, a medida que se incrementa el número de variables independientes en el modelo, \(\frac{{n - 1}}{{n - k}}\) se hace más pequeño al mismo tiempo que el \(R^2\) aumenta. Por tanto, \(\left( {1 - {R^2}} \right)\) disminuye a medida que se incluyen variables, pero el factor \(\frac{{n - 1}}{{n - k}}\) se hace más pequeño. Así, el efecto de un aumento en el número de regresores no necesariamente implica un aumento en el \({\bar R^2}\). El \({\bar R^2}\) incrementará únicamente si el aumento en el \(R^2\) es lo suficientemente grande para compensar la penalización que se hace por incluir más variables. Por lo tanto, un \({\bar R^2}\) grande será mejor que uno pequeño.

De esta manera, con el \({\bar R^2}\) se pueden comparar los modelos (4.1), (4.2) y (4.3), siempre y cuando se emplee la misma muestra para estimar los tres modelos (el mismo \(n\) y \(SST\)). El único problema que se presenta con el \({\bar R^2}\), es que este carece de una interpretación a diferencia del \(R^2\).

Existen otros problemas, que se estudiarán más adelante, que pueden “inflar” el \(R^2\) como la multicolinealidad (ver Capítulo 8). De esta manera, aunque este estadístico presenta una interpretación muy clara e intuitiva, hay que tener cuidado antes de sacar conclusiones de este estadístico y del \({\bar R^2}\).

Existen otras medidas de bondad de ajuste que permiten comparar entre modelos (como el \({\bar R^2}\)) pero que no tienen una interpretación como tal. Dos de estas medidas muy conocidas son el Criterio de Información de Akaike (por su nombre en inglés: \(AIC\)) y el Criterio de Información de Bayesiano⁴⁴ (por su nombre en inglés: \(BIC\)). Tanto el AIC como el BIC son estadísticos que fueron desarrollados en otra filosofía de estimación (Máxima Verosimilitud) y penalizan por la inclusión de más variables explicativas como el \({\bar R^2}\). Estas dos medidas adaptadas para el caso de la regresión múltiple estimada por MCO se definen como: \[\begin{equation*} AIC = n + n\log \left( {2\pi } \right) + n\log \left( {\frac{{SSR}}{n}} \right) + 2\left( {k + 1} \right) \end{equation*}\]

\[\begin{equation*} BIC = n + n\log \left( {2\pi } \right) + n\log \left( {\frac{{SSR}}{n}} \right) + \log \left( n \right)\left( {k + 1} \right). \end{equation*}\]

Cuando empleamos el \({\bar R^2}\) para comparar modelos con la misma variable dependiente (y misma muestra) se selecciona el modelo que lo maximice. En el caso del \(AIC\) y \(BIC\), se prefiere el modelo que minimice estos criterios. Por otro lado, la diferencia entre el \(AIC\) y el \(BIC\) es la forma en cómo se penaliza la inclusión de más variables explicativas, por esto es posible que el modelo seleccionado sea diferente para estos dos criterios. Está muy documentado que el \(AIC\) tiene siempre una probabilidad alta de seleccionar modelos con más regresores que el \(BIC\). Y por el otro lado, también se sabe que el \(BIC\) tiende a encontrar modelos relativamente pequeños en términos de variables explicativas. Adicionalmente, es importante anotar que estas medidas de bondad de ajuste permiten comparar modelos anidados y no anidados.

En general, para la comparación de modelos se recomienda emplear las tres medidas (o métricas) (\({\bar R^2}\), \(AIC\), \(BIC\)) de bondad de ajuste; pero los resultados de cada uno de estos criterios no necesariamente concordarán. Por ejemplo, supongamos un caso en el que \(AIC\) recomienda un modelo con 5 variables, el \(BIC\) recomienda 2 variables y el \({\bar R^2}\) sugiere 3 variables. En este caso se puede emplear los tres modelos seleccionados y comparar entre modelos empleando inferencia como veremos en la siguiente sección.

4.3 Comparación de modelos empleando inferencia

4.3.1 Modelos anidados

En el capítulo anterior discutimos la inferencia para restricciones de la forma \({R_{\left( {r} \right) \times k}}{\beta_{k \times 1}} = {C_{\left( {r} \right) \times 1}}\). Discutimos como dichas restricciones se podrán probar con una prueba F o con una prueba de Wald (ver sección 3.5.

De hecho, la comparación de modelos anidados es un caso especial de dichas pruebas. Es decir, estamos comparando modelos anidados cuando consideramos una hipótesis nula en la cual un grupo de coeficientes⁴⁵ es conjuntamente igual a cero, es decir \({H_o}:{\beta_p} = {\beta_{p + 1}} = \ldots = {\beta_l} = 0\), para \(0 < p < l \le k\) (versus la \(H_{A}\): No \(H_0\)). Esta prueba se conoce como una prueba de significancia conjunta. Es decir, la hipótesis nula es equivalente a \({H_o}:{Y_i} = {\beta_1} + {\beta_2}{X_{2i}} + \ldots + {\beta_{p - 1}}{X_{p - 1i}} + {\beta_{l + 1}}{X_{l + 1i}} + \ldots + {\beta_k}{X_{ki}} + {\varepsilon _i}\) (es decir un Modelo Restringido (R) del original) y la hipótesis alterna es \({H_A}:{Y_i} = {\beta_1} + {\beta_2}{X_{2i}} + {\beta_3}{X_{3i}} + \ldots + {\beta_k}{X_{ki}} + {\varepsilon _i}\) (el Modelo sin restricción (U). En este caso especial, el \(F_{Calculado}\) de la expresión (3.12) es equivalente a: \[\begin{equation*} F_{calculado} = \frac{\left ( SSE_R - SSE_U \right )/r}{SSE_U / (n-k)} \end{equation*}\] donde \(r = l - p\), \(SSE_R\) representa la suma de los cuadrados de los residuos estimados del modelo restringido y \(SSE_U\) representa la suma de los cuadrados de los residuos estimados del modelo sin restringir.

Así, para probar la hipótesis nula que \({H_o}:{\beta_p} = {\beta_{p + 1}} = \ldots = {\beta_l} = 0\), para \(0 < p < l \le k\) (versus la \(H_{A}\): No \(H_0\)), podemos estimar por MCO el modelo restringido, implicado por la hipótesis nula, como el modelo sin restringir y encontrar su correspondiente \(SSE\). A partir de estas cantidades se puede calcular el \(F_{Calculado}\), el cual se compara con \({F_{\alpha ,(r,(n - k))}}\). A este tipo de test se le conoce como una prueba de modelo restringido versus modelo no restringido. Como se discutió en el capítulo anterior, este tipo de restricciones también se pueden probar con una prueba de Wald (función waldtest() del paquete AER (Kleiber & Zeileis, 2008)).

4.3.2 Modelos no anidados

Formalmente, en este caso tendremos que la hipótesis nula de la prueba será: \[\begin{equation*} H_0 :{\mathbf{y}} = {\mathbf{X}}\beta + \varepsilon _0 \end{equation*}\] y la hipótesis alterna es: \[\begin{equation*} H_A :{\mathbf{y}} = Z\gamma + \varepsilon _1 . \end{equation*}\]

Para comprobar este tipo de hipótesis tenemos dos opciones: la prueba J y la prueba de Cox. Estas opciones se discuten a continuación.

4.3.2.1 Prueba J

La idea de esta prueba es bastante evidente. La prueba J considera los dos modelos en uno. Es decir, implica estimar el siguiente modelo: \[\begin{equation*} {\mathbf{y}} = \left( {1 - \lambda } \right){\mathbf{X}}\beta + \left( \lambda \right){\mathbf{Z}}\gamma + \varepsilon \end{equation*}\] La idea de la prueba es primero estimar \(\gamma\) a partir de la regresión de \({\mathbf{y}}\) en \({\mathbf{Z}}\). En un segundo paso correr una regresión de \({\mathbf{y}}\) en \({\mathbf{X}}\) y \({\mathbf{Z}}\hat \gamma\). Noten que \(H_0\) es verdad si \(\lambda = 0\). Por tanto, la prueba implica el siguiente estadístico de prueba: \[\begin{equation*} \frac{{\hat \lambda }} {{s.e.\left( {\hat \lambda } \right)}} \end{equation*}\] Esta demostrado que este estadístico sigue una distribución aproximada a la estándar normal.

4.3.2.2 Prueba de Cox

La prueba de Cox es un tipo de prueba denominada razón de máxima verosimilitud (LR). Esta prueba implica el siguiente estadístico:

\[\begin{equation*} q = \frac{{c_{01} }} {{\sqrt {\frac{{s_{\mathbf{Z}}^2 }} {{s_{{\mathbf{ZX}}}^4 {\mathbf{b}}^T {\mathbf{X}}^T {\mathbf{M}}_{\mathbf{Z}} {\mathbf{M}}_{\mathbf{X}} {\mathbf{M}}_{\mathbf{Z}} {\mathbf{Xb}}}}} }} \end{equation*}\] donde, \[\begin{equation*} c_{01} = \frac{n} {2}\ln \left[ {\frac{{s_{\mathbf{Z}}^2 }} {{s_{{\mathbf{ZX}}}^2 }}} \right] \end{equation*}\] \[\begin{equation*} s_{\mathbf{Z}}^2 = \frac{{{\mathbf{e}}_{\mathbf{Z}}^T {\mathbf{e}}_{\mathbf{Z}} }} {n} \end{equation*}\] \[\begin{equation*} s_{\mathbf{X}}^2 = \frac{{{\mathbf{e}}_{\mathbf{X}}^T {\mathbf{e}}_{\mathbf{X}} }} {n} \end{equation*}\]

\[\begin{equation*} s_{{\mathbf{ZX}}}^2 = s_{\mathbf{X}}^2 + \frac{{{\mathbf{b}}^T {\mathbf{X}}^T {\mathbf{M}}_{\mathbf{Z}} {\mathbf{Xb}}}} {n}. \end{equation*}\] Este estadístico sigue una distribución Chi-cuadrado con grados de libertad igual al número de parámetros en el numerador.

4.4 Práctica en R: Escogiendo el mejor modelo

Para mostrar como emplear R para seleccionar modelos, emplearemos unos datos simulados para una variable dependiente (\(y_i\)) y 10 variables explicativas \(X_j,i\) donde \(j=1,2,...,10\). Para cada variable se simulan 150 observaciones \(i=1,2,...,150\). Los datos están disponibles en el archivo selModel.txt⁴⁶. Carguemos los datos y constatemos que quedan bien cargados.

# Cargar datos
datos <- read.table("./Data/selModel.txt", header = TRUE, sep = ",")
# Chequear datos cargados
head(datos, 2)

##   X    x1    x2    x3    x4    x5    x6    x7    x8    x9   x10      y
## 1 1 5.925 6.512 6.883 5.157 6.420 5.833 5.356 6.626 5.730 6.108 42.799
## 2 2 4.565 3.638 4.378 5.056 3.165 4.598 4.576 4.171 4.259 4.443 30.372

class(datos)

## [1] "data.frame"

str(datos)

## 'data.frame':    150 obs. of  12 variables:
##  $ X  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ x1 : num  5.92 4.57 5.44 4.91 5.7 ...
##  $ x2 : num  6.51 3.64 5.5 5.85 6.33 ...
##  $ x3 : num  6.88 4.38 6.6 6.54 5.28 ...
##  $ x4 : num  5.16 5.06 4.48 5.71 5.02 ...
##  $ x5 : num  6.42 3.17 4.99 4.67 4.23 ...
##  $ x6 : num  5.83 4.6 5.52 5.26 5.77 ...
##  $ x7 : num  5.36 4.58 6.24 5.79 4.47 ...
##  $ x8 : num  6.63 4.17 5.4 5.07 4.77 ...
##  $ x9 : num  5.73 4.26 5.46 5.5 6.06 ...
##  $ x10: num  6.11 4.44 5.57 4.64 5.5 ...
##  $ y  : num  42.8 30.4 36.3 36.6 38.3 ...

Noten que se cargó una primera columna que corresponde al número de la observación, esto no lo necesitaremos. Eliminemos esa variable.

datos <- datos[, -1] 
head(datos, 3)

##      x1    x2    x3    x4    x5    x6    x7    x8    x9   x10      y
## 1 5.925 6.512 6.883 5.157 6.420 5.833 5.356 6.626 5.730 6.108 42.799
## 2 4.565 3.638 4.378 5.056 3.165 4.598 4.576 4.171 4.259 4.443 30.372
## 3 5.436 5.501 6.596 4.476 4.992 5.524 6.235 5.402 5.465 5.568 36.338

str(datos)

## 'data.frame':    150 obs. of  11 variables:
##  $ x1 : num  5.92 4.57 5.44 4.91 5.7 ...
##  $ x2 : num  6.51 3.64 5.5 5.85 6.33 ...
##  $ x3 : num  6.88 4.38 6.6 6.54 5.28 ...
##  $ x4 : num  5.16 5.06 4.48 5.71 5.02 ...
##  $ x5 : num  6.42 3.17 4.99 4.67 4.23 ...
##  $ x6 : num  5.83 4.6 5.52 5.26 5.77 ...
##  $ x7 : num  5.36 4.58 6.24 5.79 4.47 ...
##  $ x8 : num  6.63 4.17 5.4 5.07 4.77 ...
##  $ x9 : num  5.73 4.26 5.46 5.5 6.06 ...
##  $ x10: num  6.11 4.44 5.57 4.64 5.5 ...
##  $ y  : num  42.8 30.4 36.3 36.6 38.3 ...

Todas las variables cargadas son numéricas. Ahora consideremos los siguientes tres modelos: \[\begin{equation} {y_i} = {\beta_1} + {\beta_2}{X_{1,i}} + {\beta_3}{X_{2,i}} + {\beta_4}{X_{3,i}} + {\beta_5}{X_{6,i}} + {\beta_6}{X_{7,i}} + {\varepsilon _i} \tag{4.4} \end{equation}\] \[\begin{equation} {y_i} = {\beta_1} + {\beta_2}{X_{1,i}} + {\beta_3}{X_{2,i}} + {\beta_4}{X_{3,i}} + {\varepsilon _i} \tag{4.5} \end{equation}\] \[\begin{equation} {y_i} = {\beta_1} + {\beta_2}{X_{4,i}} + {\beta_3}{X_{5,i}} + {\beta_8}{X_{8,i}} + {\beta_9}{X_{9,i}} + {\beta_{10}}{X_{10,i}} + {\varepsilon _i} \tag{4.6} \end{equation}\] El modelo (4.5) está anidado en el modelo (4.4). El modelo (4.6) no se encuentra anidado en ninguno de los otros modelos, ni al revés.

Procedamos a estimar los tres modelos y compararlos empleando las métricas de bondad de ajuste primero y luego empleando pruebas de hipótesis. Pero antes estimemos los tres modelos descritos arriba.

res1 <- lm(y ~ x1 + x2 + x3 + x6 + x7, datos)
res2 <- lm(y ~ x1 + x2 + x3 , datos)
res3 <- lm(y ~ x4 + x5 + x8 + x9 + x10, datos)

Los resultados de estos tres modelos se presentan en el Cuadro 4.1.

Cuadro 4.1: Modelos estimados por MCO
	Model 1	Model 2	Model 3
(Intercept)	12.39^***	13.18^***	14.03^***
	(1.49)	(1.43)	(1.59)
x1	0.95^***	1.05^***
	(0.27)	(0.27)
x2	1.94^***	2.13^***
	(0.30)	(0.27)
x3	1.04^***	1.16^***
	(0.29)	(0.28)
x6	0.41
	(0.27)
x7	0.17
	(0.29)
x4			1.43^***
			(0.31)
x5			1.53^***
			(0.30)
x8			0.44
			(0.32)
x9			0.18
			(0.32)
x10			0.59
			(0.33)
R²	0.65	0.64	0.57
Adj. R²	0.64	0.64	0.55
Num. obs.	150	150	150
^*p < 0.001; ^p < 0.01; ^*p < 0.05

4.4.1 Medidas de bondad de ajuste

El \({\bar R^2}\) y los criterios de información \(AIC\) y \(BIC\) se pueden calcular fácilmente en R. Para calcular el \({\bar R^2}\) debemos emplear la función summary() y extraer del objeto que crea summary el slot (compartimiento) que contiene a \({\bar R^2}\).

R1 <- summary(res1) 
attributes(R1)

## $names
##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled" 
## 
## $class
## [1] "summary.lm"

R1$adj.r.squared

## [1] 0.6388607

Los criterios de información se calculan empleando las funciones AIC() y BIC() del paquete base de R. Las dos funciones tienen como único argumento el objeto que contiene la regresión (objeto de clase lm).

AIC(res1)

## [1] 733.9571

BIC(res1)

## [1] 755.0315

Puedes calcular estos indicadores para los otros modelos. Los resultados que se reportan en el Cuadro 4.2.

Cuadro 4.2: Medidas de bondad de ajuste para los tres modelos estimados
	Modelo 1	Modelo 2	Modelo 3
R2.ajustado	0.639	0.636	0.551
AIC	733.957	733.323	766.718
BIC	755.032	748.376	787.792
Fuente: elaboración propia.

El \({\bar R^2}\) sugiere que el mejor modelo es el 1 ((4.4), mientras que los dos criterios de información sugieren que el mejor modelo es el 2 ((4.5)). Ahora veamos que decisión tomamos al emplear pruebas estadísticas.

4.4.2 Pruebas estadísticas

4.4.2.1 Modelos anidados

Como se mencionó anteriormente, el modelo (4.5) está anidado en el modelo (4.4). Comparemos estos dos modelos empleando la función . Esta función ya la había empleado en la sección 3.6, en ese momento solo empleamos un argumento (objeto de clase lm) para obtener la tabla ANOVA. Esta misma función permite efectuar una prueba F de modelo restringido versus modelo no restringido si se emplean dos argumentos: primero el modelo restringido (objeto de clase lm) y segundo el modelo sin restringir ⁴⁷ (objeto de clase lm). Empleando esta función encontramos lo siguiente⁴⁸:

anova(res2, res1)

## Analysis of Variance Table
## 
## Model 1: y ~ x1 + x2 + x3
## Model 2: y ~ x1 + x2 + x3 + x6 + x7
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    146 1091.1                           
## 2    144 1066.9  2    24.207 1.6337 0.1988

En este caso el estadístico F es igual a 1.63 y el respectivo valor p es 0.1988. Por tanto no se puede rechazar la hipótesis nula de que el modelo restringido (el modelo 2; es decir (4.5)) es mejor que el modelo sin restringir (modelo 1; es decir (4.4)).

4.4.2.2 Modelos no anidados

Ahora comparemos el modelo 3 ((4.6)) con los otros dos modelos. La prueba J se calcula empleando la función jtest() del paquete AER () . Esta función solo tiene dos argumentos, dos objetos de clase jtest() que contienen los dos modelos no anidados. En este caso el orden de los argumentos no es importante para efectuar la prueba pero si para analizar sus resultados como veremos a continuación.

library(AER)
J.res1.3 <- jtest(res1, res3)
J.res1.3

## J test
## 
## Model 1: y ~ x1 + x2 + x3 + x6 + x7
## Model 2: y ~ x4 + x5 + x8 + x9 + x10
##                 Estimate Std. Error t value  Pr(>|t|)    
## M1 + fitted(M2)  0.48221   0.103789   4.646 7.606e-06 ***
## M2 + fitted(M1)  0.81939   0.095489   8.581 1.427e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El primer \(t_{calculado}\) (4.65) es el que permite probar la hipótesis nula de que el modelo (4.4) (este modelo fue el primer argumento de la función jtest()) es mejor que el modelo (4.6) (el segundo argumento empleado en dicha función). Este resultado permite rechazar la nula (valor p de 0), lo cuál significa que la prueba no puede concluir en favor del modelo 1. Es decir, hay suficiente evidencia para rechazar la hipótesis nula de que el modelo (4.4) (“Model 1” en la salida de R) es mejor que el modelo (4.6) (“Model 2” en la salida).

Por otro lado, el segundo \(t_{calculado}\) (8.58) es el que permite probar la hipótesis nula de que el modelo (4.6) es mejor que el modelo (4.4). En este caso también se rechaza la nula (valor p de 0). Es decir, existe suficiente evidencia para rechazar la nula de que el modelo (4.6) (“Model 2” en la salida de R) es mejor que el modelo (4.4) (“Model 1” en la salida ). Esto significa que la prueba presenta una contradicción o dicho de otra manera, esta prueba no es concluyente para este caso.

Para el caso de la comparación del modelo (4.5) y el modelo (4.6) se obtiene un resultado similar.

J.res2.3 <- jtest(res2, res3)
J.res2.3

## J test
## 
## Model 1: y ~ x1 + x2 + x3
## Model 2: y ~ x4 + x5 + x8 + x9 + x10
##                 Estimate Std. Error t value  Pr(>|t|)    
## M1 + fitted(M2)  0.46880   0.093090  5.0360 1.391e-06 ***
## M2 + fitted(M1)  0.78589   0.090552  8.6788 8.131e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En conclusión, para esta muestra la prueba J no es concluyente al momento de comparar los modelos no anidados.

Ahora empleemos la prueba de Cox. Esta prueba se puede calcular empleando la función del paquete AER (Kleiber & Zeileis, 2008)). Esta función solo tiene dos argumentos, dos objetos de clase lm que contienen los dos modelos no anidados.

Cox.res1.3 <- coxtest(res1, res3)
Cox.res1.3

## Cox test
## 
## Model 1: y ~ x1 + x2 + x3 + x6 + x7
## Model 2: y ~ x4 + x5 + x8 + x9 + x10
##                 Estimate Std. Error  z value  Pr(>|z|)    
## fitted(M1) ~ M2  -20.255     4.4843  -4.5169 6.275e-06 ***
## fitted(M2) ~ M1  -44.422     4.1659 -10.6632 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cox.res2.3 <- coxtest(res2, res3)
Cox.res2.3

## Cox test
## 
## Model 1: y ~ x1 + x2 + x3
## Model 2: y ~ x4 + x5 + x8 + x9 + x10
##                 Estimate Std. Error  z value  Pr(>|z|)    
## fitted(M1) ~ M2  -24.191     4.6851  -5.1634 2.425e-07 ***
## fitted(M2) ~ M1  -47.896     4.1226 -11.6180 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La interpretación de la salida de esta prueba en R es similar a la prueba J. En este caso, pero no necesariamente siempre, los resultados de la prueba Cox son los mismos que obtuvimos con la prueba J (no siempre coinciden los resultados de estas dos pruebas). Por tanto, no es posible concluir en favor de un modelo.

De esta manera, uniendo todos los resultados encontramos que el (4.6) no es mejor que los (4.4) y (4.5). Este resultado lo obtenemos de las métricas de bondad de ajuste, pues con las pruebas Cox y J no podemos derivar una conclusión. Y entre el modelo (4.4) y (4.5), la prueba estadística y los dos criterios de información sugieren que el modelo (4.5) es mejor modelo. En el caso del \({\bar R^2}\) no existe una diferencia muy grande entre los dos modelos. Dado que las pruebas estadísticas nos permiten tener un 99% de confianza de que el modelo (4.5) es mejor que el modelo (4.4), por eso es aconsejable seleccionar como mejor modelo (4.5). Es decir, el mejor modelo será:

\[\begin{equation*} {y_i} = {\beta_1} + {\beta_2}{X_{1,i}} + {\beta_3}{X_{2,i}} + {\beta_4}{X_{3,i}} + {\varepsilon _i} \end{equation*}\]

4.5 Ejercicios

Las respuestas a estos ejercicios se encuentran en la sección 16.2.

4.5.1 Ejercicio práctico (continuación)

Continuemos con el ejercicio del final del Capítulo 2. En ese caso se estimó un modelo con todas las variables disponibles en la base de datos como posibles variables explicativas. No obstante, diferentes científicos de datos han llegado a diferentes posibles modelos para explicar los ingresos del sector (\(I\) medidos en millones de dólares). Estos posibles modelos son: \[\begin{equation} {I_t} = {\alpha _1} + {\alpha _2}{CE_t} + {\alpha _3}{CD_t} + {\alpha _4}{LDies_t} + {\alpha _5}{LEI_t} + {\alpha _6}{V_t} + {\varepsilon _t} \tag{4.7} \end{equation}\]

\[\begin{equation} {I_t} = {\alpha _1} + {\alpha _2}LDies_t + {\alpha _3}LEI_t + {\varepsilon _t} \tag{4.8} \end{equation}\] \[\begin{equation} {I_t} = {\alpha _1} + {\alpha _2}CE_t + {\alpha _3}LEI_t + {\varepsilon _t} \tag{4.9} \end{equation}\] \[\begin{equation} {I_t} = {\alpha _1} + {\alpha _2}CE_t + {\alpha _3}CD_t + {\alpha _4}{V_t} + {\varepsilon _t} \tag{4.10} \end{equation}\] \[\begin{equation} {I_t} = {\alpha _1} + {\alpha _2}{CE_t} + {\alpha _3}{CD_t} + {\alpha _4}{LDies_t} + {\alpha _5}{LEI_t} + {\varepsilon _t} \tag{4.11} \end{equation}\]

Tu tarea es determinar cuál es el mejor modelo para explicar los ingresos del sector ferroviario e interpretar los resultados. Emplea los mismos datos que se encuentran en el archivo regmult.csv.

Referencias

Kleiber, C., & Zeileis, A. (2008). Applied econometrics with r. Springer-Verlag. https://CRAN.R-project.org/package=AER

Este criterio se conoce también como criterio de información de Schwarz y se reconocen por las siguientes siglas que vienen del inglés: SIC, SBC o SBIC.↩︎
Diferentes al término constante.↩︎
Los datos se pueden descargar de la página web del libro: https://www.icesi.edu.co/editorial/modelo-clasico/↩︎
En caso de cometer un error y asignar como primer argumento el modelo sin restringir, observaras que los grados de libertad serán negativos al igual que la suma cuadrada. Esto claramente es imposible. Nota que en todo caso el F estadístico es el correcto así como su respectivo valor p.↩︎
Cuando emplees esta función hay que tener un poco de cuidado. Por defecto, esta función le asigna el nombre de “Model 1” al primer argumento (modelo con restricción) y “Model 2” al segundo argumento (modelo sin restricción). Claramente este es un rótulo arbitrario que es asignado por la función y nada tiene que ver que la numeración le hayas asignado a tus modelos en el análisis, como es nuestro caso en este ejemplo.↩︎