Em um modelo multinível, quais são as implicações práticas de estimar versus não estimar parâmetros de correlação de efeito aleatório?

27

Em um modelo multinível, quais são as implicações práticas e relacionadas à interpretação de estimar versus não estimar parâmetros de correlação de efeito aleatório? A razão prática para perguntar isso é que, na estrutura anterior em R, não há método implementado para estimar valores de p por meio de técnicas MCMC quando são feitas estimativas no modelo de correlações entre parâmetros.

Por exemplo, olhando para este exemplo (partes citadas abaixo), quais são as implicações práticas de M2 ​​versus M3. Obviamente, em um caso, P5 não será estimado e, no outro, será.

Questões

  1. Por razões práticas (o desejo de obter um valor-p através das técnicas MCMC), pode-se querer ajustar um modelo sem correlações entre efeitos aleatórios, mesmo que P5 seja substancialmente diferente de zero. Se alguém fizer isso e depois estimar valores de p através da técnica MCMC, os resultados são interpretáveis? (Eu sei que o @Ben Bolker mencionou anteriormente que "combinar o teste de significância com o MCMC é um pouco incoerente, estatisticamente, embora eu compreenda o desejo de fazê-lo (obter intervalos de confiança é mais suportável)" , por isso, se você for dormir melhor à noite finja que disse intervalos de confiança.)
  2. Se alguém não estimar P5, é o mesmo que afirmar que é 0?
  3. Se P5 realmente for diferente de zero, de que maneira os valores estimados de P1-P4 são afetados?
  4. Se P5 realmente for diferente de zero, de que maneira as estimativas de erro para P1-P4 são afetadas?
  5. Se P5 realmente for diferente de zero, então de que maneira as interpretações de um modelo estão deixando de incluir P5 com falhas?

Tomando emprestado a resposta de @Mike Lawrence (os que têm mais conhecimento do que eu são livres para substituí-lo pela notação completa do modelo, não estou totalmente confiante de que posso fazê-lo com razoável fidelidade):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2)(Estimativas P1 - P4)

M3: V1 ~ (1+V3|V2) + V3(estimativas P1-P5)

Parâmetros que podem ser estimados:

P1 : Uma interceptação global

P2 : Efeito aleatório intercepta para V2 (ou seja, para cada nível de V2, o desvio desse intercepto em relação à interceptação global)

P3 : Uma estimativa global única para o efeito (declive) de V3

P4 : O efeito de V3 em cada nível de V2 (mais especificamente, o grau em que o efeito de V3 em um determinado nível se desvia do efeito global de V3), enquanto impõe uma correlação zero entre os desvios de interceptação e os desvios de efeito V3 nos níveis de V2.

P5 : A correlação entre os desvios de interceptação e os desvios de V3 nos níveis de V2

Respostas derivadas de uma simulação suficientemente grande e ampla, juntamente com o código que acompanha R em Lmer, seriam aceitáveis.

russellpierce
fonte
@JackTanner: Também não parece que você tenha satisfação lá. Seria ótimo se suas preocupações também fossem abordadas na resposta a esta pergunta.
22413 russellpierce
4
Dar uma resposta exata a muitas das suas perguntas - "o que acontece com _______ quando não especifico o modelo de maneira _______" - é provavelmente impossível sem aprofundar a teoria, possivelmente intratável (embora esse possa ser um caso especial em que algo é possível - I não tenho certeza). A estratégia que eu provavelmente usaria é simular dados quando a inclinação e a interceptação estiverem altamente correlacionadas, ajustar o modelo que restringe os dois a não serem correlacionados e comparar os resultados com quando o modelo for especificado corretamente (ou seja, "análise de sensibilidade").
Macro
4
Para suas perguntas, tenho 80 (mas não 100)% de certeza do seguinte: re. # 2, sim, se você não estimar a correlação, força-a a 0; de resto, se a correlação não for exatamente 0, você está especificando incorretamente a não independência de seus dados. No entanto, os betas podem ser imparciais, mas os valores de p serão desativados (e se eles são muito altos ou muito baixos dependem e podem não ser conhecidos). Assim, interpretações dos betas podem ser capazes de prosseguir normalmente, mas interpretações de "significados" serão imprecisas.
gung - Restabelece Monica
2
@ Macro: Minha esperança era que uma recompensa pudesse liberar uma boa resposta baseada na teoria e não na simulação. Com uma simulação, frequentemente ficarei preocupado em não pegar um caso de borda apropriado. Sou ótimo em executar simulações, mas sempre me sinto um pouco ... incerto de que estou executando todas as simulações corretas (embora suponha que eu possa deixar isso para os editores de periódicos decidirem). Talvez eu precise fazer outra pergunta sobre quais cenários incluir.
26413 russellpierce

