Supongamos que tengo participantes, cada uno de los cuales da una respuesta 20 veces, 10 en una condición y 10 en otra. Encajo un modelo lineal de efectos mixtos que compara en cada condición. Aquí hay un ejemplo reproducible que simula esta situación usando el paquete en :S $N$$Y$$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)

El modelo mproduce dos efectos fijos (una intersección y pendiente para la condición) y tres efectos aleatorios (una intersección aleatoria por participante, una pendiente aleatoria por participante para la condición y una correlación de pendiente de intercepción).

Me gustaría comparar estadísticamente el tamaño de la varianza de interceptación aleatoria por participante entre los grupos definidos por condition(es decir, calcular el componente de varianza resaltado en rojo por separado dentro de las condiciones de control y experimentales, luego probar si la diferencia en el tamaño de los componentes es distinto de cero) ¿Cómo haría esto (preferiblemente en R)?

---

PRIMA

Digamos que el modelo es un poco más complicado: cada participante experimenta 10 estímulos 20 veces cada uno, 10 en una condición y 10 en otra. Por lo tanto, hay dos conjuntos de efectos aleatorios cruzados: efectos aleatorios para el participante y efectos aleatorios para el estímulo. Aquí hay un ejemplo reproducible:

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)

Me gustaría comparar estadísticamente la magnitud de la varianza de interceptación aleatoria por participante entre los grupos definidos por condition. ¿Cómo haría eso, y el proceso es diferente al de la situación descrita anteriormente?

---

EDITAR

Para ser un poco más específico sobre lo que estoy buscando, quiero saber:

1. Es la pregunta, "son las respuestas medias condicionales dentro de cada condición (es decir, valores de intercepción aleatoria en cada condición) sustancialmente diferentes entre sí, más allá de lo que esperaríamos debido al error de muestreo", una pregunta bien definida (es decir, es esta pregunta incluso teóricamente responsable)? ¿Si no, porque no?
2. Si la respuesta a la pregunta (1) es sí, ¿cómo la respondería? Preferiría una Rimplementación, pero no estoy atado al lme4paquete; por ejemplo, parece que el OpenMxpaquete tiene la capacidad de acomodar análisis multigrupo y multinivel ( https: //openmx.ssri.psu. edu / openmx-features ), y este parece ser el tipo de pregunta que debería responderse en un marco SEM.

Hay más de una forma de probar esta hipótesis. Por ejemplo, el procedimiento descrito por @amoeba debería funcionar. Pero me parece que la forma más simple y más conveniente de probarlo es usar una buena prueba de razón de probabilidad antigua que compara dos modelos anidados. La única parte potencialmente complicada de este enfoque es saber cómo configurar el par de modelos para que la eliminación de un solo parámetro pruebe de manera clara la hipótesis deseada de variaciones desiguales. A continuación explico cómo hacer eso.

Respuesta corta

Cambie a la codificación de contraste (suma a cero) para su variable independiente y luego haga una prueba de razón de probabilidad comparando su modelo completo con un modelo que obligue a que la correlación entre pendientes aleatorias e intercepciones aleatorias sea 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)

Explicación visual / intuición

Para que esta respuesta tenga sentido, debe tener una comprensión intuitiva de lo que implican los diferentes valores del parámetro de correlación para los datos observados. Considere las líneas de regresión específicas del sujeto (que varían aleatoriamente). Básicamente, el parámetro de correlación controla si las líneas de regresión del participante "se despliegan a la derecha" (correlación positiva) o "se despliegan a la izquierda" (correlación negativa) en relación con el punto $X=0$ , donde X es su código independiente de contraste variable. Cualquiera de estos implica una variación desigual en las respuestas medias condicionales de los participantes. Esto se ilustra a continuación:

En esta gráfica, ignoramos las múltiples observaciones que tenemos para cada sujeto en cada condición y, en su lugar, graficamos las dos medias aleatorias de cada sujeto, con una línea que las conecta, que representa la pendiente aleatoria de ese sujeto. (Esto se compone de datos de 10 sujetos hipotéticos, no de los datos publicados en el OP).

En la columna de la izquierda, donde hay una fuerte correlación negativa de pendiente-intersección, las líneas de regresión se extienden hacia la izquierda en relación con el punto $X=0$ . Como puede ver claramente en la figura, esto conduce a una mayor variación en las medias aleatorias de los sujetos en la condición $X=-1$ que en la condición $X=1$ .

