Użyj funkcji rmvnorm (), wymaga 3 argumentów: macierzy kowariancji wariancji, średnich i liczby wierszy.
Sigma będzie miała 3 * 5 = 15 wierszy i kolumn. Jeden na każdą obserwację każdej zmiennej. Istnieje wiele sposobów ustawiania tych 15 ^ 2 parametrów (ar, dwustronna symetria, nieustrukturyzowany ...). Jednak wypełniając tę macierz, należy pamiętać o założeniach, zwłaszcza gdy ustawisz korelację / kowariancję na zero lub gdy ustawisz dwie wariancje na równe. Na początek macierz sigma może wyglądać mniej więcej tak:
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)
Zatem sigma [1,12] wynosi .2, a to oznacza, że kowariancja między pierwszą obserwacją Y1 a drugą obserwacją Y3 wynosi .2, zależnie od wszystkich pozostałych 13 zmiennych. Nie wszystkie rzędy po przekątnej muszą mieć tę samą liczbę: to moje uproszczone założenie. Czasami ma to sens, a czasem nie. Zasadniczo oznacza to, że korelacja między trzecią obserwacją a czwartą jest taka sama jak korelacja między pierwszą a drugą.
Potrzebujesz także środków. To może być tak proste jak
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)
Tutaj pierwsze 5 jest średnimi dla 5 obserwacji Y1, ... ostatnich 5 to obserwacje Y3
uzyskaj 2000 obserwacji swoich danych za pomocą:
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")
Gdzie Y11 jest pierwszą obserwacją Y1, ..., Y15 jest piątą obsadą Y1 ...
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
. Podobne podejście wygeneruje drugi przykład. Mają jednak wspólny problem: utraciłeś kowariancje wśród