Respostas:

16

Considere os dados do estudo do sono, incluídos no lme4. Bates discute isso em seu livro on-line sobre lme4. No capítulo 3, ele considera dois modelos para os dados.

M0 0:Reação1+Dias+(1|Sujeito)+(0 0+Dias|Sujeito)

e

MUMA:Reação1+Dias+(Dias|Sujeito)

O estudo envolveu 18 indivíduos, estudados por um período de 10 dias privados de sono. Os tempos de reação foram calculados na linha de base e nos dias subsequentes. Há um efeito claro entre o tempo de reação e a duração da privação do sono. Também existem diferenças significativas entre os sujeitos. O modelo A permite a possibilidade de uma interação entre os efeitos de interceptação aleatória e inclinação: imagine, digamos, que pessoas com tempos de reação ruins sofram mais agudamente dos efeitos da privação do sono. Isso implicaria uma correlação positiva nos efeitos aleatórios.

No exemplo de Bates, não houve correlação aparente do gráfico Lattice e nenhuma diferença significativa entre os modelos. No entanto, para investigar a questão colocada acima, decidi pegar os valores ajustados do estudo do sono, aumentar a correlação e analisar o desempenho dos dois modelos.

Como você pode ver na imagem, longos tempos de reação estão associados a uma maior perda de desempenho. A correlação utilizada para a simulação foi de 0,58

insira a descrição da imagem aqui

Simulei 1000 amostras, usando o método simule no lme4, com base nos valores ajustados dos meus dados artificiais. Encaixei M0 e Ma em cada uma e observei os resultados. O conjunto de dados original teve 180 observações (10 para cada um dos 18 indivíduos) e os dados simulados têm a mesma estrutura.

A linha inferior é que há muito pouca diferença.

  1. Os parâmetros fixos têm exatamente os mesmos valores nos dois modelos.
  2. Os efeitos aleatórios são ligeiramente diferentes. Existem 18 efeitos aleatórios de interceptação e 18 de inclinação para cada amostra simulada. Para cada amostra, esses efeitos são forçados a adicionar a 0, o que significa que a diferença média entre os dois modelos é (artificialmente) 0. Mas as variações e covariâncias diferem. A covariância mediana em MA foi de 104, contra 84 em M0 (valor real, 112). As variações de declives e interceptações foram maiores sob M0 que MA, presumivelmente para obter espaço extra de manobra necessário na ausência de um parâmetro de covariância livre.
  3. O método ANOVA para lmer fornece uma estatística F para comparar o modelo Slope com um modelo com apenas uma interceptação aleatória (nenhum efeito devido à privação do sono). Claramente, esse valor era muito grande nos dois modelos, mas era tipicamente (mas nem sempre) maior no MA (média 62 vs média de 55).
  4. A covariância e variância dos efeitos fixos são diferentes.
  5. Cerca da metade do tempo, sabe que o MA está correto. O valor p mediano para comparar M0 a MA é 0,0442. Apesar da presença de uma correlação significativa e 180 observações equilibradas, o modelo correto seria escolhido apenas cerca da metade do tempo.
  6. Os valores previstos diferem nos dois modelos, mas muito ligeiramente. A diferença média entre as previsões é 0, com dp de 2,7. O sd dos valores previstos é 60,9

