使用 R 中的 delta 方法估计线性回归的截距标准误差答案

【问题标题】：Estimation of the standard error of intercept of the linnear regression using the delta method in R使用 R 中的 delta 方法估计线性回归的截距标准误差
【发布时间】：2019-11-13 16:58:37
【问题描述】：

我尝试通过“car”包中的“deltaMethod”函数使用 delta 方法从线性回归中估计转换后的截距 (10^a) 的 SE。但是 deltaMethod 总是为截距的 SE 或其任何变换给出零！例如，对于未转换的截距：

> x <- 1:100; y <- rnorm(100, 4*x, 5)
> m1<-lm(y~x) 
> library(car)
> deltaMethod(m1, "Intercept")
          Estimate      SE   2.5 %  97.5 %
Intercept  -1.9561  0.0000 -1.9561 -1.9561

根据“deltaMethod”函数，截距是0.000，但应该是0.9034276：

> (summary(m1))$ coefficients [1,2]
[1] 0.9034276

“msm”包中的“deltamethod”函数给出截距 SE 的正确值：

> library(msm)
> deltamethod (~ x1, coef(m1), vcov(m1))
[1] 0.9034276

但是我在“car”包中的“deltaMethod”函数代码有什么问题？（出于某种原因，我更喜欢使用“car”包而不是“msm”）。谢谢！

【问题讨论】：

这对我有用（我从来没有得到 0）。重启 R 后你试过了吗？
只是为了让你知道我也遇到了同样的问题。我一直无法找到问题根源的答案。实际上，我最初对评估 1/(Intercept) 的标准误差很感兴趣，但是 deltaMethod 报告了零标准误差，然后我尝试获取 (Intercept) 的标准误差，也得到了零标准误差。因此，我使用 msm 包中的 deltamethod 以及解析解决方案获得了 1/(Intercept) 的标准误差，并且两者都同意并且当然不同于零。让我担心的是另一个 est

标签： r linear-regression intercept

【解决方案1】：

我似乎无法重现您的问题。我们可以测试这一点的一种方法是编写一个重复创建数据、拟合模型并返回 SE 的函数，就像您在上面所做的那样。我们可以将它传递给replicate 函数并多次重复该过程。这是使用您的代码进行调查的一种方法：

set.seed(42)

library(car)

delta_sim <- function(){
  x <- 1:100; y <- rnorm(100, 4*x, 5)
  m1<-lm(y~x) 
  a <- deltaMethod(m1, "Intercept")
  list(a)
}

out <- replicate(1000, delta_sim())

out <- do.call(rbind, out)

hist(out$SE)

这将产生以下结果：

可能值得重新启动您的 R 会话以确保您没有任何本地变量潜入您的函数。

【讨论】：

【解决方案2】：

看来，问题出在我的 R 或 PC 上，而不是“汽车”包中。是的，我重新启动了 R。我什至将 R 和我的所有软件包都更新到了最新版本。唉，错误仍然存在：“汽车”总是显示正确的截距（或其转换），但总是给它 SE=0。看起来我必须使用“msm”。总的来说，这不是那么重要，但我只是好奇这是什么错误。

【讨论】：

【解决方案3】：

我遇到了同样的问题。使用 MDEWITT 提供的代码在我的机器上产生了全零。我追踪到car::deltaMethod.default的以下部分

    para <- object         
    para.names <- names(para)
    para.names[1] <- gsub("\\(Intercept\\)", "Intercept", para.names[1])
    g. <- parse(text = g.)
    q <- length(para)

    envir <- new.env(parent=envir)
    for (i in 1:q) {
        assign(para.names[i], para[i], envir)
    }
    if(!missing(constants)){
     for (i in seq_along(constants)) assign(names(constants[i]), constants[[i]], envir)}
    est <- eval(g., envir)
    names(est) <- NULL
    gd <- rep(0, q)
    for (i in 1:q) {
        gd[i] <- eval(D(g., names(para)[i]), envir)
    }

注意使用names(para)[i] 代替para.names[i]。前者在截距上有括号，但后者将其清除，可能应该改用。我可以通过使用deltaMethod 的parameterNames 参数来解决这个问题，例如


    set.seed(42)
    
    library(car)
    
    delta_sim <- function(){
      x <- 1:100; y <- rnorm(100, 4*x, 5)
      m1<-lm(y~x) 
      a <- deltaMethod(m1, "b0", parameterNames=c("b0", "b1"))
      list(a)
    }
    
    out <- replicate(1000, delta_sim())
    
    out <- do.call(rbind, out)
    
    hist(out$SE)

产生了所需的输出。

【讨论】：