如何从MCMCglmm的后验密度模拟感兴趣的量? [英] How can one simulate quantities of interest from the posterior density in MCMCglmm?

问题描述

我想用Zelig软件包的方式或多或少地模拟使用MCMCglmm估计的模型中的感兴趣量.在Zelig中，可以为独立值设置所需的值，然后软件将计算结果变量的结果(预期值，概率等).一个例子:

I would like to simulate quantities of interest from a model estimated with MCMCglmm more or less the way Zelig package does. In Zelig you can set the values you want for the independent values and software calculates the result for the outcome variable (expected value, probability, etc). An example:

# Creating a dataset:
set.seed(666)
df <- data.frame(y=rnorm(100,20,20),z=rnorm(100,50,70))

# Loading Zelig
library(Zelig)

# Model
m1.zelig <- zelig(y~z, data=df, model="ls")
summary(m1.zelig)

# Simulating z = 10
s1 <- setx(m1.zelig, z = 10)
simulation <- sim(m1.zelig, x = s1)
summary(simulation)

我们可以看到，如果z = 10 y大约为17.

As we can see, if z = 10 y is approximately 17.

# Same model with MCMCglmm
library(MCMCglmm)
m1.mcmc <- MCMCglmm(y~z, data=df, family = "gaussian", verbose = FALSE)
summary(m1.mcmc)

有什么方法可以用MCMCglmm估计的后验分布模拟z = 10并获得y的期望值?非常感谢你！

Is there any way to simulate z = 10 with the posterior distribution estimated by MCMCglmm and get the expected value of y? Thank you very much!

推荐答案

您可以模拟，但不像Zelig那样容易.您需要进一步了解要拟合的模型的结构以及参数在MCMCglmm对象中的存储方式.

You can simulate, but not as easily as in Zelig. You have to know a little more about the structure of the model you're fitting and the way the parameters are stored in the MCMCglmm object.

设置数据并拟合:

set.seed(666)
df <- data.frame(y=rnorm(100,20,20),z=rnorm(100,50,70))
library(MCMCglmm)
m1.mcmc <- MCMCglmm(y~z, data=df, family = "gaussian", verbose=FALSE)

R中用于预测和仿真的最常见协议是建立一个具有与原始数据相同结构的新数据帧:

The most common protocol in R for prediction and simulation is to set up a new data frame with the same structure as the original data:

predframe <- data.frame(z=10)

构造线性模型的模型矩阵:

Construct the model matrix for the linear model:

X <- model.matrix(~z,data=predframe)

现在使用存储在MCMCglmm对象的Sol(解决方案")组件中的系数链；为了方便起见，我将其设置为矩阵计算.

Now use the chains of coefficients stored in the Sol ("solution") component of the MCMCglmm object; for convenience, I've set this up as a matrix calculation.

predframe$y <- X %*% t(m1.mcmc$Sol)

如果要模拟更复杂的模型(线性或广义线性混合模型)，则需要加倍努力(适当处理随机效应和反向链接函数)...

If you want to simulate for more complicated models (linear or generalized linear mixed models) then you'll need to work a little harder (handle random effects and inverse-link functions appropriately) ...

这篇关于如何从MCMCglmm的后验密度模拟感兴趣的量?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！