La columna de la derecha muestra la imagen inversa de este patrón. En este caso, existe una mayor varianza en las medias aleatorias de los sujetos en la condición $X=1$ que en la condición $X=-1$ .

La columna en el medio muestra lo que sucede cuando las pendientes aleatorias y las intersecciones aleatorias no están correlacionadas. Esto significa que las líneas de regresión se despliegan a la izquierda exactamente tanto como a la derecha, en relación con el punto $X=0$ . Esto implica que las varianzas de las medias de los sujetos en las dos condiciones son iguales.

Aquí es crucial que hayamos utilizado un esquema de codificación de contraste de suma a cero, no códigos ficticios (es decir, no establecer los grupos en $X=0$ vs. $X=1$ ). Es solo bajo el esquema de codificación de contraste que tenemos esta relación en la que las variaciones son iguales si y solo si la correlación pendiente-intersección es 0. La siguiente figura intenta construir esa intuición:

Lo que muestra esta figura es el mismo conjunto de datos exacto en ambas columnas, pero con la variable independiente codificada de dos maneras diferentes. En la columna de la izquierda usamos códigos de contraste: esta es exactamente la situación de la primera figura. En la columna de la derecha usamos códigos ficticios. Esto altera el significado de las intercepciones: ahora las intercepciones representan las respuestas predichas de los sujetos en el grupo de control. El panel inferior muestra la consecuencia de este cambio, a saber, que la correlación pendiente-intersección ya no está cerca de 0, a pesar de que los datos son los mismos en un sentido profundo y las variaciones condicionales son iguales en ambos casos. Si esto todavía no parece tener mucho sentido, estudiar esta respuesta anterior mía donde hablo más sobre este fenómeno puede ayudar.

Prueba

Sea $y_{ijk}$ la respuesta $j$ del sujeto $i$ bajo la condición $k$ . (Aquí solo tenemos dos condiciones, entonces $k$ es solo 1 o 2.) Entonces el modelo mixto puede escribirse

$$y_{ijk} = \alpha_i + \beta_ix_k + e_{ijk},$$ donde $\alpha_i$ son los sujetos ' intercepta al azar y tiene varianza $\sigma^2_\alpha$ , $\beta_i$son la pendiente aleatoria de los sujetos y tienen una varianza $\sigma^2_\beta$ , $e_{ijk}$ es el término de error de nivel de observación, y $\text{cov}(\alpha_i, \beta_i)=\sigma_{\alpha\beta}$ .

Queremos mostrar que

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

Comenzando con el lado izquierdo de esta implicación, tenemos

$$\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}$$

Los códigos de contraste de suma a cero implican que $x_1 + x_2 = 0$ y $x_1^2 = x_2^2 = x^2$ . Entonces podemos reducir aún más la última línea de lo anterior a

$$\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 es lo que queríamos probar. (Para establecer la otra dirección de la implicación, podemos seguir estos mismos pasos a la inversa).

Para reiterar, esto muestra que si la variable independiente está codificada por contraste (suma a cero) , entonces las variaciones de las medias aleatorias de los sujetos en cada condición son iguales si y solo si la correlación entre pendientes aleatorias e interceptaciones aleatorias es 0. La clave El punto clave de todo esto es que probar la hipótesis nula de que $\sigma_{\alpha\beta} = 0$ probará la hipótesis nula de varianzas iguales descritas por el OP.

Esto NO funciona si la variable independiente es, por ejemplo, un código ficticio. Específicamente, si conectamos los valores $x_1=0$ y $x_2=1$ en las ecuaciones anteriores, encontramos que

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

Puede probar la importancia de los parámetros del modelo con la ayuda de intervalos de confianza estimados para los que el paquete lme4 tiene la confint.merModfunción.

bootstrapping (ver por ejemplo Intervalo de confianza de bootstrap )

> 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 probabilidad (ver por ejemplo ¿Cuál es la relación entre la probabilidad de perfil y los intervalos de confianza? )

> 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

---

• También hay un método, 'Wald'pero esto se aplica solo a efectos fijos.

