2

我想从一个MCMCglmm或多或少地以Zelig包的方式估计的模型中模拟感兴趣的数量。Zelig您可以为独立值设置所需的值，软件计算结果变量的结果（预期值、概率等）。一个例子：

# 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。

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

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

4

1 回答 1

3

您可以模拟，但不像 Zelig 那样容易。您必须更多地了解您正在拟合的模型的结构以及参数在MCMCglmm对象中的存储方式。

设置数据并拟合：

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 中用于预测和模拟的最常见协议是建立一个与原始数据具有相同结构的新数据帧：

predframe <- data.frame(z=10)

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

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

Sol现在使用存储在对象（“解决方案”）组件中的系数链MCMCglmm；为方便起见，我将其设置为矩阵计算。

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

如果您想模拟更复杂的模型（线性或广义线性混合模型），那么您需要更加努力（适当地处理随机效应和反向链接函数）......

于 2014-01-11T02:14:20.353 回答