14

Suponha que eu tenha $N$ participantes, cada um deles dando uma resposta $Y$ 20 vezes, 10 em uma condição e 10 em outra. Ajustei um modelo linear de efeitos mistos comparando em cada condição. Aqui está um exemplo reproduzível simulando essa situação usando o pacote em : $Y$ lme4R

library(lme4)
fml <- "~ condition + (condition | participant_id)"
d <- expand.grid(participant_id=1:40, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))

set.seed(23432)
d <- cbind(d, simulate(formula(fml), 
                       newparams=list(beta=c(0, .5), 
                                      theta=c(.5, 0, 0), 
                                      sigma=1), 
                       family=gaussian, 
                       newdata=d))

m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)

O modelo mproduz dois efeitos fixos (uma interceptação e inclinação para condição) e três efeitos aleatórios (uma interceptação aleatória por participante, uma inclinação aleatória por participante para condição e uma correlação intercepto).

Gostaria de comparar estatisticamente o tamanho da variação de interceptação aleatória por participante nos grupos definidos por condition(ou seja, computar o componente de variação destacado em vermelho separadamente nas condições de controle e experimentais, em seguida, testar se a diferença no tamanho dos componentes é diferente de zero). Como eu faria isso (de preferência em R)?

BÔNUS

Digamos que o modelo seja um pouco mais complicado: cada participante experimenta 10 estímulos 20 vezes cada, 10 em uma condição e 10 em outra. Assim, existem dois conjuntos de efeitos aleatórios cruzados: efeitos aleatórios para participante e efeitos aleatórios para estímulo. Aqui está um exemplo reproduzível:

library(lme4)
fml <- "~ condition + (condition | participant_id) + (condition | stimulus_id)"
d <- expand.grid(participant_id=1:40, stimulus_id=1:10, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))

set.seed(23432)
d <- cbind(d, simulate(formula(fml), 
                       newparams=list(beta=c(0, .5), 
                                      theta=c(.5, 0, 0, .5, 0, 0), 
                                      sigma=1), 
                       family=gaussian, 
                       newdata=d))

m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)

Eu gostaria de comparar estatisticamente a magnitude da variação aleatória de interceptação por participante entre os grupos definidos por condition. Como eu faria isso, e o processo é diferente do da situação descrita acima?

EDITAR

Para ser um pouco mais específico sobre o que estou procurando, quero saber:

