1

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

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.rqwhich 对所有模型运​​行 Wald 测试，并为每个协变量输出一个 pvalue，告诉我每个协变量的影响是否在我的结果分布中显着变化。

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无法对自举估计执行 anova (Wald) 检验。软件包上的信息文件quantreg明确指出，“应该对要在 anova.rq 中使用的方法进行扩展”关于 boostrapping 方法。

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

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

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

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

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

4

2 回答 2

0

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

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

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

请注意，如果您在每个帖子中仅包含 1 个问题，您可能会得到更好的问题答复。

于 2015-08-17T17:08:44.983 回答
0

咨询了一位同事，他确认这不太可能，Hinv并且J可能是从自举估计中“逆向”计算的。然而，我们认为可以使用 Wald 检验比较来自不同 taus 的估计值，如下所示。

从rqs产生的对象

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

在这种情况下，您提取感兴趣变量的引导 Beta 值，每个 tauVAR的第一个协变量FML

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)

您还想验证引导的 Beta 是否是正态分布的，如果您正在运行许多 taus，检查所有这些 QQ 图可能很麻烦，因此只需逐行求和

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

这似乎有效，如果有人能想到我的解决方案有什么问题，请分享

于 2015-08-17T23:29:23.343 回答