最小二乘回归系数非线性函数的标准误差和置信区间答案

【问题标题】：Standard error and confidence interval for nonlinear function of least squares regression coefficients最小二乘回归系数非线性函数的标准误差和置信区间
【发布时间】：2017-05-08 08:07:01
【问题描述】：

我在 R 中运行 OLS 回归，从中得到几个系数。以下是部分代码：

Attacks <- Treat.Terr.Dataset$Attacks[2:30]
Attackslag <- Treat.Terr.Dataset$Attacks[1:29]
TreatmentEffect <- Treat.Terr.Dataset$TreatmentEffect[2:30]
TreatmentEffectlag <- Treat.Terr.Dataset$TreatmentEffect[1:29]

olsreg <- lm(TreatmentEffect ~ TreatmentEffectlag + Attacks + Attackslag)
coeffs<-olsreg$coefficients

然后我需要计算：(Attacks + Attackslag) / (1 - TreatmentEffectlag)。问题是我可以使用(coeffs[3] + coeffs[4]) / (1 - coeffs[2]) 在 R 上执行此操作，但结果是一个没有任何 p 值或置信区间的固定数字，就像计算器会显示给我一样。

有谁知道我是否可以使用任何函数来计算置信区间？

编者注

如果目标量是回归系数的线性函数，则问题会简化为一般线性假设检验，其中可以进行精确推断。

【问题讨论】：

欢迎来到 StackOverflow！请阅读有关how to ask a good question 的信息以及如何提供reproducible example。这将使其他人更容易帮助您。
你打算用这些做什么？这将确定最佳响应。

标签： r regression linear-regression lm coefficients

【解决方案1】：

## variance-covariance of relevant coefficients
V <- vcov(olsreg)[2:4, 2:4]
## point estimate (mean) of relevant coefficients
mu <- coef(olsreg)[2:4]

## From theory of OLS, coefficients are normally distributed: `N(mu, V)`
## We now draw 2000 samples from this multivariate distribution
beta <- MASS::mvrnorm(n = 2000, mu, V)

## With those 2000 samples, you can get 2000 samples for your target quantity
z <- (beta[, 2] + beta[, 3]) / (1 - beta[, 1])

## You can get Monte Carlo standard error, and Monte Carlo Confidence Interval
mean(z)
sd(z)
quantile(z, prob = c(0.025, 0.975))

## You can of course increase sample size from 2000 to 5000

【讨论】：

不错。我打算通过 delta 方法来做到这一点。有关此处发生的情况的说明，您可以将 OP 参考ms.mcmaster.ca/~bolker/emdbook/chap7A.pdf 的第 5 节（特别是 5.3）

【解决方案2】：

这是一个使用 'car' 包中的 delta 方法的独立示例：

# Simulate data
dat <- data.frame(Attacks = rnorm(30), Trt=rnorm(30))
dat <- transform(dat, AttacksLag = lag(Attacks), TrtLag = lag(Trt))
dat <- dat[2:30,]

# Fit linear model
m1 <- lm(Trt ~  TrtLag + Attacks + AttacksLag, data=dat)

# Use delta method
require("car")
del1 <- deltaMethod(m1, "(Attacks + AttacksLag) / (1 - TrtLag)")

# Simple Wald-type conf int
del1$Est +  c(-1,1) * del1$SE * qt(1-.1/2, nrow(dat)-length(coef(m1)))
# [1] -0.2921529  0.6723991

【讨论】：