É a questão ", são as respostas médias condicionais dentro de cada condição (ou seja, valores de interceptação aleatória em cada condição) substancialmente diferentes umas das outras, além do que seria de esperar devido ao erro de amostragem" uma questão bem definida (ou seja, essa é a questão teoricamente responsável)? Se não, por que não?
Se a resposta à pergunta (1) for afirmativa, como eu a responderia? Eu preferiria uma Rimplementação, mas não estou vinculado ao lme4pacote - por exemplo, parece que o OpenMxpacote tem capacidade para acomodar análises de vários grupos e de vários níveis ( https: //openmx.ssri.psu. edu / openmx-features ), e esse parece ser o tipo de pergunta que deve ser respondida em uma estrutura SEM.

r mixed-model lme4-nlme Patrick S. Forscher
fonte

1

@ MarkWhite, atualizei a pergunta em resposta aos seus comentários. Quero dizer que quero comparar o desvio padrão das interceptações dos participantes quando eles dão respostas na condição de controle versus quando eles dão respostas na condição experimental. Eu quero fazer isso estatisticamente, ou seja, testar se a diferença no desvio padrão dos intercepta é diferente de 0.

Patrick S. Forscher

2

Eu escrevi uma resposta, mas vou dormir nela porque não tenho certeza de que seja muito útil. A questão se resume a que eu não acho que alguém possa fazer o que você está perguntando. O efeito aleatório da interceptação é a variação na média dos participantes quando eles estão na condição de controle. Portanto, não se pode olhar para a variação daqueles para observações na condição experimental. As interceptações são definidas no nível da pessoa e a condição está no nível de observação. Se você está tentando comparar variações entre condições, eu pensaria em modelos condicionalmente heterocedásticos.

Mark

2

Estou trabalhando na revisão e reenvio de um artigo em que tenho participantes que dão respostas a conjuntos de estímulos. Cada participante é exposto a várias condições e cada estímulo recebe uma resposta em várias condições - em outras palavras, meu estudo simula a configuração que eu descrevo na descrição "BÔNUS". Em um dos meus gráficos, parece que a resposta média dos participantes tem maior variabilidade em uma das condições do que nas outras. Um revisor me pediu para testar se isso é verdade.

Patrick S. Forscher

2

Consulte aqui stats.stackexchange.com/questions/322213 para saber como configurar um modelo lme4 com diferentes parâmetros de variação para cada nível de uma variável de agrupamento. Não sei ao certo como fazer um teste de hipótese sobre se dois parâmetros de variação são iguais; pessoalmente, eu sempre preferiria iniciar bootstrap sobre assuntos e estímulos para obter um intervalo de confiança, ou talvez configurar algum tipo de teste de hipótese tipo permutação (baseado em reamostragem).

Ameba diz Reinstate Monica

3

Eu concordo com o comentário de @MarkWhite de que a questão "são as variações aleatórias de interceptação substancialmente diferentes umas das outras ..." é, na melhor das hipóteses, obscura e, na pior das hipóteses, sem sentido, porque a interceptação se refere necessariamente a valores Y em um grupo específico (o grupo atribuído o valor de 0), portanto, comparar "intercepta" entre grupos estritamente falando não faz sentido. Eu acho que uma maneira melhor de reformular sua pergunta, como eu a entendo, seria algo como: "as variações das respostas médias condicionais dos participantes na condição A vs. condição B são desiguais?"

Jake Westfall

6

Há mais de uma maneira de testar essa hipótese. Por exemplo, o procedimento descrito por @amoeba deve funcionar. Mas parece-me que a maneira mais simples e conveniente de testá-lo é usando um bom e antigo teste de razão de verossimilhança comparando dois modelos aninhados. A única parte potencialmente complicada dessa abordagem é saber como configurar o par de modelos para que o abandono de um único parâmetro teste claramente a hipótese desejada de variações desiguais. Abaixo, explico como fazer isso.

Resposta curta

Alterne para a codificação de contraste (soma a zero) para sua variável independente e faça um teste de razão de probabilidade comparando seu modelo completo com um modelo que força a correlação entre inclinações aleatórias e interceptações aleatórias a 0:

# switch to numeric (not factor) contrast codes
d$contrast <- 2*(d$condition == 'experimental') - 1

# reduced model without correlation parameter
mod1 <- lmer(sim_1 ~ contrast + (contrast || participant_id), data=d)

# full model with correlation parameter
mod2 <- lmer(sim_1 ~ contrast + (contrast | participant_id), data=d)

# likelihood ratio test
anova(mod1, mod2)

Explicação visual / intuição

Para que essa resposta faça sentido, você precisa ter uma compreensão intuitiva do que valores diferentes do parâmetro de correlação implicam para os dados observados. Considere as linhas de regressão específicas do sujeito (variando aleatoriamente). Basicamente, o parâmetro de correlação controla se as linhas de regressão do participante "se expandem para a direita" (correlação positiva) ou "se expandem para a esquerda" (correlação negativa) em relação ao ponto $X=0$ , onde X é seu código independente de contraste variável. Qualquer um desses implica variação desigual nas respostas médias condicionais dos participantes. Isso é ilustrado abaixo:

Neste gráfico, ignoramos as múltiplas observações que temos para cada sujeito em cada condição e, em vez disso, apenas plotamos as duas médias aleatórias de cada sujeito, com uma linha conectando-as, representando a inclinação aleatória desse sujeito. (São dados de 10 assuntos hipotéticos, não os dados publicados no OP.)

Na coluna da esquerda, onde há uma forte correlação negativa de interceptação de inclinação, as linhas de regressão se espalham para a esquerda em relação ao ponto $X=0$ . Como você pode ver claramente na figura, isso leva a uma maior variação nas médias aleatórias dos sujeitos na condição $X=-1$ que na condição $X=1$ .

A coluna à direita mostra a imagem invertida desse padrão. Nesse caso, há uma maior variação nas médias aleatórias dos sujeitos na condição $X=1$ que na condição $X=-1$ .

A coluna no meio mostra o que acontece quando as inclinações e interceptações aleatórias não são correlacionadas. Isso significa que as linhas de regressão se espalham para a esquerda exatamente como se espalham para a direita, em relação ao ponto $X=0$ . Isso implica que as variações das médias dos sujeitos nas duas condições são iguais.

É crucial aqui que usamos um esquema de codificação de contraste de soma para zero, e não códigos fictícios (ou seja, não definimos os grupos em $X=0$ vs. $X=1$ ). É somente sob o esquema de codificação de contraste que temos essa relação em que as variações são iguais se e somente se a correlação intercepto-inclinação é 0. A figura abaixo tenta construir essa intuição:

O que esta figura mostra é o mesmo conjunto de dados exato nas duas colunas, mas com a variável independente codificada de duas maneiras diferentes. Na coluna da esquerda, usamos códigos de contraste - essa é exatamente a situação da primeira figura. Na coluna à direita, usamos códigos fictícios. Isso altera o significado das interceptações - agora as interceptações representam as respostas previstas dos sujeitos no grupo de controle. O painel inferior mostra a conseqüência dessa mudança, a saber, que a correlação entre intercepto e inclinação não está mais perto de 0, mesmo que os dados sejam os mesmos em um sentido profundo e as variações condicionais sejam iguais nos dois casos. Se isso ainda não parece fazer muito sentido, estudar essa resposta anterior, em que falo mais sobre esse fenômeno, pode ajudar.

Prova

Seja $y_{ijk}$ a $j$ ésima resposta do $i$ sujeito na condição $k$ . (Temos apenas duas condições aqui, então $k$ é apenas 1 ou 2.) Então o modelo misto pode ser escrito

y_{i j k} = α_{i} + β_{i} x_{k} + e_{i j k},

$y_{ijk} = \alpha_i + \beta_ix_k + e_{ijk},$ onde

são os assuntos interceptações aleatórias e tem variância

,

são a inclinação aleatória dos sujeitos e têm variação

,

α_{i}

$\alpha_i$

σ_{α}^{2}

$\sigma^2_\alpha$

β_{i}

$\beta_i$

σ_{β}^{2}

$\sigma^2_\beta$

e_{i j k}

$e_{ijk}$ é o termo de erro no nível de observação e

cov (α_{i}, β_{i}) = σ_{α β}

$\text{cov}(\alpha_i, \beta_i)=\sigma_{\alpha\beta}$ .

Queremos mostrar que

var (α_{i} + β_{i} x_{1}) = var (α_{i} + β_{i} x_{2}) \Leftrightarrow σ_{α β} = 0.

$\text{var}(\alpha_i + \beta_ix_1) = \text{var}(\alpha_i + \beta_ix_2) \Leftrightarrow \sigma_{\alpha\beta}=0.$

Começando com o lado esquerdo desta implicação, temos

\begin{aligned} var (α_{i} + β_{i} x_{1}) & = var (α_{i} + β_{i} x_{2}) \\ σ_{α}^{2} + x_{1}^{2} σ_{β}^{2} + 2 x_{1} σ_{α β} & = σ_{α}^{2} + x_{2}^{2} σ_{β}^{2} + 2 x_{2} σ_{α β} \\ σ_{β}^{2} (x_{1}^{2} - x_{2}^{2}) + 2 σ_{α β} (x_{1} - x_{2}) & = 0. \end{aligned}

$\begin{aligned} \text{var}(\alpha_i + \beta_ix_1) &= \text{var}(\alpha_i + \beta_ix_2) \\ \sigma^2_\alpha + x^2_1\sigma^2_\beta + 2x_1\sigma_{\alpha\beta} &= \sigma^2_\alpha + x^2_2\sigma^2_\beta + 2x_2\sigma_{\alpha\beta} \\ \sigma^2_\beta(x_1^2 - x_2^2) + 2\sigma_{\alpha\beta}(x_1 - x_2) &= 0. \end{aligned}$

Códigos de contraste soma a zero implica que $x_1 + x_2 = 0$ e $x_1^2 = x_2^2 = x^2$ . Em seguida, podemos reduzir ainda mais a última linha acima para

\begin{aligned} σ_{β}^{2} (x^{2} - x^{2}) + 2 σ_{α β} (x_{1} + x_{1}) & = 0 \\ σ_{α β} & = 0, \end{aligned}

$\begin{aligned} \sigma^2_\beta(x^2 - x^2) + 2\sigma_{\alpha\beta}(x_1 + x_1) &= 0 \\ \sigma_{\alpha\beta} &= 0, \end{aligned}$ que é o que queríamos provar. (Para estabelecer a outra direção da implicação, podemos apenas seguir os mesmos passos ao contrário.)

Para reiterar, isso mostra que, se a variável independente é codificada por contraste (soma a zero) , as variações das médias aleatórias dos sujeitos em cada condição são iguais se e somente se a correlação entre inclinações aleatórias e interceptações aleatórias for 0. A chave O ponto de partida de tudo isso é que testar a hipótese nula de que $\sigma_{\alpha\beta} = 0$ testará a hipótese nula de variâncias iguais descritas pelo OP.

Isso NÃO funciona se a variável independente for, digamos, codificada como fictícia. Especificamente, se ligue os valores $x_1=0$ e $x_2=1$ nas equações acima, descobrimos que

var (α_{i}) = var (α_{i} + β_{i}) \Leftrightarrow σ_{α β} = - \frac{σ_{β}^{2}}{2} .

$\text{var}(\alpha_i) = \text{var}(\alpha_i + \beta_i) \Leftrightarrow \sigma_{\alpha\beta} = -\frac{\sigma^2_\beta}{2}.$

Jake Westfall
fonte

Esta já é uma resposta fantástica, obrigado! Acho que isso está mais próximo de responder à minha pergunta, então estou aceitando e dando a você a recompensa (está prestes a expirar), mas eu adoraria ver uma justificativa algébrica se você tiver tempo e energia para isso.

Patrick S. Forscher 19/08/19

1

@ PatrickS.Forscher Acabei de adicionar uma prova

Jake Westfall

1

@JakeWestfall No meu exemplo de brinquedo, os sujeitos lançaram respostas nas duas condições. Se um sujeito tem resposta

em condição A e

em estado B, então o que seria o valor BLUP de intercepção aleatória para este assunto quando usamos modelo? Eu acho que só pode ser 0. Se todos os sujeitos tiverem BLUPs iguais a zero, a variação da interceptação aleatória também será zero. Portanto, este modelo é incapaz de se encaixar nesse exemplo de brinquedo. Por outro lado, o modelo definido acima via terá dois BLUPs para cada sujeito, e eles podem ser facilmente

e

. Estou faltando alguma coisa aqui?

a

$a$

- a

$-a$ (1 | subject)dummy

a

$a$

- a

$-a$

Ameba diz Reinstate Monica

1

Vejo agora que você está certo @amoeba, obrigado por explicar. Vou editar minha resposta de acordo.

Jake Westfall

1

@amoeba Você está certo que é possível que os BLUPs possam ser correlacionados mesmo sem um parâmetro de correlação no modelo. Mas acredito que, para fins de teste, o procedimento ainda funcione conforme o planejado (por exemplo, possui a taxa de erro nominal do tipo 1) porque apenas o modelo com o parâmetro de correlação é capaz de incorporá-lo à função de probabilidade e, assim, "receber crédito" por isso. . Ou seja, mesmo que os BLUPs sejam correlacionados no modelo mais simples, ainda é como se os efeitos não fossem correlacionados no que diz respeito à probabilidade total, para que o teste de LR funcione. Eu acho que :)

