在R语言中，MCMCglmm多项式模型

Question

我将尝试在R中使用MCMCglmm软件包创建一个模型。

数据结构如下，其中dyad、focal、other都是随机效应，predict1-2是预测变量，response 1-5是捕捉不同子类型观察行为数量的结果变量：

 dyad focal other r    present  village  resp1 resp2 resp3 resp4 resp5 
 1    10101 14302 0.5  3        1        0     0     4     0     5
 2    10405 11301 0.0  5        0        0     0     1     0     1
 …

因此，只有一个结果（教学）的模型如下：

 prior_overdisp_i <- list(R=list(V=diag(2),nu=0.08,fix=2), 
 G=list(G1=list(V=1,nu=0.08), G2=list(V=1,nu=0.08), G3=list(V=1,nu=0.08), G4=list(V=1,nu=0.08)))

 m1 <- MCMCglmm(teaching ~ trait-1 + at.level(trait,1):r + at.level(trait,1):present, 
 random= ~idh(at.level(trait,1)):focal + idh(at.level(trait,1)):other + 
 idh(at.level(trait,1)):X + idh(at.level(trait,1)):village, 
 rcov=~idh(trait):units, family = "zipoisson", prior=prior_overdisp_i, 
 data = data, nitt = nitt.1, thin = 50, burnin = 15000, pr = TRUE, pl = TRUE, verbose = TRUE, DIC  = TRUE)

Hadfield的课堂笔记（第五章）提供了一个多项式模型的示例，该模型仅使用一个具有3个水平的结果变量（3种羊角类型）。类似的处理可以在这里找到：http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/。这不完全适用于我正在做的事情，但包含了有用的背景信息。

另一个参考资料(Hadfield 2010)给出了一个多响应MCMCglmm的示例，遵循相同的格式，但使用cbind()来预测一个响应向量，而不是一个单一的结果。具有多个响应的相同模型将如下所示：

 m1 <- MCMCglmm(cbind(resp1, resp2, resp3, resp4, resp5) ~ trait-1 + 
 at.level(trait,1):r + at.level(trait,1):present, 
 random= ~idh(at.level(trait,1)):focal + idh(at.level(trait,1)):other + 
 idh(at.level(trait,1)):X + idh(at.level(trait,1)):village, 
 rcov=~idh(trait):units, 
 family = cbind("zipoisson","zipoisson","zipoisson","zipoisson","zipoisson"), 
 prior=prior_overdisp_i, 
 data = data, nitt = nitt.1, thin = 50, burnin = 15000, pr = TRUE, pl = TRUE, verbose = TRUE, DIC  = TRUE)

我这里有两个编程问题：

1. 我该如何为这个模型指定先验概率？我看了这篇帖子中提到的材料，但还是无法理解。

2. 我运行了一个类似的版本，只有两个反应变量，但我只得到了一个斜率——我以为每个resp变量都应该有不同的斜率。我错在哪里了？还是我误解了这个模型？

- M.A.Kline

---

你有没有检查过在 R=list(V=diag(2),nu=0.08,fix=2) 中 fix = 2 是否真的有意义？在我对 MCMCglmm 先验规范的理解中，fix 应该被视为一个布尔值：fix = 0 是不将方差固定为 V 的默认值，而 fix = 1 表示“将方差固定为 V 的值”。因此，fix = 2（或类似的值）在我的看法中应该没有任何意义。（但是在 Hadfield 的课程笔记第103页中使用了这种规范：ftp://cran.r-project.org/pub/R/web/packages/MCMCglmm/vignettes/CourseNotes.pdf） - Qaswed

@Qaswed 我几年后回到这些数据，并再次查看这些模型。我的理解是，“修复”组件与先验模型的哪个部分有关...因为有一个分类组件（预测零）和一个连续组件（预测非零值）。这是特定于zipoisson模型的，它们在技术上属于多项式模型。注意：我可能会感到困惑！ - M.A.Kline

Answer 1

根据HLP帖子和一位同事/统计顾问的帮助，回答我的第一个问题：

# values for prior
k <- 5 # originally: length(levels(dative$SemanticClass)), so k = # of outcomes for SemanticClass     aka categorical outcomes 
I <- diag(k-1) #should make matrix of 0's with diagonal of 1's, dimensions k-1 rows and k-1 columns
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1)) # should make k-1 x k-1 matrix of 1's 

对于我的模型，使用multinomial5家族和5个结果变量，先验为:

prior = list(
             R = list(fix=1, V=0.5 * (I + J), n = 4),
             G = list(
               G1 = list(V = diag(4), n = 4))

对于我的第二个问题，我需要在这个模型的固定效应中添加一个交互项：

 m <- MCMCglmm(cbind(Resp1, Resp2...) ~ -1 + trait*predictorvariable,
 ...

结果显示响应变量的主效应以及响应/预测交互的后验估计值（预测变量对每个响应变量的影响）。

- M.A.Kline