如何在 R 中模拟重复测量的多变量结果?

机器算法验证 r 重复测量 模拟
2022-03-13 03:17:40

@whuber 演示了如何模拟一个时间点的多变量结果)。y1y2y3

众所周知,纵向数据经常出现在医学研究中。我的问题是如何在 R 中模拟重复测量的多元结果?例如,我们在 5 个不同的时间点对两个不同的治疗组y1y2y3

2个回答

使用 rmvnorm() 函数,它需要 3 个参数:方差协方差矩阵、均值和行数。

sigma 将有 3*5=15 行和列。每个变量的每个观察值一个。有很多方法可以设置这 15^2 个参数(ar、双边对称、非结构化......)。但是,您填写此矩阵时请注意假设,特别是当您将相关/协方差设置为零时,或者当您将两个方差设置为相等时。作为起点,sigma 矩阵可能看起来像这样:

 sigma=matrix(c(
    #y1             y2             y3 
    3 ,.5, 0, 0, 0, 0, 0, 0, 0, 0,.5,.2, 0, 0, 0,
    .5, 3,.5, 0, 0, 0, 0, 0, 0, 0,.2,.5,.2, 0, 0,
    0 ,.5, 3,.5, 0, 0, 0, 0, 0, 0, 0,.2,.5,.2, 0,
    0 , 0,.5, 3,.5, 0, 0, 0, 0, 0, 0, 0,.2,.5,.2,
    0 , 0, 0,.5, 3, 0, 0, 0, 0, 0, 0, 0, 0,.2,.5,
    0 ,0 ,0 ,0 , 0, 3,.5, 0, 0, 0, 0, 0, 0, 0, 0,
    0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0, 0, 0, 0, 0, 0,
    0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0, 0, 0, 0, 0,
    0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0, 0, 0, 0,
    0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3, 0, 0, 0, 0, 0,
    .5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 , 0, 3,.5, 0, 0, 0,
    .2,.5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0,
    0 ,.2,.5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0,
    0 ,0 ,.2,.5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5,
    0 ,0 ,0 ,.2,.5,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3

    ),15,15)

所以 sigma[1,12] 是 0.2,这意味着 Y1 的第一次观察和 Y3 的第二次观察之间的协方差是 0.2,以所有其他 13 个变量为条件。对角线行不必都是相同的数字:这是我所做的简化假设。有时它是有道理的,有时它没有。一般来说,这意味着第 3 次观察和第 4 次观察之间的相关性与第 1 次和第 2 次观察之间的相关性相同。

你也需要手段。它可能很简单

 meanTreat=c(1:5,51:55,101:105)
 meanControl=c(1,1,1,1,1,50,50,50,50,50,100,100,100,100,100)

这里前 5 个是 Y1 的 5 个观测值的均值,...,后 5 个是 Y3 的观测值

然后通过以下方式对您的数据进行 2000 次观察:

sampleT=rmvnorm(1000,meanTreat,sigma)
sampleC=rmvnorm(1000,meanControl,sigma)
 sample=data.frame(cbind(sampleT,sampleC) )
  sample$group=c(rep("Treat",1000),rep("Control",1000) )

colnames(sample)=c("Y11","Y12","Y13","Y14","Y15",
                   "Y21","Y22","Y23","Y24","Y25",
                   "Y31","Y32","Y33","Y34","Y35")

其中 Y11 是 Y1 的第一个观测值,...,Y15 是 Y1 的第 5 个观测值...

要生成具有指定相关结构的多元正态数据,您需要构建方差协方差矩阵并使用该chol函数计算其 Cholesky 分解。所需 vcov 矩阵的 Cholesky 分解与观察的独立随机法线向量的乘积将产生具有该方差协方差矩阵的随机法线数据。

v <- matrix(c(2,.3,.3,2), 2)
cv <- chol(v)

o <- replicate(1000, {
  y <- cv %*% matrix(rnorm(100),2)

  v1 <- var(y[1,])
  v2 <- var(y[2,])
  v3 <- cov(y[1,], y[2,])

  return(c(v1,v2,v3))
})

## MCMC means should estimate components of v
rowMeans(o)