2017-06-14 77 views
2

我要适合在这excellent example所提供的数据模型如何计算该响应的95%置信区间,进行逻辑回归后:predic.glm返回什么标准错误(...,type =“response”,se.fit = TRUE)?

foo <- mtcars[,c("mpg","vs")]; names(foo) <- c("x","y") 
mod <- glm(y ~ x, data = foo, family = binomial) 
preddata <- with(foo, data.frame(x = seq(min(x), max(x), length = 100))) 
preds <- predict(mod, newdata = preddata, type = "link", se.fit = TRUE) 

critval <- 1.96 ## approx 95% CI 
upr <- preds$fit + (critval * preds$se.fit) 
lwr <- preds$fit - (critval * preds$se.fit) 
fit <- preds$fit 

fit2 <- mod$family$linkinv(fit) 
upr2 <- mod$family$linkinv(upr) 
lwr2 <- mod$family$linkinv(lwr) 

现在,我的问题来自于一个事实你可以直接得到的回应,只是通过询问

predict(..., type = 'response', se.fit = TRUE) 

没有转化

predict(..., type = 'link', se.fit = TRUE) 

但是,这会产生不同的标准错误。 这些错误是什么?它们可以直接用于计算拟合响应值的置信区间吗?

predict.glm的代码,相关的是这样的:

switch(type, response = { 
    se.fit <- se.fit * abs(family(object)$mu.eta(fit)) 
    fit <- family(object)$linkinv(fit) 
}, link = , terms =) 

和比较输出:

preds_2 <- predict(mod, newdata = preddata, type = "response", se.fit = TRUE) 


> head(preds_2$fit) 
     1   2   3   4   5   6 
0.01265744 0.01399994 0.01548261 0.01711957 0.01892627 0.02091959 

> head(preds_2$se.fit) 
     1   2   3   4   5   6 
0.01944681 0.02098841 0.02263473 0.02439022 0.02625902 0.02824491 

它似乎并不十分明显如何从上面去:

> head(fit2) 
     1   2   3   4   5   6 
0.01265744 0.01399994 0.01548261 0.01711957 0.01892627 0.02091959 
> head(upr2) 
     1   2   3   4   5   6 
0.2130169 0.2184891 0.2240952 0.2298393 0.2357256 0.2417589 
> head(lwr2) 
      1   2   3   4   5   6 
0.0006067975 0.0007205942 0.0008555502 0.0010155472 0.0012051633 0.0014297930 

回答

3

如果你看看?family,你会看到$mu.etamu的衍生物,而不是eta(即反向链接函数的导数)。因此,se.fit = TRUE对于type = "response"给出了真实标准误差的一阶近似值。这被称为“delta方法”。

当然,由于反向链接函数通常是非线性的,所以置信区间关于拟合值不是对称的。因此,在链接规模上获得置信区间,然后将其映射回响应规模即可。

+1

谢谢,我希望他们会从你的第一段开始写一行,并放在'?predict.glm'文件中。 – Alex

相关问题