【发布时间】:2017-04-13 18:33:13
【问题描述】:
我想将使用MCMCglmm 估计的响应变量的建模(“拟合”)值与零膨胀泊松分布与数据中的观察值进行比较。任何人都可以就从这样的模型中提取响应变量的条件平均值的机制提出建议吗?本质上,我只是想使用观察到的数据和模型中响应变量的预测值来计算残差,但 MCMCglmm 中的 zipoisson 系列模型的数据规模尚未实现 predict。
require(MCMCglmm)
example = data.frame(response=rbinom(10000,1,0.05), predictorA=rnorm(10000,100,10),
predictorB=rnorm(10000,50,5), predictorC=rnorm(10000,1000,100),
predictorD=rnorm(10000,10,1), predictorE=rnorm(10000,10000,1),
randomA=runif(10000,1,10), randomB=runif(10000,75,90),
randomC=runif(10000,800,10000))
pois_example = round(example,0)
gen_lin_mix_mod = MCMCglmm(fixed = response ~ predictorA + predictorB + predictorC +
predictorD + predictorE,
random = ~ randomA + randomB + randomC,
family = "zipoisson", data = pois_example, nitt = 10000,
burnin = 1000, rcov = ~ idh(trait):units, DIC = TRUE)
更新
在 Jarrod Hadfield 的课程笔记文件中,似乎有一种用于计算预测值的“手动”方法:
在对 MCMCglmm 的调用中,我们指定了 saveX=TRUE 和 saveZ=TRUE 表明我们想要保存设计矩阵。我们可以结合 将这些矩阵放入设计矩阵 W 并乘以参数 向量 θ 得到预测(见公式 2.9):
W.1
prediction.1
xyplot(weight+prediction.1@x~Time|Chick, data=ChickWeight)
但是,当我尝试实现此功能时,出现以下错误:
design_mat = cBind(gen_lin_mix_mod$X, gen_lin_mix_mod$Z)
fitted = design_mat %*% posterior.mode(gen_lin_mix_mod$Sol)
observed = example
> design_mat = cBind(gen_lin_mix_mod$X, gen_lin_mix_mod$Z)
> fitted = design_mat %*% posterior.mode(gen_lin_mix_mod$Sol)
Error in design_mat %*% posterior.mode(gen_lin_mix_mod$Sol) :
Cholmod error 'X and/or Y have wrong dimensions' at file ../MatrixOps/cholmod_sdmult.c,
line 90
【问题讨论】:
-
你能给我们一个可重现的例子吗? (您确定要比较均值而不是条件均值吗?)
-
添加了示例。而且,是的,好点;条件手段确实是我想要的。谢谢你帮我澄清。
标签: r mixed-models