Teste t pareado como um caso especial de modelagem linear de efeitos mistos

Sabemos que um teste t emparelhado é apenas um caso especial de ANOVA unidirecional de medidas repetidas (ou dentro do sujeito), bem como modelo linear de efeitos mistos, que pode ser demonstrado com a função lme (), o pacote nlme em R como mostrado abaixo.

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

Quando executo o seguinte teste t emparelhado:

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

Eu obtive este resultado (você obterá um resultado diferente por causa do gerador aleatório):

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

Com a abordagem ANOVA, podemos obter o mesmo resultado:

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

Agora posso obter o mesmo resultado no lme com o modelo a seguir, assumindo uma matriz de correlação simétrica definida positiva para as duas condições:

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

Ou outro modelo, assumindo uma simetria composta para a matriz de correlação das duas condições:

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

Com o teste t emparelhado e a ANOVA de medidas repetidas unidirecionais, posso escrever o modelo tradicional de média de células como

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

onde i indexa a condição, j indexa o sujeito, Y _ij é a variável de resposta, μ é constante para o efeito fixo para a média geral, α _i é o efeito fixo para a condição, β _j é o efeito aleatório do sujeito após N (0, σ _p ² ) (σ _p ² é variação da população) e ε _ij é residual após N (0, σ ² ) (σ ² é a variação dentro do sujeito).

Eu pensei que o modelo de média de células acima não seria apropriado para os modelos lme, mas o problema é que não consigo criar um modelo razoável para as duas abordagens lme () com a suposição da estrutura de correlação. O motivo é que o modelo lme parece ter mais parâmetros para os componentes aleatórios do que o modelo de média de células acima oferece. Pelo menos o modelo lme fornece exatamente o mesmo valor F, graus de liberdade e valor p, o que os gls não podem. Mais especificamente, gls fornece DFs incorretos devido ao fato de que ele não explica o fato de que cada sujeito tem duas observações, levando a DFs muito inflados. O modelo lme provavelmente é superparâmetro na especificação dos efeitos aleatórios, mas não sei qual é o modelo e quais são os parâmetros. Portanto, o problema ainda não foi resolvido para mim.

A equivalência dos modelos pode ser observada calculando a correlação entre duas observações do mesmo indivíduo, como segue:

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

Observe que os modelos, no entanto, não são exatamente equivalentes, pois o modelo de efeito aleatório força a correlação a ser positiva. O modelo CS e o modelo t-test / anova não.

EDIT: Existem outras duas diferenças também. Primeiro, os modelos CS e efeito aleatório assumem normalidade para o efeito aleatório, mas o modelo t-test / anova não. Em segundo lugar, os modelos CS e efeito aleatório são adequados usando a máxima verossimilhança, enquanto a anova é adequada usando quadrados médios; quando tudo está equilibrado, eles concordam, mas não necessariamente em situações mais complexas. Finalmente, eu seria cauteloso ao usar valores F / df / p dos vários ajustes como medidas de quanto os modelos concordam; veja a famosa mesa de Doug Bates em df's para mais detalhes. (END EDIT)

O problema com o seu Rcódigo é que você não está especificando a estrutura de correlação corretamente. Você precisa usar glscom a corCompSymmestrutura de correlação.

Gere dados para que haja um efeito de assunto:

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

A seguir, veja como você ajustaria os efeitos aleatórios e os modelos de simetria composta.

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

Os erros padrão do modelo de efeitos aleatórios são:

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

E a correlação e variação residual do modelo CS é:

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

E eles são iguais ao esperado:

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

Outras estruturas de correlação geralmente não se ajustam a efeitos aleatórios, mas simplesmente especificando a estrutura desejada; Uma exceção comum é o modelo de efeito aleatório AR (1) +, que tem um efeito aleatório e correlação AR (1) entre observações no mesmo efeito aleatório.

EDIT2: Quando encaixo as três opções, obtenho exatamente os mesmos resultados, exceto que o gls não tenta adivinhar o df para o termo de interesse.

> 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 **

(A interceptação é diferente aqui porque, com a codificação padrão, não é a média de todos os assuntos, mas a média do primeiro assunto.)

Também é interessante notar que o lme4pacote mais recente fornece os mesmos resultados, mas nem tenta calcular um valor-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

Você também pode considerar usar a função mixedno pacote afexpara retornar valores de p com a aproximação de Kenward-Roger df, que retorna valores de p idênticos como um teste t emparelhado:

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

Ou

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