• También existe algún tipo de expresión anova (razón de probabilidad) en el paquete lmerTestque se nombra ranova. Pero parece que esto no tiene sentido. La distribución de las diferencias en logLikelihood, cuando la hipótesis nula (varianza cero para el efecto aleatorio) es verdadera, no está distribuida por chi-cuadrado (posiblemente cuando el número de participantes y ensayos es alto, la prueba de razón de probabilidad podría tener sentido).

---

Varianza en grupos específicos

Para obtener resultados de varianza en grupos específicos, puede volver a parametrizar

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

Donde agregamos dos columnas al marco de datos ^{(esto solo es necesario si desea evaluar 'control' y 'experimental' no correlacionados, la función (0 + condition || participant_id)no conduciría a la evaluación de los diferentes factores en condición como no 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")

Ahora lmerdará variación para los 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 

Y puede aplicar los métodos de perfil a estos. Por ejemplo, ahora confint da intervalos de confianza para el control y la varianza ejercimental.

> 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

---

Sencillez

Podría usar la función de probabilidad para obtener comparaciones más avanzadas, pero hay muchas maneras de hacer aproximaciones a lo largo del camino (por ejemplo, podría hacer una prueba conservadora anova / lrt, pero ¿es eso lo que quiere?).

En este punto, me pregunto cuál es el punto de esta comparación (no tan común) entre las variaciones. Me pregunto si comienza a volverse demasiado sofisticado. ¿Por qué la diferencia entre las variaciones en lugar de la relación entre las variaciones (que se relaciona con la distribución F clásica)? ¿Por qué no solo informar intervalos de confianza? Necesitamos dar un paso atrás y aclarar los datos y la historia que se supone que debe contar, antes de entrar en vías avanzadas que pueden ser superfluas y perder el contacto con el asunto estadístico y las consideraciones estadísticas que en realidad son el tema principal.

Me pregunto si uno debería hacer mucho más que simplemente declarar los intervalos de confianza (que en realidad pueden decir mucho más que una prueba de hipótesis. Una prueba de hipótesis da un sí, no una respuesta, pero no hay información sobre la propagación real de la población. Con suficientes datos puede hacer una pequeña diferencia para ser reportado como una diferencia significativa). Para profundizar en el asunto (para cualquier propósito), se requiere, creo, una pregunta de investigación más específica (definida de manera más estricta) para guiar la maquinaria matemática para hacer las simplificaciones adecuadas (incluso cuando un cálculo exacto sea factible o cuando se puede aproximar mediante simulaciones / bootstrapping, incluso en algunos casos aún requiere una interpretación adecuada). Compare con la prueba exacta de Fisher para resolver una pregunta (particular) (sobre tablas de contingencia) exactamente,

Ejemplo simple

Para proporcionar un ejemplo de la simplicidad que es posible, muestro a continuación una comparación (mediante simulaciones) con una evaluación simple de la diferencia entre las dos variaciones grupales basada en una prueba F realizada al comparar las variaciones en las respuestas medias individuales y al comparar El modelo mixto derivó variaciones.

$j$

$$\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$

Puede ver esto en la simulación del gráfico a continuación, donde, aparte de la puntuación F basada en la muestra, se calcula una puntuación F basada en las variaciones predichas (o sumas de error al cuadrado) del modelo.

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

Puedes ver que hay alguna diferencia. Esta diferencia puede deberse al hecho de que el modelo lineal de efectos mixtos está obteniendo las sumas de error al cuadrado (para el efecto aleatorio) de una manera diferente. Y estos términos de error al cuadrado ya no se expresan (ya) como una distribución simple de Chi-cuadrado, pero aún están estrechamente relacionados y pueden 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$}

Entonces, el modelo basado en los medios es muy exacto. Pero es menos poderoso. Esto muestra que la estrategia correcta depende de lo que quiere / necesita.

^{En el ejemplo anterior, cuando establece los límites de la cola derecha en 2.1 y 3.1, obtiene aproximadamente el 1% de la población en el caso de varianza igual (resp 103 y 104 de los 10 000 casos) pero en el caso de varianza desigual, estos límites difieren mucho (dando 5334 y 6716 de los 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)))

Una forma relativamente directa podría ser utilizar pruebas de razón de probabilidad anovasegún se describe en las lme4preguntas frecuentes .

