obtendo graus de liberdade do lmer

11

Eu ajustei um modelo mais recente com o seguinte (embora composto):

Random effects:
 Groups        Name        Std.Dev.
 day:sample (Intercept)    0.09
 sample        (Intercept) 0.42
 Residual                  0.023 

Eu realmente gostaria de criar um intervalo de confiança para cada efeito usando a seguinte fórmula:

(n1)s2χα/2,n12,(n1)s2χ1α/2,n12

Existe uma maneira de obter convenientemente os graus de liberdade?

user1357015
fonte
1
Eu acho que você deseja verificar o lmerTest . Existem várias aproximações para aproximar o df em um modelo de efeitos mistos para os efeitos fixos (por exemplo , Satterthwaite , Kenward-Roger, etc.) Para o seu caso, parece-me que você complicou demais sua vida. Você assume que cada efeito é gaussiano. Basta usar o desvio padrão para obter o intervalo de confiança de sua escolha.
usar o seguinte comando
3
@ usεr11852 Em um modelo de efeitos mistos, você assume que cada efeito é gaussiano, mas o parâmetro é a variação da distribuição gaussiana, não a média. Portanto, a distribuição de seu estimador será muito distorcida e o intervalo de confiança normal dos desvios padrão de ± ~ 2 não será apropriado.
Karl Ove Hufthammer
1
@KarlOveHufthammer: Você tem razão em apontar isso; Entendo o que você (e provavelmente o OP) quer dizer. Eu pensei que ele estava preocupado com os meios e / ou as realizações dos efeitos aleatórios ao mencionar graus de liberdade.
precisa saber é o seguinte
graus de liberdade são "problemáticos" para modelos mistos, consulte: stat.ethz.ch/pipermail/r-help/2006-May/094765.html e stats.stackexchange.com/questions/84268/…
Tim

Respostas:

17

Em vez disso, eu apenas criaria intervalos de confiança de probabilidade de perfil . Eles são confiáveis ​​e muito fáceis de calcular usando o pacote 'lme4'. Exemplo:

> library(lme4)
> fm = lmer(Reaction ~ Days + (Days | Subject),
            data=sleepstudy)
> summary(fm)
[]
Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.09   24.740       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       

Agora você pode calcular os intervalos de confiança da probabilidade do perfil com a confint()função:

> confint(fm, oldNames=FALSE)
Computing profile confidence intervals ...
                               2.5 %  97.5 %
sd_(Intercept)|Subject        14.381  37.716
cor_Days.(Intercept)|Subject  -0.482   0.685
sd_Days|Subject                3.801   8.753
sigma                         22.898  28.858
(Intercept)                  237.681 265.130
Days                           7.359  13.576

Você também pode usar a auto-inicialização paramétrica para calcular intervalos de confiança. Aqui está a sintaxe R (usando o parmargumento para restringir para quais parâmetros queremos intervalos de confiança):

> confint(fm, method="boot", nsim=1000, parm=1:3)
Computing bootstrap confidence intervals ...
                              2.5 % 97.5 %
sd_(Intercept)|Subject       11.886 35.390
cor_Days.(Intercept)|Subject -0.504  0.929
sd_Days|Subject               3.347  8.283

Os resultados variam naturalmente um pouco para cada execução. Você pode aumentar nsimpara diminuir essa variação, mas isso também aumentará o tempo necessário para estimar os intervalos de confiança.

Karl Ove Hufthammer
fonte
1
Boa resposta (+1). Eu também mencionaria o fato de que também é possível obter CIs de inicialização paramétrica neste caso. Este tópico contém uma discussão muito interessante sobre o assunto.
precisa saber é o seguinte
@ usεr11852 Obrigado pela sugestão. Agora adicionei um exemplo usando a inicialização paramétrica.
Karl Ove Hufthammer
6

Graus de liberdade para modelos mistos são "problemáticos". Para ler mais sobre isso, você pode conferir os valores-p, lmer e tudo o que postou por Douglas Bates. Além disso , a FAQ do r-sig-mixed-models resume as razões pelas quais é incômoda:

  • Em geral, não está claro que a distribuição nula da razão computada de somas de quadrados seja realmente uma distribuição F, para qualquer escolha de graus de liberdade do denominador. Embora isso seja verdade para casos especiais que correspondem a desenhos experimentais clássicos (aninhados, plotagem dividida, bloco aleatório etc.), aparentemente não é verdade para desenhos mais complexos (desequilibrado, GLMMs, correlação temporal ou espacial etc.).
  • Para cada receita simples de graus de liberdade que foi sugerida (traço da matriz do chapéu etc.), parece haver pelo menos um contraexemplo bastante simples em que a receita falha muito.
  • Outros esquemas de aproximação df sugeridos (Satterthwaite, Kenward-Roger, etc.) aparentemente seriam bastante difíceis de implementar no lme4 / nlme,
    (...)
  • Como os autores principais do lme4 não estão convencidos da utilidade da abordagem geral de teste com referência a uma distribuição nula aproximada, e por causa da sobrecarga de qualquer outra pessoa que está cavando o código para habilitar a funcionalidade relevante (como um patch ou um complemento) -em), é improvável que esta situação mude no futuro.

O FAQ também oferece algumas alternativas

  • use MASS :: glmmPQL (usa regras nlme antigas aproximadamente equivalentes às regras 'internas-externas' do SAS) para GLMMs ou (n) lme para LMMs
  • Adivinhe o denominador df a partir de regras padrão (para projetos padrão) e aplique-os a testes t ou F
  • Execute o modelo em lme (se possível) e use o denominador df relatado lá (que segue uma regra simples 'interna-externa' que deve corresponder à resposta canônica para projetos simples / ortogonais), aplicada a testes t ou F. Para a especificação explícita das regras que o lme usa, consulte a página 91 de Pinheiro e Bates - esta página está disponível no Google Livros
  • usar SAS, Genstat (AS-REML), Stata?
  • Suponha que o denominador infinito df (ou seja, teste Z / qui-quadrado ao invés de t / F) se o número de grupos for grande (> 45?) Foram colocadas várias regras práticas para o tamanho "aproximadamente infinito", incluindo [em Angrist e Pischke's '' Econometria principalmente inofensiva ''], 42 (em homenagem a Douglas Adams)

Mas se você estiver interessado em intervalos de confiança, existem abordagens melhores, por exemplo, baseadas no bootstrap, conforme sugerido por Karl Ove Hufthammer em sua resposta, ou as propostas no FAQ.

Tim
fonte
"Adivinhe o denominador df a partir das regras padrão (para projetos padrão) e aplique-os aos testes t ou F"; Eu realmente gostaria que alguém pudesse elaborar isso. Por exemplo, para um projeto aninhado (por exemplo, pacientes x controles, várias amostras por sujeito; com o ID do sujeito sendo o efeito aleatório), como obtemos os graus de liberdade para esse projeto?
Arnaud A