r - 如何从 MCMCregress 计算模型残差

Question

我正在使用贝叶斯推理做类作业。为此，我使用来自 MCMCpack 的 MCMCregress 函数。当我想要获得残差时，问题就出现了，因为该函数不提供残差，所以我必须“手动”计算它们(在 R 中)。

我的模型是:

Model <- MCMCregress(Y~X1+X2+X3+X4+X5, data=DATA)

其中 X1 和 X5 是连续的，而 X2、X3 和 X4 > 是二分法。模型输出为我提供了每个变量的估计:

(Intercept) = 1.90, 
X1 = -0.02, 
X2 = -0.05, 
X3 = 0.32,
X4 = 0.61,
X5 = -0.003,
sigma2 = 0.54

我想我必须这样做:

1.90 - 0.02*X1 - 0.05*X2 + 0.32*X3 ...

但我知道我在 R 代码中遗漏了一些重要的东西，所以我想知道哪个是正确的 R 代码以获得残差。

<小时/>

这是一个可重现的示例(尽管它与原始数据不对应):

Y  <- c(0.2,0.8,1,4.3,5,3.5,3.2,1,3.3,1,2,4,3.6,3,5,4.3,3.2,4,1,2)
X1 <- c(17,13,15,NA,12,24,15,NA,12,22,14,12,18,NA,10,13,12,11,26,10)
X2 <- c(0,0,0,1,0,0,1,1,1,1,0,1,0,0,0,0,0,1,NA,NA)
X3 <- c(0,0,1,1,1,0,1,0,0,0,0,1,0,NA,0,NA,NA,0,1,0)
X4 <- c(1,0,1,0,0,1,0,0,NA,0,NA,NA,0,0,1,0,1,1,0,1)
X5 <- c(2.46,4.56,32.1,NA,NA,NA,NA,NA,NA,3.76,5.67,4.56,NA,
        17.32,12.2,4.56,7.2,1.2,NA,9.2)

X2 <- ifelse(X2 > 0, c("Yes"), c("No")) 
X3 <- ifelse(X3 > 0, c("Yes"), c("No")) 
X4 <- ifelse(X4 > 0, c("Yes"), c("No")) 

DATA <- data.frame(Y,X1,X2,X3,X4,X5)


library(MCMCpack)
MCMC <- MCMCregress(Y~X1+X2+X3+X4+X5, data=DATA)
summary(MCMC)

如何获得残差？

Answer 1

我发现 MCMC 链在接近尾声时发生了一些可怕的事情运行 - 查看 xyplot(MCMC,layout=c(2,4))。我不知道什么问题是，你必须解决这个问题，但与此同时我将使用中位数而不是平均值作为系数估计。

coefs <- apply(MCMC,2,median)[1:6]
X <- model.matrix(~X1+X2+X3+X4+X5, data=DATA)
pred0 <- X %*% coefs

最复杂的部分是处理 NA 值:这是内置的大多数现有的residuals()方法...

pred <- napredict(attr(na.exclude(DATA),"na.action"),pred0)
resids <- pred-DATA$Y

请注意，在这种情况下，一旦删除 NA 值，您就会将 5 参数回归模型拟合到 6 个完整的观测值。您可能需要对插补进行一些合理的处理才能获得任何合理的答案......