# Reconciling the survey adjusted AIC for multiple regression models in the R survey package version 4.1-1

I am unable to reconcile the new survey adjusted AIC calculation (dAIC) used for multiple regression in version 4.1-1 of the survey package on CRAN for modelling data with complex survey weights. I have reached out to the author, Thomas Lumley, but haven't received a response.

The survey design adjusted dAIC is given in the paper of Lumley and Scott (2015)

$$d\text{AIC} = -2n \widehat \ell (\widehat{\boldsymbol \theta}) + 2 p \widehat{\overline \delta}$$

In the most recent version of survey, the dAIC is computed for Multiple regression directly. I have derived the equations from the accompanying paper but cannot determine why there are two extra operations done in the survey package.

1. Why is the final result divided by $$\widehat \sigma^2$$ in https://github.com/cran/survey/blob/9addea4be893748c60974732b01ef7a393019217/R/dAIC.R#L121? This affects the estimates of the AIC, $$\overline \delta$$ and the design adjusted effective df.
2. Why are the intercept elements removed from the covariance matrices that are used to estimate $$\overline \delta$$ in lines https://github.com/cran/survey/blob/9addea4be893748c60974732b01ef7a393019217/R/dAIC.R#L100-L103?

# Derivation of result for Multiple regression model

Assuming $$\boldsymbol\theta = (\beta_0, \beta_1, \ldots,\beta_p, \sigma^2)^T = (\boldsymbol{\beta}^T, \sigma^2)^T$$ with the model