Então porque isso acontece? @gung imaginou, razoavelmente, que a falha em incluir a possibilidade de uma correlação força os efeitos aleatórios a não serem correlacionados. Talvez deva; mas nessa implementação, os efeitos aleatórios podem ser correlacionados, o que significa que os dados são capazes de puxar os parâmetros na direção certa, independentemente do modelo. O erro do modelo errado aparece na probabilidade, razão pela qual você pode (às vezes) distinguir os dois modelos nesse nível. O modelo de efeitos mistos está basicamente ajustando regressões lineares para cada sujeito, influenciado pelo que o modelo pensa que deveria ser. O modelo errado força o ajuste de valores menos plausíveis do que o modelo correto. Mas os parâmetros, no final do dia, são governados pelo ajuste aos dados reais.

insira a descrição da imagem aqui

Aqui está o meu código um pouco desajeitado. A idéia era ajustar os dados do estudo do sono e, em seguida, construir um conjunto de dados simulados com os mesmos parâmetros, mas uma correlação maior para os efeitos aleatórios. Esse conjunto de dados foi alimentado para simular.lmer () para simular 1000 amostras, cada uma das quais se encaixava nos dois sentidos. Depois de emparelhar os objetos ajustados, eu poderia extrair diferentes recursos do ajuste e compará-los, usando testes t ou qualquer outra coisa.

    # Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2) 
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01) 
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study. 
yhat <- XX.f %*%  fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10) 
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
    sim.data <- vector(mode="list",)
    tempReaction <- simulate(fm.sim, nsim=n.sim)
    tempdata <- model.frame(fm.sim)
    for (i in 1:n.sim){
        tempdata$Reaction <- tempReaction[,i]
			output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
			output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
			temp <- anova(output0,output1)
			pval <- temp$`Pr(>Chisq)`[2]
        sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
    }
Placidia
fonte
1
Esse é um trabalho interessante. Obrigado. Quero ver que outros comentários serão apresentados nos próximos dias e como as coisas se generalizam para outros casos antes de aceitar a resposta. Você também consideraria incluir o código R relevante em sua resposta, além de especificar a versão domer usada? Seria interessante inserir os mesmos casos simulados no PROC MIXED para ver como ele lida com a correlação de efeitos aleatórios não especificada.
russellpierce
1
@rpierce Adicionei o exemplo de código conforme solicitado. Eu o escrevi originalmente no LaTeX / Sweave, então as linhas de código foram entrelaçadas com meus comentários para mim mesmo. Eu usei a versão 1.1-6 do lme4, que é a versão atual em junho de 2014.
Placidia
@ Ben quando você diz "O modelo A permite" no segundo parágrafo, não deveria ser MO?
Nzcoops
Eu acho que o texto está correto (tudo o que fiz para esta questão era embelezar a fórmula um pouco)
Ben Bolker
+6. Excelente resposta, obrigado por sua atenção para perguntas antigas, mas valiosas.
Ameba diz Reinstate Monica
4

O Placidia já forneceu uma resposta completa usando dados simulados com base no sleepstudyconjunto de dados. Aqui está outra resposta (menos rigorosa) que também usa os sleepstudydados.

Vemos que é possível afetar a correlação estimada entre a interceptação aleatória e a inclinação aleatória "deslocando" a variável preditora aleatória. Veja os resultados dos modelos fm1e fm2abaixo:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

A partir da saída do modelo, vemos que a correlação de variação aleatória mudou. No entanto, as inclinações (fixas e aleatórias) permaneceram as mesmas, assim como a estimativa da variância residual. As estimativas de interceptações (fixas e aleatórias) foram alteradas em resposta à variável deslocada.

A covariância aleatória de intercepto-inclinação para LMMs é discutida nas notas de aula do Dr. Jack Weiss aqui . Weiss observa que reduzir a correlação de variância dessa maneira às vezes pode ajudar na convergência do modelo, entre outras coisas.

O exemplo acima varia a correlação aleatória (parâmetro "P5"). Parcialmente abordando o Q3 do OP, vemos pela saída acima que:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected
kakarot
fonte
Obrigado por adicionar sinal a esta pergunta de longa data!
russellpierce
Nota: todas as excelentes palestras de Jack Weiss e exercícios de classe / notas estão ligados de este post
theforestecologist
Como devemos então interpretar os dados em questão? Qual é a correlação "verdadeira"? Aquele do primeiro ou do segundo modelo? Ou os dos BLUPs?
User33268