Comparar estadísticos de diversos modelos de regresión en R

Esta entrada forma parte de una función elaborada en el marco del curso de econometría computacional, del centro Empírica, y su objetivo es doble: i) crear una herramienta que permita comparar fácilmente los estadísticos asociados a los modelos de regresión estimados; e, ii) ilustrar como utilizar R para desarrollar funciones adaptables a nuestras necesidades (en el ejemplo se intenta explicar la lógica usada, aunque no pretende mostrar como crear estas desde cero, solo dar un ejemplo de uso, debido a que hay bastante material en línea para estos fines).

La función creada (ModeloStat) recibe como argumentos la lista de estimaciones que resultan de estimar un modelo de regresión en R (list<-lm(y~x)), especificado en puntos suspensivos dentro de la función (function(...){) para brindar flexibilidad en su uso y evitar se precise de un número fijos de modelos. Usando estos modelos, la función utiliza la función vectorizada para lapply para aplicar la función summary sobre cada uno de los elementos incluidos como argumentos en la función. En otras palabras, cada una de las salidas de regresión se guarda en una sola lista (list(…)) y luego se crea una lista llamada modelos1 donde se guarda el summary alcanzado por los modelos, entonces tendremos dos objetos tipo lista, donde estarán todos los modelos colocados en el argumento y desde donde se puedan extraer los estadísticos de regresión:

list(...)

lapply(list(...), summary)

> library(wooldridge)

> result1 <-lm(wage~educ+exper, data=wage1)

> result1

Call:

lm(formula = wage ~ educ + exper, data = wage1)

Coefficients:

(Intercept)         educ        exper 

    -3.3905       0.6443       0.0701 

summary(result1)

Call:

lm(formula = wage ~ educ + exper, data = wage1)

Residuals:

    Min      1Q  Median      3Q     Max

-5.5532 -1.9801 -0.7071  1.2030 15.8370

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept) -3.39054    0.76657  -4.423 1.18e-05 ***

educ         0.64427    0.05381  11.974  < 2e-16 ***

exper        0.07010    0.01098   6.385 3.78e-10 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.257 on 523 degrees of freedom

Multiple R-squared:  0.2252,         Adjusted R-squared:  0.2222

F-statistic: 75.99 on 2 and 523 DF,  p-value: < 2.2e-16

Posteriormente se utiliza la indexación para extraer elementos de cada modelo en la lista anteriores. En el caso de la lista creada con todos los resultados de las regresiones estimadas (list(result1,result2)) guardamos los resultados de los modelos estimados, luego se usa la función sapply para para aplicar la función AIC sobre cada uno de los objetos que componen la lista (es lo mismo decir, para cada uno de los modelos colocados en el argumento).

Llanamente, se guardan todos los modelos en una lista y se usa la función sapply para aplicar la función deseada sobre cada uno de estos:

1 –  

library(wooldridge)

result1 <-lm(wage~educ+exper, data=wage1)

result2 <-lm(log(wage)~educ+exper, data=wage1)

 

r_list <- list(result1,result2)

r_list

[[1]]

 

Call:

lm(formula = wage ~ educ + exper, data = wage1)

Coefficients:

(Intercept)         educ        exper 

    -3.3905       0.6443       0.0701 

[[2]]

Call:

lm(formula = log(wage) ~ educ + exper, data = wage1)

Coefficients:

(Intercept)         educ        exper 

    0.21685      0.09794      0.01035 

> sapply(r_list, FUN=AIC)

[1] 2739.9375  684.0187

2 –

En el caso de la lista donde se guarda el summary de cada uno de los modelos estimados (lapply(list(...), summary)), la indexación se crea mediante una función propia, que permita acceder a los elementos utilizando el operador $ (no una función como se hizo anteriormente). Primero guardamos el summary de cada modelo en la lista llamada modelos1.

modelos1<- lapply(r_list, summary)

modelos1

[[1]]

Call:

lm(formula = wage ~ educ + exper, data = wage1)

Residuals:

    Min      1Q  Median      3Q     Max

-5.5532 -1.9801 -0.7071  1.2030 15.8370

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept) -3.39054    0.76657  -4.423 1.18e-05 ***

educ         0.64427    0.05381  11.974  < 2e-16 ***

exper        0.07010    0.01098   6.385 3.78e-10 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.257 on 523 degrees of freedom

Multiple R-squared:  0.2252, Adjusted R-squared:  0.2222

F-statistic: 75.99 on 2 and 523 DF,  p-value: < 2.2e-16

 

[[2]]

Call:

lm(formula = log(wage) ~ educ + exper, data = wage1)

Residuals:

     Min       1Q   Median       3Q      Max

-2.05800 -0.30136 -0.04539  0.30601  1.44425

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept) 0.216854   0.108595   1.997   0.0464 * 

educ        0.097936   0.007622  12.848  < 2e-16 ***

exper       0.010347   0.001555   6.653 7.24e-11 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4614 on 523 degrees of freedom

Multiple R-squared:  0.2493, Adjusted R-squared:  0.2465

F-statistic: 86.86 on 2 and 523 DF,  p-value: < 2.2e-16

Ahora, para acceder a los componentes de esta segunda lista creamos una función propia que llame el estadístico deseado o que realice el cálculo en función de alguno de los elementos reservados en las listas de listas que hemos creado. Debajo un ejemplo con el caso del R2. Pero note que como esta función es personalizada, podemos asumir cálculos personalizados de estadísticos que no estén disponibles originalmente en las colecciones de listas guardadas.

sapply(modelos1, function(x){x$r.squared})

[1] 0.2251622 0.2493434

Ahora colocamos el código complete de la función.

ModeloStat <- function(...){

 modelos1<- lapply(list(...), summary)

tabla<-sapply(modelos1,

         function(x){

           rbind(length(x$residuals),

                 x$fstatistic[1],

                 x$r.squared,

                 x$adj.r.squared,

                 var(x$residuals))

           }

    )

 

tabla<-rbind(tabla,sapply(list(...), AIC))

tabla<-rbind(tabla,sapply(list(...), BIC))

 

rownames(tabla) <- c("Obs.","Fstat","R2","R2adj","ds_e","AIC","BIC")

colnames(tabla) <- c(paste("modelo ", 1:length(list(...)), sep=""))

 

return(tabla)

}

Finalmente, dejamos un ejemplo sobre el uso de la función, donde solo es necesario colocar dentro el nombre de los modelos que hemos estimados:

library(wooldridge)

result1 <-lm(wage~educ+exper, data=wage1)

result2 <-lm(log(wage)~educ+exper, data=wage1)

result3 <-lm(log(wage)~log(1+educ)+log(1+exper), data=wage1)

Entradas populares de este blog

Valores perdidos (NA) en R: identificación y tratamiento (I)

Ejemplos de tablas en Stata

Métodos y técnicas de Investigación Económica