MCMCglmm 中的可交换相关性和恒定方差？

问题描述

我正试图用 MCMCglmm 拟合一个非常简单的模型，但我陷入了困境。

假设一个班级（30 名学生）在整个学期获得两篇论文的成绩，其中论文作业完全相同（我们不想模拟论文之间的平均分数差异，没有“学习效果” ，我们可以假设成绩的方差是相同的。）

让 $i = 1...30$ 索引该学生，$y_{i1}$ 和 $y_{i2}$ 索引该学生第一篇和第二篇论文的分数。

建模此数据的一种方法是对学生分数使用随机截距，以说明每个学生分数之间的相关性。令 $\mu_i$ 为学生截距，$sigma$ 为残差 sd，$\sigma_{\mu}$ 为截距的 sd。然后我们在 $f(y_{ij}|\mu_i) = Normal(\mu_i,\sigma)$ 和 $f(\mu_i) = Normal(\mu,\sigma_{\亩)$。

编写此模型的另一种方法是更明确地对残差相关结构进行建模。也就是说，我们会写出 ${y_{i1},y_{i2}}$ 有一个多元正态分布，均值 ${\mu,\mu}$ 方差 $\tau = \sigma^2 + \sigma_{\ mu}^2$ 和相关性 $\rho = \sigma_{\mu}^2 / (\sigma^2 + \sigma_{\mu}^2)$.

需要明确的是，这些模型在数学上是等价的，但统计软件通常会为每个模型提供特定的实现。例如，我们可以用 nlme 分别拟合这两种方法：

library(nlme)
library(tidyverse)
library(MCMCglmm)
df <-
tibble(id = factor(rep(1:100,each = 20))) %>%
mutate(paper = 1:n()) %>%
group_by(id) %>%
mutate(mu = rnorm(1),y = mu + rnorm(n(),3))

gls(data = df,model = y~1,correlation = corCompSymm(form = ~ 1 | id))

lme(data = df,fixed = y ~ 1,random = ~1|id)

似乎 MCMCglmm 可以很好地拟合模型的第一个参数化（随机截距）。

MCMCglmm(data = df,random = ~id,nitt = 1000,burnin = 0,thin = 1)

但是，我没有看到实现第二种方法的方法。我的最佳尝试包括“扩大”数据框并拟合多响应模型。

df.wide <- df %>% select(- paper) %>%
pivot_wider(values_from = "y",names_from = "obs",names_prefix = "paper") %>%
    as.data.frame

MCMCglmm(fixed = cbind(paper1,paper2) ~ 1,rcov = ~us(trait):units,data = df.wide)

但是，（1）我不确定我是否正确地拟合了这个模型，（2）我不确定如何解释拟合值（特别是因为我的后验平均协方差似乎太小了）和（3）那里似乎不是在特征之间获得恒定差异的方法。

附言我很感激不被告知只适合随机截距模型。我正在写一些课程材料，希望学生能够更直接地将可交换相关模型与我们在有两个以上观察（即 AR、Toeplitz 等）时可能使用的其他类型的相关结构进行比较，以及我希望我的学生能够自己进行两种参数化的比较，就像我使用 nlme 时所做的那样。

后续：我目前正在尝试使用 BRMS 拟合模型，但仍然对 MCMCglmm 中的任何“黑客”持开放态度。

model1 <- brms::brm(data = df,formula = y ~ 1 + cosy(gr = id,time = obs),family =  "gaussian",chains = 4,thin = 1,iter = 5000,warmup = 100)

解决方法

可交换性 + 等方差是否与我所说的复合对称性相同？ （我想是的，因为您在 corCompSymm() 中使用了 nlme）... 据我所知这是不可能的（我不能排除有一些方法可以用可用的方差结构来破解它，但这远非显而易见......）来自?MCMCglmm：

目前，唯一可用的‘variance.functions’是‘idv’， “idh”、“us”、“cor[]”和“ante[]”。 ‘idv’适合一个常数 “公式”中所有组件的差异。 “idh”和 “我们”在每个组件中拟合不同的差异 “公式”，但“我们”也将适合协方差。 ‘corg’ 将沿对角线的方差固定为 1 和“corgh” 将沿对角线的方差固定为中指定的那些 之前的。 ‘cors’允许相关子矩阵。 '赌注[]' 适合不同顺序的前依存结构（例如 ante1,ante2)，并且数字可以以“c”为前缀 保持所有相同阶次的回归系数相等。这 数字也可以后缀“v”来保存所有创新 方差相等（例如‘antec2v’有3个参数）。

通过使用 us()（非结构化，nlme 将 pdSymm 称为“正定对称”）结构，我相信您不是将相关参数限制为完全相同（即违反可交换性）。

就其价值而言，想要明确指定复合对称相关矩阵而不是通过组合组级和个人级随机效应的总和的一个原因（教学法除外）是，如果您想建模负复合对称性（随机效应总和方法只能模拟rho>0）。

我的猜测是您也只能使用 MCMCglmm 来回答，但是如果“一些贝叶斯 MCMC 方法”足够好，那么您可以通过 brms 或（有点模糊，有点）glmmTMB + tmbstan（尽管这种组合目前不使用信息先验！）