Supposons que j'ai $N$ participants, chacun donnant une réponse $Y$ 20 fois, 10 dans une condition et 10 dans une autre. J'adapte un modèle linéaire à effets mixtes comparant $Y$ dans chaque condition. Voici un exemple reproductible simulant cette situation en utilisant le lme4package dans R:

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)

Le modèle mproduit deux effets fixes (une interception et une pente pour la condition) et trois effets aléatoires (une interception aléatoire par participant, une pente aléatoire par participant pour la condition et une corrélation intercept-pente).

Je voudrais comparer statistiquement la taille de la variance d'interception aléatoire par participant parmi les groupes définis par condition(c.-à-d. Calculer la composante de variance surlignée en rouge séparément dans les conditions de contrôle et expérimentales, puis tester si la différence de taille des composantes est différent de zéro). Comment pourrais-je faire cela (de préférence en R)?

---

PRIME

Disons que le modèle est légèrement plus compliqué: les participants éprouvent chacun 10 stimuli 20 fois chacun, 10 dans une condition et 10 dans une autre. Ainsi, il existe deux ensembles d'effets aléatoires croisés: les effets aléatoires pour le participant et les effets aléatoires pour le stimulus. Voici un exemple reproductible:

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)

Je voudrais comparer statistiquement l'ampleur de la variance d'interception aléatoire par participant à travers les groupes définis par condition. Comment pourrais-je faire cela, et le processus est-il différent de celui dans la situation décrite ci-dessus?

---

ÉDITER

Pour être un peu plus précis sur ce que je recherche, je veux savoir:

1. La question «est-ce que les réponses moyennes conditionnelles à l'intérieur de chaque condition (c.-à-d. Les valeurs d'interception aléatoire dans chaque condition) sont sensiblement différentes les unes des autres, au-delà de ce à quoi nous nous attendions en raison d'une erreur d'échantillonnage»? Une question bien définie (c.-à-d., Cette question même théoriquement responsable)? Sinon, pourquoi pas?
2. Si la réponse à la question (1) est oui, comment pourrais-je y répondre? Je préférerais une Rimplémentation, mais je ne suis pas lié au lme4package - par exemple, il semble que le OpenMxpackage ait la capacité d'accueillir des analyses multi-groupes et multi-niveaux ( https: //openmx.ssri.psu. edu / openmx-features ), et cela semble être le genre de question à laquelle il faudrait répondre dans un cadre SEM.

Il existe plusieurs façons de tester cette hypothèse. Par exemple, la procédure décrite par @amoeba devrait fonctionner. Mais il me semble que le moyen le plus simple et le plus rapide de le tester est d'utiliser un bon vieux test de rapport de vraisemblance comparant deux modèles imbriqués. La seule partie potentiellement délicate de cette approche consiste à savoir comment configurer la paire de modèles afin que l'abandon d'un seul paramètre teste proprement l'hypothèse souhaitée de variances inégales. Ci-dessous, j'explique comment faire cela.

Réponse courte

Passez au codage de contraste (somme à zéro) pour votre variable indépendante, puis effectuez un test de rapport de vraisemblance comparant votre modèle complet à un modèle qui force la corrélation entre les pentes aléatoires et les intersections aléatoires à 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)

Explication visuelle / intuition

Pour que cette réponse ait un sens, vous devez avoir une compréhension intuitive de ce que les différentes valeurs du paramètre de corrélation impliquent pour les données observées. Considérez les lignes de régression spécifiques au sujet (variant au hasard). Fondamentalement, le paramètre de corrélation contrôle si les lignes de régression des participants "s'étalent vers la droite" (corrélation positive) ou "s'étalent vers la gauche" (corrélation négative) par rapport au point $X=0$ , où X est votre indépendant codé par contraste. variable. L'un ou l'autre de ces facteurs implique une variance inégale dans les réponses moyennes conditionnelles des participants. Ceci est illustré ci-dessous:

