【问题标题】:MCMCglmm: Posterior mean of modeled response variableMCMCglmm:建模响应变量的后验平均值
【发布时间】: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


【解决方案1】:

除了saveX = TRUEsaveZ = TRUE 选项之外,您是否添加了pr = TRUE?如果没有,gen_lin_mix_mod$Sol 只包含固定效果,不包含随机效果。这就是为什么尺寸不对。

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2021-11-10
    • 2021-10-08
    相关资源
    最近更新 更多