I get different AIC values when using svyglm()$aic as opposed to using stats::AIC(model). I understand that topic has been discussed here before. However, there was no detailed explanation for the difference—”only” a suggestion to use one method over the other. I suspect it has something to do with this, but unfortunately my statistics is not strong enough to fully understand and draw conclusions.
I would appreciate it if somebody could give me a hand, telling me which is the correct way (if there is one) and why. As shown in the example below, the two methods can yield quite different results, even changing the interpretation of which model seems to have a better fit.
Minimal example:
library(survey)
df <- data.frame(
y = c(1.2, 2.4, 3.1, 4.5, 5.6, 6.3, 8.6),
a = c(0.5, 1.3, 1.7, 2.2, 3.1, 3.8, 4.0),
b = c(2.1, 2.7, 3.3, 3.8, 4.2, 4.5, 5.1),
weights = c(0.7, 1.4, 0.9, 1.1, 1.0, 0.8, 1.2)
)
weighted <- svydesign(ids = ~1, data = df, weights = ~weights)
A <- svyglm(y ~ a + b, design = weighted, family = "gaussian")
B <- svyglm(y ~ a, design = weighted, family = "gaussian")
AIC(A,B)[,2]
[1] 19.29119 17.03959
> c(A$aic,B$aic)
[1] 15.06302 17.45406
I am using R 4.4.1 and survey 4.4-2
$endgroup$
The short answer is that you should always* use accessor methods such as AIC() rather than extracting elements from objects using $ or (especially) @. There is no guarantee that an element labeled $aic is actually an AIC value, or that the definition will stay stable across package releases. You should assume that the package authors know what they’re doing and have provided the best possible version of the AIC (or whatever) in the method that comes with the package. (stats::AIC() will call the specific method provided in the package for objects of a particular class, if it exists.)
The help page for ?AIC.svyglm reads:
‘AIC’ is defined using the Rao-Scott approximation to the weighted
loglikelihood (Lumley and Scott, 2015). It replaces the usual
penalty term p, which is the null expectation of the log
likelihood ratio, by the trace of the generalised design effect
matrix, which is the expectation under complex sampling. For
computational reasons everything is scaled so the weights sum to
the sample size.
If you want to dig down and see exactly what AIC() is doing in this case you will eventually get to the code of survey:::extractAIC_svylm.
To understand exactly what the $aic slot contains you’d have to dig around in the survey package code (or ask the maintainer), but you are not supposed to have to know what this quantity represents, as it is an internal component of the fitted model object.
One possibly confusing ‘feature’ is that when printing a svyglm object (which has class c("svyglm", "glm", "lm")), survey-specific machinery is used to print out the details of the survey design etc., but R then falls through to the standard print.glm() method, which prints the $aic component (!!). This seems like an ‘infelicity’ (== arguably a bug) in the package machinery, as the results of AIC() should be more reliable. It may be worth contacting the package maintainer to see if this can be changed (I don’t see an official bug-reports channel listed on the package’s CRAN page)
* Happy to hear counterexamples in the comment section.
$endgroup$
2
You should not use $aic — as @BenBolker says, you shouldn’t even know it’s there. It’s a byproduct of using glm to fit the model, and it is meaningless — the model isn’t fitted by maximum likelihood, so standard approaches don’t even let you define AIC meaningfully. If you use the usual AIC formula, the result depends on the scale of the weights, and can be massively misleading.
You should use AIC(), which calls extractAIC.svyglm(), which returns a version of AIC that Alastair Scott and I invented for models using weighted estimation in survey data, and which is documented.
I note that stargazer doesn’t use methods for anything, really, and relies on being able to pull the right bits of information from the undocumented internals of the objects. They do that even when there are very standard accessors like coef and vcov dating from well before R. The problem with doing that is you need to understand what you’re extracting, and in the case of the $aic component of a svyglm it looks like they didn’t understand.
In the next version of the package, the $aic component will be NA to try to stop people using it.
$endgroup$
2