Timeline for What is the lme4::lmer equivalent of a three-way repeated measures ANOVA?
Current License: CC BY-SA 3.0
14 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Dec 14, 2019 at 14:45 | comment | added | statmerkur |
@amoebasaysReinstateMonica (continued) Because variance components are constrained to be non-negative in the optimization of lme4, $\sigma^2_1$, which represents the covariance of a subject in the 1st and 2nd condition (which is negative in our data), has to be estimated as zero. This is not the case for m1 where $\Sigma$ looks like this $\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}$ (given the sum-coding) where $\sigma^2_1$ and $\sigma^2_2$ are the variance parameters for the intercept and the contrast, respectively.
|
|
| Dec 14, 2019 at 14:30 | comment | added | statmerkur |
@amoebasaysReinstateMonica (continued) The subject's intercept variance in m2 is estimated as zero and m1 fits the data better than m2 (log-likelihoods -1182.1 vs. -1187.7). For m2, this is the covariance matrix of the random vector for a subject in the two conditions: $\Sigma = \begin{bmatrix}\sigma^2_1 +\sigma^2_2 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1 + \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.
|
|
| Dec 14, 2019 at 14:08 | comment | added | statmerkur |
@amoebasaysReinstateMonica I think models like mod1 and mod3 are forced to estimate some of the variance components as zero in some situations. Consider a simpler data set with only two conditions set.seed(123); d <- MASS::mvrnorm(100, c(0, 10), cbind(c(1, -0.5), c(-0.5, 1))); d <- data.frame(y = rep(c(d[, 1], d[, 2]), 4) + rnorm(400), c = rep(rep(c(-1, 1), each = 100), 4), subj = gl(100, 1, 400)); d$f <-factor(d$c) where the 1st and the 2nd condition are negatively correlated. We define m1 <- lmer(y ~ f + (1|subj) + (0 + c|subj), d) and m2 <- lmer(y ~ f + (1|subj) + (1|subj:f), d).
|
|
| Jun 28, 2017 at 22:33 | comment | added | amoeba |
I am still not quite sure about the case when e.g. a has $k=3$ levels. Then instead of your numeric variable A, you would have A1 and A2. How would you modify the (0+A|subject) term then, given that it should (?) still have only 1 variance parameter in order to correspond to the (1|a:subject) term which in turn corresponds to the classical RM-ANOVA's Error(a:subject) term? And regarding two variance parametrizations... hmm, I have to think about it. I don't get why it's mod2/4 variances that match EMS() output when it's mod1/3 that mimic Error() structure from aov.
|
|
| Jun 28, 2017 at 15:35 | comment | added | Jake Westfall |
@amoeba That's a good insight about the zero random effect variances in mod3 vs. mod4. I agree that's probably finally the source of the difference in log-likelihoods! Perhaps in an unconstrained optimization that allowed negative variance components the two models would always return equivalent results, but in the constrained optimization of lme4 where variance components must be non-negative, different (but theoretically statistically equivalent) parameterizations can return different results when one of the parameterizations bumps into the zero boundary.
|
|
| Jun 28, 2017 at 15:30 | comment | added | Jake Westfall |
@amoeba The variances differ because the random effects are parameterized differently. In mod2 and mod4 they are parameterized the same way as the fixed effects: in terms of intercepts (=group means) and slopes (=group differences (divided by 2)). In mod1 and mod3 they are parameterized in terms of intercepts and cell means; i.e. (for $n$ subjects) the means of the $2n$ a cells (collapsing over b) and the $2n$ b cells (collapsing over a). Finally, if a or b had $k>2$ levels, I would manually replace them with $k-1$ numeric columns of the desired contrast codes.
|
|
| Jun 28, 2017 at 11:17 | comment | added | amoeba |
Getting back to this issue... I noticed that for the two-factor case where two lmer calls produce identical anova() output, the random effect variances are nevertheless quite different: see VarCorr(mod1) and VarCorr(mod2). I don't quite understand why this happens; do you? For mod3 and mod4, one can see that four out of seven variances for mod3 are actually equal to zero (for mod4 all seven are non-zero); this "singularity" in mod3 is probably why the anova tables differ. Apart from that, how would you use your "preferred way" if a and b had more than two levels?
|
|
| Oct 12, 2016 at 15:59 | comment | added | amoeba |
I keep investigating this issue and realized the following: as soon as repeated-measures factors have more than two levels, these two methods become entirely different! If A has $k$ levels then (1|A:sub) estimates only one variance parameter, whereas (0+A|sub) estimates $k-1$ variance parameters and $k(k-1)/2$ correlations between them. My understanding is that classical ANOVA estimates only one variance parameter and so should be equivalent to the first method, not the second. Right? But for $k=2$ both methods estimate one parameter, and I am still not quite sure why they disagree.
|
|
| Oct 12, 2016 at 9:43 | comment | added | amoeba | Thanks for your reply. I am rather unconvinced: I tried to fiddle with optimizer settings and could not change anything in the outcomes; my impression is that both models are well converged. I might ask this is a separate question at some moment. | |
| Oct 10, 2016 at 21:18 | comment | added | Jake Westfall | @amoeba My current thinking is still pretty much the same as then: AFAIK, the two models are statistically equivalent (in the sense that they make the same predictions about the data and imply the same standard errors), even though the random effects are parameterized differently. I think the observed differences -- which seem to only occur sometimes -- are just due to computational issues. In particular, I suspect that you could fiddle around with optimizer settings (like varying starting points or using more strict convergence criteria) until the two models returned exactly the same answer. | |
| Oct 10, 2016 at 19:05 | comment | added | amoeba |
+1. Re I am not sure why these two methods give different results, as again I think they are "in principle" equivalent, but maybe it is for some numerical/computational reasons - do you perhaps understand this better now (two years later)? I tried to figure out what is the difference, but don't get it either...
|
|
| Nov 6, 2014 at 19:22 | history | bounty ended | Henrik | ||
| Nov 4, 2014 at 18:41 | vote | accept | Henrik | ||
| Nov 4, 2014 at 18:23 | history | answered | Jake Westfall | CC BY-SA 3.0 |