Sparowany test t jako specjalny przypadek liniowego modelowania z mieszanym efektem

Wiemy, że sparowany test t jest tylko specjalnym przypadkiem jednokierunkowej ANOVA z powtarzanymi pomiarami (lub wewnątrz podmiotu), a także liniowym modelem mieszanego efektu, który można zademonstrować za pomocą funkcji lme () pakiet nlme w języku R jak pokazano niżej.

#response data from 10 subjects under two conditions
x1<-rnorm(10)
x2<-1+rnorm(10)

# Now create a dataframe for lme
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Kiedy uruchamiam następujący sparowany test t:

t.test(x1, x2, paired = TRUE)

Mam ten wynik (otrzymasz inny wynik z powodu losowego generatora):

t = -2.3056, df = 9, p-value = 0.04657

Dzięki podejściu ANOVA możemy uzyskać ten sam wynik:

summary(aov(y ~ x + Error(subj/x), myDat))

# the F-value below is just the square of the t-value from paired t-test:
          Df  F value Pr(>F)
x          1  5.3158  0.04657

Teraz mogę uzyskać ten sam wynik w lme przy użyciu następującego modelu, zakładając dodatnio określoną symetryczną macierz korelacji dla dwóch warunków:

summary(fm1 <- lme(y ~ x, random=list(subj=pdSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.3142115  9 -0.7918878  0.4488
# xx2          1.3325786 0.5779727  9  2.3056084  0.0466

Lub inny model, zakładając złożoną symetrię dla macierzy korelacji dwóch warunków:

summary(fm2 <- lme(y ~ x, random=list(subj=pdCompSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.4023431  9 -0.618428  0.5516
# xx2          1.3325786 0.5779727  9  2.305608  0.0466

Za pomocą sparowanego testu t i jednokierunkowej ANOVA z powtarzanymi pomiarami mogę zapisać tradycyjny model średnich komórek jako

Yij = μ + αi + βj + εij, i = 1, 2; j = 1, ..., 10

gdzie i indeksuje warunek, j indeksuje podmiot, Y _ij jest zmienną odpowiedzi, μ jest stały dla ustalonego efektu dla średniej ogólnej, α _i jest ustalonym efektem dla warunku, β _j jest efektem losowym dla badanego po N (0, σ _p ² ) (σ _p ² to wariancja populacyjna), a ε _ij jest resztkowe po N (0, σ ² ) (σ ² jest wariancją wewnątrzosobniczą).

Myślałem, że powyższy model średniej komórkowej nie byłby odpowiedni dla modeli lme, ale problem polega na tym, że nie mogę wymyślić rozsądnego modelu dla dwóch podejść lme () przy założeniu struktury korelacji. Powodem jest to, że model lme wydaje się mieć więcej parametrów dla składników losowych niż model średniej komórki powyżej oferuje. Przynajmniej model lme zapewnia dokładnie tę samą wartość F, stopnie swobody, a także wartość p, czego gls nie może. Mówiąc dokładniej, gls podaje nieprawidłowe DF ze względu na fakt, że nie uwzględnia faktu, że każdy badany ma dwie obserwacje, co prowadzi do znacznie zawyżonego DF. Model lme najprawdopodobniej jest sparametryzowany w określaniu losowych efektów, ale nie wiem, co to jest model i jakie są parametry. Tak więc problem nadal jest dla mnie nierozwiązany.

Równoważność modeli można zaobserwować, obliczając korelację między dwoma obserwacjami tej samej osoby, w następujący sposób:

$Y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}$$\beta_j \sim N(0, \sigma_p^2)$$\epsilon_{ij} \sim N(0, \sigma^2)$$Cov(y_{ik}, y_{jk}) = Cov(\mu + \alpha_i + \beta_k + \epsilon_{ik}, \mu + \alpha_j + \beta_k + \epsilon_{jk}) = Cov(\beta_k, \beta_k) = \sigma_p^2$$Var(y_{ik}) = Var(y_{jk}) = \sigma_p^2 + \sigma^2$$\sigma_p^2/(\sigma_p^2 + \sigma^2)$

Zauważ, że modele nie są jednak całkowicie równoważne, ponieważ model efektu losowego wymusza korelację na dodatnią. Model CS i model t-test / anova nie.

EDYCJA: Istnieją również dwie inne różnice. Po pierwsze, modele CS i losowe zakładają normalność efektu losowego, ale model t-test / anova nie. Po drugie, modele CS i losowe są dopasowane przy użyciu maksymalnego prawdopodobieństwa, podczas gdy anova jest dopasowana przy użyciu średnich kwadratów; kiedy wszystko jest zrównoważone, zgodzą się, ale niekoniecznie w bardziej złożonych sytuacjach. Wreszcie, będę ostrożny w stosowaniu wartości F / df / p z różnych dopasowań jako miary tego, jak bardzo modele się zgadzają; zobacz słynny stół Douga Batesa na df po więcej szczegółów. (EDYCJA KOŃCOWA)

Problem z twoim Rkodem polega na tym, że nie określasz poprawnie struktury korelacji. Musisz używać glsze corCompSymmstrukturą korelacji.

Wygeneruj dane, aby uzyskać efekt podmiotowy:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Oto jak dopasujesz losowe efekty i złożone modele symetrii.

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

Standardowe błędy z modelu efektów losowych to:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

Korelacja i wariancja rezydualna z modelu CS jest następująca:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

I są równe oczekiwanym:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Inne struktury korelacji zwykle nie pasują do efektów losowych, ale po prostu przez określenie pożądanej struktury; jednym powszechnym wyjątkiem jest model efektu losowego AR (1) +, który ma efekt losowy i korelację AR (1) między obserwacjami tego samego efektu losowego.

EDYCJA 2: Kiedy dopasuję trzy opcje, otrzymuję dokładnie takie same wyniki, z tym wyjątkiem, że gls nie próbuje zgadnąć df dla interesującego nas terminu.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(Punkt przecięcia jest tutaj inny, ponieważ przy domyślnym kodowaniu nie jest to średnia z wszystkich przedmiotów, ale średnia z pierwszego przedmiotu).

Warto również zauważyć, że nowszy lme4pakiet daje takie same wyniki, ale nawet nie próbuje obliczyć wartości p.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283

Możesz również rozważyć użycie funkcji mixedw pakiecie afexdo zwrócenia wartości p z aproksymacją df Kenwarda-Rogera, która zwraca identyczne wartości p jako sparowany test t:

library(afex)
mixed(y ~ x + (1|subj), type=3,method="KR",data=myDat) 

Lub

library(lmerTest)
options(contrasts=c('contr.sum', 'contr.poly'))
anova(lmer(y ~ x + (1|subj),data=myDat),ddf="Kenward-Roger")