Comment simuler des résultats multivariés de mesures répétées dans R?

@whuber a montré comment simuler des résultats multivariés ( , y 2 et$y_1$$y_2$ ) pour un point dans le temps.$y_3$

Comme nous le savons, les données longitudinales se produisent souvent dans les études médicales. Ma question est de savoir comment simuler des résultats multivariés de mesures répétées dans R? Par exemple, nous mesurons à plusieurs reprises , y 2 et y $y_1$$y_2$$y_3$ à 5 moments différents pour deux groupes de traitement différents.

Pour générer des données normales multivariées avec une structure de corrélation spécifiée, vous devez construire la matrice de covariance de variance et calculer sa décomposition de Cholesky en utilisant la cholfonction. Le produit de la décomposition de Cholesky de la matrice vcov souhaitée et des vecteurs normaux aléatoires indépendants des observations produira des données normales aléatoires avec cette matrice de covariance de variance.

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)

Utilisez la fonction rmvnorm (), elle prend 3 arguments: la matrice de covariance de variance, les moyennes et le nombre de lignes.

Le sigma aura 3 * 5 = 15 lignes et colonnes. Un pour chaque observation de chaque variable. Il existe de nombreuses façons de régler ces 15 ^ 2 paramètres (ar, symétrie bilatérale, non structurée ...). Quelle que soit la manière dont vous remplissez cette matrice, tenez compte des hypothèses, en particulier lorsque vous définissez une corrélation / covariance sur zéro, ou lorsque vous définissez deux variances pour qu'elles soient égales. Pour un point de départ, une matrice sigma pourrait ressembler à ceci:

 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)

Le sigma [1,12] est donc de 0,2 et cela signifie que la covariance entre la première observation de Y1 et la deuxième observation de Y3 est de 0,2, conditionnelle à toutes les 13 autres variables.Les lignes diagonales ne doivent pas toutes être du même nombre: c'est une hypothèse simplificatrice que j'ai faite. Parfois, cela a du sens, parfois non. En général, cela signifie que la corrélation entre une 3ème observation et une 4ème est la même que la corrélation entre une 1ère et une seconde.

Vous avez également besoin de moyens. Cela pourrait être aussi simple que

 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)

Ici les 5 premiers sont les moyennes des 5 observations de Y1, ..., les 5 derniers sont les observations de Y3

puis obtenez 2000 observations de vos données avec:

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")

Où Y11 est la 1ère observation de Y1, ..., Y15 est la 5ème obs de Y1 ...