As stated in my question, I would like to access the Fisher weights used in PIRLS model fitting for GLMMs from my glmer() fit in the R package lme4. For such a simple task, I was surprised that I couldn't find any information in the documentation or on the internet at all. By looking at the structure of the glmer fit, I found two possible quantities that might correspond to what I want (but I have no way to know):
glmmfit@resp$sqrtWrkWt() and glmmfit@pp$Xwts. They seem to be the same thing (up to very small numerical error). Are these the Fisher weights, or the square root of them? Or is this something else entirely?
P.S.: Could someone also confirm that glmmfit@resp$wrkResp() gives the working responses z=G(y-mu)+eta (sometimes called pseudodata), where G is a matrix containing the derivaties of the link function? Unexpectedly, it turns out that when I do GLMM_model@resp$eta+GLMM_model@resp$wrkResids()-GLMM_model@resp$wrkResp(), having added an offset of 4 to the model, I get a vector of fours, not zeros as I would expect..
This is a harder question to answer than it should be, but let's try.
There is a draft paper describing the implementation of glmer (an unpublished sequel to the Bates et al. JSS paper available via vignette("lmer", package = "lme4") in this directory (PDF here), but — although it is useful as background reading — it doesn't make a direct connection with the code.
double glmResp::updateWts() {
d_sqrtrwt = (d_weights.array() / variance()).sqrt();
d_sqrtXwt = muEta() * d_sqrtrwt.array();
return updateWrss();
}
i.e., d_sqrtrtwt is the square root of the working weights [on the linear predictor or link scale] (to be honest I'm not sure what the r signifies); d_sqrtXwt is those weights transformed back on to the response/data scale (by multiplying by dmu/deta, the derivative of the inverse-link function).
From here, sqrtWrkWt is the same as the d_sqrtXwt value computed in updateWts.
Here we can see that the weights(., type = "working") returns object@pp$Xwts^2, and we can even see the comment that
the working weights available through pp$Xwts should be equivalent to: object@resp$weights*(object@resp$muEta()^2)/object@resp$variance()
However, the unit tests in tests/glmmWeights.R suggest that this equivalence is approximate. This may be fine, however, if the discrepancy is due to another instance of the general problem of reference class fields not being updated at the optimum, then this could cause real problems. see for example: https://github.com/lme4/lme4/issues/166
Here we see that wrkResp is defined as (d_eta - d_offset).array() + wrkResids(); and here that wrkResids() is (d_y - d_mu).array() / muEta();
Hopefully you should be able to access all the pieces you need without poking around in the guts this way ... e.g. weights(., "working") should give you the weights; family(.)$mu.eta should give you the derivative of the inverse-link function; residuals(., "working") should give you the working residuals.
The clue to why your "PS" is not working is that, as you can see from the code listed above, the $eta component of the @resp slot does not include the offset ... another reason it's best to try to work with accessor methods whenever possible instead of digging around ...