Jake Westfall

6

Você pode testar a significância dos parâmetros do modelo com a ajuda de intervalos de confiança estimados para os quais o pacote lme4 tem a confint.merModfunção.

inicialização (consulte, por exemplo, Intervalo de confiança na inicialização )

> confint(m, method="boot", nsim=500, oldNames= FALSE)
Computing bootstrap confidence intervals ...
                                                           2.5 %     97.5 %
sd_(Intercept)|participant_id                         0.32764600 0.64763277
cor_conditionexperimental.(Intercept)|participant_id -1.00000000 1.00000000
sd_conditionexperimental|participant_id               0.02249989 0.46871800
sigma                                                 0.97933979 1.08314696
(Intercept)                                          -0.29669088 0.06169473
conditionexperimental                                 0.26539992 0.60940435

perfil de probabilidade (ver, por exemplo, qual é a relação entre a probabilidade do perfil e os intervalos de confiança? )

> confint(m, method="profile", oldNames= FALSE)
Computing profile confidence intervals ...
                                                          2.5 %     97.5 %
sd_(Intercept)|participant_id                         0.3490878 0.66714551
cor_conditionexperimental.(Intercept)|participant_id -1.0000000 1.00000000
sd_conditionexperimental|participant_id               0.0000000 0.49076950
sigma                                                 0.9759407 1.08217870
(Intercept)                                          -0.2999380 0.07194055
conditionexperimental                                 0.2707319 0.60727448

Há também um método, 'Wald'mas isso é aplicado apenas a efeitos fixos.
Também existe algum tipo de expressão anova (razão de verossimilhança) no pacote lmerTestque é chamado ranova. Mas eu não consigo entender isso. A distribuição das diferenças no logLikelihood, quando a hipótese nula (variação zero para o efeito aleatório) for verdadeira é não distribuída por qui-quadrado (possivelmente quando o número de participantes e tentativas for alto, o teste da razão de verossimilhança poderá fazer sentido).

Variação em grupos específicos

Para obter resultados de variação em grupos específicos, você pode remameter

# different model with alternative parameterization (and also correlation taken out) 
fml1 <- "~ condition + (0 + control + experimental || participant_id) "

Onde adicionamos duas colunas ao quadro de dados ^{(isso é necessário apenas se você deseja avaliar 'controle' e 'experimental' não correlacionados, a função (0 + condition || participant_id)não levaria à avaliação dos diferentes fatores na condição de não correlacionados)}

#adding extra columns for control and experimental
d <- cbind(d,as.numeric(d$condition=='control'))
d <- cbind(d,1-as.numeric(d$condition=='control'))
names(d)[c(4,5)] <- c("control","experimental")

Agora lmer dará variação para os diferentes grupos

> m <- lmer(paste("sim_1 ", fml1), data=d)
> m
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: paste("sim_1 ", fml1)
   Data: d
