Por que é preciso usar REML (em vez de ML) para escolher entre modelos var-covar aninhados?

16

Várias descrições na seleção de modelos sobre efeitos aleatórios de Modelos Mistos Lineares instruem a usar REML. Conheço a diferença entre REML e ML em algum nível, mas não entendo por que REML deve ser usado porque ML é tendencioso. Por exemplo, é errado realizar um LRT em um parâmetro de variação de um modelo de distribuição normal usando ML (veja o código abaixo)? Não entendo por que é mais importante ser imparcial do que ser ML, na seleção de modelos. Acho que a resposta final deve ser "porque a seleção de modelos funciona melhor com REML do que com ML", mas eu gostaria de saber um pouco mais do que isso. Não li as derivações de LRT e AIC (não sou bom o suficiente para entendê-las completamente), mas se REML for usado explicitamente nas derivações, apenas sabendo que será realmente suficiente (por exemplo,

n <- 100
a <- 10
b <- 1
alpha <- 5
beta <- 1
x <- runif(n,0,10)
y <- rnorm(n,a+b*x,alpha+beta*x)

loglik1 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
  -sum(dnorm(y,a+b*x,alpha,log=T))
}

loglik2 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
   beta <- p[4]
  -sum(dnorm(y,a+b*x,alpha+beta*x,log=T))
}

m1 <- optim(c(a,b,alpha),loglik1,x=x,y=y)$value
m2 <- optim(c(a,b,alpha,beta),loglik2,x=x,y=y)$value
D <- 2*(m1-m2)
1-pchisq(D,df=1) # p-value
queixa
fonte
1
Sobre REML e AIC, você deve dar uma olhada nesta pergunta .
Elvis

Respostas:

13

Uma resposta muito curta: o REML é um ML; portanto, o teste baseado no REML está correto de qualquer maneira. Como a estimativa dos parâmetros de variação com REML é melhor, é natural usá-lo.

Por que REML é um ML? Considere, por exemplo, um modelo com X R n × p , Z R n × q e β R p é o vetor dos efeitos fixos, u N ( 0 , τ I q ) é o vetor de efeitos aleatórios e e N ( 0 , σ 2 I n )

Y=Xβ+Zu+e
XRn×pZRn×qβRpuN(0,τIq)eN(0,σ2In). A Probabilidade restrita pode ser obtida considerando -se os contrastes para "remover" os efeitos fixos. Mais precisamente, deixe C R ( n - p ) × n , de modo que C X = 0 e C C = I n - p (ou seja, as colunas de C são uma base ortonormal do espaço vetorial ortogonal ao espaço gerado pelas colunas de X ); então C Y = C Z u +npCR(np)×nCX=0CC=InpCX com ε ~ N ( 0 , σ 2 I n - p ) , e a probabilidade para τ , σ 2 dado C Y é a probabilidade restrita.
CY=CZu+ϵ
ϵN(0,σ2Inp)τ,σ2CY
Elvis
fonte
Boa resposta (+1), estou correto em dizer que a matriz é dependente do modelo para a média? Então você só pode comparar estimativas REML para a mesma matriz C ? CC
Sim, depende de X (editarei a resposta em um minuto para deixar claro), portanto seus modelos aninhados precisam ter as mesmas variáveis ​​com efeitos fixos. CX
Elvis
REML não é um ML! O ML é definido exclusivamente para um determinado modelo de probabilidade, mas o REML depende da parametrização de efeitos fixos. Veja, por exemplo, este comentário de Doug Bates (bem como muitos históricos sobre modelos mistos R-SIG).
Livius
1
@Livius Acho que minha resposta afirma com clareza suficiente como a probabilidade restrita é construída. Ele é uma probabilidade, não é apenas a possibilidade dada ao observado no modelo de escrita na primeira equação apresentada, mas dado o vector projectado C Y no modelo escrito na segunda equação apresentada. REML é o ML obtido com essa probabilidade. YCY
Elvis
2
Eu acho que esse é o ponto dos protestos de DBates sobre esse assunto: é um modelo diferente e é um modelo para o qual as comparações são difíceis porque o modelo e a parametrização estão entrelaçados. Portanto, você não está computando o ML para o seu modelo original, mas o ML para um modelo diferente, resultante de uma parametrização específica do seu modelo original. Portanto, os modelos equipados com REML com estruturas de efeitos fixos aninhados não são mais modelos aninhados (como você mencionou acima). Mas os modelos equipados com ML ainda estão aninhados, porque você está maximizando a probabilidade no modelo especificado.
Livius
9

