covariance matrix of residuals from a fitted model to decorrelate residuals I fit a geeglm model with clustered data and now I would like to decorrelate the residuals of the model in order to run model diagnostics.
I read that if I can obtain the covariance matrix of the residuals I can use the Cholesky decomposition to transform the residuals and thus decorrelate them.
How can I obtain this covariance matrix of residuals?
Furthermore, I understand that doing this involves a decent understanding of linear algebra in R. Can someone point me in the right direction on this topic?
 A: The standard residual diagnostics of a linear model do not apply to GEE for several reasons. In this case, a strong rationale is needed to know what you actually need to do in this case.
With a standard linear model, you might plot residuals against fitted on the chance of catching either a) a mean model misspecification (real trend different from model) or b) heteroscedasticity (non-constant mean-variance relation). "Funnel shaped residuals" would indicate B and trending residuals indicate A. The point of a GEE is to use the linear model as a rule of thumb, so for instance if the trend is in fact curvilinear, the linear parameter estimates provide the "expected slope" of the response over a range of the predictor.
The GEE estimates the model parameters using a generalized least squares approach in addition to a sandwich variance estimator. More precisely, the estimating equation is given by:
$$ U(\beta) = D^T V^{-1} (Y - g^{-1}(X^T \beta)$$
Where $V$ is the working covariance matrix, and $D = \frac{\partial}{\partial \beta} g^{-1}(X^T \beta)$. The value of $\beta$ that provides a root to the above equation is the GEE estimate. The $A$ (bread) matrix is given by $A = \frac{\partial^2}{\partial^2 \beta} g^{-1}(X^T \beta)$ and the $B$ (meat) matrix $B = U(\beta)U(\beta)^T$. The sandwich variance estimate is $A^{-1}BA^{-1}$. This is also called the heteroscedasticity consistent variance estimate or "HC". Basically point B from the residual diagnostic exercise is completely deprecated with a robust enough sample (small sample performance of the sandwich is quite good, a number of 40 is often quoted). It turns out that model misspecification can conceptually be considered a form of heteroscedasticity. In each case, a misstep in the GEE specification should lead the standard errors to being exploded to accommodate the problem. In that regard, getting the model "right" only matters insofar as you can get the "optimal GEE" in which case $A=B$ and the covariance estimate is $A^{-1}$.
To inspect these issues, I suggest you fit the generalized least squares model, and obtain the Pearsonized residuals. With those you can perform a standard inspection of funnel and trend in the plot of Pearson residuals versus mean. But "optimality" is rarely a concern in its own right.
The advantage of a GEE is that you do not need assumptions related to residual diagnostics to apply to claim that the model is interpretable or valid.
