MCMCglmm?中的r-可交换相关性和常方差



我试图用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_id(=Normal(\mu,\sigma\mu($处写下(简写(我们的随机截距模型。

编写该模型的另一种方法是更明确地对残差相关结构进行建模。也就是说,我们将写${y_{i1},y_{i2}}$具有多元正态分布,其均值${\mu,\mu}$方差$\tau=\sigma_2+\sigma_{\mu}^2$和相关性$\rho=\sigmon_{\mu}^2/(\sigma\^2+\sigmon\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(), 0, 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, 
fixed = y ~ 1, 
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)

可交换性+等方差与我所说的复合对称性相同吗?(我想是的,因为您在nlme中使用corCompSymm()(。。。据我所知,这是不可能的(我不能排除有一些方法可以用可用的方差结构破解它,但这还远不明显…(来自?MCMCglmm:

目前,唯一可用的"方差函数"是"idv","idh"、"us"、"cor[]"one_answers"ante[]"idv'符合常数"公式"中所有成分的差异。"idh"和我们在"formula",但"us"也将拟合协变量corg'将对角线上的方差固定为1和"corgh"将对角线上的方差固定为中指定的方差"cors允许关联子矩阵ante[]'适合不同阶数的事前依赖结构(例如ante1,ante2(,并且该数字可以以"c"为前缀保持相同阶数的所有回归系数相等。这个数字后面也可以加一个"v"来表示所有创新方差相等(例如"antec2v"有3个参数(。

通过使用us()(非结构化,nlme将其称为"正定对称"的pdSymm(结构,我相信您不是将相关参数约束为完全相同(即,违反可交换性(。

值得一提的是,如果你想对复合对称性进行建模(随机效应之和方法只能对rho>0进行建模(,那么想要明确指定复合对称相关矩阵,而不是通过组合组级和个人级随机效应的总和。

我的猜测是,你也被限制为使用MCMCglmm的答案,但如果";一些贝叶斯MCMC方法";已经足够好了,那么您可以通过brms(有点模糊,有点像(glmmTMB + tmbstan(尽管这种组合目前没有使用信息先验!(

最新更新