REML criterion at convergence: 2408.186
Random effects:
 Groups           Name         Std.Dev.
 participant_id   control      0.4963  
 participant_id.1 experimental 0.4554  
 Residual                      1.0268  
Number of obs: 800, groups:  participant_id, 40
Fixed Effects:
          (Intercept)  conditionexperimental  
               -0.114                  0.439

E você pode aplicar os métodos de perfil a eles. Por exemplo, agora o confint fornece intervalos de confiança para o controle e a variação do esforço.

> confint(m, method="profile", oldNames= FALSE)
Computing profile confidence intervals ...
                                    2.5 %     97.5 %
sd_control|participant_id       0.3490873 0.66714568
sd_experimental|participant_id  0.3106425 0.61975534
sigma                           0.9759407 1.08217872
(Intercept)                    -0.2999382 0.07194076
conditionexperimental           0.1865125 0.69149396

Simplicidade

Você pode usar a função de probabilidade para obter comparações mais avançadas, mas existem muitas maneiras de fazer aproximações ao longo da estrada (por exemplo, você pode fazer um teste conservador de anova / lrt, mas é isso que você deseja?).

Nesse ponto, fico me perguntando qual é realmente o ponto dessa comparação (não tão comum) entre as variações. Eu me pergunto se isso começa a se tornar muito sofisticado. Por que a diferença entre as variações em vez da proporção entre as variações (que se refere à distribuição F clássica)? Por que não apenas relatar intervalos de confiança? Precisamos dar um passo atrás e esclarecer os dados e a história que devemos contar, antes de seguirmos caminhos avançados que podem ser supérfluos e com pouco contato com o assunto estatístico e as considerações estatísticas que são realmente o tópico principal.

Eu me pergunto se alguém deveria fazer muito mais do que simplesmente declarar os intervalos de confiança (o que pode realmente dizer muito mais do que um teste de hipóteses. Um teste de hipóteses dá uma resposta sim não, mas não há informações sobre a expansão real da população. Dados suficientes que você possa faça qualquer pequena diferença para ser relatada como uma diferença significativa). Para aprofundar o assunto (para qualquer finalidade), requer, acredito, uma questão de pesquisa mais específica (definida de maneira restrita), a fim de guiar o mecanismo matemático para fazer as devidas simplificações (mesmo quando um cálculo exato for viável ou quando pode ser aproximado por simulações / bootstrapping, mesmo que em algumas configurações ainda exija alguma interpretação apropriada). Compare com o teste exato de Fisher para resolver exatamente uma pergunta (específica) (sobre tabelas de contingência),

Exemplo simples

Para fornecer um exemplo da simplicidade possível, mostro abaixo uma comparação (por simulações) com uma avaliação simples da diferença entre as duas variações de grupo com base em um teste F realizado pela comparação de variações nas respostas médias individuais e realizado comparando as variações derivadas do modelo misto.

$j$

{\hat{Y}}_{Eu, j} \sim N (μ_{j}, σ_{j}^{2} + \frac{σ_{ϵ}^{2}}{10})

$\hat{Y}_{i,j} \sim N(\mu_j, \sigma_j^2 + \frac{\sigma_{\epsilon}^2}{10})$

$\sigma_\epsilon$ $\sigma_{j}$ $j = \lbrace 1,2 \rbrace$

Você pode ver isso na simulação do gráfico abaixo, onde, além da pontuação F com base na amostra, uma pontuação F é calculada com base nas variações previstas (ou somas de erro ao quadrado) do modelo.

^{$\sigma_{j=1} = \sigma_{j=2} = 0.5$ $\sigma_\epsilon=1$}

Você pode ver que há alguma diferença. Essa diferença pode dever-se ao fato de que o modelo linear de efeitos mistos está obtendo as somas de erro ao quadrado (para o efeito aleatório) de uma maneira diferente. E esses termos de erro ao quadrado não são (mais) bem expressos como uma distribuição qui-quadrado simples, mas ainda estão intimamente relacionados e podem ser aproximados.

$\sigma_{j=1} \neq \sigma_{j=2}$ $\hat{Y}_{i,j}$ $\sigma_j$ $\sigma_\epsilon$

^{$\sigma_{j=1} = 0.5$ $\sigma_{j=2} = 0.25$ $\sigma_\epsilon=1$}

Portanto, o modelo baseado nos meios é muito exato. Mas é menos poderoso. Isso mostra que a estratégia correta depende do que você deseja / precisa.

^{No exemplo acima, quando você define os limites da cauda direita em 2,1 e 3,1, obtém aproximadamente 1% da população no caso de variância igual (respectivamente 103 e 104 dos 10.000 casos), mas no caso de variância desigual, esses limites diferem muito (dando 5334 e 6716 dos casos)}

código:

set.seed(23432)

# different model with alternative parameterization (and also correlation taken out)
fml1 <- "~ condition + (0 + control + experimental || participant_id) "
fml <- "~ condition + (condition | participant_id)"

n <- 10000

theta_m <- matrix(rep(0,n*2),n)
theta_f <- matrix(rep(0,n*2),n)

# initial data frame later changed into d by adding a sixth sim_1 column
ds <- expand.grid(participant_id=1:40, trial_num=1:10)
ds <- rbind(cbind(ds, condition="control"), cbind(ds, condition="experimental"))
  #adding extra columns for control and experimental
  ds <- cbind(ds,as.numeric(ds$condition=='control'))
  ds <- cbind(ds,1-as.numeric(ds$condition=='control'))
  names(ds)[c(4,5)] <- c("control","experimental")

# defining variances for the population of individual means
stdevs <- c(0.5,0.5) # c(control,experimental)

pb <- txtProgressBar(title = "progress bar", min = 0,
                    max = n, style=3)
for (i in 1:n) {

  indv_means <- c(rep(0,40)+rnorm(40,0,stdevs[1]),rep(0.5,40)+rnorm(40,0,stdevs[2]))
  fill <- indv_means[d[,1]+d[,5]*40]+rnorm(80*10,0,sqrt(1)) #using a different way to make the data because the simulate is not creating independent data in the two groups 
  #fill <- suppressMessages(simulate(formula(fml), 
  #                     newparams=list(beta=c(0, .5), 
  #                                    theta=c(.5, 0, 0), 
  #                                    sigma=1), 
  #                     family=gaussian, 
  #                     newdata=ds))
  d <- cbind(ds, fill)
  names(d)[6] <- c("sim_1")


  m <- lmer(paste("sim_1 ", fml1), data=d)
  m
  theta_m[i,] <- m@theta^2

  imeans <- aggregate(d[, 6], list(d[,c(1)],d[,c(3)]), mean)
  theta_f[i,1] <- var(imeans[c(1:40),3])
  theta_f[i,2] <- var(imeans[c(41:80),3])

  setTxtProgressBar(pb, i)
}
close(pb)

