I have two linear models created with lm that I would like to compare with a table in the stargazer package. For the most part, I like the results I'm getting. But the Akaike Information Criterion is not showing. The docs say I can pass "aic" in the keep.stat argument to include it. But it's not there. No error messages.
stargazer(model1, model2, type="text", report="vc", header=FALSE,
title="Linear Models Predicting Forest Land",
keep.stat=c("aic", "rsq", "n"), omit.table.layout="n")
Linear Models Predicting Forest Land
==========================================
Dependent variable:
--------------------
forest
(1) (2)
------------------------------------------
log.MS.MIL.XPND.GD.ZS -11.948 -12.557
log.TX.VAL.AGRI.ZS.UN 2.310 2.299
log.NY.GDP.MKTP.CD 0.505
Constant 40.857 28.365
------------------------------------------
Observations 183 183
R2 0.142 0.146
==========================================
I don't see any reason why it wouldn't be able to include it. Calling the global AIC function on these models works fine.
> AIC(model1)
[1] 1586.17
> AIC(model2)
[1] 1587.208
The problem is given by the .AIC function defined inside stargazer:::.stargazer.wrap.
As one can see, this function does not calculate AIC for lm models:
.AIC <- function(object.name) {
model.name <- .get.model.name(object.name)
if (model.name %in% c("coeftest")) {
return(NA)
}
if (model.name %in% c("lmer", "lme", "nlme", "glmer",
"nlmer", "ergm", "gls", "Gls", "lagsarlm", "errorsarlm",
"", "Arima")) {
return(as.vector(AIC(object.name)))
}
if (model.name %in% c("censReg")) {
return(as.vector(AIC(object.name)[1]))
}
if (model.name %in% c("fGARCH")) {
return(object.name@fit$ics["AIC"])
}
if (model.name %in% c("maBina")) {
return(as.vector(object.name$w$aic))
}
if (model.name %in% c("arima")) {
return(as.vector(object.name$aic))
}
else if (!is.null(.summary.object$aic)) {
return(as.vector(.summary.object$aic))
}
else if (!is.null(object.name$AIC)) {
return(as.vector(object.name$AIC))
}
return(NA)
}
The .get.model.name function in .AIC calls .model.identify. If the component call of the model is lm(), then .model.identify returns ls:
if (object.name$call[1] == "lm()") {
return("ls")
}
Solution 1: Use add.lines.
set.seed(12345)
n <- 100
df <- data.frame(y=rnorm(n), x1=rnorm(n), x2=rnorm(n))
model1 <- lm(y ~ x1, data=df)
model2 <- lm(y ~ x2, data=df)
library(stargazer)
stargazer(model1, model2, type="text", report="vc", header=FALSE,
title="Linear Models Predicting Forest Land",
keep.stat=c("rsq", "n"), omit.table.layout="n",
add.lines=list(c("AIC", round(AIC(model1),1), round(AIC(model2),1))))
and the output is:
Linear Models Predicting Forest Land
=================================
Dependent variable:
--------------------
y
(1) (2)
---------------------------------
x1 0.115
x2 -0.052
Constant 0.240 0.243
---------------------------------
AIC 309.4 310.3
Observations 100 100
R2 0.011 0.002
=================================
Solution 2: Add the component AIC to model objects.
model1 <- lm(y ~ x1, data=df)
model2 <- lm(y ~ x2, data=df)
model1$AIC <- AIC(model1)
model2$AIC <- AIC(model2)
stargazer(model1, model2, type="text", report="vc", header=FALSE,
title="Linear Models Predicting Forest Land",
keep.stat=c("aic", "rsq", "n"), omit.table.layout="n")
and the output is
Linear Models Predicting Forest Land
======================================
Dependent variable:
--------------------
y
(1) (2)
--------------------------------------
x1 0.115
x2 -0.052
Constant 0.240 0.243
--------------------------------------
Observations 100 100
R2 0.011 0.002
Akaike Inf. Crit. 309.413 310.318
======================================