LM的覆盖概率计算 这是 10 万次重复后的样子

编程问答 2022-06-01

问题描述

我正在尝试计算我在回归的截距和斜率上生成的一组剩余引导复制的覆盖概率。谁能告诉我如何计算置信区间的覆盖概率？非常感谢。

请注意，我使用 Qr 分解手动运行回归，但如果更容易，您可以使用 lm()。我只是认为手动执行会更快。

set.seed(42)  ## for sake of reproducibility
n <- 100
x <- rnorm(n)
e <- rnorm(n)
y <- as.numeric(50 + 25*x + e)
dd <- data.frame(id=1:n,x=x,y=y)

mo <- lm(y ~ x,data=dd)

# Manual Residual Bootstrap
resi <- residuals(mo)
fit <- fitted(mo)
ressampy <- function() fit + sample(resi,length(resi),replace=TRUE)
# Sample y values:
head(ressampy())
# Qr decomposition of X values
qrX <- qr(cbind(Intercept=1,dd[,"x",drop=FALSE]),LAPACK=TRUE)
# faster than LM
qr.coef(qrX,"y"])
# One Bootstrap replication
boot1 <- qr.coef(qrX,ressampy())
# 1000 bootstrap replications
boot <- t(replicate(1000,qr.coef(qrX,ressampy())))

编辑 结合jay.sf的回答，我重写了运行lm()方法的代码，比较了jay.sf分享的链接中计算覆盖概率的第一种方法和第二种方法：

library(lmtest);library(sandwich)
ci <- confint(coeftest(mo,vcov.=vcovHC(mo,type="HC3")))
ci

FUNInter <- function() {
  X <- model.matrix(mo)
  ressampy.2 <- fit + sample(resi,replace = TRUE)
  bootmod <- lm(ressampy.2 ~ X-1)
  confint(bootmod,"X(Intercept)",level = 0.95)
}

FUNBeta <- function() {
  X <- model.matrix(mo)
  ressampy.2 <- fit + sample(resi,"Xx",level = 0.95)
}

set.seed(42)
R <- 1000
Interres <- replicate(R,FUNInter(),simplify=FALSE)
Betares <- replicate(R,FUNBeta(),simplify=FALSE)
ciinter <- t(sapply(Interres,function(x,y) x[grep(y,rownames(x)),],"X\\(Intercept\\)"))
cibeta <- t(sapply(Betares,"Xx"))

#second approach of calculating CP
sum(ciinter[,"2.5 %"] <=50 & 50 <= ciinter[,"97.5 %"])/R
[1] 0.842
sum(cibeta[,"2.5 %"] <=25 & 25 <= cibeta[,"97.5 %"])/R
[1] 0.945
#first approach of calculating CP
sum(apply(ciinter,1,function(x) {
  all(data.table::between(x,ci[1,1],2]))
}))/R
[1] 0.076
sum(apply(cibeta,ci[2,2]))
}))/R
[1] 0.405

解决方法

根据 Morris et. al 2019,Table 6，覆盖概率被定义为真实 theta 在自举置信区间 (CI) 内的概率（即该模型基于实际数据应用于许多样本，或者——换句话说——新的实验）：

enter image description here

因此，我们希望根据 OP 提出的 i.i.d 计算 CI。 bootstrap declare module '.*mp3'; 次并计算 theta 在这些 CI 中出现或不出现的频率的比率。

首先，我们使用实际数据估计我们的模型 R。

mo

为了避免重复中不必要的解包拟合值mo <- lm(y ~ x) 、残差yhat、模型矩阵u和系数X，我们预先提取它们。

coef0

在引导函数 yhat <- mo$fitted.values u <- as.matrix(mo$residuals) X <- model.matrix(mo) theta <- c(50,25) ## known from data generating process of simulation 中，我们将想要执行的所有步骤包装在一个复制中。为了应用非常快的 FUN，我们必须手动计算白色标准误差（与 .lm.fit 相同）。

lmtest::coeftest(fit,vcov.=sandwich::vcovHC(fit,type="HC1"))

现在我们使用 FUN <- function() { ## resampling residuals y.star <- yhat + sample(u,length(u),replace=TRUE) ## refit model fit <- .lm.fit(X,y.star) coef <- fit$coefficients[sort.list(fit$pivot)] ## alternatively using QR,but `.lm.fit` is slightly faster # qrX <- qr(X,LAPACK=TRUE) # coef <- qr.coef(qrX,y.star) ## white standard errors v.cov <- chol2inv(chol(t(X) %*% X)) meat <- t(X) %*% diag(diag(u %*% t(u))) %*% X ## degrees of freedom adjust (HC1) d <- dim(X) dfa <- d[1] / (d[1] - d[2]) white.se <- sqrt(diag(v.cov %*% meat %*% v.cov)*dfa) ## 95% CIs ci <- coef + qt(1 - .025,d[1] - d[2])*white.se %*% t(c(-1,1)) ## coverage c(intercept=theta[1] >= ci[1,1] & theta[1] <= ci[1,2],x=theta[2] >= ci[2,1] & theta[2] <= ci[2,2]) } 执行引导程序。

replicate

两列中的 R <- 5e3 set.seed(42) system.time(res <- t(replicate(R,FUN()))) # user system elapsed # 71.19 28.25 100.28 head(res,3) # intercept x # [1,] TRUE TRUE # [2,] FALSE TRUE # [3,] TRUE TRUE 的 mean 同时跨行，或分别在每列中，给出了我们正在寻找的覆盖概率。

TRUE

这是 10 万次重复后的样子

enter image description here

数据：

(cp.t <- mean(apply(v <- res,1,all)))  ## coverage probability total  
(cp.i <- colMeans(res))  ## coverage probability individual coefs
(cp <- c(total=cp.t,cp.i))
#  total intercept         x
# 0.8954    0.9478    0.9444

## values with other R:
#   total intercept         x
# 0.90700   0.95200   0.95200  ## R ==   1k
# 0.89950   0.95000   0.94700  ## R ==   2k
# 0.89540   0.94780   0.94440  ## R ==   5k
# 0.89530   0.94570   0.94680  ## R ==  10k
# 0.89722   0.94694   0.94777  ## R == 100k
confidence-intervalrrregressionresampling