英文:
How to access the Fisher Weight matrix W from glmer() fit?
问题
根据我的问题,我想从R包lme4的glmer()拟合中访问用于GLMMs的PIRLS模型拟合中使用的Fisher权重。对于这样一个简单的任务,我很惊讶在文档或互联网上找不到任何信息。通过查看glmer拟合的结构,我找到了两个可能与我想要的内容对应的量(但我无法确定):
glmmfit@resp$sqrtWrkWt()
和 glmmfit@pp$Xwts
。它们似乎是相同的(除了非常小的数值误差)。这些是Fisher权重还是它们的平方根?还是完全不同的东西?
附注:有人可以确认一下 glmmfit@resp$wrkResp()
是否给出工作响应 z=G(y-mu)+eta(有时称为伪数据),其中 G 是包含链接函数导数的矩阵吗?出乎意料的是,当我执行 GLMM_model@resp$eta+GLMM_model@resp$wrkResids()-GLMM_model@resp$wrkResp()
时,在模型上添加了一个偏移量为4后,我得到的是一个四的向量,而不是我所期望的零。
英文:
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..
答案1
得分: 2
这是一个比应该更难回答的问题,但让我们尝试一下。
有一篇草稿论文描述了glmer
的实现(这是Bates等人的未发表续篇,可通过vignette("lmer", package = "lme4")
在此目录中获取(PDF在此),但尽管它作为背景阅读很有用,它并没有与代码直接联系起来。
- 权重在这里更新,使用以下代码:
double glmResp::updateWts() {
d_sqrtrwt = (d_weights.array() / variance()).sqrt();
d_sqrtXwt = muEta() * d_sqrtrwt.array();
return updateWrss();
}
即d_sqrtrtwt
是工作权重的平方根[在线性预测器或链接比例上](老实说,我不确定r
代表什么);d_sqrtXwt
是将这些权重转换回响应/数据比例的权重(通过乘以dmu/deta
,反链接函数的导数)。
从这里可以看到,sqrtWrkWt
与updateWts
中计算的d_sqrtXwt
值相同。
这里我们可以看到weights(., type = "working")
返回object@pp$Xwts^2
,甚至可以看到注释:
可通过pp$Xwts获得的工作权重应该等同于:object@resp$weights*(object@resp$muEta()^2)/object@resp$variance()
。
然而,tests/glmmWeights.R中的单元测试表明,这种等价性是近似的。如果差异是由于引用类字段在最优时未更新的另一个实例造成的,则可能会导致实际问题。请参阅例如:https://github.com/lme4/lme4/issues/166
这里我们看到wrkResp
被定义为(d_eta - d_offset).array() + wrkResids();
,这里则显示wrkResids()
为(d_y - d_mu).array() / muEta();
希望您能够访问所有需要的部分,而无需在内部细节中查找……例如,weights(., "working")
应该会给您权重;family(.)$mu.eta
应该会给您反链接函数的导数;residuals(., "working")
应该会给您工作残差。
为什么您的“PS”不起作用的线索在于,正如您从上面列出的代码中所看到的,@resp
槽的$eta
部分不包括偏移...这也是尽量在可能的情况下使用访问器方法而不是深入挖掘的原因之一...
英文:
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.
- The weights are updated here via
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() <br>
> 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 ...
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论