p1 <- hist(theta_f[,1]/theta_f[,2], breaks = seq(0,6,0.06))       
fr <- theta_m[,1]/theta_m[,2]
fr <- fr[which(fr<30)]
p2 <- hist(fr, breaks = seq(0,30,0.06))



plot(-100,-100, xlim=c(0,6), ylim=c(0,800), 
     xlab="F-score", ylab = "counts [n out of 10 000]")
plot( p1, col=rgb(0,0,1,1/4), xlim=c(0,6), ylim=c(0,800), add=T)  # means based F-score
plot( p2, col=rgb(1,0,0,1/4), xlim=c(0,6), ylim=c(0,800), add=T)  # model based F-score
fr <- seq(0, 4, 0.01)
lines(fr,df(fr,39,39)*n*0.06,col=1)
legend(2, 800, c("means based F-score","mixed regression based F-score"), 
       fill=c(rgb(0,0,1,1/4),rgb(1,0,0,1/4)),box.col =NA, bg = NA)
legend(2, 760, c("F(39,39) distribution"), 
       lty=c(1),box.col = NA,bg = NA)
title(expression(paste(sigma[1]==0.5, " , ", sigma[2]==0.5, " and ", sigma[epsilon]==1)))

Sextus Empiricus
fonte

Isso é útil, mas parece não abordar a questão de como comparar variações em duas condições.

Ameba diz Reinstate Monica

@amoeba Descobri que essa resposta fornece o núcleo do problema (sobre como testar os componentes de variação aleatória). O que o OP deseja com precisão é difícil de ler no texto inteiro. A que "as variações de interceptação aleatória" se referem? (o plural em relação à interceptação me confunde) Um caso possível pode ser o uso do modelo. sim_1 ~ condition + (0 + condition | participant_id)"Nesse caso, você obtém uma parametrização em dois parâmetros (um para cada grupo) em vez de dois parâmetros, um para a interceptação e outro para o efeito (que precisam ser combinados para os grupos).

Sextus Empiricus

Cada indivíduo tem alguma resposta média na condição A e alguns resposta média na condição B. A questão é se a variação entre os indivíduos em A é diferente da variância entre os indivíduos em B.

ameba diz Reintegrar Monica

Isso não conclui a tarefa apresentada no título "Comparar o componente de variação aleatória entre os níveis de uma variável de agrupamento". Percebi que havia um erro de digitação confuso no corpo da pergunta, que corrigi. Também tentei esclarecer melhor a redação da pergunta.

Patrick S. Forscher 13/08/19

Pode ser possível responder à pergunta usando car::linearHypothesisTest( math.furman.edu/~dcs/courses/math47/R/library/car/html/… ), que permite ao usuário testar hipóteses arbitrárias com um modelo ajustado. No entanto, eu teria que usar o método de @ amoeba para obter as duas interceptações aleatórias no mesmo modelo ajustado para que elas possam ser comparadas com esta função. Também estou um pouco incerto quanto à validade do método.

Patrick S. Forscher

5

Uma maneira relativamente direta pode ser o uso de testes de razão de verossimilhança, anovaconforme descrito nas lme4Perguntas frequentes .

Começamos com um modelo completo no qual as variações são irrestritas (ou seja, duas variações diferentes são permitidas) e, em seguida, ajustamos um modelo restrito no qual as duas variações são consideradas iguais. Nós simplesmente compará-los com anova()(note que eu definir REML = FALSEembora REML = TRUEcom anova(..., refit = FALSE)é completamente viável ).

m_full <- lmer(sim_1 ~ condition + (condition | participant_id), data=d, REML = FALSE)
summary(m_full)$varcor
 # Groups         Name                  Std.Dev. Corr  
 # participant_id (Intercept)           0.48741        
 #                conditionexperimental 0.26468  -0.419
 # Residual                             1.02677     

