Forum Replies Created
March 10, 2022 at 09:06 GMT+0000 in reply to: using multcomp to update glmm emmeans #434
Solving specific problems without a reproducible example is generally difficult. At least for the data oin the vignette, it still works:
library("afex") # needed for mixed() and attaches lme4 automatically. library("emmeans") # emmeans is needed for follow-up tests library("multcomp") # for advanced control for multiple testing/Type 1 errors. data("fhch2010") # load fhch <- droplevels(fhch2010[ fhch2010$correct,]) # remove errors str(fhch2010) # structure of the data m9s <- mixed(log_rt ~ task*stimulus*density*frequency + (stimulus+frequency||id)+ (task|item), fhch, control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE) emm_options(lmer.df = "asymptotic") emm_i1 <- emmeans(m9s, "frequency", by = c("stimulus", "task")) update(pairs(emm_i1), by = NULL, adjust = "holm") summary(as.glht(update(pairs(emm_i1), by = NULL)), test = adjusted("free"))
Note: adjust = "tukey" was changed to "sidak" because "tukey" is only appropriate for one set of pairwise comparisons Simultaneous Tests for General Linear Hypotheses Linear Hypotheses: Estimate Std. Error z value Pr(>|z|) low - high, word, naming == 0 0.05877 0.01491 3.941 0.000162 *** low - high, nonword, naming == 0 -0.20699 0.01497 -13.823 < 2e-16 *** low - high, word, lexdec == 0 0.06353 0.01669 3.807 0.000162 *** low - high, nonword, lexdec == 0 0.11086 0.01675 6.619 1.08e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- free method)June 12, 2020 at 14:52 GMT+0000 in reply to: Contrast coding, LMMS and interactions #422
If you do want the sequential type tests in which each test does not the effects of a higher level, you can simply use
type = 2, this is exactly what it does. In
mixedthe implementation of Type II depends on the specific method, but the general approach is the same across methods.
So in this case I would run Type II tests (i.e.,
type = 2) and then look at the highest order effect that is significant. Then, refit the model that only includes those effects and use this model to inspect the effects. However, if there are interactions with continus variables, I would suggest using
emmeans::emtrendsfor this instead of looking at the parameter estimates directly (i.e.,
emmeans::emtrendsfor looking at conditional slopes of interactions of continuous with categorical/continuous variables).June 11, 2020 at 19:25 GMT+0000 in reply to: Contrast coding, LMMS and interactions #419
If you have interactions involving categorical variables, I do not see why you would use treatment contrasts. In my eyes sum-to-zero just makes more sense. So I would not spend too much time thinking about what treatment contrasts would mean in that case.
The answer to you question only depends on the balance of the two categorical variables and is essentially the same as the Type III versus Type II sums of squares question. In the balanced case it does not matter, both models would lead to the same outcome. In the case of imbalance, removing the three-way interaction is the same as using Type II tests.
In case of imbalance, refitting without the three-way interaction does not adjust (or “control”) the two-way interactions for the three-way interaction. In case of observational data that is probably preferred. In case of experimental data, I would use the full model with the three-way interaction.May 30, 2020 at 14:36 GMT+0000 in reply to: trouble when re-running (G)LMM #415
Without additional information it is difficult to say what the issue is here. It might be that in the original fit you used an older version of
lme4and they have changed some optimizer settings. But it seems unlikely that an older result was more reliable than one from the current version. In other words, if with the recent version you get different results I would probably trust those recent results more.
You might try to set
all_fit = TRUEto see if that helps. And in the second model you might want to remove the correlation as you do in the first model.
One more comment, usually one only enters variables one wants to adjust for (i.e., gender and age) only as main effects. I am not sure if the interaction of those with group is a good idea.
I do not feel compotent to answer this question. More specifically, I have no experience in changing the
nAGQsetting not am I aware of any study evaluating its effect on the error rates of a procedure.
However, if the results are consistent across different settings of
nAGQand random-effects structures that would encourage confidence in the results.May 15, 2020 at 17:58 GMT+0000 in reply to: Compute effect sizes for mixed() objects #410
The best is always to report specific contrasts from
emmeans(i.e., the estimate from the contrast). They particularly answer your research question and can easily be interpreted.
I do not understand really understand your further suggestions. There is no such thing as one
bin the case of more than two levels. Consider for example the case of the
Machinesdata with three levels.
library("afex") library("emmeans") data("Machines", package = "MEMSS") emm_options(lmer.df = "asymptotic") # simple model with random-slopes for repeated-measures factor m1 <- mixed(score ~ Machine + (Machine|Worker), data=Machines) em1 <- emmeans(m1, "Machine") em1 # Machine emmean SE df lower.CL upper.CL # A 52.4 1.68 5 48.0 56.7 # B 60.3 3.53 5 51.3 69.4 # C 66.3 1.81 5 61.6 70.9 # # Degrees-of-freedom method: kenward-roger # Confidence level used: 0.95
The parameter that is actually estimated from the model is the difference from the intercept (i.e., unweighted grand mean) as the following code shows:
round(summary(m1)$coefficients, 3) # Estimate Std. Error df t value Pr(>|t|) # (Intercept) 59.650 2.145 4.999 27.811 0.000 # Machine1 -7.294 1.079 4.999 -6.760 0.001 # Machine2 0.672 1.539 4.999 0.437 0.681 summary(em1)$emmean - mean(summary(em1)$emmean) #  -7.2944444 0.6722222 6.6222222
We see that the parameters (the two
Machine2are the first two differences. The third difference is the negative sum of the two. Thus, the average of the three differences is:
mean(summary(em1)$emmean - mean(summary(em1)$emmean)) #  0
So try to figure out which contrasts are of interest to you or answer your research question and report those.May 1, 2020 at 17:12 GMT+0000 in reply to: Correct Split-plot Mixed Model Specification #407
The warning does not seem to be emitted from
afex, but from the package performing the tests (
car), so I cannot really explain it. Try using a different
method = "S") and hopefully this works.May 1, 2020 at 11:03 GMT+0000 in reply to: Correct Split-plot Mixed Model Specification #405
scale()does not return a vector. It is safer to calculate directly. One of the following should work:
df$blockID.c <- df$blockID - 2.5
df$blockID.c <- df$blockID - mean(df$blockID)May 1, 2020 at 09:16 GMT+0000 in reply to: Correct Split-plot Mixed Model Specification #403
As far as I understand your design, the
blockIDinteraction only varies across participants. If this is indeed the case, it cannot receive random slopes. In this case,
(blockID + SM.blockDiff|VP)is indeed the maximal model justified by the design for the by-participant random intercept. For more on this I recommend our chapter:
I also do not see how centering
blockIDwould do any harm. That is usually a good idea for conintuous covariates. Again, I recommend reading:
Dalal, D. K., & Zickar, M. J. (2012). Some Common Myths About Centering Predictor Variables in Moderated Multiple Regression and Polynomial Regression. Organizational Research Methods, 15(3), 339–362. https://doi.org/10.1177/1094428111430540May 1, 2020 at 06:11 GMT+0000 in reply to: Correct Split-plot Mixed Model Specification #399
mixedshould only emit a warning if you enter a continuous covariate, not an error. This warning tells you that all model terms
blockIDinteracts with are now evaluated when
blockIDis zero. In your case this refers to all terms of your
m1(i.e., main effects of
SM.blockDiffas well as their interaction), they are now all simple effects when
blockID= 0. That is why centering makes sense as they are then evaluated at the mean of
blockID(if centered at the mean). For more on this issue see:
Dalal, D. K., & Zickar, M. J. (2012). Some Common Myths About Centering Predictor Variables in Moderated Multiple Regression and Polynomial Regression. Organizational Research Methods, 15(3), 339–362. https://doi.org/10.1177/1094428111430540
You can avoid the issue of centering by only adding
blockIDas a main effect:
quest ~ blockID + condition * SM.blockDiff + (blockID + SM.blockDiff|VP). However, I am generally not a big fan of “conrolling for” factors in this sense. It is often unclear to me what kind of substantial question are answered by that if indded shared variance is removed then.
Note that in any case, adding
blockIDonly accounts for linear effects of time. If you want to take the idea of exhaustion or attention seriously, you need to consider more complicated models that allow for non-linear effects. See:
Baayen, H., Vasishth, S., Kliegl, R., & Bates, D. (2017). The cave of shadows: Addressing the human factor with generalized additive mixed models. Journal of Memory and Language, 94, 206–234. https://doi.org/10.1016/j.jml.2016.11.006January 16, 2020 at 06:58 GMT+0000 in reply to: cl causes an error when passing weights variable #394
Sorry for the slow reply. Can you try using
weights = performance.expt2$wt?
If this still gives an error you might have to export
clusterExport(cl, "performance.expt2")November 1, 2019 at 17:15 GMT+0000 in reply to: Mixed-design experiment with non-binary non-ordinal response variable #382
If the DV is categorical with five levels, you cannot use such a mixed model. Binomial models only support categorical data with two categories. For more categories you need to use multinomial logistic models. I find these model quite advanced and their results are not easy to interpret.
You can however split your categories in what is called nested dichotomies and analyse each of these with a binomial model. The theory is described in the books by John Fox.
Let’s say you only have three categories, A, B, C. Then you first pick one category of interest, such as A. You then make a binary variable with A versus the rest (i.e., B and C). You analyse this with a binomial model. In the next step you discard all observations in which the response is A and only analyse the remaining data with a new binomial variable B versus C. We did exactly this in one paper once: http://singmann.org/download/publications/Douven_Elqayam_Singmann_Wijnbergen-Huitink-HIT_in_press.pdfOctober 31, 2019 at 14:41 GMT+0000 in reply to: Mixed-design experiment with non-binary non-ordinal response variable #376
I am not sure I fully understand the question, but I guess the simple answer is yes. The question of how to map aspects of your design onto the random-effects part of the model (i.e., for which grouping-factors to you want to have random intercepts and for which fixed factors random slopes) is independent of the choice of the family or link function. So you can easily do this for a response variable that is assumed to have a conditional normal distribution.
For an introduction to this topics see our chapter:
May 29, 2019 at 13:44 GMT+0000 in reply to: Compute effect sizes for mixed() objects #369
The idea of unstandardised effects is to simply report the effect on the response scale. For example, if the dependent variable is response time in seconds you report the effect in time (e.g., “the difference between condition was 0.1 seconds” or “the estimated slope was 1.5 seconds”).May 27, 2019 at 14:02 GMT+0000 in reply to: Default Sphericity correction method: GG? #364
The reason I have not implemented the approach you suggest is because I believe it is not statistically fully appropriate. Specifically, the Mauchly test, as any significance test, has a specific power and Type I error rate (the latter of which is controlled by the alpha/significance level). I do not believe it is worthwhile to incorporate these additional error probabilities into an analysis. Instead, I believe that the small loss in power by always using Greenhouse-Geisser is overall a more appropriate strategy and more conservative.
Consider the case in which the Mauchly test fails to correctly detect a violation because the power is too low in a specific case. If this happens, you get inflated Type I error rates. I consider this a potentially larger drawback than occasionally committing a Type II error in case Greenhouse-Geisser is too strict.
I believe that in general, significance tests for testing assumptions is not a good idea. For some arguments why see:
Nevertheless, the current development version of
afexnow makes it easier to get assumption tests, see:
https://github.com/singmann/afex/pull/69April 24, 2019 at 15:46 GMT+0000 in reply to: Non-significant explanatory variables #358
It does not matter how many levels each of your factor has. What matters if each participants sees both levels of each factor (i.e., it varies within-participants). If this is the case, corresponding random slopes should be included.April 24, 2019 at 11:29 GMT+0000 in reply to: What is happening here? #354
That can happen due to hierarchical shrinkage, which nudges the individual-level effects to follow a normal distribution. If the individual-level effects show a normal distribution around their mean, which is one of the assumptions of the mixed-model framework, this should not have too dramatic effects. Your plot suggests that this assumption is violated here. In any case, it suggests also that the specific pattern across levels of threat is not very strong.
It is difficult to say more without additional details (at least the corresponding standard errors). But removing random-slopes should only be done mildly. Try removing the correlations among slopes first. Please have a look at this discussion in my chapter: http://singmann.org/download/publications/singmann_kellen-introduction-mixed-models.pdfApril 24, 2019 at 11:18 GMT+0000 in reply to: Non-significant explanatory variables #353
I forgot to say one more thing. If any of your factors or additional variables varies within-participants, it should be added as a random slope. Random-intercept only models, as the one you fit, seem generally quite dubious. Please see my discussion on this issue at: http://singmann.org/download/publications/singmann_kellen-introduction-mixed-models.pdfApril 24, 2019 at 11:12 GMT+0000 in reply to: Non-significant explanatory variables #352
I have some answer to your questions, but it makes sense to read something more on variable selection in regression modeling as this is probably the problem your are struggling with and it is not a trivial problem. A good starting point may be ‘Regression Modeling Strategies’ by Frank Harrel (2015): https://www.springer.com/gb/book/9781441929181
Note that this book deals a lot with simple regression models (i.e., not mixed models), but the problem you have holds in the same manner.
Model selection using AIC and hypothesis testing using significance test do different things. Reasons why they can diverge can be because they use different cut-offs. However, usually a variable that dramatically improves AIC should also be significant in a significance test. Because what this test essentially does is fit a restricted model in which this variable is withheld and then compare the fit using something akin a likelihood-ratio test that compares the fit of two nested models.
If a variable alone improves AIC but not jointly is probably a case of multicollinearity: https://en.wikipedia.org/wiki/Multicollinearity
Your additional predictor variables are probably itself correlated with each other. So the predictive ability by one can then often be covered by other variables that are still correlated with it.
To answer your question directly: A variable should not really be non-significant but improve the model considerably in terms of AIC. At least not with all the other variables are in there but only the one is withhold.
I repeat again that variable selection is generally a difficult problem and many complicated procedures (not implemented in afex or lme4) exist to tackle this in a more principled manner.April 24, 2019 at 10:51 GMT+0000 in reply to: Including a covariate #351
Sorry for the slow reply, your post somehow went through the cracks. You need to set
factorize = FALSEin the call to use numerical covariates.
For the example data that would be:
data(obk.long, package = "afex") aov_ez("id", "value", obk.long, between = c("treatment", "gender"), within = c("phase", "hour"), covariate = "age", observed = c("gender", "age"), factorize = FALSE)
In your case this would be:
fit.grat2 <-aov_car(grat ~ tgrat+condition*time + Error(PPT/time), data=mf_long2, covariate=centgrat, factorize = FALSE)
I should probably add a corresponding note to the documentation.March 5, 2019 at 16:51 GMT+0000 in reply to: Mixed model specification and centering of predictor #343
Treating something as a random-effects grouping factor (i.e., to the right side of
|) leads to shrinkage. Levels of this factor (e.g., specific electrodes) for which the effect differs quite strongly from the mean are shrunk towards the overall mean. If this is an assumption you are willing to make (i.e., that specific electrodes for which the effect diverges relatively strongly from the others is probably a measurement error), then it makes sense to treat it as a random effect.
The benefit of treating a factor as random-effects grouping factor is that under this assumption (i.e., that the effect across the different electrodes is randomly distributed around the grand mean effect) the overall prediction will be better.
There are two downsides: (1) if the assumption is false, you miss out on “true” differences. (2) You are unable to easily check for differences between electrodes.
Having said that,
m2looks somewhat reasonable. However, the same caveat as before regarding the continuous covariates holds. The main effects are tested when the covariates are 0. So 0 needs to be meaningful given the data (and one way to achieve this is by centering them).
I only wonder why not have the
Distressinteraction as random slopes:
... + (FrB * Cluster | PatID)November 7, 2018 at 13:12 GMT+0000 in reply to: Trouble with ordered contrasts and lmer_alt #328
Thanks for the bug report. I have fixed this in the development version on github, which you can get via:
It might take some time until it gets on CRAN as there are some things I still want to fix before that.October 26, 2018 at 17:37 GMT+0000 in reply to: Logistic models using mixed() #326
I seem to be getting convergence errors when using afex that don’t occur when just using lme4 and glmer(). Does this mean that the results of glmer() shouldn’t be trusted?
No, this is not a good reason. Note that
glmer. So if the model uses the same parameterization (which can be achieved by running afex::set_sum_contrasts()
before fitting withglmer
) then the (full) model should be identical. In this case, runningsummary
on both models (the one fitted withafex
and the one fitted withglmer`) will reval them being identical.
The problem is really that Wald tests for generalized models are not particularly trustworthy. In my opinion, LRTs are way better then. However, as said above, if computationally feasible parametric bootstrap is the better choice. But of course, if fitting the model takes very long, then this is not a good option (as parametric bootstrap fits the model several times, preferably 1000 times or more).
Note that the convergence warnings can be false positives, some more on that:
https://biologyforfun.wordpress.com/2018/04/09/help-i-have-convergence-warnings/ (note that this blog seems somewhat too alarmist for my taste)October 25, 2018 at 18:49 GMT+0000 in reply to: Logistic models using mixed() #323
There are several questions here:
1. Is afex the best option for me?
If you want to stay in the frequentist realm, than
afexis probably as easy as it gets (but I might of course be biased).
If you can consider going Bayesian then both packages,
brmsare pretty good options. However, note that if you use those packages you have to make sure to use the correct coding (e.g.,
afex::set_sum_contrasts()) when running the model.
afexdoes so automatically, but not the other packages.
2. Should I use parametric bootstrapping?
If computationally feasible that would be great. If not, LRTs is your only remaining option.
3. However, due to the nonlinear nature of most link functions, the interpretations of most model predictions, specifically of lower-order effects in factorial designs, can be quite challenging.
What we mean here is that due to the non-linear nature of the link function, the lower order effects might not faithfully represent your lower-order effect in the data. So it is worth checking whether or not the lower order effect actually represent pattern that are in the data and not an artifact. So compare the marginal estimates on the response scale with those in the data. If this makes sense, you should be mostly fine.September 12, 2018 at 13:35 GMT+0000 in reply to: how to pass 'weights' to mixed() #311
Yeah, that should probably be better documented. You need to pass the whole vector:
> mixed(pCorrect ~ touchType + cued + (1 |exptPID), weights=touchCompare$nTotal, family=binomial, data = touchCompare)