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)