Comenzamos con un modelo completo en el que las variaciones no están restringidas (es decir, se permiten dos variaciones diferentes) y luego ajustamos un modelo restringido en el que se supone que las dos variaciones son iguales. Simplemente los comparamos con anova()(tenga en cuenta que configuré REML = FALSEaunque REML = TRUEcon anova(..., refit = FALSE)es completamente factible ).

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

Sin embargo, esta prueba es probablemente conservadora . Por ejemplo, las preguntas frecuentes dicen:

Tenga en cuenta que las pruebas de hipótesis nula basadas en LRT son conservadoras cuando el valor nulo (como σ2 = 0) está en el límite del espacio factible; en el caso más simple (variación de efecto aleatorio simple), el valor p es aproximadamente el doble de lo que debería ser (Pinheiro y Bates 2000).

Hay varias alternativas:

1. Cree una distribución de prueba adecuada, que generalmente consiste en una mezcla de $\chi^2$distribuciones Ver, por ejemplo,
   Self, SG y Liang, K.-Y. (1987) Propiedades asintóticas de estimadores de máxima verosimilitud y pruebas de razón de verosimilitud en condiciones no estándar. Revista de la Asociación Americana de Estadística, 82 (398), 605. https://doi.org/10.2307/2289471 Sin embargo, esto es bastante complicado.

2. Simule la distribución correcta usando RLRsim(como también se describe en las preguntas frecuentes).

Demostraré la segunda opción a continuación:

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, el resultado sugiere que con REML = TRUEnosotros habríamos obtenido resultados exactos. Pero esto se deja como un ejercicio para el lector.

Con respecto a la bonificación, no estoy seguro si RLRsimpermite la prueba simultánea de múltiples componentes, pero si es así, esto se puede hacer de la misma manera.

---

Respuesta al comentario:

Entonces, es cierto que, en general, la pendiente aleatoria $\theta_X$ permite la intercepción aleatoria $\theta_0$ para variar a través de los niveles de $X$?

No estoy seguro de que esta pregunta pueda recibir una respuesta razonable.

• Una intercepción aleatoria permite una diferencia idiosincrásica en el nivel general para cada nivel del factor de agrupación. Por ejemplo, si la variable dependiente es el tiempo de respuesta, algunos participantes son más rápidos y otros más lentos.
• Una pendiente aleatoria permite a cada nivel del factor de agrupación un efecto idiosincrásico del factor para el cual se estiman las pendientes aleatorias. Por ejemplo, si el factor es la congruencia, algunos participantes pueden tener un efecto de congruencia más alto que otros.

Entonces, ¿las pendientes aleatorias afectan la intersección aleatoria? En cierto sentido, esto podría tener sentido, ya que permiten que cada nivel del factor de agrupación tenga un efecto completamente idiosincrásico para cada condición. Al final, estimamos dos parámetros idiosincrásicos para dos condiciones. Sin embargo, creo que la distinción entre el nivel general capturado por la intersección y el efecto específico de la condición capturado por la pendiente aleatoria es importante y, a continuación, la pendiente aleatoria no puede afectar realmente la intersección aleatoria. Sin embargo, todavía permite que cada nivel del factor de agrupación sea idiosincrásico por separado para cada nivel de la condición.

Sin embargo, mi prueba sigue haciendo lo que quiere la pregunta original. Comprueba si la diferencia en las variaciones entre las dos condiciones es cero. Si es cero, entonces las variaciones en ambas condiciones son iguales. En otras palabras, solo si no hay necesidad de una pendiente aleatoria, la varianza en ambas condiciones es idéntica. Espero que tenga sentido.

Su modelo

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

ya permite que la varianza entre sujetos en la condición de control difiera de la varianza entre sujetos en la condición experimental. Esto puede hacerse más explícito mediante una re-parametrización equivalente:

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

La matriz de covarianza aleatoria ahora tiene una interpretación más simple:

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

Aquí las dos variaciones son precisamente las dos variaciones que le interesan: la variación [entre sujetos] de las respuestas medias condicionales en la condición de control y lo mismo en la condición experimental. En su conjunto de datos simulado, son 0.25 y 0.21. La diferencia viene dada por

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

y es igual a 0.039. Desea probar si es significativamente diferente de cero.

EDITAR: me di cuenta de que la prueba de permutación que describo a continuación es incorrecta; no funcionará según lo previsto si los medios en control / condición experimental no son los mismos (porque entonces las observaciones no son intercambiables bajo nulo). Puede ser una mejor idea arrancar temas (o temas / elementos en el caso de Bonus) y obtener el intervalo de confianza para delta.

