I am analysing an experiment where groups received either medication or placebo treatment. I want to find out whether the different treatment types have an impact on test scores over time (compared to placebo) and whether the participant's sex also has an influence. Scores are measured 6 times over several days.

So I have 6 groups defined by sex and treatment (with t being the verum treatment, p being placebo):

  • female, p
  • female, t1
  • female, t2
  • male, p
  • male, t1
  • male, t2

How would I go for the mixed model? I was thinking of: score ~ session * sex * treatment. However, when used with lme, I have no idea what is being used as the reference as several groups (in my case all females and t1) disappear from the output. Is it possible to create a new factor and use that in the model? The factor would be derived from sex and treatment. E. g. score ~ fac * session

  • female_p
  • female_t1
  • female_t2
  • male_p
  • male_t1
  • male_t2

Does that still count as an interaction that tells me whether sex in relation to treatment has an impact over time? What alternatives are there?


2 Answers 2


Yes, in your case females are the default human being so you only get parameter estimates telling you how males differ from females. You do have two treatments and one control? 1. If you are primarily interested in the relationship between test score and treatment group, you want treatment group as a fixed effect whatever else you do. 2. You could handle sex as a random effect, and have a random effect on the intercept and on the treatment effect parameter. (1|sex) or (treatment|sex). 3. Or you could handle sex as a fixed effect, and have an interaction as treatment + sex + treatment:sex. By default, the way R works you will have contrasts such that there is a default subject (female, not treated), and you will only see parameter estimates for the non-default (male, treated, and hence two interactions male:treatment 1 and male:treatment 2). You can change the contrasts: https://rcompanion.org/rcompanion/h_01.html 4. The fun with adding a time element into this is that you might invalidate all your assumptions about errors being independent. Again, you could handle time as a random effect OR as a fixed effect with interactions. 5. The reason I suggested the longwinded notation is that by 4., you might prefer treatment + sex + time + treatment:sex + treatment:time so that you are controlling which interactions you want (rightly or wrongly, I have omitted a sex:time interaction and the three way interaction).


I'll answer just a slice of the question. (Actually, the slice that is really about R programming and not statistics).

I have no idea what is being used as the reference as several groups (in my case all females and t1) disappear from the output.

What you are seeing is likely a result of using the summary function in R after fitting a model. For this type of output, yes, some levels of the independent category variables will not be shown. To see why this makes sense, start with a simple linear model with a categorical independent variable, and use the reported estimates for the intercept and dummy variables to predict the dependent variable for a given set of independent values. You will see that it would make no sense to include an estimate for the reference levels in this kind of output.

This summary output is useful in some situations. In other situations, people are expecting an anova-like table. For lm and glm models, this can be reported with e.g. library(car); Anova(model). For lme models, you might use the lmerTest package. Comparisons among levels can be made with the emmeans package.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.