m_red <- lmer(sim_1 ~ condition + (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

anova(m_full, m_red)
# Data: d
# Models:
# m_red: sim_1 ~ condition + (1 | participant_id)
# m_full: sim_1 ~ condition + (condition | participant_id)
#        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
# m_red   4 2396.6 2415.3 -1194.3   2388.6                         
# m_full  6 2398.7 2426.8 -1193.3   2386.7 1.9037      2      0.386

No entanto, este teste é provavelmente conservador . Por exemplo, o FAQ diz:

Lembre-se de que os testes de hipótese nula baseados em LRT são conservadores quando o valor nulo (como σ2 = 0) está no limite do espaço viável; no caso mais simples (variância de efeito aleatório único), o valor de p é aproximadamente duas vezes maior do que deveria (Pinheiro e Bates 2000).

Existem várias alternativas:

$\chi^2$
Simule a distribuição correta usando RLRsim(como também descrito nas Perguntas frequentes).

Vou demonstrar a segunda opção da seguinte maneira:

library("RLRsim")
## reparametrize model so we can get one parameter that we want to be zero:
afex::set_sum_contrasts() ## warning, changes contrasts globally
d <- cbind(d, difference = model.matrix(~condition, d)[,"condition1"])

m_full2 <- lmer(sim_1 ~ condition + (difference | participant_id), data=d, REML = FALSE)
all.equal(deviance(m_full), deviance(m_full2))  ## both full models are identical

## however, we need the full model without correlation!
m_full2b <- lmer(sim_1 ~ condition + (1| participant_id) + 
                   (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_full2b)$varcor
 # Groups           Name        Std.Dev.
 # participant_id   (Intercept) 0.44837 
 # participant_id.1 difference  0.13234 
 # Residual                     1.02677 

## model that only has random effect to be tested
m_red <- update(m_full2b,  . ~ . - (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name       Std.Dev.
 # participant_id difference 0.083262
 # Residual                  1.125116

## Null model 
m_null <- update(m_full2b,  . ~ . - (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_null)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

exactRLRT(m_red, m_full2b, m_null)
# Using restricted likelihood evaluated at ML estimators.
# Refit with method="REML" for exact results.
# 
#   simulated finite sample distribution of RLRT.
#   
#   (p-value based on 10000 simulated values)
# 
# data:  
# RLRT = 1.9698, p-value = 0.0719

Como podemos ver, o resultado sugere que REML = TRUE, se tivéssemos obtido resultados exatos. Mas isso é deixado como um exercício para o leitor.

Em relação ao bônus, não tenho certeza se RLRsimpermite testes simultâneos de vários componentes, mas, se houver, isso pode ser feito da mesma maneira.

Resposta ao comentário:

$\theta_X$ $\theta_0$ $X$

Não tenho certeza de que esta pergunta possa receber uma resposta razoável.

Uma interceptação aleatória permite uma diferença idiossincrática no nível geral para cada nível do fator de agrupamento. Por exemplo, se a variável dependente é o tempo de resposta, alguns participantes são mais rápidos e outros mais lentos.
Uma inclinação aleatória permite a cada nível do fator de agrupamento um efeito idiossincrático do fator para o qual as inclinações aleatórias são estimadas. Por exemplo, se o fator for congruência, alguns participantes podem ter um efeito de congruência mais alto que outros.

Então, inclinações aleatórias afetam a interceptação aleatória? Em certo sentido, isso pode fazer sentido, pois permite a cada nível do fator de agrupamento um efeito completamente idiossincrático para cada condição. No final, estimamos dois parâmetros idiossincráticos para duas condições. No entanto, acho que a distinção entre o nível geral capturado pela interceptação e o efeito específico da condição capturado pela inclinação aleatória é importante e, em seguida, a inclinação aleatória não pode realmente afetar a interceptação aleatória. No entanto, ele ainda permite que cada nível do fator de agrupamento seja idiossincrático separadamente para cada nível da condição.

No entanto, meu teste ainda faz o que a pergunta original quer. Ele testa se a diferença nas variações entre as duas condições é zero. Se for zero, as variações nas duas condições são iguais. Em outras palavras, somente se não houver necessidade de uma inclinação aleatória é que a variação nas duas condições é idêntica. Espero que faça sentido.

Henrik
fonte

1

Você utiliza contrastes de tratamento ( contr.treatment) para os quais a condição de controle é a referência (ou seja, para a qual a interceptação aleatória é calculada). A parametrização que proponho usar contrastes de soma (ie contr.sum) e a interceptação é a grande média. Sinto que faz mais sentido testar se a diferença é nula quando a interceptação é a grande média em vez da condição de controle (mas escrever isso sugere que pode ser relativamente inconseqüente). Você pode ler as páginas 24 a 26 de: singmann.org/download/publications/…

Henrik

1

Obrigado! Minhas perguntas são um pouco diferentes, no entanto: (1) Sua resposta parece implicar que minha pergunta se reduz a "é a inclinação aleatória da condição diferente de 0". Isso é verdade? (2) Se a resposta a (1) for "sim", isso sugere outra interpretação da inclinação aleatória para condition: permite que a interceptação aleatória varie entre os níveis de condition. Isso é verdade?

Patrick S. Forscher

2

Meu contraexemplo de 2 ¢: @amoeba ao procedimento proposto por Henrik está correto. Henrik está quase correto, mas ele compara o par errado de modelos. A comparação de modelos que a pergunta de Patrick de resposta é a comparação entre os modelos Henrik chamados m_fullvs. m_full2b. Ou seja: as variações das respostas médias condicionais dos participantes em A vs. B são desiguais se a correlação aleatória interceptar-inclinação for diferente de zero - importante, sob a parametrização da codificação de contraste soma-zero . Testar a variação aleatória da inclinação não é necessário. Tentando pensar como explicar isso de forma sucinta ...

Jake Westfall

2

Esta não é realmente uma explicação adequada, mas estudar minha resposta aqui pode lançar um pouco de luz sobre o assunto. Basicamente, o parâmetro de correlação controla se as linhas de regressão do participante "se espalham para a direita" (corr. Positiva) ou "se espalham para a esquerda" (corr. Negativa). Qualquer um destes implica variação desigual nas respostas médias condicionais dos participantes. Soma-to-zero codificação em seguida, garante que nós estamos olhando para a correlação no ponto certo em X

Jake Westfall

2

Vou considerar postar uma resposta com fotos, se eu encontrar tempo ...

Jake Westfall

5

Seu modelo

m = lmer(sim_1 ~ condition + (condition | participant_id), data=d)

já permite que a variação entre sujeitos na condição de controle seja diferente da variação entre sujeitos na condição experimental. Isso pode ser explicitado por uma parametrização equivalente:

m = lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=d)

A matriz de covariância aleatória agora tem uma interpretação mais simples:

Random effects:
 Groups         Name                  Variance Std.Dev. Corr
 participant_id conditioncontrol      0.2464   0.4963       
                conditionexperimental 0.2074   0.4554   0.83

Aqui, as duas variações são precisamente as duas variações nas quais você está interessado: a variação [entre sujeitos] das respostas médias condicionais na condição de controle e a mesma na condição experimental. No seu conjunto de dados simulado, eles são 0,25 e 0,21. A diferença é dada por

delta = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]

e é igual a 0,039. Você deseja testar se é significativamente diferente de zero.

EDIT: percebi que o teste de permutação que descrevi abaixo está incorreto; não funcionará como pretendido se os meios na condição controle / experimental não forem os mesmos (porque as observações não podem ser trocadas sob o valor nulo). Pode ser uma idéia melhor inicializar assuntos (ou assuntos / itens no caso Bonus) e obter o intervalo de confiança para delta.

Vou tentar corrigir o código abaixo para fazer isso.

Sugestão original baseada em permutação (incorreta)

Costumo achar que alguém pode se salvar de muitos problemas fazendo um teste de permutação. De fato, neste caso, é muito fácil de configurar. Permitamos controle / condições experimentais para cada sujeito separadamente; então qualquer diferença nas variações deve ser eliminada. Repetir isso muitas vezes produzirá a distribuição nula para as diferenças.

(Eu não programa em R; todo mundo sinta-se à vontade para reescrever o seguinte em um estilo R melhor.)

set.seed(42)
nrep = 100
v = matrix(nrow=nrep, ncol=1)
for (i in 1:nrep)
{
   dp = d
   for (s in unique(d$participant_id)){             
     if (rbinom(1,1,.5)==1){
       dp[p$participant_id==s & d$condition=='control',]$condition = 'experimental'
       dp[p$participant_id==s & d$condition=='experimental',]$condition = 'control'
     }
   }
  m <- lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=dp)
  v[i,] = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]
}
pvalue = sum(abs(v) >= abs(delta)) / nrep

Executar isso gera o valor p $p=0.7$ . Pode-se aumentar nreppara 1000 ou mais.

Exatamente a mesma lógica pode ser aplicada no seu caso de bônus.

ameba diz restabelecer Monica
fonte

Super interessante, obrigado! Vou ter que pensar mais sobre por que sua reparameterização funciona, pois essa parece ser a principal ideia dessa resposta.

Patrick S. Forscher

Estranhamente, os valores de interceptação por grupo em sua resposta parecem diferir daqueles da resposta de @MartijnWeterings.

Patrick S. Forscher

@ PatrickS.Forscher Isso porque, eu acho, gera um conjunto de dados diferente. Posso usar a sim_1 ~ 0 + condition + (0 + dummy(condition, "control") + dummy(condition, "experimental") | participant_id)formulação e obter o mesmo resultado da minha resposta.

Ameba diz Reinstate Monica

1

@ PatrickS.Forscher Não, usei os dados gerados pelo seu código (com sua semente). Defino a semente para 42 somente ao realizar o teste de permutação. Foi Martijn quem mudou o conjunto de dados, não eu.

Ameba diz Reinstate Monica

1

Esta proposta é definitivamente sólida. Como eu acho que você já experimentou, a configuração de testes de permutação para dados em vários níveis não é totalmente simples. Uma abordagem semelhante que seria um pouco mais fácil de implementar seria o bootstrapping paramétrico, o que é bastante simples de fazer com o lme4, usando o método simulate () dos objetos lmer ajustados, ou seja, chame simulate (m) muitas vezes para criar o bootstrap distribuição. Apenas uma idéia para brincar.

Jake Westfall

0

Olhando para esse problema de uma perspectiva um pouco diferente e partindo da forma "geral" do modelo linear misto, temos

y_{Eu j k} = μ + α_{j} + d_{Eu j} + e_{Eu j k}, d_{Eu} \sim N (0 0, Σ), e_{Eu j k} \sim N (0 0, σ^{2})

$y_{ijk} = \mu + \alpha_j + d_{ij} + e_{ijk}, \quad d_{i} \sim N(0, \Sigma), \quad e_{ijk} \sim N(0, \sigma^2)$ Onde

α_{j}

$\alpha_j$ é o efeito fixo do

j

$j$ condição e

d_{i} = (d_{i 1}, \dots, d_{i J})^{⊤}

$d_i = (d_{i1}, \ldots, d_{iJ})^\top$ é um vetor aleatório (acho que alguns chamam de efeito aleatório com valor vetorial) para o

i

$i$ participante do

j

$j$ condição.
No seu exemplo, temos duas condições

y_{i 1 k}

$y_{i1k}$ e

y_{i 2 k}

$y_{i2k}$ que vou denotar como

A

$A$ e

B

$B$ no que se segue. Portanto, a matriz de covariância do vetor aleatório bidimensional

d_{i}

$d_{i}$ é da forma geral

Σ = [\begin{matrix} σ_{UMA}^{2} & σ_{UMA B} \\ σ_{UMA B} & σ_{B}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_A & \sigma_{AB}\\ \sigma_{AB} & \sigma^2_{B} \end{bmatrix}$

com não negativo $\sigma^2_A$ e $\sigma^2_{B}$ .

Vamos primeiro ver como a versão re-parametrizada do $\Sigma$ parece quando usamos contrastes de soma.

A variação da interceptação, que corresponde à média geral, é

σ_{1}^{2} : = Var (média geral) = Var (\frac{1}{2} \cdot (UMA + B)) = \frac{1}{4} \cdot (Var (UMA) + Var (B) + 2 \cdot Cov (UMA, B)) .

$\sigma^2_{1} := \text{Var(grand mean)} = \text{Var}(\frac{1}{2} \cdot (A+B)) = \frac{1}{4} \cdot (\text{Var}(A) + \text{Var}(B) + 2 \cdot \text{Cov}(A,B)).$

A variação do contraste é

σ_{2}^{2} : = Var (contraste) = Var (\frac{1}{2} \cdot (UMA - B)) = \frac{1}{4} \cdot (Var (UMA) + Var (B) - 2 \cdot Cov (UMA, B)) .

$\sigma^2_{2} := \text{Var(contrast)} = \text{Var}(\frac{1}{2} \cdot (A-B)) = \frac{1}{4} \cdot (\text{Var}(A) + \text{Var}(B) - 2 \cdot \text{Cov}(A,B)).$

E a covariância entre o intercepto e o contraste é

σ_{12} : = Cov (grande média, contraste) = Cov (\frac{1}{2} \cdot (UMA + B), \frac{1}{2} \cdot (UMA - B)) = \frac{1}{4} \cdot (Var (UMA) - Var (B)) .

$\sigma_{12} := \text{Cov}(\text{grand mean, contrast}) = \text{Cov}(\frac{1}{2} \cdot (A+B), \frac{1}{2} \cdot (A-B)) \\ = \frac{1}{4} \cdot (\text{Var}(A) - \text{Var}(B)).$

Assim, o re-parametrizado $\Sigma$ é

Σ = [\begin{matrix} σ_{1}^{2} + σ_{2}^{2} + 2 σ_{12} & σ_{1}^{2} - σ_{2}^{2} \\ σ_{1}^{2} - σ_{2}^{2} & σ_{1}^{2} + σ_{2}^{2} - 2 σ_{12} \end{matrix}] = [\begin{matrix} σ_{UMA}^{2} & σ_{UMA B} \\ σ_{UMA B} & σ_{B}^{2} \end{matrix}] .

$\Sigma = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 + 2\sigma_{12} & \sigma^2_1 - \sigma^2_2\\ \sigma^2_1 - \sigma^2_2 & \sigma^2_1 + \sigma^2_2 - 2\sigma_{12} \end{bmatrix} = \begin{bmatrix}\sigma^2_A & \sigma_{AB}\\ \sigma_{AB} & \sigma^2_{B} \end{bmatrix}.$

$\Sigma$ pode ser decomposto em

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} σ_{2}^{2} & - σ_{2}^{2} \\ - σ_{2}^{2} & σ_{2}^{2} \end{matrix}] + 2 [\begin{matrix} σ_{12} & 0 0 \\ 0 0 & - σ_{12} \end{matrix}] .

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & -\sigma^2_2\\ -\sigma^2_2 & \sigma^2_2\end{bmatrix} + 2\begin{bmatrix}\sigma_{12} & 0\\ 0 & -\sigma_{12}\end{bmatrix}.$

Configurando o parâmetro de covariância $\sigma_{12}$ a zero chegamos

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} σ_{2}^{2} & - σ_{2}^{2} \\ - σ_{2}^{2} & σ_{2}^{2} \end{matrix}] = [\begin{matrix} σ_{1}^{2} + σ_{2}^{2} & σ_{1}^{2} - σ_{2}^{2} \\ σ_{1}^{2} - σ_{2}^{2} & σ_{1}^{2} + σ_{2}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & -\sigma^2_2\\ -\sigma^2_2 & \sigma^2_2\end{bmatrix} = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 & \sigma^2_1 - \sigma^2_2 \\ \sigma^2_1 - \sigma^2_2 & \sigma^2_1 + \sigma^2_2\end{bmatrix}$

que, como @Jake Westfall derivou um pouco diferente, testa a hipótese de variâncias iguais quando comparamos um modelo sem esse parâmetro de covariância com um modelo em que o parâmetro de covariância ainda está incluído / não está definido como zero.

Notavelmente, a introdução de outro fator de agrupamento aleatório cruzado (como estímulos) não altera a comparação do modelo que deve ser feita, ou seja, anova(mod1, mod2)(opcionalmente com o argumento refit = FALSEao usar a estimativa REML) onde mod1e mod2são definidos como @Jake Westfall.

Tirando $\sigma_{12}$ e o componente de variação para o contraste $\sigma^2_2$ (o que @Henrik sugere) resulta em

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix}$

