如何从glmer()拟合中访问Fisher权重矩阵W?

huangapple go评论62阅读模式
英文:

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,反链接函数的导数)。

这里可以看到,sqrtWrkWtupdateWts中计算的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(., &quot;working&quot;) should give you the weights; family(.)$mu.eta should give you the derivative of the inverse-link function; residuals(., &quot;working&quot;) 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 ...

huangapple
  • 本文由 发表于 2023年6月1日 22:39:42
  • 转载请务必保留本文链接:https://go.coder-hub.com/76383075.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定