Dependent data can be modeled using covariance structures like compound symmetry, spherical, AR-1, and other. Using generalized least squares, inference can be made on the regression coefficients using model based standard errors. While the covariance structure accommodates some forms of heteroscedasticity, if we still have model misspecification a robust variance estimate should technically allow us to calculate correct 95% CIs for the trend. In R, the sandwich package implements this type of error calculation for generalized linear models, but not generalized least squares. How can this be done?
$\begingroup$
$\endgroup$
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
Generalized least squares jointly models fixed effects and a covariance structure in data to yield Gauss-Markov optimal estimators of the fixed effects. We can also refer to the fixed effects as the "trend", borrowing heavily from semiparametric inference and theory.
If the variance structure in a GLS model is misspecified, the estimate is an unbiased, but inefficient estimate of the trend, just like with OLS and independent data. The sandwich estimator, also called the heteroscedasticity-consistent estimator, provides asymptotically correct 95% confidence intervals even when the model is misspecified, so tests of trend are of the correct 0.05 level and achieve relatively good power.
The GLS model solves the estimating equation:
$$0 = \sum_{i=1}^n S_i (\hat{\beta}) = \mathbf{X}^T W^{-1} \left( Y - \mathbf{X}^T\hat{\beta} \right)$$
Where $W$ is the inverse covariance matrix estimated from a number of possible ways. As an interesting note, most GEE estimators use method-of-moments estimators for covariance parameters, but I believe the GLS which is an exact likelihood based procedure uses the EM algorithm. I do not know if this is a problem with my proposed solution here.
The robust sandwich variance estimator of $\hat{\beta}$ is given by:
$$\text{var}(\hat{\beta}) = \frac{1}{n}\mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-1}$$
Where $A = \dfrac{\partial S(\hat{\beta})}{\partial \hat{\beta}}$ and $B = S(\hat{\beta})S(\hat{\beta})^T/n$. Putting it all together, we should be able to construct a sandwich estimate from a GLS model by hand fairly easily.
Using the Ovary example from GLS, the model fm1 has a non-zero unweighted sum of residuals (sum(residuals(fm1)) = 1.65). Here is an example of checking the weighted score to verify that, indeed, it sums to zero thus solves the estimating equation
X <- model.matrix(fm1, data=Ovary)
Mares <- Ovary$Mare
W <- lapply(lapply(levels(Mares), getVarCov, obj=fm1), solve)
Xsplit <- lapply(split(as.data.frame(X), Mares), as.matrix)
Resids <- split(residuals(fm1), Mares)
Scores <- mapply(function(x, w, r) {
t(x)%*%w%*%r}, Xsplit, W, Resids)
rowSums(Scores)
gives
> rowSums(Scores)
[1] 1.318667e-13 -1.706968e-14 4.593981e-14
okay, this forms the basis of generating the model based variance-covariance matrix for the trend.
xwx <- Reduce(`+`, mapply(function(x,w) t(x)%*%w%*%x, Xsplit, W, SIMPLIFY=F))
solve(xwx)
vcov(fm1)
> solve(xwx)
(Intercept) sin(2 * pi * Time) cos(2 * pi * Time)
(Intercept) 4.417512e-01 -8.395204e-08 -1.363831e-01
sin(2 * pi * Time) -8.395204e-08 4.160866e-01 7.475767e-08
cos(2 * pi * Time) -1.363831e-01 7.475767e-08 4.865597e-01
> vcov(fm1)
(Intercept) sin(2 * pi * Time) cos(2 * pi * Time)
(Intercept) 4.417512e-01 -8.395204e-08 -1.363831e-01
sin(2 * pi * Time) -8.395204e-08 4.160866e-01 7.475767e-08
cos(2 * pi * Time) -1.363831e-01 7.475767e-08 4.865597e-01
Now the sandwich is obtained by using the residuals as plug-in estimators and obtaining the output:
Amat <- xwx
Bmat <- mapply(function(x, w, r) t(x)%*%w%*%diag(c(r)^2)%*%w%*%x, Xsplit, W, Resids, SIMPLIFY=F)
Bmat <- Reduce(`+`, Bmat)
SandCov <- solve(Amat) %*% Bmat %*% solve(Amat)
And these robust covariance estimates can be used in place of the model-based estimates for heteroscedasticity robust, consistent inference.