que testa a hipótese de que as variações nas duas condições são iguais e que são iguais à covariância (positiva) entre as duas condições.

Quando temos duas condições, um modelo que se ajusta a uma matriz de covariância com dois parâmetros em uma estrutura simétrica composta (positiva) também pode ser escrito como

# code snippet from Jake Westfall
d$contrast <- 2*(d$condition == 'experimental') - 1

# new model
mod3 <- lmer(sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id), 
             data = d, REML = FALSE)

ou (usando a variável / fator categórico condition)

mod4 <- lmer(sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id), 
             data = d, REML = FALSE)

com

Σ = [\begin{matrix} σ_{1}^{2} + σ_{2}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} + σ_{2}^{2} \end{matrix}] = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} σ_{2}^{2} & 0 0 \\ 0 0 & σ_{2}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1 + \sigma^2_2 \end{bmatrix} = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & 0\\ 0 & \sigma^2_2\end{bmatrix}$

Onde $\sigma^2_1$ e $\sigma^2_2$ são os parâmetros de variação para o participante e a interceptação participante-condição-combinação, respectivamente. Note que este $\Sigma$ possui um parâmetro de covariância não negativo.

Abaixo vemos que mod1, mod3e mod4deu acessos equivalentes:

# code snippet from Jake Westfall
d$contrast <- 2*(d$condition == 'experimental') - 1