Dans ce graphique, nous ignorons les multiples observations que nous avons pour chaque sujet dans chaque condition et au lieu de cela, nous traçons simplement les deux moyennes aléatoires de chaque sujet, avec une ligne les reliant, représentant la pente aléatoire de ce sujet. (Il s'agit de données provenant de 10 sujets hypothétiques, et non des données publiées dans le PO.)

Dans la colonne de gauche, où il existe une forte corrélation négative d'interception de pente, les lignes de régression s'étendent vers la gauche par rapport au point $X=0$ . Comme vous pouvez le voir clairement sur la figure, cela conduit à une plus grande variance des moyennes aléatoires des sujets dans la condition $X=-1$ que dans la condition $X=1$ .

La colonne de droite montre l'image miroir inversée de ce motif. Dans ce cas, il y a une plus grande variance dans les moyennes aléatoires des sujets dans la condition $X=1$ que dans la condition $X=-1$ .

La colonne du milieu montre ce qui se passe lorsque les pentes aléatoires et les intersections aléatoires ne sont pas corrélées. Cela signifie que les lignes de régression se déploient vers la gauche exactement autant qu'elles se déploient vers la droite, par rapport au point $X=0$ . Cela implique que les variances des moyennes des sujets dans les deux conditions sont égales.

Il est crucial ici que nous ayons utilisé un schéma de codage de contraste de somme à zéro, pas des codes fictifs (c'est-à-dire ne définissant pas les groupes à $X=0$ contre $X=1$ ). C'est seulement sous le schéma de codage de contraste que nous avons cette relation dans laquelle les variances sont égales si et seulement si la corrélation pente-ordonnée à l'origine est 0. La figure ci-dessous essaie de construire cette intuition:

Ce que cette figure montre est le même ensemble de données exact dans les deux colonnes, mais avec la variable indépendante codée de deux manières différentes. Dans la colonne de gauche, nous utilisons des codes de contraste - c'est exactement la situation de la première figure. Dans la colonne de droite, nous utilisons des codes fictifs. Cela modifie le sens des interceptions - maintenant les interceptions représentent les réponses prédites des sujets dans le groupe témoin. Le panneau du bas montre la conséquence de ce changement, à savoir que la corrélation interception de pente n'est plus nulle part proche de 0, même si les données sont les mêmes dans un sens profond et les variances conditionnelles sont égales dans les deux cas. Si cela ne semble toujours pas très logique, étudier ma réponse précédente où je parle davantage de ce phénomène peut aider.

Preuve

Soit $y_{ijk}$ la $j$ ème réponse du $i$ ème sujet sous la condition $k$ . (Nous n'avons ici que deux conditions, donc $k$ est juste soit 1 soit 2.) Alors le modèle mixte peut être écrit

$$y_{ijk} = \alpha_i + \beta_ix_k + e_{ijk},$$ où $\alpha_i$ sont les sujets ' interceptions aléatoires et ont la variance $\sigma^2_\alpha$ , $\beta_i$sont la pente aléatoire des sujets et ont une variance $\sigma^2_\beta$ , $e_{ijk}$ est le terme d'erreur au niveau de l'observation, et $\text{cov}(\alpha_i, \beta_i)=\sigma_{\alpha\beta}$ .

Nous souhaitons montrer que

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

En commençant par le côté gauche de cette implication, nous avons

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

Les codes de contraste de somme à zéro impliquent que $x_1 + x_2 = 0$ et $x_1^2 = x_2^2 = x^2$ . Ensuite, nous pouvons encore réduire la dernière ligne de ce qui précède à

$$\begin{aligned} \sigma^2_\beta(x^2 - x^2) + 2\sigma_{\alpha\beta}(x_1 + x_1) &= 0 \\ \sigma_{\alpha\beta} &= 0, \end{aligned}$$ c'est ce que nous voulions prouver. (Pour établir l'autre sens de l'implication, nous pouvons simplement suivre ces mêmes étapes à l'envers.)

Pour le répéter, cela montre que si la variable indépendante est codée en contraste (somme à zéro) , alors les variances des moyennes aléatoires des sujets dans chaque condition sont égales si et seulement si la corrélation entre les pentes aléatoires et les intersections aléatoires est de 0. La clé Le point à retenir de tout cela est que tester l'hypothèse nulle selon laquelle $\sigma_{\alpha\beta} = 0$ testera l'hypothèse nulle des variances égales décrites par l'OP.

Cela ne fonctionne PAS si la variable indépendante est, disons, codée fictivement. Plus précisément, si nous introduisons les valeurs $x_1=0$ et $x_2=1$ dans les équations ci-dessus, nous constatons que

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

Vous pouvez tester la signification des paramètres du modèle à l'aide d' intervalles de confiance estimés pour lesquels le package lme4 a leconfint.merMod fonction.

bootstrapping (voir par exemple Intervalle de confiance du 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 

profil de vraisemblance (voir par exemple Quelle est la relation entre la probabilité de profil et les intervalles de confiance? )

> 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

---

• Il y a aussi une méthode 'Wald' mais elle n'est appliquée qu'aux effets fixes.

• Il existe également une sorte d'expression de type anova (rapport de vraisemblance) dans le package lmerTestnommé ranova. Mais je n'arrive pas à donner un sens à cela. La distribution des différences de logLik vraisemblance, lorsque l'hypothèse nulle (variance nulle pour l'effet aléatoire) est vraie, n'est pas distribuée en chi carré (peut-être lorsque le nombre de participants et d'essais est élevé, le test du rapport de vraisemblance pourrait avoir un sens).

---

Variance dans des groupes spécifiques

Pour obtenir des résultats de variance dans des groupes spécifiques, vous pouvez reparamétrer

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

Lorsque nous avons ajouté deux colonnes à la base de données ^{(cela n'est nécessaire que si vous souhaitez évaluer un `` contrôle '' et un `` expérimental '' non corrélés, la fonction (0 + condition || participant_id)ne conduirait pas à l'évaluation des différents facteurs en condition de non corrélés)}

#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")

Donnera maintenant la lmervariance pour les différents groupes

> 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 

Et vous pouvez leur appliquer les méthodes de profil. Par exemple maintenant confint donne des intervalles de confiance pour le contrôle et la variance exerimentale.

> 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

---

Simplicité

Vous pouvez utiliser la fonction de vraisemblance pour obtenir des comparaisons plus avancées, mais il existe de nombreuses façons de faire des approximations le long de la route (par exemple, vous pouvez faire un test anova / lrt conservateur, mais est-ce ce que vous voulez?).

À ce stade, cela me fait me demander quel est réellement le point de cette comparaison (pas si courante) entre les variances. Je me demande si ça commence à devenir trop sophistiqué. Pourquoi la différence entre les variances au lieu du ratio entre les variances (qui se rapporte à la distribution F classique)? Pourquoi ne pas simplement signaler les intervalles de confiance? Nous devons prendre du recul et clarifier les données et l'histoire qu'elles sont censées raconter, avant de s'engager dans des voies avancées qui peuvent être superflues et en contact lâche avec la question statistique et les considérations statistiques qui sont en fait le sujet principal.

Je me demande si l'on devrait faire bien plus que simplement énoncer les intervalles de confiance (ce qui peut en fait en dire bien plus qu'un test d'hypothèse. Un test d'hypothèse donne oui non, mais aucune information sur la répartition réelle de la population. Avec suffisamment de données, vous pouvez faire une légère différence à signaler comme différence significative). Pour approfondir la question (à quelque fin que ce soit), je crois qu'une question de recherche plus spécifique (étroitement définie) doit guider la machine mathématique à effectuer les simplifications appropriées (même lorsqu'un calcul exact pourrait être faisable ou lorsque il pourrait être approximé par des simulations / bootstrapping, même dans certains cas, cela nécessite encore une interprétation appropriée). Comparez avec le test exact de Fisher pour résoudre une question (particulière) (sur les tableaux de contingence) exactement,

