Angenommen , ich habe Teilnehmer, von denen jeder eine Antwort gibt mal 20, 10 in einem Zustand und 10 in einem anderen. Ich passe ein lineares Mischeffektmodell an, das in jeder Bedingung vergleicht. Hier ist ein reproduzierbares Beispiel, das diese Situation anhand des Pakets in simuliert :Y $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)

Das Modell mliefert zwei feste Effekte (einen Schnittpunkt und eine Steigung für die Bedingung) und drei zufällige Effekte (einen zufälligen Schnittpunkt für jeden Teilnehmer, eine zufällige Steigung für jeden Teilnehmer für die Bedingung und eine Korrelation zwischen Schnittpunkt und Steigung).

Ich möchte die Größe der zufälligen Intercept-Varianz nach Teilnehmern statistisch über die durch definierten Gruppen vergleichen condition(dh die rot hervorgehobene Varianzkomponente innerhalb der Kontroll- und Versuchsbedingungen getrennt berechnen und dann testen, ob der Unterschied in der Größe der Komponenten besteht ist nicht Null). Wie würde ich das machen (am besten in R)?

---

BONUS

Nehmen wir an, das Modell ist etwas komplizierter: Die Teilnehmer erleben jeweils 10 Stimuli 20-mal, 10 in einer Bedingung und 10 in einer anderen. Somit gibt es zwei Sätze gekreuzter zufälliger Effekte: zufällige Effekte für Teilnehmer und zufällige Effekte für Stimulus. Hier ist ein reproduzierbares Beispiel:

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)

Ich möchte die Größe der zufälligen Abfangvarianz nach Teilnehmern über die durch definierten Gruppen statistisch vergleichen condition. Wie würde ich das machen und unterscheidet sich der Prozess von dem in der oben beschriebenen Situation?

---

BEARBEITEN

Um etwas genauer zu wissen, wonach ich suche, möchte ich Folgendes wissen:

1. Ist die Frage "Sind die bedingten mittleren Antworten innerhalb jeder Bedingung (dh zufällige Schnittwerte in jeder Bedingung) wesentlich voneinander verschieden, über das hinaus, was wir aufgrund von Stichprobenfehlern erwarten würden", eine genau definierte Frage (dh ist diese Frage)? sogar theoretisch beantwortbar)? Wenn nein, warum nicht?
2. Wenn die Antwort auf Frage (1) Ja lautet, wie würde ich sie beantworten? Ich würde eine RImplementierung vorziehen , aber ich bin nicht an das lme4Paket gebunden - zum Beispiel scheint es, als ob das OpenMxPaket die Fähigkeit hat, Analysen auf mehreren Gruppen und Ebenen ( https: //openmx.ssri.psu) durchzuführen. edu / openmx-features ), und dies scheint die Art von Frage zu sein, die in einem SEM-Framework beantwortet werden sollte.

Es gibt mehr als einen Weg, diese Hypothese zu testen. Zum Beispiel sollte das von @amoeba beschriebene Verfahren funktionieren. Mir scheint jedoch, dass die einfachste und zweckmäßigste Methode zum Testen die Verwendung eines guten alten Likelihood-Ratio-Tests ist, bei dem zwei verschachtelte Modelle verglichen werden. Der einzige potenziell knifflige Teil dieses Ansatzes besteht darin, zu wissen, wie das Modellpaar so eingerichtet wird, dass das Herausfallen eines einzelnen Parameters die gewünschte Hypothese ungleicher Varianzen sauber testet. Im Folgenden erkläre ich, wie das geht.

Kurze Antwort

Wechseln Sie zur Kontrastcodierung (Summe zu Null) für Ihre unabhängige Variable und führen Sie dann einen Likelihood-Ratio-Test durch, bei dem Sie Ihr vollständiges Modell mit einem Modell vergleichen, bei dem die Korrelation zwischen zufälligen Steigungen und zufälligen Abschnitten 0 ist:

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

Visuelle Erklärung / Intuition

Damit diese Antwort sinnvoll ist, müssen Sie intuitiv verstehen, welche unterschiedlichen Werte des Korrelationsparameters für die beobachteten Daten gelten. Betrachten Sie die (zufällig variierenden) fachspezifischen Regressionslinien. Grundsätzlich steuert der Korrelationsparameter, ob sich die Regressionslinien der Teilnehmer relativ zum Punkt $X=0$ "nach rechts auffächern" (positive Korrelation) oder "nach links auffächern" (negative Korrelation) , wobei X Ihre kontrastcodierte Unabhängigkeit ist Variable. Beides impliziert eine ungleiche Varianz der bedingten Mittelwerte der Teilnehmer. Dies ist unten dargestellt:

In diesem Diagramm ignorieren wir die mehreren Beobachtungen, die wir für jedes Subjekt in jeder Bedingung haben, und zeichnen stattdessen nur die zwei zufälligen Mittelwerte jedes Subjekts auf, wobei eine Linie sie verbindet und die zufällige Steigung dieses Subjekts darstellt. (Dies sind Daten von 10 hypothetischen Personen, nicht die im OP veröffentlichten Daten.)

In der linken Spalte, in der eine starke negative Korrelation zwischen Steigung und Achsenabschnitt besteht, werden die Regressionslinien relativ zum Punkt $X=0$ nach links aufgefächert . Wie Sie in der Abbildung deutlich sehen können, führt dies bei Bedingung $X=-1$ zu einer größeren Varianz der Zufallsmittelwerte der Probanden als bei Bedingung $X=1$ .

Die rechte Spalte zeigt das umgekehrte Spiegelbild dieses Musters. In diesem Fall gibt es größere Varianz in den statistischen Mittel der Probanden in Zustand $X=1$ als im Zustand $X=-1$ .

Die Spalte in der Mitte zeigt, was passiert, wenn die zufälligen Steigungen und Abschnitte nicht korreliert sind. Dies bedeutet, dass sich die Regressionslinien genau so nach links auffächern, wie sie sich nach rechts auffächern, bezogen auf den Punkt $X=0$ . Dies impliziert, dass die Varianzen der Mittelwerte der Probanden unter den beiden Bedingungen gleich sind.

Es ist hier von entscheidender Bedeutung, dass wir ein Kontrastcodierungsschema von Summe zu Null verwendet haben, nicht Dummy-Codes (dh, dass die Gruppen nicht auf $X=0$ vs. $X=1$ ). Es ist nur unter dem Kontrast - Codierschema , dass wir diese Beziehung haben , wobei die Varianzen gleich sind , wenn und nur wenn die Steigungsabschnitt-Korrelation ist 0. Die nachfolgende Abbildung versucht aufzubauen , die Intuition:

Was diese Abbildung zeigt, ist derselbe exakte Datensatz in beiden Spalten, wobei die unabhängige Variable jedoch auf zwei verschiedene Arten codiert ist. In der linken Spalte verwenden wir Kontrastcodes - das ist genau die Situation aus der ersten Abbildung. In der rechten Spalte verwenden wir Dummy-Codes. Dies ändert die Bedeutung der Abschnitte - jetzt stellen die Abschnitte die vorhergesagten Antworten der Probanden in der Kontrollgruppe dar. Das untere Feld zeigt die Konsequenz dieser Änderung, nämlich, dass die Korrelation zwischen Steigung und Achsenabschnitt nicht mehr nahe bei 0 liegt, obwohl die Daten im tiefen Sinne gleich sind und die bedingten Varianzen in beiden Fällen gleich sind. Wenn dies immer noch nicht viel Sinn macht, kann es hilfreich sein, meine vorherige Antwort zu studieren, in der ich mehr über dieses Phänomen spreche.

Beweis

Sei $y_{ijk}$ die $j$ te Antwort des $i$ ten Subjekts unter der Bedingung $k$ . (Wir haben hier nur zwei Bedingungen, $k$ ist also entweder 1 oder 2.) Dann kann das gemischte Modell geschrieben werden:

$$y_{ijk} = \alpha_i + \beta_ix_k + e_{ijk},$$ wobei $\alpha_i$ die Subjekte sind. zufällige Abschnitte und haben Varianz $\sigma^2_\alpha$ , $\beta_i$sind die zufällige Steigung der Versuchspersonen und haben eine Varianz $\sigma^2_\beta$ , $e_{ijk}$ ist der Beobachtungsebenen-Fehlerterm und $\text{cov}(\alpha_i, \beta_i)=\sigma_{\alpha\beta}$ .

Wir wollen zeigen, dass

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

Beginnend mit der linken Seite dieser Implikation haben wir

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

Kontrastcodes von Summe zu Null implizieren, dass $x_1 + x_2 = 0$ und $x_1^2 = x_2^2 = x^2$ . Dann können wir weiter die letzte Zeile der obigen reduzieren , um

$$\begin{aligned} \sigma^2_\beta(x^2 - x^2) + 2\sigma_{\alpha\beta}(x_1 + x_1) &= 0 \\ \sigma_{\alpha\beta} &= 0, \end{aligned}$$ das wollten wir beweisen. (Um die andere Richtung der Implikation zu bestimmen, können wir genau diese Schritte in umgekehrter Reihenfolge ausführen.)

Um es noch einmal zu wiederholen, dies zeigt, dass, wenn die unabhängige Variable kontrastcodiert (Summe zu Null) ist , die Varianzen der zufälligen Mittelwerte der Subjekte in jeder Bedingung nur dann gleich sind, wenn die Korrelation zwischen zufälligen Steigungen und zufälligen Abschnitten 0 ist. Der Schlüssel mitnehmen Punkt von all dieser ist , dass die Prüfung der Nullhypothese , dass $\sigma_{\alpha\beta} = 0$ wird die Nullhypothese von gleichen Varianzen Test von dem OP beschrieben ist .

Dies funktioniert NICHT, wenn die unabhängige Variable beispielsweise Dummy-codiert ist. Insbesondere, wenn wir die Werte $x_1=0$ und $x_2=1$ in die obigen Gleichungen einfügen, finden wir, dass

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

Sie können die Signifikanz von Modellparametern mit Hilfe von geschätzten Konfidenzintervallen testen, für die das lme4-Paket die confint.merModFunktion hat.

Bootstrapping (siehe zum Beispiel Konfidenzintervall vom 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 

Wahrscheinlichkeitsprofil (siehe z. B. Welche Beziehung besteht zwischen der Wahrscheinlichkeit des Profils und den Konfidenzintervallen? )

> 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

---

• Es gibt auch eine Methode, die 'Wald'jedoch nur auf feste Effekte angewendet wird.

• Es gibt auch eine Art von Anova-Ausdruck (Likelihood-Verhältnis) in dem Paket, lmerTestdas benannt ist ranova. Aber ich kann nicht scheinen, einen Sinn daraus zu machen. Die Verteilung der Unterschiede in logLikelihood ist, wenn die Nullhypothese (Nullvarianz für den Zufallseffekt) wahr ist, nicht Chi-Quadrat-verteilt (möglicherweise ist der Likelihood-Ratio-Test sinnvoll, wenn die Anzahl der Teilnehmer und Versuche hoch ist).

---

Varianz in bestimmten Gruppen

Um Ergebnisse für die Varianz in bestimmten Gruppen zu erhalten, können Sie diese neu parametrisieren

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

Wenn wir dem Datenrahmen zwei Spalten hinzugefügt haben ^{(dies ist nur erforderlich, wenn Sie die nicht korrelierte "Kontrolle" und "experimentell" bewerten möchten. Die Funktion (0 + condition || participant_id)würde nicht zur Bewertung der verschiedenen Faktoren in der Bedingung als nicht korreliert führen.)}

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

Nun lmerwird die Varianz für die verschiedenen Gruppen angegeben

> 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 

Und Sie können die Profilmethoden auf diese anwenden. Zum Beispiel gibt Confint jetzt Konfidenzintervalle für die Kontrolle und die experimentelle Varianz an.

> 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

---

Einfachheit

Sie könnten die Likelihood-Funktion verwenden, um genauere Vergleiche zu erhalten, aber es gibt viele Möglichkeiten, Annäherungen auf der Straße vorzunehmen (z. B. könnten Sie einen konservativen Anova- / Lrt-Test durchführen, aber ist das das, was Sie wollen?).

An dieser Stelle frage ich mich, worum es eigentlich bei diesem (nicht so häufigen) Vergleich von Varianzen geht. Ich frage mich, ob es zu raffiniert wird. Warum der Unterschied zwischen Varianzen anstelle des Verhältnisses zwischen Varianzen (was sich auf die klassische F-Verteilung bezieht)? Warum nicht einfach Konfidenzintervalle melden? Wir müssen einen Schritt zurücktreten und die Daten und die Geschichte, die sie erzählen sollen, klären, bevor wir auf fortgeschrittene Pfade eingehen, die überflüssig sind und den Kontakt mit der statistischen Materie und den statistischen Überlegungen verlieren, die eigentlich das Hauptthema sind.

Ich frage mich, ob man viel mehr tun sollte, als nur die Konfidenzintervalle anzugeben (die tatsächlich viel mehr aussagen als einen Hypothesentest. Ein Hypothesentest gibt eine Ja-Nein-Antwort, aber keine Information über die tatsächliche Ausbreitung der Population. Wenn Sie genügend Daten haben, können Sie einen geringfügigen Unterschied machen, der als signifikanter Unterschied ausgewiesen wird). Um tiefer in die Materie einzusteigen (für welchen Zweck auch immer), bedarf es meines Erachtens einer spezifischeren (eng definierten) Forschungsfrage, um den mathematischen Maschinen die richtigen Vereinfachungen zu geben (auch wenn eine genaue Berechnung möglich ist oder wann) es könnte durch Simulationen / Bootstrapping angenähert werden, selbst dann bedarf es in einigen Einstellungen noch einer angemessenen Interpretation). Vergleichen Sie mit dem genauen Test von Fisher, um eine (bestimmte) Frage (über Kontingenztabellen) genau zu lösen.

Einfaches Beispiel

Um ein Beispiel für die Einfachheit zu geben, die möglich ist, zeige ich im Folgenden einen Vergleich (durch Simulationen) mit einer einfachen Bewertung des Unterschieds zwischen den beiden Gruppenvarianzen auf der Grundlage eines F-Tests, der durch Vergleichen der Varianzen in den einzelnen mittleren Antworten und durch Vergleichen durchgeführt wird das gemischte Modell abgeleitet Varianzen.

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

Sie können dies in der Simulation der folgenden Grafik sehen, in der neben dem auf der Stichprobe basierenden F-Score ein F-Score berechnet wird, der auf den vorhergesagten Varianzen (oder Quadratsummen) des Modells basiert.

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

Sie können sehen, dass es einen Unterschied gibt. Dieser Unterschied kann auf die Tatsache zurückzuführen sein, dass das lineare Modell mit gemischten Effekten die Quadratsummen des Fehlers (für den Zufallseffekt) auf andere Weise ermittelt. Und diese quadratischen Fehlerausdrücke werden (nicht mehr) gut als einfache Chi-Quadrat-Verteilung ausgedrückt, sind aber immer noch eng miteinander verwandt und können angenähert werden.

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

Das auf den Mitteln basierende Modell ist also sehr genau. Aber es ist weniger mächtig. Dies zeigt, dass die richtige Strategie davon abhängt, was Sie wollen / brauchen.

^{Wenn Sie im obigen Beispiel die rechten Endgrenzen auf 2,1 und 3,1 setzen, erhalten Sie bei gleicher Varianz ungefähr 1% der Bevölkerung (bzw. 103 und 104 der 10 000 Fälle), bei ungleicher Varianz unterscheiden sich diese Grenzen jedoch viel (mit 5334 und 6716 Fällen)}

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

Ein relativ einfacher Weg könnte sein, Likelihood-Ratio-Tests über anovawie in den lme4FAQ beschrieben zu verwenden .

Wir beginnen mit einem vollständigen Modell, in dem die Varianzen nicht eingeschränkt sind (dh zwei verschiedene Varianzen sind zulässig), und passen dann zu einem eingeschränkten Modell, in dem angenommen wird, dass die beiden Varianzen gleich sind. Wir vergleichen sie einfach mit anova()(beachte, dass ich setze, REML = FALSEobwohl es REML = TRUEmit anova(..., refit = FALSE)durchaus machbar ist ).

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

Dieser Test ist jedoch wahrscheinlich konservativ . In den FAQ heißt es beispielsweise:

Beachten Sie, dass LRT-basierte Nullhypothesentests konservativ sind, wenn sich der Nullwert (z. B. σ2 = 0) an der Grenze des zulässigen Raums befindet. Im einfachsten Fall (Varianz einzelner zufälliger Effekte) ist der p-Wert ungefähr doppelt so groß wie er sein sollte (Pinheiro und Bates 2000).

Es gibt verschiedene Alternativen:

1. Erstellen Sie eine geeignete Testverteilung, die in der Regel aus einer Mischung von besteht $\chi^2$Verteilungen. Siehe z. B.
   Self, SG & Liang, K.-Y. (1987). Asymptotische Eigenschaften von Maximum-Likelihood-Schätzern und Likelihood-Ratio-Tests unter nicht standardisierten Bedingungen. Journal of the American Statistical Association, 82 (398), 605. https://doi.org/10.2307/2289471 Dies ist jedoch recht kompliziert.

2. Simulieren Sie die korrekte Verteilung mit RLRsim(wie auch in den FAQ beschrieben).

Ich werde die zweite Option im Folgenden demonstrieren:

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

Wie wir sehen können, deutet die Ausgabe darauf hin, dass REML = TRUEwir mit genaue Ergebnisse erzielt hätten. Dies ist jedoch eine Übung für den Leser.

In Bezug auf den Bonus bin ich mir nicht sicher, ob das RLRsimgleichzeitige Testen mehrerer Komponenten möglich ist, aber in diesem Fall kann dies auf die gleiche Weise erfolgen.

---

Antwort auf Kommentar:

Es ist also wahr, dass im Allgemeinen die zufällige Steigung $\theta_X$ erlaubt das zufällige Abfangen $\theta_0$ über Ebenen von variieren $X$?

Ich bin nicht sicher, ob diese Frage eine vernünftige Antwort erhalten kann.

• Ein zufälliger Schnittpunkt ermöglicht einen eigenwilligen Unterschied in der Gesamtstufe für jede Stufe des Gruppierungsfaktors. Wenn die abhängige Variable beispielsweise die Antwortzeit ist, sind einige Teilnehmer schneller und andere langsamer.
• Eine zufällige Steigung ermöglicht jeder Ebene des Gruppierungsfaktors einen eigenwilligen Effekt des Faktors, für den zufällige Steigungen geschätzt werden. Wenn der Faktor beispielsweise Kongruenz ist, können einige Teilnehmer einen höheren Kongruenzeffekt haben als andere.

Beeinflussen Zufallssteigungen also den Zufallsabschnitt? In gewissem Sinne kann dies sinnvoll sein, da sie jeder Ebene des Gruppierungsfaktors für jede Bedingung einen völlig eigenwilligen Effekt verleihen. Am Ende schätzen wir zwei idiosynkratische Parameter für zwei Zustände. Ich denke jedoch, dass die Unterscheidung zwischen dem durch den Schnittpunkt erfassten Gesamtpegel und dem durch die zufällige Steigung erfassten bedingungsspezifischen Effekt wichtig ist, und dass die zufällige Steigung dann den zufälligen Schnittpunkt nicht wirklich beeinflussen kann. Es lässt jedoch weiterhin zu, dass jede Ebene des Gruppierungsfaktors für jede Ebene der Bedingung eine eigene Identität aufweist.

Trotzdem macht mein Test immer noch, was die ursprüngliche Frage will. Es wird geprüft, ob der Unterschied in den Abweichungen zwischen den beiden Bedingungen Null ist. Wenn es Null ist, sind die Varianzen in beiden Zuständen gleich. Mit anderen Worten, nur wenn keine Zufallssteigung erforderlich ist, ist die Varianz unter beiden Bedingungen identisch. Ich hoffe das ergibt Sinn.

Dein Modell

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

ermöglicht bereits, dass sich die subjektübergreifende Varianz in der Kontrollbedingung von der subjektübergreifenden Varianz in der Versuchsbedingung unterscheidet. Dies kann durch eine äquivalente Neuparametrisierung deutlicher gemacht werden:

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

Die zufällige Kovarianzmatrix ist jetzt einfacher zu interpretieren:

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

Hierbei sind die beiden Varianzen genau die beiden Varianzen, an denen Sie interessiert sind: die Varianz der bedingten Mittelwerte unter den Kontrollbedingungen und dieselbe unter den experimentellen Bedingungen. In Ihrem simulierten Datensatz sind dies 0,25 und 0,21. Der Unterschied ist gegeben durch

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

und ist gleich 0,039. Sie möchten testen, ob es sich erheblich von Null unterscheidet.

EDIT: Ich habe festgestellt, dass der unten beschriebene Permutationstest falsch ist. es funktioniert nicht wie beabsichtigt, wenn die Mittelwerte im Kontroll- / Versuchszustand nicht die gleichen sind (da dann die Beobachtungen nicht unter der Null austauschbar sind). Es ist möglicherweise eine bessere Idee, Themen (oder Themen / Gegenstände im Bonusfall) zu booten und das Konfidenzintervall für zu erhalten delta.

Ich werde versuchen, den folgenden Code zu korrigieren, um das zu tun.

---

Ursprünglicher permutationsbasierter Vorschlag (falsch)

Ich stelle oft fest, dass man sich durch einen Permutationstest viel Ärger ersparen kann. In diesem Fall ist die Einrichtung sehr einfach. Lassen Sie uns die Kontroll- / Versuchsbedingungen für jedes Subjekt separat permutieren. dann sollte jeder Unterschied in den Abweichungen beseitigt werden. Wenn Sie dies mehrmals wiederholen, erhalten Sie die Nullverteilung für die Differenzen.

(Ich programmiere nicht in R; jeder kann das Folgende in einem besseren R-Stil umschreiben.)

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

Wenn Sie dies ausführen, erhalten Sie den p-Wert $p=0.7$. Man kann nrepauf 1000 oder so erhöhen .

Genau die gleiche Logik kann in Ihrem Bonusfall angewendet werden.

Betrachtet man dieses Problem aus einer etwas anderen Perspektive und geht man von der "allgemeinen" Form des linearen Mischmodells aus

$$y_{ijk} = \mu + \alpha_j + d_{ij} + e_{ijk}, \quad d_{i} \sim N(0, \Sigma), \quad e_{ijk} \sim N(0, \sigma^2)$$ wo $\alpha_j$ ist die feste Wirkung der $j$'th Bedingung und $d_i = (d_{i1}, \ldots, d_{iJ})^\top$ ist ein Zufallsvektor (manche nennen ihn vektorwertiger Zufallseffekt, glaube ich) für die $i$Teilnehmer an der $j$'th Bedingung.
In Ihrem Beispiel haben wir zwei Bedingungen$y_{i1k}$ und $y_{i2k}$ was ich als bezeichnen werde $A$ und $B$im folgenden. Also die Kovarianzmatrix des zweidimensionalen Zufallsvektors$d_{i}$ ist von der allgemeinen Form

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

mit nicht negativ $\sigma^2_A$ und $\sigma^2_{B}$.

Schauen wir uns zunächst an, wie die Version von umparametriert wird $\Sigma$ sieht aus, wenn wir Summenkontraste verwenden.

Die Varianz des Abschnitts, die dem großen Mittelwert entspricht, ist

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

Die Varianz des Kontrastes ist

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

Und die Kovarianz zwischen dem Schnittpunkt und dem Kontrast ist

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