mod1 <- lmer(sim_1 ~ contrast + (contrast || participant_id),
             data = d, REML = FALSE)

mod2 <- lmer(sim_1 ~ contrast + (contrast | participant_id),
             data = d, REML = FALSE)

# new models 
mod3 <- lmer(sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id), 
             data = d, REML = FALSE) 

mod4 <- lmer(sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id), 
             data = d, REML = FALSE)

anova(mod3, mod1)
# Data: d
# Models:
# mod3: sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id)
# mod1: sim_1 ~ contrast + ((1 | participant_id) + (0 + contrast | participant_id))
#      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
# mod3  5 2396.9 2420.3 -1193.5   2386.9                        
# mod1  5 2396.9 2420.3 -1193.5   2386.9     0      0          1

anova(mod4, mod3)
# Data: d
# Models:
# mod4: sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id)
# mod3: sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id)
#      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
# mod4  5 2396.9 2420.3 -1193.5   2386.9                        
# mod3  5 2396.9 2420.3 -1193.5   2386.9     0      0          1

Com os contrastes do tratamento (o padrão em R ), o parametrizado $\Sigma$ é

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} + σ_{12} \\ σ_{1}^{2} + σ_{12} & σ_{1}^{2} + σ_{2}^{2} + 2 σ_{12} \end{matrix}] = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} 0 0 & 0 0 \\ 0 0 & σ_{2}^{2} \end{matrix}] + [\begin{matrix} 0 0 & σ_{12} \\ σ_{12} & 2 σ_{12} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1 + \sigma_{12}\\ \sigma^2_1 + \sigma_{12} & \sigma^2_1 + \sigma^2_2 + 2\sigma_{12} \end{bmatrix} = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}0 & 0\\ 0 & \sigma^2_2\end{bmatrix} + \begin{bmatrix}0 & \sigma_{12}\\ \sigma_{12} & 2\sigma_{12}\end{bmatrix}$

Onde $\sigma^2_1$ é o parâmetro de variação para a interceptação (condição $A$ ), $\sigma^2_2$ o parâmetro de variação para o contraste ( $A - B$ ) e $\sigma_{12}$ o parâmetro de covariância correspondente.

Podemos ver que nem a configuração $\sigma_{12}$ para zero nem configuração $\sigma^2_2$ zerar (apenas) a hipótese de variâncias iguais.

No entanto, como mostrado acima, ainda podemos mod4testar a hipótese, pois a alteração dos contrastes não afeta a parametrização de $\Sigma$ para este modelo.

statmerkur
fonte

Modelo de efeitos mistos: compare o componente de variação aleatória entre os níveis de uma variável de agrupamento

Respostas:

Resposta curta

Explicação visual / intuição

Prova