Exemple simple

Pour donner un exemple de la simplicité possible, je montre ci-dessous une comparaison (par simulations) avec une évaluation simple de la différence entre les deux variances de groupe basée sur un test F fait en comparant les variances dans les réponses moyennes individuelles et fait en comparant les variances dérivées du modèle mixte.

$j$

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

si la variance de l'erreur de mesure $\sigma_\epsilon$$\sigma_{j}$$j = \lbrace 1,2 \rbrace$

Vous pouvez le voir dans la simulation du graphique ci-dessous où, à part le score F basé sur l'échantillon, un score F est calculé en fonction des variances prévues (ou sommes d'erreur quadratique) du modèle.

^{L'image est modélisée avec 10 000 répétitions en utilisant $\sigma_{j=1} = \sigma_{j=2} = 0.5$$\sigma_\epsilon=1$}

Vous pouvez voir qu'il y a une différence. Cette différence peut être due au fait que le modèle linéaire à effets mixtes obtient les sommes d'erreur quadratique (pour l'effet aléatoire) d'une manière différente. Et ces termes d'erreur au carré ne sont pas (plus) bien exprimés sous la forme d'une simple distribution du Chi au carré, mais restent étroitement liés et peuvent être approximés.

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

Le modèle basé sur les moyennes est donc très exact. Mais c'est moins puissant. Cela montre que la bonne stratégie dépend de ce que vous voulez / avez besoin.

^{Dans l'exemple ci-dessus, lorsque vous définissez les limites de la queue droite à 2,1 et 3,1, vous obtenez environ 1% de la population en cas de variance égale (respectivement 103 et 104 des 10000 cas), mais en cas de variance inégale, ces limites diffèrent beaucoup (donnant 5334 et 6716 des cas)}

code:

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

Un moyen relativement simple pourrait être d'utiliser des tests de rapport de vraisemblance via anovacomme décrit dans la lme4FAQ .

Nous commençons avec un modèle complet dans lequel les variances ne sont pas contraintes (c'est-à-dire que deux variances différentes sont autorisées), puis nous ajustons un modèle contraint dans lequel les deux variances sont supposées égales. Nous les comparons simplement avec anova()(notez que je les fixe REML = FALSEbien REML = TRUEqu'avec anova(..., refit = FALSE)soit complètement réalisable ).

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

Cependant, ce test est probablement conservateur . Par exemple, la FAQ dit:

Gardez à l'esprit que les tests d'hypothèse nulle basés sur LRT sont conservateurs lorsque la valeur nulle (telle que σ2 = 0) se trouve à la limite de l'espace réalisable; dans le cas le plus simple (variance à effet aléatoire unique), la valeur de p est environ deux fois plus grande qu'elle devrait l'être (Pinheiro et Bates 2000).

Il existe plusieurs alternatives:

1. Créez une distribution de test appropriée, qui consiste généralement en un mélange de $\chi^2$distributions. Voir, par exemple,
   Self, SG et Liang, K.-Y. (1987). Propriétés asymptotiques des estimateurs du maximum de vraisemblance et tests de rapport de vraisemblance dans des conditions non standard. Journal de l'American Statistical Association, 82 (398), 605. https://doi.org/10.2307/2289471 Cependant, c'est assez compliqué.

2. Simulez la distribution correcte en utilisant RLRsim(comme également décrit dans la FAQ).

Je vais démontrer la deuxième option dans ce qui suit:

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

Comme nous pouvons le voir, la sortie suggère qu'avec REML = TRUE nous, nous aurions obtenu des résultats exacts. Mais cela est laissé au lecteur comme exercice.

Concernant le bonus, je ne sais pas si RLRsim permet le test simultané de plusieurs composants, mais si c'est le cas, cela peut être fait de la même manière.

---

Réponse au commentaire:

Il est donc vrai qu'en général la pente aléatoire $\theta_X$ permet l'interception aléatoire $\theta_0$ varier selon les niveaux de $X$?

Je ne suis pas sûr que cette question puisse recevoir une réponse raisonnable.

• Une interception aléatoire permet une différence idiosyncratique dans le niveau global pour chaque niveau du facteur de regroupement. Par exemple, si la variable dépendante est le temps de réponse, certains participants sont plus rapides et d'autres plus lents.
• Une pente aléatoire permet à chaque niveau du facteur de regroupement un effet idiosyncrasique du facteur pour lequel les pentes aléatoires sont estimées. Par exemple, si le facteur est la congruence, certains participants peuvent avoir un effet de congruence plus élevé que d'autres.

Les pentes aléatoires affectent-elles donc l'interception aléatoire? Dans un certain sens, cela pourrait avoir un sens, car ils permettent à chaque niveau du facteur de regroupement un effet complètement idiosyncrasique pour chaque condition. Au final, nous estimons deux paramètres idiosyncratiques pour deux conditions. Cependant, je pense que la distinction entre le niveau global capturé par l'interception et l'effet spécifique à la condition capturé par la pente aléatoire est importante et que la pente aléatoire ne peut pas vraiment affecter l'interception aléatoire. Cependant, il permet toujours à chaque niveau du facteur de regroupement d'être idiosyncratique séparément pour chaque niveau de la condition.

Néanmoins, mon test fait toujours ce que la question d'origine veut. Il teste si la différence de variances entre les deux conditions est nulle. S'il est nul, alors les variances dans les deux conditions sont égales. En d'autres termes, ce n'est que s'il n'y a pas besoin d'une pente aléatoire que la variance dans les deux conditions est identique. J'espère que cela à du sens.

Votre modèle

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

déjàpermet à la variance entre sujets dans la condition de contrôle de différer de la variance entre sujets dans la condition expérimentale. Cela peut être rendu plus explicite par une re-paramétrisation équivalente:

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

La matrice de covariance aléatoire a maintenant une interprétation plus simple:

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

Ici, les deux variances sont précisément les deux variances qui vous intéressent: la variance [entre sujets] des réponses moyennes conditionnelles dans la condition de contrôle et la même dans la condition expérimentale. Dans votre jeu de données simulé, ils sont de 0,25 et 0,21. La différence est donnée par

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

et est égal à 0,039. Vous voulez tester s'il est significativement différent de zéro.

ÉDITER: J'ai réalisé que le test de permutation que je décris ci-dessous est incorrect; cela ne fonctionnera pas comme prévu si les moyens en condition contrôle / expérimentale ne sont pas les mêmes (car alors les observations ne sont pas échangeables sous le nul). Il serait peut-être préférable de démarrer des sujets (ou des sujets / objets dans le cas Bonus) et d'obtenir l'intervalle de confiance pourdelta .

Je vais essayer de corriger le code ci-dessous pour ce faire.

---

Suggestion originale basée sur la permutation (fausse)

Je trouve souvent que l'on peut se sauver beaucoup de problèmes en faisant un test de permutation. En effet, dans ce cas, il est très facile à mettre en place. Permutons les conditions de contrôle / expérimentales pour chaque sujet séparément; alors toute différence de variance doit être éliminée. Répéter cela plusieurs fois donnera la distribution nulle pour les différences.

(Je ne programme pas en R; tout le monde, n'hésitez pas à réécrire ce qui suit dans un meilleur style 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

L'exécution de ceci donne la valeur de p $p=0.7$. On peut augmenter nrepà 1000 environ.

Exactement la même logique peut être appliquée dans votre cas Bonus.

En examinant ce problème sous un angle légèrement différent et en partant de la forme «générale» du modèle mixte linéaire, nous avons

$$y_{ijk} = \mu + \alpha_j + d_{ij} + e_{ijk}, \quad d_{i} \sim N(0, \Sigma), \quad e_{ijk} \sim N(0, \sigma^2)$$ where $\alpha_j$ is the fixed effect of the $j$'th condition and $d_i = (d_{i1}, \ldots, d_{iJ})^\top$ is a random vector (some call it vector-valued random effect, I think) for the $i$'th participant in the $j$'th condition.
In your example we have two conditions $y_{i1k}$ and $y_{i2k}$ which I'll denote as $A$ and $B$ in what follows. So the covariance matrix of the two-dimensional random vector $d_{i}$ is of the general form

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

with non-negative $\sigma^2_A$ and $\sigma^2_{B}$.

Let's first see how the re-parameterized version of $\Sigma$ looks when we use sum contrasts.

The variance of the intercept, which corresponds to the grand mean, is

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

The variance of the contrast is

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

And the covariance between the intercept and the contrast is

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

which, as @Jake Westfall derived slightly differently, tests the hypothesis of equal variances when we compare a model without this covariance parameter to a model where the covariance parameter is still included/not set to zero.

Notably, introducing another crossed random grouping factor (such as stimuli) does not change the model comparison that has to be done, i.e., anova(mod1, mod2) (optionally with the argument refit = FALSE when you use REML estimation) where mod1 and mod2 are defined as @Jake Westfall did.

Taking out $\sigma_{12}$ and the variance component for the contrast $\sigma^2_2$ (what @Henrik suggests) results in

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

which tests the hypothesis that the variances in the two conditions are equal and that they are equal to the (positive) covariance between the two conditions.

---

When we have two conditions, a model that fits a covariance matrix with two parameters in a (positive) compound symmetric structure can also be written as

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

or (using the categorical variable/factor condition)

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

with

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