Wie simuliert man multivariate Ergebnisse mit wiederholten Messungen in R?

@whuber hat gezeigt, wie multivariate Ergebnisse ( , und y_3 ) für einen Zeitpunkt simuliert werden können.$y_1$$y_2$$y_3$

Wie wir wissen, treten in medizinischen Studien häufig Längsschnittdaten auf. Meine Frage ist, wie man multivariate Ergebnisse mit wiederholten Messungen in R simuliert. Zum Beispiel messen wir wiederholt $y_1$ , $y_2$ und $y_3$ zu 5 verschiedenen Zeitpunkten für zwei verschiedene Behandlungsgruppen.

Um multivariate Normaldaten mit einer bestimmten Korrelationsstruktur zu generieren, müssen Sie die Varianz-Kovarianz-Matrix erstellen und ihre Cholesky-Zerlegung mithilfe der cholFunktion berechnen . Das Produkt der Cholesky-Zerlegung der gewünschten vcov-Matrix und unabhängiger zufälliger normaler Beobachtungsvektoren ergibt zufällige normale Daten mit dieser Varianz-Kovarianzmatrix.

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)

Verwenden Sie die Funktion rmvnorm (). Es werden 3 Argumente benötigt: die Varianz-Kovarianz-Matrix, die Mittelwerte und die Anzahl der Zeilen.

Das Sigma hat 3 * 5 = 15 Zeilen und Spalten. Eine für jede Beobachtung jeder Variablen. Es gibt viele Möglichkeiten, diese 15 ^ 2 Parameter einzustellen (ar, bilaterale Symmetrie, unstrukturiert ...). Wenn Sie diese Matrix ausfüllen, müssen Sie sich der Annahmen bewusst sein, insbesondere wenn Sie eine Korrelation / Kovarianz auf Null setzen oder wenn Sie zwei Varianzen gleich setzen. Als Ausgangspunkt könnte eine Sigma-Matrix ungefähr so ​​aussehen:

 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)

Das Sigma [1,12] ist also 0,2 und das bedeutet, dass die Kovarianz zwischen der ersten Beobachtung von Y1 und der zweiten Beobachtung von Y3 0,2 beträgt, abhängig von allen anderen 13 Variablen. Die diagonalen Reihen müssen nicht alle dieselbe Zahl sein: Das ist eine vereinfachende Annahme, die ich gemacht habe. Manchmal macht es Sinn, manchmal nicht. Im Allgemeinen bedeutet dies, dass die Korrelation zwischen einer dritten und einer vierten Beobachtung dieselbe ist wie die Korrelation zwischen einer ersten und einer zweiten.

Sie brauchen auch Mittel. Es könnte so einfach sein wie

 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)

Hier sind die ersten 5 die Mittelwerte für die 5 Beobachtungen von Y1, ..., die letzten 5 sind die Beobachtungen von Y3

Dann erhalten Sie 2000 Beobachtungen Ihrer Daten mit:

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

Wo Y11 die 1. Beobachtung von Y1 ist, ..., ist Y15 die 5. Beobachtung von Y1 ...