Intentaré arreglar el siguiente código para hacerlo.

---

Sugerencia original basada en permutación (incorrecta)

A menudo encuentro que uno puede ahorrarse muchos problemas haciendo una prueba de permutación. De hecho, en este caso es muy fácil de configurar. Permutamos las condiciones de control / experimentales para cada sujeto por separado; entonces cualquier diferencia en las variaciones debe ser eliminada. Repetir esto muchas veces dará como resultado una distribución nula para las diferencias.

(No programo en R; todos pueden reescribir lo siguiente en un mejor estilo R).

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

Ejecutar esto produce el valor p $p=0.7$. Uno puede aumentarnrep a 1000 más o menos.

Se puede aplicar exactamente la misma lógica en su caso de Bonus.

Mirando este problema desde una perspectiva ligeramente diferente y comenzando desde la forma "general" del modelo lineal mixto, tenemos

$$y_{ijk} = \mu + \alpha_j + d_{ij} + e_{ijk}, \quad d_{i} \sim N(0, \Sigma), \quad e_{ijk} \sim N(0, \sigma^2)$$ dónde $\alpha_j$ es el efecto fijo de la $j$condición y $d_i = (d_{i1}, \ldots, d_{iJ})^\top$ es un vector aleatorio (creo que algunos lo llaman efecto aleatorio con valor vectorial) para el $i$'th participante en el $j$Condición.
En tu ejemplo tenemos dos condiciones$y_{i1k}$ y $y_{i2k}$ que denotaré como $A$ y $B$en lo que sigue. Entonces, la matriz de covarianza del vector aleatorio bidimensional$d_{i}$ es de la forma general

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

con no negativo $\sigma^2_A$ y $\sigma^2_{B}$.

Primero veamos cómo la versión re-parametrizada de $\Sigma$ se ve cuando usamos contrastes de suma.

La varianza de la intersección, que corresponde a la gran media, es

$$\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)).$$

La varianza del contraste es

$$\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)).$$

Y la covarianza entre la intersección y el contraste es

$$\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)).$$

Thus, the re-parameterized $\Sigma$ is

$$\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$ can be decomposed into

$$\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}.$$

Setting the covariance parameter $\sigma_{12}$ to zero we get

$$\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 dedujo de manera ligeramente diferente, prueba la hipótesis de varianzas iguales cuando comparamos un modelo sin este parámetro de covarianza con un modelo en el que el parámetro de covarianza todavía se incluye / no se establece en cero.

En particular, la introducción de otro factor de agrupación aleatorio cruzado (como los estímulos) no cambia la comparación del modelo que debe hacerse, es decir, anova(mod1, mod2)(opcionalmente con el argumento refit = FALSEcuando usa la estimación REML) donde mod1y mod2se definen como @Jake Westfall.

Sacando $\sigma_{12}$ y el componente de varianza para el contraste $\sigma^2_2$ (lo que sugiere @Henrik) da como resultado

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

que prueba la hipótesis de que las varianzas en las dos condiciones son iguales y que son iguales a la covarianza (positiva) entre las dos condiciones.

---

Cuando tenemos dos condiciones, un modelo que se ajusta a una matriz de covarianza con dos parámetros en una estructura simétrica compuesta (positiva) también se puede escribir 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) 

o (usando la variable / factor categórico condition)

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

con

$$\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}$$

where $\sigma^2_1$ and $\sigma^2_2$ are the variance parameters for the participant and the participant-condition-combination intercepts, respectively. Note that this $\Sigma$ has a non-negative covariance parameter.

Below we see that mod1, mod3, and mod4 yield equivalent fits:

# 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

---

With treatment contrasts (the default in R) the re-parameterized $\Sigma$ is

$$\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}$$

where $\sigma^2_1$ is the variance parameter for the intercept (condition $A$), $\sigma^2_2$ the variance parameter for the contrast ($A - B$), and $\sigma_{12}$ the corresponding covariance parameter.

We can see that neither setting $\sigma_{12}$ to zero nor setting $\sigma^2_2$ to zero tests (only) the hypothesis of equal variances.

However, as shown above, we can still use mod4 to test the hypothesis as changing the contrasts has no impact on the parameterization of $\Sigma$ for this model.