Como simular medidas repetidas resultados multivariados em R?

9

O @whuber demonstrou como simular resultados multivariados ( , e ) por um ponto no tempo.y 2y1y2y3

Como sabemos, dados longitudinais geralmente ocorrem em estudos médicos. Minha pergunta é como simular resultados multivariados de medidas repetidas em R? Por exemplo, medimos repetidamente , e em 5 pontos no tempo para dois grupos de tratamento diferentes.y 2 y 3y1y2y3

Tu.2
fonte

Respostas:

2

Para gerar dados normais multivariados com uma estrutura de correlação especificada, é necessário construir a matriz de covariância de variância e calcular sua decomposição de Cholesky usando a cholfunção O produto da decomposição de Cholesky da matriz vcov desejada e vetores normais aleatórios independentes de observações produzirão dados normais aleatórios com essa matriz de covariância de variância.

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)
AdamO
fonte
2

Use a função rmvnorm (). São necessários 3 argumentos: a matriz de covariância das variações, as médias e o número de linhas.

O sigma terá 3 * 5 = 15 linhas e colunas. Um para cada observação de cada variável. Existem muitas maneiras de definir esses 15 ^ 2 parâmetros (ar, simetria bilateral, não estruturada ...). No entanto, você preenche essa matriz, esteja ciente das premissas, principalmente quando define uma correlação / covariância como zero ou quando define duas variações para serem iguais. Para um ponto de partida, uma matriz sigma pode ser algo como isto:

 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)

Portanto, o sigma [1,12] é 0,2 e isso significa que a covariância entre a primeira observação de Y1 e a segunda observação de Y3 é 0,2, condicional a todas as outras 13 variáveis. Nem todas as linhas diagonais precisam ter o mesmo número: essa é uma suposição simplificadora que eu fiz. Às vezes faz sentido, às vezes não. Em geral, significa que a correlação entre uma 3ª observação e uma 4ª é igual à correlação entre uma 1ª e uma segunda.

Você também precisa de meios. Poderia ser tão simples quanto

 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)

Aqui os 5 primeiros são os meios para as 5 observações de Y1, ..., os 5 últimos são as observações de Y3

obtenha 2000 observações de seus dados com:

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

Onde Y11 é a 1ª observação de Y1, ..., Y15 é o 5º obs de Y1 ...

Seth
fonte
11
Para criar matrizes como no primeiro exemplo, Seth, tente o seguinte: n <- 3*5; sigma <- diag(1, nrow=n, ncol=n); sigma[rbind(cbind(1:n-1,1:n),cbind(1:n,1:n-1))] <- 1/2. Uma abordagem semelhante irá gerar o segundo exemplo. No entanto, eles têm um problema comum: você perdeu as covariâncias entre os 's durante cada período - essas matrizes não refletem uma estrutura de medidas repetidas. y
whuber
@whuber, sua sintaxe é útil, mas é diferente do que escrevi. Eu acho que a diferença é um pouco importante. Penso no que escrevi como AR (1) e você tem uma entrada na correlação cruzada entre a última observação de uma variável e a primeira observação da próxima variável. Em outras palavras, eu acho sigma [5,6] deve ser 0.
Seth
Ah, agora eu vejo o que você está fazendo: você está criando três séries AR (1) no primeiro exemplo. Perdi isso porque acredito que o OP também se preocupa com correlações entre as séries: é isso que se entende por "resultados multivariados". A partir desse ponto de vista, parece que você deseja ver estas matrizes de correlação como sendo por matrizes de bloco com cada entrada de um por matriz, em vez de como um por matriz de bloco com por blocos. 5 3 3 3 3 5 555333355
whuber
Eu pensei que meu segundo sigma era um exemplo simples de permitir que a variação entre Y1 e Y3 fosse positiva. Editei a resposta um pouco para deixar mais claro que a matriz existe para ser configurada, dependendo do processo de geração de dados. Definitivamente, existem muitas maneiras de esfolar esse gato.
Seth
É justo, mas sua abordagem cria dificuldades, porque não é trivial combinar a correlação multivariada entre o e um modelo de AR. Por exemplo, você sabia que a segunda matriz falha em ser positiva definida? (A "variação" de c (-102, 177, -204, 177, -102, 0, 0, 0, 0, 0, 102, -177, 204, -177, 102) é negativa.)yi
whuber