对分位数回归的自举估计执行方差分析

Question

所以我使用 R 中的 quantreg 包进行分位数回归分析，以测试我的预测变量的影响如何随结果分布而变化。

FML <- as.formula(outcome ~ VAR + c1 + c2 + c3)
quantiles <- c(0.25, 0.5, 0.75)

q.Result <- list()
for (i in quantiles){
    i.no <- which(quantiles==i)
    q.Result[[i.no]] <- rq(FML, tau=i, data, method="fn", na.action=na.omit)
}

然后我调用 anova.rq 对所有模型运行 Wald 测试并输出每个协变量的 p 值告诉我每个协变量的影响是否在我的结果分布中有显着差异。

anova.Result <- anova(q.Result[[1]], q.Result[[2]], q.Result[[3]], joint=FALSE)

这很好用。然而，对于我的特定数据（以及一般情况下？），引导我的估计和他们的错误是可取的。我对上面的代码稍作修改。

q.Result <- rqs(FML, tau=quantiles, data, method="fn", na.action=na.omit)
q.Summary <- summary(Q.mod, se="boot", R=10000, bsmethod="mcmb",
                            covariance=TRUE)

这是我卡住的地方。 quantreg 目前无法对 boostrapped 估计执行方差分析 (Wald) 检验。 quantreg 包上的信息文件明确指出 "extensions of the methods to be used in anova.rq should be made" 关于 boostrapping 方法。

正在查看 anova.rq 方法的详细信息。我可以看到它需要 2 个在引导时分位数模型中不存在的组件。

1) Hinv（逆海森矩阵）。包信息文件明确指出“请注意，对于 se = "boot"，无法将估计的协方差矩阵拆分为其三明治组成部分。”

2) J 根据信息文件，它是“如果 cov=TRUE 和 se != "iid" 则返回梯度矩阵的未缩放外积。Huber 三明治是cov = tau(1-tau) Hinv %*% J %*% Hinv。至于Hinv分量，se == "boot"时没有J分量。（注意做Huber三明治需要加tau (1-tau) 自己做蛋黄酱。）"

我可以根据自举估计值计算或估计 Hinv 和 J 吗？如果不是，最好的方法是什么？

非常感谢对此的任何帮助。这是我第一次在这里发布问题，尽管我在过去对其他人的问题的回答中受益匪浅。

Answer 1

对于问题2：您可以使用R =进行重采样。例如：

anova(object, ..., test = "Wald", joint = TRUE, score =
"tau", se = "nid", R = 10000, trim = NULL)

其中 R 是检验的 anowar 形式的重采样重复次数，用于估计检验统计量的参考分布。

请注意，如果每个 post 只包含 1 个问题，您的问题可能会得到更好的答复。

Answer 2

咨询了一位同事，他确认 Hinv 和 J 不太可能 'reverse' 根据自举估计计算得出。然而，我们解决了可以使用 Wald 测试比较不同 taus 的估计值，如下所示。

来自

生成的对象 rqs

q.Summary <- summary(Q.mod, se="boot", R=10000, bsmethod="mcmb", covariance=TRUE)

在这种情况下，您提取感兴趣变量的自举 Beta 值 VAR，FML 中每个 tau

的第一个协变量

boot.Bs <- sapply(q.Summary, function (x) x[["B"]][,2])
B0 <- coef(summary(lm(FML, data)))[2,1] # Extract liner estimate data linear estimate

然后计算 wald 统计量并获得具有自由度分位数的 pvalue

Wald <- sum(apply(boot.Bs, 2, function (x) ((mean(x)-B0)^2)/var(x)))
Pvalue <- pchisq(Wald, ncol(boot.Bs), lower=FALSE)

您还想验证自举 Betas 是否呈正态分布，如果您运行许多 taus，检查所有这些 QQ 图可能很麻烦，因此只需按行对它们求和

qqnorm(apply(boot.Bs, 1, sum))
qqline(apply(boot.Bs, 1, sum), col = 2)

这似乎有效，如果有人认为我的解决方案有任何问题，请分享

对分位数回归的自举估计执行方差分析

Performing Anova on Bootstrapped Estimates from Quantile Regression

r

regression

anova

statistics-bootstrap

quantreg