Os testes de razão de verossimilhança são testes estatísticos de hipóteses baseados em uma razão de duas verossimilhanças. Suas propriedades estão ligadas à estimativa de máxima verossimilhança (MLE). (ver, por exemplo, estimativa de máxima verossimilhança (MLE) em termos leigos ).

No seu caso (consulte a pergunta), você deseja '' escolher '' entre dois modelos var-covar aninhados, digamos que você queira escolher entre um modelo em que o var-covar é e um modelo em que o var-covar é Σ s onde o segundo (modelo simples) é um caso especial do primeiro (o geral). ΣgΣs

O teste baseia-se na probabilidade proporção . Onde Σ s e Σ g são o estimadores de probabilidade máxima.LR=2(log(Ls(Σ^s))log(Lg(Σ^g))Σ^sΣ^g

A estatística é, assintoticamente (!) × 2 . LR χ2

Sabe-se que os estimadores de probabilidade máxima são consistentes; no entanto, em muitos casos, eles são tendenciosos. Este é o caso para os estimadores MLE para a e Σ g, pode ser mostrar que eles são tendenciosos. Isso ocorre porque eles são calculados usando uma média derivada dos dados, de modo que a dispersão em torno dessa "média estimada" é menor que a dispersão em torno da média verdadeira (consulte Por exemplo,explicação intuitiva para dividir porn-1ao calcular o desvio padrão ?)Σ^sΣ^gn1

A estatística acima é χ 2LRχ2 em amostras de grandes dimensões, isto é só por causa do facto de que, em grandes e Σ g convergem para os verdadeiros valores (MLE são consistentes). (Nota: no link acima, para amostras muito grandes, dividindo por n ou por (n-1), não fará diferença)Σ^sΣ^g

Para amostras mais pequenas, a MLE estima de Σ s e Σ g serão polarizados e, por conseguinte, a distribuição de L R irá desviar-se do χ 2Σ^sΣ^gLRχ2 , enquanto que as estimativas REML vai dar estimativas imparciais para e Σ g , de modo que se estiver a utilizar , para a selecção do modelo var-covar, o REML calcula, em seguida, o L I vai para amostras menores ser melhor aproximada pela χ 2 .ΣsΣgLRχ2

Observe que REML deve ser usado apenas para escolher entre estruturas var-covar aninhadas de modelos com a mesma média; para modelos com médias diferentes, a REML não é apropriada; para modelos com diferentes formas, deve-se usar ML.


fonte
A afirmação "A estatística LR é, assintoticamente (!) Χ2" não é verdadeira neste caso. Isso ocorre porque se está aninhado em Σ g , então Σ s está no limite de Σ g . Nesse caso, a distribuição χ 2 não se mantém. Por exemplo, veja aquiΣsΣgΣsΣgχ2
Cliff AB
@Cliff AB, é isso que é explicado abaixo dessa declaração e é o motivo pelo qual você deve usar o REML.
-4

Eu tenho uma resposta que tem mais a ver com senso comum do que com estatística. Se você der uma olhada no PROC MIXED no SAS, a estimativa poderá ser realizada com seis métodos:

http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect008.htm

mas REML é o padrão. Por quê? Aparentemente, a experiência prática mostrou que tem o melhor desempenho (por exemplo, a menor chance de problemas de convergência). Portanto, se seu objetivo for atingível com o REML, faz sentido usar o REML em oposição aos outros cinco métodos.

James
fonte
2
Isso tem a ver com a "teoria das grandes amostras" e a parcialidade das estimativas do MLE, veja minha resposta.
1
"É o padrão no SAS" não é uma resposta aceitável para uma pergunta "por que" neste site.
Paul
Os valores p para modelos mistos fornecidos pelo SAS por padrão não estão disponíveis por design na biblioteca lme4 para R porque não são confiáveis ​​( stat.ethz.ch/pipermail/r-help/2006-May/094765.html ). Portanto, o "SAS padrão" pode estar errado.
Tim