VGAMs for survey data: theory

The VGAM package’s vglm() function, like the survey package’s svymle() function, allows for maximum likelihood fitting where linear predictors are added to one or more parameters of a distribution — but vglm() is a lot faster and has many distributions already built in. This is how we make svyVGAM handle complex sampling.

I will write β for the regression parameters, θ for the base parameters of the response distribution, and η for the linear predictors. So, for example, in a classical linear model there would be two parameters θ: the mean (θ1) and variance (θ2). The mean would have a set of regression parameters and the variance would have a single parameter. Collectively, these would be β (maybe β11β1p and β21), and the two combinations that are plugged in as θ would be called η1 and η2. The big advantage of VGAM is that it does a lot of the work for the user: while the user can add new families, there are many pre-prepared ones, and there are built-in ways to constrain parameters to be equal or related in some other way.

To provide survey versions of vglm(), we need to (a) get design-consistent point estimates out of vglm(), and (b) construct design-based standard errors for the fit. The first is easy: vglm() accepts frequency weights, which are equivalent to sampling weights for point estimation with independent observations.

The second can be done in two ways: resampling (which is straightforward, if potentially slow), and linearisation. Linearisation requires computing the influence functions of the parameters $$h_i(\beta) = -\widehat{\cal I}^{-1}_w U_i(\beta)$$ where $\widehat{\cal I}_w$ is the weighted estimate of the population Fisher information, Ui = ∂βi(β) is the loglikelihood contribution of the ith observation, and wi is its weight. The influence functions have the property β̂ − β0 = ∑iwihi(β0) + op(∥β̂ − β0∥) so that the variance of β̂ is asymptotically the variance of the population total of the influence functions. The survey package provides a function svyrecvar() to estimate standard errors given the influence functions, or (a bit less efficiently) you can just call svytotal().

Resampling

A design object of class svyrep.design contains sets of replicate weights analogous to jackknife or bootstrap replicates. We need to call vglm() with each set of weights. It should be helpful to specify the full-sample estimates as starting values.

One complication is that sets of replicate weights will typically include some zeroes, which vglm() does not allow (eg, a bootstrap or jackknife resample will omit some observations). We set these to 10−9 times the maximum weight, which has the desired effect that the observations are present in the fit but with (effectively) zero weight.

Linearisation

The cov.unscaled slot of a summary.vglm object contains the inverse of the estimated population Fisher information, $\widehat{\cal I}^{-1}_w$.

The vglm object provides ηi(η) for the base parameters θ, and also the model matrices that specify βη, so we can construct Ui. We need to take care with the constraints, which can cause a coefficient β to appear in more than one linear predictor.

Suppose βx appears in both η1 and η2, with x values x1 and x2. The chain rule tells us βxi = ∂η1iβxη1 + ∂η2iβxη2 = (∂η1i)x1i + (∂η2i)x2i We might have x1 ≡ x2 ( = x), but that just means βxi = (∂η1i)xi + (∂η2i)xi

The constraint matrix C consists of 1s and 0s. If there are p parameters in M equations the matrix is M × p, with Cjk = 1 if parameter k is in linear predictor j. In the default, unconstrained setup, the constraint matrix consists of an M × M identity matrix for each parameter, pasted columnwise to give a M × pM matrix. In the proportional odds model, as another example, there are separate intercepts for each linear predictor but the other parameters all appear in every linear predictor. The first M × M block is the identity, and the rest of the matrix is a column of 1s for each predictor variable. Another way to say this is that Cjk = ∂(βkxk)ηj

So, if we want βi, the chain rule says

There is one further complication. The model.matrix method for vglm objects returns a model matrix of dimension Mn × p rather than n × p, so we need to sum over the rows for each observation, which can be identified from the row names, and then rescale. The right standardisation appears to come from defining $$\tilde C_{kj}=\frac{C_{kj}}{\sum_k C_{kj}}$$ and then $$\frac{\partial \ell_i}{\partial \beta_j}=\sum_k \frac{\partial \ell_i}{\partial \eta_k} \tilde C_{kj}x_{ikj}.$$