$$\boldsymbol{y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol\varepsilon \qquad \text{where}\qquad \boldsymbol\varepsilon \sim \mathcal N(\boldsymbol 0, \sigma^2 I)$$ or for a single observation

$$y_i = \boldsymbol x_i^T \boldsymbol \beta + \varepsilon_i$$ which has likelihood $$f_{\boldsymbol\theta} (y|\boldsymbol x) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left\{ - \frac{(y - \boldsymbol x^T \boldsymbol\beta)^2}{2\sigma^2} \right\}.$$

Then the weighted log-likelihood is

$$\widehat \ell(\boldsymbol \theta) = \frac{1}{N} \sum_{i \in s} w_i \ell_i(\boldsymbol \theta)$$

where the weights in the paper are scaled so that $$N = \sum_{i \in s}w_i$$ and the log-likelihood is $$\ell_i(\boldsymbol \theta) = \log f_{\boldsymbol\theta} (y_i|\boldsymbol x_i)$$. For convenience, define the diagonal matrix $$\boldsymbol W = \text{diag}(w_i)$$. Then this leads to, $$\sum_{i \in s} w_i \ell_i(\boldsymbol \theta) = - \frac{N}{2} \log(2\pi) - \frac{N}{2} \log(\sigma^2) - \frac{\boldsymbol \varepsilon^T \boldsymbol W \boldsymbol \varepsilon}{2\sigma^2},$$ with score vector $$\widehat{\boldsymbol U} = \frac{\partial \widehat \ell (\boldsymbol\theta)}{\partial \boldsymbol \theta} = \begin{pmatrix} \frac{1}{\sigma^2} \boldsymbol X^T \boldsymbol W \boldsymbol \varepsilon \\ - \frac{N}{2\sigma^2} + \frac{\boldsymbol \varepsilon^T \boldsymbol W \boldsymbol \varepsilon}{2 (\sigma^2)^2}\end{pmatrix}$$ Solving $$\widehat{\boldsymbol U} = \boldsymbol 0$$ gives the expected estimates $$\widehat{\boldsymbol \beta} = (\boldsymbol X^T \boldsymbol W \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol W \boldsymbol Y$$ and $$\widehat \sigma^2 = \frac{1}{N}(\boldsymbol Y - \boldsymbol X \widehat{\boldsymbol \beta})^T \boldsymbol W (\boldsymbol Y - \boldsymbol X \widehat{\boldsymbol \beta})$$

In the survey package the weights are further chosen such that $$w_i \mapsto \frac{w_i}{\overline w} \Rightarrow N = n$$ which has the consequence that $$-2n\widehat \ell(\widehat{\boldsymbol \theta}) =\sum_{i \in s} w_i \ell_i(\boldsymbol \theta)= N + N \log(2\pi) + N \log(\widehat \sigma^2)$$ which validates minus2ellhat on https://github.com/cran/survey/blob/9addea4be893748c60974732b01ef7a393019217/R/dAIC.R#L96:w

However, the final line in that function divides minus2ellhat (and the other terms) by sigmahat2. Why?

The design adjusted penalty also has the step of removing the intercept. The penalty term of the dAIC is, $$2 p \widehat{\overline \delta}$$ where $$p$$ is the number of parameters in the model and $$\widehat{\overline{\delta}}$$ denotes the average design affect. The average design effect is defined $$\overline \delta = tr(\boldsymbol \Delta)/p$$ where $$\boldsymbol \Delta = \boldsymbol{\mathcal I}(\boldsymbol \theta) \boldsymbol V(\boldsymbol \theta)$$ is known as the design effect, the ratio of covariance matrices with and without sampling weights ($$\boldsymbol{\mathcal I}(\boldsymbol \theta)$$ is the Fisher information matrix for the unweighted model and $$\boldsymbol V(\boldsymbol \theta)$$ being the asymptotic covariance matrix of the parameter estimates in the weighted model).

The Fisher information is defined $$\boldsymbol{\mathcal I}(\boldsymbol \theta) = \begin{pmatrix} \frac{1}{\sigma^2} \boldsymbol X^T\boldsymbol X & \boldsymbol 0\\ \boldsymbol 0^T & \frac{n}{2(\sigma^2)^2}\end{pmatrix}$$

The design effect of the slope parameters are extracted from the two appropriate covariance matrices fitted using svyglm but there is a step where the intercept components are removed https://github.com/cran/survey/blob/9addea4be893748c60974732b01ef7a393019217/R/dAIC.R#L100-L103. If the null model is the intercept model itself, this will cause an error since the matrices in this step will be empty and cannot be inverted.

The other calculations seem correct and the design effect of the $$\sigma^2$$ component is estimated correctly using the second moment of the score vector assuming the survey weighted probabilities and using the Fisher information above in the unweighted case.

# References

[1] Lumley, T., and A. Scott. “AIC and BIC for Modeling with Complex Survey Data.” Journal of Survey Statistics and Methodology 3, no. 1 (2015): 1–18. https://doi.org/10.1093/jssam/smu021.

• Looking at this with fresh eyes. The V0 matrix defined in the R code is the naive.cov element from the svyglm model which itself is the cov.unscaled from a standard glm. That is, it is $(X^TX)^{-1}$ and missing the $\widehat{\sigma}^2$ factor. So this division by sigma2 is necessary for the regression slope parameter part (including intercept) for the calculation of $dAIC$. It shouldn't be done everywhere as is currently done in survey though. Also unable to verify why the intercept elements (first row & col) is removed from both V and V0 in the new survey implementation too.
– JRW
Commented Feb 20, 2023 at 2:11

The $$d$$AIC metric for multiple regression was calculated directly in the survey package in version 4.1-1 here.

Note in the implementation in survey the rescaled weights $$w_i^*$$ are used so that $$\sum_{i \in s} w_i^* = N^* = n$$. That is, $$w_i^* = \frac{w_i}{\overline w}$$ where $$\overline w = \frac{1}{n}\sum_{i \in s} w_i$$ which has the implication that $$N = N^* = n$$.

Recall that the $$d$$AIC has two components, the likelihood component and the penalty component.

$$$$d\text{AIC} = -2n \widehat{\ell}(\widehat{\boldsymbol{\theta}}) + 2p\widehat{\overline{\delta}}$$$$

There seem to be two issues with the implementation in survey. First, the penalty component is not correctly computed for the $$\sigma^2$$ contribution. Second, both the likelihood and penalty components are divided by $$\widehat{\sigma^2}$$ which is not necessary.

1. The penalty component is computed in the extractAIC_svylm function at this point. However, the sigma component in the code this seems to be the reciprocal of the theory. See \eqref{deltabar} and the equations following it. In particular, reconcile $$\mathcal I_{\sigma^2}$$ and $$V_{\sigma^2}$$ in the theory with Isigma2 and Hsigma2 in the code. In addition to this, the penalty component (and likelihood component) is further modified with a division by $$\widehat{\sigma^2}$$ on the final line.
2. The likelihood component is correctly computed in survey at this point of the code in the body of the extractAIC_svylm function as it agrees with \eqref{ellhat}. However, the return statement on the last line divides the AIC by $$\widehat \sigma^2$$.

The theory to justify this is below and a suggested code fix below that.

# Weighted likelihood for multiple regression under sampling design

In the complex survey design approach, the observations $$\left\{ (y_i, \boldsymbol{x}_i): i \in s\right\}$$ are assumed to be derived from a sample $$s$$ of $$n$$ units drawn from a finite population of size $$N$$ using some probability sampling design. The weights associated with the $$i^{\mathrm{th}}$$ unit from the population are defined to be $$w_i$$.

The data from the finite population are assumed to be generated from some distribution with joint density, $$g(y, \boldsymbol{x})$$. The desired modelling approach for the data is assumed to be a parametric model, $$\left\{ f_{\boldsymbol{\theta}}(y|\boldsymbol{x}), \boldsymbol{\theta} \in \boldsymbol{\Theta}\right\}$$. Generally speaking the parametric family induced by $$f_{\boldsymbol{\theta}}(y|x)$$ does not need to contain the true distribution in $$g(y, x)$$ or $$g(y|x)$$. The parametric model would be considered misspecified in this case. The likelihood approach can nevertheless still be used and the objective function being maximised is the quantity, $$$$\ell(\boldsymbol{\theta}) = \mathbb{E}_g \left( \log f_{\boldsymbol{\theta}}(y|\boldsymbol{x})\right).$$$$

This is motivated by the using the Kullback-Leibler divergence of the candidate parametric model away from the true model, $$$$\mathrm{KL}(f_{\boldsymbol{\theta}}, g) = \mathbb{E}_g \left( \log g(y| \boldsymbol{x})\right) - \ell(\boldsymbol{\theta}).$$$$

The target parameter is $$\boldsymbol{\theta}^*$$ is the parametric value that is the 'least false' which is the value that minimizes the Kullback-Leibler divergence above (or equivalently, maximizer of $$\ell(\boldsymbol{\theta}))$$. In terms of the classical likelihood equations, this occurs at the new score equation,

$$$$\mathcal U(\boldsymbol{\theta}^*) = \mathbb{E}_g \left( \left.\frac{\partial \log f_{\boldsymbol{\theta}}(y|\boldsymbol{x})}{\partial \boldsymbol{\theta}}\right|_{\boldsymbol{\theta} = \boldsymbol{\theta}^*}\right) = 0.$$$$

The likelihood equation, $$\ell(\boldsymbol{\theta})$$, being a population mean parameter is estimated using the weighted estimator, $$$$\widehat{\ell}(\boldsymbol{\theta}) = \frac{1}{N} \sum_{i \in s} w_i \ell_i(\boldsymbol{\theta})$$$$ where $$\ell_i(\boldsymbol{\theta}) = \log f_{\boldsymbol{\theta}}(y_i|\boldsymbol{x})$$ and the weights are assumed to aggregate to the population total, $$\sum_{i \in s}w_i = N$$. For the multiple regression model, each observation is assumed to follow the log-likelihood, $$$$\ell_i(\boldsymbol{\theta}) = \log f_{\boldsymbol{\theta}}(y_i|\boldsymbol{x}_i) = - \frac{1}{2} \log(2 \pi) - \frac{1}{2}\log(\sigma^2) - \frac{(y_i - \boldsymbol{x}_i^T \boldsymbol{\beta})^2}{2 \sigma^2},$$$$ where the parameter space is defined $$\boldsymbol{\theta} = \left( \boldsymbol{\beta}^T, \sigma^2\right) \in \mathbb{R}^{m+1} \cup \left\{\mathbb{R}_+\setminus 0\right\}$$ and $$\boldsymbol{x}_i$$ is a vector of covariates for observation $$i$$. To estimate the maximizer of $$\ell(\boldsymbol{\theta})$$, use an estimator that solves the score equation, $$$$\widehat{\mathcal U}(\boldsymbol{\theta}) = \frac{1}{N} \sum_{i \in s} w_i \frac{\partial}{\partial \boldsymbol{\theta}}\ell_i(\boldsymbol{\theta}).$$$$

Then $$\boldsymbol{\theta}^*$$ is consistently estimated by $$\widehat{\boldsymbol{\theta}}$$ that satisfies the weighted pseudo-score equations, $$\widehat{\mathcal U}(\widehat{\boldsymbol{\theta}}) = 0$$. This is shown later using results of Fuller (2009).

Before computing the required terms in this analysis, the link between $$\widehat{\boldsymbol{\theta}}$$ and the AIC strategy of model assessment is covered. The model that matches the population likelihood the best will minimize the Kullback-Leibler divergence above. The AIC strategy is to estimate this difference by estimating $$Q_n = \mathbb{E}_g( \ell(\widehat{\boldsymbol{\theta}}))$$ and selecting the model with the largest value.

Theorem 1 of Lumley and Scott (2015) showed the statistical properties of the naive estimator, $$\widehat{\ell}(\widehat{\boldsymbol{\theta}})$$, with, $$$$\mathbb{E}_g(\widehat{\ell}(\widehat{\boldsymbol{\theta}})) = Q_n + \frac{1}{n} \mathrm{tr} \left\{ \boldsymbol{\Delta} \right\} + o_p(n^{-1})$$$$ where $$\boldsymbol{\Delta} = \mathcal I(\boldsymbol{\theta}^*)\boldsymbol{V}(\boldsymbol{\theta}^*)$$ and $$\mathcal I(\boldsymbol{\theta}) = \mathbb{E}\widehat{\mathcal J}(\boldsymbol{\theta})$$ and $$\boldsymbol{V}(\boldsymbol{\theta}^*)$$ being the asymptotic covariance matrix of $$\sqrt{n}\widehat{\boldsymbol{\theta}}$$.

In particular, their proof utilizes the regularity conditions and results of Theorem 1.3.9 Fuller (2009) under the typical framework in survey estimation. That is, it is supposed that there are a sequence of finite populations indexed by $$\nu$$ and a sequence of samples of size $$n_\nu$$ drawn from the $$N_\nu$$ units in the $$\nu^{\mathrm th}$$ population with a well-defined probability scheme. Then asymptotic normality and consistency is established with

1. $$\widehat{\boldsymbol{\theta}}_\nu$$ is asymptotically normal with $$\sqrt{n} (\widehat{\boldsymbol{\theta}}_{\nu}- \boldsymbol{\theta}^*) \stackrel{d}{\longrightarrow} \mathcal N(\boldsymbol{0}, \boldsymbol{V}(\boldsymbol{\theta}^*))$$.
2. If a sequence $$\left\{\boldsymbol{\theta}_\nu \right\}$$ are consistent estimators of $$\boldsymbol{\theta}^*$$, that is, $$\boldsymbol{\theta}_\nu \stackrel{p}{\longrightarrow} \boldsymbol{\theta}^*$$, then $$\widehat{\mathcal J}_\nu(\boldsymbol{\theta}_\nu) \stackrel{p}{\longrightarrow} \mathcal I(\boldsymbol{\theta})$$.

Consider now the precise form of the likelihood equations in this weighted sampling design context. For convenience define $$\boldsymbol{W} = \mathrm{diag}(w_i)_{i \in s}$$ as the diagonal matrix of weights.

The log-likelihood takes the form, \begin{aligned} \widehat{\ell}(\boldsymbol{\theta}) &= - \frac{1}{2} \log(2 \pi) - \frac{1}{2} \log \sigma^2 - \frac{1}{N}\frac{(\boldsymbol{Y} - \boldsymbol{X}\boldsymbol{\beta})^T\boldsymbol{W}(\boldsymbol{Y} - \boldsymbol{X}\boldsymbol{\beta})}{2\sigma^2}\\ &= - \frac{1}{2} \log(2 \pi) - \frac{1}{2} \log \sigma^2 - \frac{1}{N}\frac{\boldsymbol{\varepsilon}^T\boldsymbol{W}\boldsymbol{\varepsilon}}{2\sigma^2} \end{aligned}\tag{1}\label{ellhat} Maximising the log likelihood function in a similar manner to the unweighted case gives the weighted score vector $$\widehat{\mathcal{U}} = \frac{\partial \widehat{\ell}}{\partial \boldsymbol{\theta}} = \frac{1}{N} \sum_{i \in s} w_i \frac{\partial l_i(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}$$, and corresponding negative Hessian $$\widehat{\mathcal J}(\boldsymbol{\theta}) = - \frac{\partial^2 \widehat{\ell}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^T} = - \frac{1}{N} \sum_{i \in s} w_i \frac{\partial^2\ell_i(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}\partial\boldsymbol{\theta}^T}$$. Evaluating these leads to the equations, $$$$\widehat{\mathcal{U}}(\boldsymbol{\theta}) = \begin{pmatrix} \boldsymbol{X}^T\boldsymbol{W} \boldsymbol{\varepsilon}/ (N\sigma^2)\\ - \frac{1}{2N\sigma^2} (N - \boldsymbol{\varepsilon}^T\boldsymbol{W}\boldsymbol{\varepsilon}/(\sigma^2))\end{pmatrix} = \frac{1}{N} \sum_{i \in s} w_i \begin{pmatrix} \boldsymbol{x}_i^T \frac{(y_i - \boldsymbol{x}_i^T\boldsymbol{\beta})}{2\sigma^2}\\ - \frac{1}{2\sigma^2} + \frac{(y_i - \boldsymbol{x}_i^T \boldsymbol{\beta})^2}{2 (\sigma^2)^2}\end{pmatrix} =: \begin{pmatrix} U_{\boldsymbol{\beta}}\\ U_{\sigma^2}\end{pmatrix}.$$$$ $$$$\widehat{\mathcal J}(\boldsymbol{\theta}) = \begin{pmatrix} \frac{\boldsymbol{X}^T \boldsymbol{W} \boldsymbol{X}}{N\sigma^2} & \frac{\boldsymbol{X}^T \boldsymbol{W} \boldsymbol{\varepsilon}}{N(\sigma^2)^2} \\ \frac{\boldsymbol{\varepsilon}^T \boldsymbol{W} \boldsymbol{X}}{N(\sigma^2)^2} & -\frac{1}{2(\sigma^2)^2} + \frac{\boldsymbol{\varepsilon}^T \boldsymbol{W}\boldsymbol{\varepsilon}}{N(\sigma^2)^3} \end{pmatrix}, \qquad \mathcal J_N(\boldsymbol{\theta}) = \begin{pmatrix} \frac{\boldsymbol{X}^T \boldsymbol{X}}{N\sigma^2} & \frac{\boldsymbol{X}^T \boldsymbol{\varepsilon}}{N(\sigma^2)^2} \\ \frac{\boldsymbol{\varepsilon}^T \boldsymbol{X}}{N(\sigma^2)^2} & -\frac{1}{2(\sigma^2)^2} + \frac{\boldsymbol{\varepsilon}^T \boldsymbol{\varepsilon}}{N(\sigma^2)^3} \end{pmatrix};$$$$ where $$\mathcal J_N(\boldsymbol{\theta})$$ is the negative Hessian without the weights applied which is used in later results.

Solving the pseudo-score equation above yields the pseudo maximum likelihood parameter estimates, $$$$\widehat{\boldsymbol{\beta}} = \left(\boldsymbol{X}^T\boldsymbol{W}\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{W}\boldsymbol{Y}, \qquad \widehat{\sigma}^2 = \frac{1}{N} \boldsymbol{\varepsilon}^T\boldsymbol{W}\boldsymbol{\varepsilon}.$$$$ From (1.3.78) in Fuller (2009) $$\mathcal I(\boldsymbol{\theta}) = \mathbb{E}\widehat{\mathcal J}(\boldsymbol{\theta})$$ can be estimated by $$\mathcal J(\boldsymbol{\theta})$$, $$$$\widehat{\mathcal J}(\boldsymbol{\theta}) = -\frac{1}{N}\sum_{i \in s} w_i \frac{\partial \ell_i(\boldsymbol{\theta})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^T} = -\frac{1}{N}\sum_{i \in s} \frac{\partial \ell_i(\boldsymbol{\theta})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^T} + O_p(n^{-1/2}) =: \mathcal J_N(\boldsymbol{\theta}) + O_p(n^{-1/2}) \stackrel{a.s.}{\longrightarrow} \mathcal J(\boldsymbol{\theta})$$$$ and $$\mathcal J(\boldsymbol{\theta})$$ is non-singular.

## Estimating the design effect matrix, $$\boldsymbol{\boldsymbol{\Delta}}$$

Recall the design effect matrix term $$\boldsymbol{\boldsymbol{\Delta}} = \mathcal I(\boldsymbol{\theta}^*) \boldsymbol{V}(\boldsymbol{\theta}^*)$$ where $$\boldsymbol{V}(\boldsymbol{\theta}^*)$$ denotes the asymptotic covariance matrix of $$\sqrt{n}\widehat{\boldsymbol{\theta}}$$. The estimate of $$\mathcal I(\boldsymbol{\theta})$$ is considered first. $$$$\mathcal I(\boldsymbol{\theta}) = \mathbb{E}\widehat{\boldsymbol{\mathcal J}}(\boldsymbol{\theta}) \Rightarrow \widehat{\mathcal I}(\widehat{\boldsymbol{\theta}}) = \begin{pmatrix} \frac{\boldsymbol{X}^T\boldsymbol{X}}{N\widehat{\sigma}^2} & \boldsymbol{0} \\ \boldsymbol{0}^T & \frac{1}{2(\widehat{\sigma}^2)^2}\end{pmatrix} =: \begin{pmatrix} \mathcal I_{\boldsymbol{\beta}} & \boldsymbol{0}\\ \boldsymbol{0}^T & \mathcal I_{\sigma^2} \end{pmatrix},$$$$ where $$\boldsymbol{0}$$ is a column vector of zeros of dimension $$m +1$$. Note that the regression parameters, $$\boldsymbol{\beta}$$, and the variance nuisance parameter, $$\sigma^2$$, are orthogonal. That is, in the sense that the Information matrix is block diagonal.

The design effect adjustment is computed using $$\boldsymbol{\boldsymbol{\Delta}}$$ with, $$$$\overline{\delta} = \mathrm{trace}(\boldsymbol{\boldsymbol{\Delta}}) / p.$$$$ The blockwise diagonal form of $$\mathcal I(\boldsymbol{\theta})$$ allows the trace to be computed separately for the regression parameters and the variance nuisance parameter. Indeed, due to the form of the negative Hessian $$\widehat{\mathcal J}(\boldsymbol{\theta})$$, it follows that $$\boldsymbol{V}(\boldsymbol{\theta})$$ will be blockwise diagonal with, $$$$\boldsymbol{V}(\boldsymbol{\theta}) = \begin{pmatrix} V_{\boldsymbol{\beta}} & \boldsymbol{0} \\ \boldsymbol{0}^T & V_{\sigma^2} \end{pmatrix} \Rightarrow \mathrm{trace}(\boldsymbol{\boldsymbol{\Delta}}) = \mathrm{trace}(\mathcal I_{\boldsymbol{\beta}}\boldsymbol{V}_{\boldsymbol{\beta}}) + \mathrm{trace}(\mathcal I_{\boldsymbol{\sigma^2}}V_{\sigma^2}) = \mathrm{trace}(\mathcal I_{\boldsymbol{\beta}}\boldsymbol{V}_{\boldsymbol{\beta}}) +\mathcal I_{\boldsymbol{\sigma^2}}V_{\sigma^2}\tag{2}\label{deltabar}$$$$ The trace calculation is split into the two components. One for the regression parameters and one for the variance nuisance parameter. The trace of the regression parameters is estimated routinely using the standard variance covariance matrices of the regression parameters in $$\mathcal I_{\boldsymbol{\beta}}$$ and $$V_{\boldsymbol{\beta}}$$. That is, define $$\boldsymbol{V}_{\boldsymbol{\beta}, 0}$$ to be the variance covariance matrix of $$\boldsymbol{\beta}$$ under simple random design and the variance covariance matrix under sampling weights to be $$\boldsymbol{V}_{\boldsymbol{\beta}}$$. Then simply $$\mathcal I_{\boldsymbol{\beta}}\boldsymbol{V}_{\boldsymbol{\beta}} = \boldsymbol{V}_{\boldsymbol{\beta}, 0}^{-1}\boldsymbol{V}(\boldsymbol{\beta})$$

The $$\sigma^2$$ component of the calculation involves estimating the scalars $$\mathcal I_{\sigma^2}$$ and $$V_{\sigma^2}$$. The $$\mathcal I_{\sigma^2}$$ component can be estimated directly with $$\widehat{\mathcal{I}}_{\sigma^2} = 1 /(2(\widehat{\sigma}^2)^2)$$. To estimate $$V_{\sigma^2}$$, it can be done by estimating the $$\sigma^2$$ element in the Information matrix under the sampling weights and then taking the reciporical. To estimate the $$\sigma^2$$ element in the Information matrix, use the standard estimation equation as the variance of the score equation. Since the score equation is zero at the MLE it only requires computing the second moment of the score vector. The contribution of each observation to the $$\sigma^2$$ element in the score equation is defined, $$$$U_{\sigma^2, i}(\boldsymbol{\theta}) := \frac{\partial l_i(\boldsymbol{\theta})}{\partial \sigma^2} = - \frac{1}{2\sigma^2} + \frac{(y_i - \boldsymbol{x}_i^T \boldsymbol{\beta})^2}{2(\sigma^2)^2}.$$$$ Therefore an estimate of the second moment and consequently variance of this component since the score vector is zero at the MLE is given by, $$$$\widehat{\mathbb{V}\text{ar}}U_{\sigma^2} = \frac{1}{N}\sum_{i \in s} w_i U_{\sigma^2, i}^2(\widehat{\boldsymbol{\theta}}) = \frac{1}{N}\sum_{i \in s} w_i \left( - \frac{1}{2\widehat{\sigma}^2} + \frac{(y_i - \boldsymbol{x}_i^T \widehat{\boldsymbol{\beta}})^2}{2(\widehat{\sigma}^2)^2}\right)^2$$$$ Thus the estimate of the covariance element $$V_{\sigma^2}$$ can be estimated with, $$$$\widehat V_{\sigma^2} = \left(\widehat{\mathbb{V}\text{ar}} U_{\sigma^2}\right)^{-1}.$$$$ Finally $$\widehat{\overline{\delta}} = \widehat{\delta}_{\boldsymbol{\beta}} + \widehat{\delta}_{\sigma^2} = \mathrm{trace}(\widehat{\boldsymbol{V}}_{\boldsymbol{\beta}, 0}^{-1} \widehat{\boldsymbol{V}}_{\boldsymbol{\beta}}) + \widehat{\mathcal{I}}_{\sigma^2}\widehat{V}_{\sigma^2}$$.

# Code fix

Below includes a proposal fix using the theory and comparison in the other answer by @JRW.

The following calculations are setup assuming survey version 4.2-1 which has the first CRAN implementation for the direct dAIC calculation including the effects for estimating $$\sigma^2$$ as well as the $$\beta$$ coefficients in a survey weighted linear multiple regression model.

set.seed(2023)


Load VGAM to use their TIC calculation

library(survey, warn.conflicts = FALSE, quietly = TRUE)
library(VGAM, warn.conflicts = FALSE, quietly = TRUE)
library(knitr)


To compute old survey output and alternative proposal

require(testthat, warn.conflicts = FALSE, quietly = TRUE)
stopifnot(packageVersion('survey') == "4.2.1")


proposed function

extractSvyLmAIC <- function(fit, k = 2) {
# Extract weights used in survey
w <- fit[["prior.weights"]]
n.hat <- sum(w)
y <- fit[["y"]]
mu.hat <- fit[["linear.predictors"]]
eps <- y - mu.hat
# Use ML estimate of sigma2
sigma2.hat <- sum(eps^2 * w) / n.hat
# Compute the likelihood component of the AIC
minus.2.ell.hat <- n.hat * (log(sigma2.hat) + 1 + log(2 * pi))
# Construct the beta contribution to the design effect matrix
v.beta.zero <- fit[["naive.cov"]] * sigma2.hat
v.beta <- vcov(fit)
# Compute the delta matrix for the regression coefficients
delta.beta.matrix <- solve(v.beta.zero, v.beta)
# Compute the sigma2 component of the delta matrix
## Information matrix for sigma2 in unweighted case
i.sigma2 <- 1 / (2 * sigma2.hat^2)
## Estimate the covariance of sigma2 under sampling weights
## Use the score equation estimator
u.sigma2.i <- -1 / (2 * sigma2.hat) + eps^2 / (2 * sigma2.hat^2)
var.sigma2 <-  1 / mean(w * u.sigma2.i^2)
delta.sigma2 <- i.sigma2 * var.sigma2
# Compute the overall design effects
delta.beta <- diag(delta.beta.matrix)
delta.bar <- mean(c(delta.beta, delta.sigma2))
eff.p <- sum(delta.beta, delta.sigma2)
# Compute the dAIC = -2LL + k * p * deltabar
aic <- minus.2.ell.hat + k * eff.p
c(eff.p = eff.p, AIC = aic, deltabar = delta.bar)
}


Number of decimal places in output

ndp <- 2L


Basic check with unit weights For a first test of the dAIC code, a simple linear regression model is considered where the outcome and predictors have no relationship whatsoever and are pure gaussian noise with correlation matrix and matrix scatterplot shown below.

n <- 100L
dat <- replicate(3L, rnorm(n), simplify = FALSE) |>
as.data.frame() |>
setNames(c("y", "x1", "x2"))


Inspect the various AIC, dAIC and TIC calculations by first fitting the appropriate models Shown below is the various AIC calculations, including the standard AIC generated for regular lm/glm models, along with the TIC for the vglm models. The $$d$$AIC calculations proposed by survey 4.0and 4.2-1 (current as of 11/07/23) are shown along with the newly proposed correction in extractSvyLmAIC function included in this document. To generate the [email protected] values, the code uses the survey 4.2-1 package intercepting the check that the model is a gaussian multiple regression model forcing it to compute the dAIC via the glm framework

#" as was done in [email protected].
numPredictors <- function(model) {
if (inherits(model, "vglm"))
model <- model@misc
attr(terms(model$formula), "term.labels") |> length() } dfOrder <- function(df) { desired.order <- c("value", "predictors", "calculation", "df", "delta") df <- df[match(colnames(df), desired.order)] }  Intercept the is.svylm function to force the [email protected] calculation s40AIC <- function(model) { aic <- testthat::with_mock(is.svylm = function(x) FALSE, AIC(model), .env = "survey") aic <- data.frame(value = aic[["AIC"]], df = aic[["eff.p"]], delta = aic[["deltabar"]]) aic[["calculation"]] <- "survey v4.0 (dAIC)" aic[["predictors"]] <- numPredictors(model) aic }  Current survey calc sAIC <- function(model) { aic <- AIC(model) aic <- data.frame(value = aic[["AIC"]], df = aic[["eff.p"]], delta = aic[["deltabar"]]) aic[["calculation"]] <- "survey (dAIC)" aic[["predictors"]] <- numPredictors(model) aic }  Proposed survey calc dAIC <- function(model) { aic <- extractSvyLmAIC(model) aic <- data.frame(value = aic[["AIC"]], df = aic[["eff.p"]], delta = aic[["deltabar"]]) aic[["calculation"]] <- "proposal (dAIC)" aic[["predictors"]] <- numPredictors(model) aic } calcTIC <- function(model) { aic <- data.frame(value = TICvlm(model)) aic[["df"]] <- NA aic[["delta"]] <- NA aic[["calculation"]] <- "vglm (TIC)" aic[["predictors"]] <- numPredictors(model) aic } glmAIC <- function(model) { aic <- AIC(model) n.preds <- numPredictors(model) aic <- data.frame(value = aic[[1]], df = n.preds + 2) aic[["delta"]] <- NA aic[["calculation"]] <- "glm (AIC)" aic[["predictors"]] <- n.preds aic } emptyDF <- function(calculation) { data.frame(value = NA, df = NA, delta = NA, calculation = calculation, predictors = NA) } aicCalcs <- function(models) { output <- lapply(models, function(x) { if (inherits(x, "svyglm")) { proposeddAIC <- dAIC(x) s40 <- tryCatch(s40AIC(x), error = function(e) emptyDF("[email protected] (dAIC)")) s41 <- tryCatch(sAIC(x), error = function(e) emptyDF("survey (dAIC)")) return(rbind.data.frame(proposeddAIC, s40, s41)) } if (inherits(x, "vglm")) { return(calcTIC(x)) } glmAIC(x) }) output <- do.call(rbind, output) output <- dfOrder(output) output } unit.design <- svydesign(ids = ~1, weights = ~1, data = dat)  Fitting survey lm models with unit weights m0 <- svyglm(y ~ x1, design = unit.design) summary(m0) #> #> Call: #> svyglm(formula = y ~ x1, design = unit.design) #> #> Survey design: #> svydesign(ids = ~1, weights = ~1, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 0.07124 0.09692 0.735 0.464 #> x1 -0.01786 0.09958 -0.179 0.858 #> #> (Dispersion parameter for gaussian family taken to be 0.9857939) #> #> Number of Fisher Scoring iterations: 2 m1 = svyglm(y ~ x1 + x2, design = unit.design) summary(m1) #> #> Call: #> svyglm(formula = y ~ x1 + x2, design = unit.design) #> #> Survey design: #> svydesign(ids = ~1, weights = ~1, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 0.069979 0.098790 0.708 0.480 #> x1 -0.017611 0.100177 -0.176 0.861 #> x2 -0.007479 0.101843 -0.073 0.942 #> #> (Dispersion parameter for gaussian family taken to be 0.9857391) #> #> Number of Fisher Scoring iterations: 2 ## base R regular glm models with glm o0 <- glm(y ~ x1, data = dat) summary(o0) #> #> Call: #> glm(formula = y ~ x1, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 0.07124 0.10226 0.697 0.488 #> x1 -0.01786 0.10710 -0.167 0.868 #> #> (Dispersion parameter for gaussian family taken to be 0.995853) #> #> Null deviance: 97.621 on 99 degrees of freedom #> Residual deviance: 97.594 on 98 degrees of freedom #> AIC: 287.35 #> #> Number of Fisher Scoring iterations: 2 o1 <- glm(y ~ x1 + x2, data = dat) summary(o1) #> #> Call: #> glm(formula = y ~ x1 + x2, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 0.069979 0.104201 0.672 0.503 #> x1 -0.017611 0.107701 -0.164 0.870 #> x2 -0.007479 0.101869 -0.073 0.942 #> #> (Dispersion parameter for gaussian family taken to be 1.006064) #> #> Null deviance: 97.621 on 99 degrees of freedom #> Residual deviance: 97.588 on 97 degrees of freedom #> AIC: 289.35 #> #> Number of Fisher Scoring iterations: 2 ## VGAM computations with vglm to compute TIC vm0 <- vglm(y ~ x1, data = dat, family = uninormal) summary(vm0) #> #> Call: #> vglm(formula = y ~ x1, family = uninormal, data = dat) #> #> Coefficients: #> Estimate Std. Error z value Pr(>|z|) #> (Intercept):1 0.07124 0.10124 0.704 0.482 #> (Intercept):2 -0.01218 0.07071 -0.172 0.863 #> x1 -0.01786 0.10603 -0.168 0.866 #> #> Names of linear predictors: mean, loglink(sd) #> #> Log-likelihood: -140.6759 on 197 degrees of freedom #> #> Number of Fisher scoring iterations: 3 #> #> No Hauck-Donner effect found in any of the estimates vm1 <- vglm(y ~ 1 + x1 + x2, data = dat, family = uninormal) summary(vm1) #> #> Call: #> vglm(formula = y ~ 1 + x1 + x2, family = uninormal, data = dat) #> #> Coefficients: #> Estimate Std. Error z value Pr(>|z|) #> (Intercept):1 0.069979 0.102634 0.682 0.495 #> (Intercept):2 -0.012207 0.070711 -0.173 0.863 #> x1 -0.017611 0.106081 -0.166 0.868 #> x2 -0.007479 0.100337 -0.075 0.941 #> #> Names of linear predictors: mean, loglink(sd) #> #> Log-likelihood: -140.6732 on 196 degrees of freedom #> #> Number of Fisher scoring iterations: 3 #> #> No Hauck-Donner effect found in any of the estimates models <- list(o0, o1, vm0, vm1, m0, m1) unit.aics <- aicCalcs(models)  Unit weighted model comparisons. In the case of unit weights, the dAIC calculations should agree with the TIC as a special case along with the standard AIC computation in glm models. pack.index <- table(unit.aics[["calculation"]]) unit.aics <- unit.aics[order(unit.aics[["calculation"]]), ] numeric.cols <- vapply(unit.aics, is.numeric, logical(1L)) unit.aics[numeric.cols] <- lapply(unit.aics[numeric.cols], round, ndp) unit.aics |> kable(row.names = FALSE, format = "pipe")  value calculation predictors delta df 287.35 glm (AIC) 1 NA 3.00 289.35 glm (AIC) 2 NA 4.00 287.21 proposal (dAIC) 1 0.98 2.93 289.26 proposal (dAIC) 2 0.99 3.96 290.09 survey (dAIC) 1 0.92 1.83 292.15 survey (dAIC) 2 0.95 2.86 99.32 survey v4.0 (dAIC) 1 0.86 0.86 101.32 survey v4.0 (dAIC) 2 0.93 1.87 287.03 vglm (TIC) 1 NA NA 289.07 vglm (TIC) 2 NA NA Scenario 2 Consider now when x2 is significant in explaining y. It should be included in the model and reflected in the AIC calculations. n <- 100L dat <- replicate(3, rnorm(n), simplify = FALSE) |> as.data.frame() |> setNames(c("y", "x1", "x2")) dat[["y"]] <- dat[["x2"]] + rnorm(n, sd = 0.5) unit.design <- svydesign(ids = ~1, weights = ~1, data = dat) m0 <- svyglm(y ~ x1, design = unit.design) summary(m0) #> #> Call: #> svyglm(formula = y ~ x1, design = unit.design) #> #> Survey design: #> svydesign(ids = ~1, weights = ~1, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) -0.04666 0.11948 -0.391 0.697 #> x1 0.06003 0.12155 0.494 0.623 #> #> (Dispersion parameter for gaussian family taken to be 1.401867) #> #> Number of Fisher Scoring iterations: 2 m12 <- svyglm(y ~ x1 + x2, design = unit.design) summary(m12) #> #> Call: #> svyglm(formula = y ~ x1 + x2, design = unit.design) #> #> Survey design: #> svydesign(ids = ~1, weights = ~1, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) -0.03446 0.05347 -0.644 0.521 #> x1 0.06865 0.05903 1.163 0.248 #> x2 1.02685 0.04582 22.413 <2e-16 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> (Dispersion parameter for gaussian family taken to be 0.2792745) #> #> Number of Fisher Scoring iterations: 2  glm versions o0 <- glm(y ~ x1, data = dat) summary(o0) #> #> Call: #> glm(formula = y ~ x1, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) -0.04666 0.11969 -0.390 0.697 #> x1 0.06003 0.12639 0.475 0.636 #> #> (Dispersion parameter for gaussian family taken to be 1.416172) #> #> Null deviance: 139.10 on 99 degrees of freedom #> Residual deviance: 138.78 on 98 degrees of freedom #> AIC: 322.56 #> #> Number of Fisher Scoring iterations: 2 o12 <- glm(y ~ x1 + x2, data = dat) summary(o12) #> #> Call: #> glm(formula = y ~ x1 + x2, data = dat) #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) -0.03446 0.05370 -0.642 0.523 #> x1 0.06865 0.05670 1.211 0.229 #> x2 1.02685 0.05200 19.746 <2e-16 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> (Dispersion parameter for gaussian family taken to be 0.2850327) #> #> Null deviance: 139.104 on 99 degrees of freedom #> Residual deviance: 27.648 on 97 degrees of freedom #> AIC: 163.23 #> #> Number of Fisher Scoring iterations: 2  VGAM versions in vglm for TIC vm0 <- vglm(y ~ x1, data = dat, family = uninormal) summary(vm0) #> #> Call: #> vglm(formula = y ~ x1, family = uninormal, data = dat) #> #> Coefficients: #> Estimate Std. Error z value Pr(>|z|) #> (Intercept):1 -0.04666 0.11849 -0.394 0.6937 #> (Intercept):2 0.16388 0.07071 2.318 0.0205 * #> x1 0.06003 0.12512 0.480 0.6314 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Names of linear predictors: mean, loglink(sd) #> #> Log-likelihood: -158.2816 on 197 degrees of freedom #> #> Number of Fisher scoring iterations: 3 #> #> No Hauck-Donner effect found in any of the estimates vm12 <- vglm(y ~ 1 + x1 + x2, data = dat, family = uninormal) summary(vm12) #> #> Call: #> vglm(formula = y ~ 1 + x1 + x2, family = uninormal, data = dat) #> #> Coefficients: #> Estimate Std. Error z value Pr(>|z|) #> (Intercept):1 -0.03446 0.05289 -0.652 0.515 #> (Intercept):2 -0.64281 0.07071 -9.091 <2e-16 *** #> x1 0.06865 0.05585 1.229 0.219 #> x2 1.02685 0.05122 20.049 <2e-16 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Names of linear predictors: mean, loglink(sd) #> #> Log-likelihood: -77.6133 on 196 degrees of freedom #> #> Number of Fisher scoring iterations: 7 #> #> No Hauck-Donner effect found in any of the estimates sig.models <- list(o0, o12, vm0, vm12, m0, m12) x2.aics <- aicCalcs(sig.models)  Below are the calculated criterions for this example and the model with x2 should be the clear superior model as it includes the most informative predictor against y. Again the proposed output agrees with the TIC and standard AIC calculation. The survey calculations vary wildly, the survey v4.0 has the correct conclusion but the values are not in agreement with TIC. Also the design effect adjustment in deltabar shows a very small deltabar for the 2 predictor model. On the other hand, survey >= 4.1-1 has the completely incorrect conclusion with an $$d$$AIC three times larger in the 2 predictor model and the design effect has blown up to 2 in deltabar. The proposal is stable with design effects near 1 for both models which is consistent with the earlier results. pack.index <- table(x2.aics[["calculation"]]) x2.aics <- x2.aics[order(x2.aics[["calculation"]]), ] numeric.cols <- vapply(x2.aics, is.numeric, logical(1L)) x2.aics[numeric.cols] <- lapply(x2.aics[numeric.cols], round, ndp) x2.aics |> kable(row.names = FALSE, format = "pipe")  value calculation predictors delta df 322.56 glm (AIC) 1 NA 3.00 163.23 glm (AIC) 2 NA 4.00 322.84 proposal (dAIC) 1 1.05 3.14 163.50 proposal (dAIC) 2 1.03 4.14 229.58 survey (dAIC) 1 0.76 1.52 569.02 survey (dAIC) 2 1.60 4.80 141.40 survey v4.0 (dAIC) 1 1.31 1.31 28.71 survey v4.0 (dAIC) 2 0.27 0.53 322.12 vglm (TIC) 1 NA NA 162.68 vglm (TIC) 2 NA NA Created on 2023-07-11 with reprex v2.0.2 • Interesting, thx for your effort. Wonder if the author will respond. Question: Tried both, original and you proposal; why are all AICs different from the AICs stored in the "svyglm" objects, i.e. fit$aic? Using 4.2-1. Commented Feb 3 at 11:58
• Author got in contact and implemented the proposal in survey 4.3. As for the fit$aic inside the svyglm models. I haven't checked to confirm but I believe survey weights the data and then calls glm under the hood. So it'll get the standard aic calculation from the glm code and it's parked in fit$aic
– JRW
Commented May 4 at 10:16
• Thanks for letting know, appreciated. Commented May 4 at 10:18