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?

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.

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.

e

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

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.

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

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