Overview of lme4 and linear effects model?

Question

I'm trying to wrap my head around the lme4 package and was having a bit of trouble. I was hoping to get a bit of help understanding a specific scenario I'm dealing with.

So I'm working on a repeated measure study to see if Drug B is significantly different than Drug A. I wanted to see if either the Drug, day, or interaction between the two is significant. Subjects in the trial were also taking outside pain meds (advil etc.) that I want to hold constant for.

Here's some simulated data:

sim <- data.frame(id=factor(rep(1:10, 3)),
                  Day=factor(rep(c("Day1", "Day2", "Day3"), each=10)),
                  Pain=round(c(runif(10, min=3, max=10), runif(10, min=1, max=7),
                               runif(10, min=0, max=5)),2), 
                  Drug=factor(rep(rep(c("A", "B"), each=5), 3)),
                  Pills=round(runif(30, min=0, max=3), 0))

Now from my understanding, an adequate model using lme4 would look like this

model<- lmer(Pain ~ Pills + Drug * Day + (1 | id), data=sim, REML=F)

Now I just have a couple of general questions about this stuff:

1. Is the model I wrote acceptable for this circumstance?
2. How does one usually look for significance in the summary of the model (I noticed that p-values are not calculated in lmer)?
3. Can lmer models handle unequal groups? If so is there something special I have to do with the model?

Any help would be greatly appreciated! Sorry if I'm asking simple questions, but I'm sort of new to this area.

Answer 1

Is the model I wrote acceptable for this circumstance?

Since you're asking about lme4 and mixed models, I'll only comment on what the mixed component in the model is accounting for (and what it isn't). Model building is a complex endeavor, and there are few clear answers on what is "acceptable". Here's some thoughts for consideration.

The (1|id) term is accounting for the variation in intercepts from the different subjects in the study. This means that you're accounting for different pain levels between individuals on Day 1 using Drug A and taking 0 Pills (or in other words, your baseline where coefficients are all 0), and you're assuming that the effect of Drug, Day, Drug*Day Interaction, and Pills are all equivalent between individuals. Put another way, the slopes are equivalent for each of the individuals. The command ranef(model) indicates you're only varying intercept for the individuals.

However, this isn't always the case. A nice tutorial is here, starting on page 11. Say you believed (this is best done a priori, before looking at the results) that the effect of the Pills was different for each person. Then you would vary the slope of Pills for each subject as:

model2 <- lmer(Pain ~ Pills + Drug * Day + (Pills | id), data= sim, REML =F)

Now compare results of model and model2:

Random Effects:

> ranef(model)
$id
   (Intercept)
1  -0.35394161
2  -0.01989078
3  -0.65715060
4  -0.22364587
5   1.25462885
6  -0.01051134
7  -0.20663670
8   0.69305108
9  -0.71553186
10  0.23962882
> ranef(model2)
$id
   (Intercept)       Pills
1   -0.6280053  0.17073199
2   -0.2478412  0.06737907
3   -1.0326309  0.28073508
4   -0.3661452  0.09954168
5    2.3461824 -0.63784233
6    0.0580454 -0.01578045
7   -0.3237023  0.08800297
8    1.3540668 -0.36812189
9   -1.1082416  0.30129090
10  -0.0517280  0.01406298

This estimates the effect of Pills for each individual, and could mean that controlling for Pills would also mean controlling for the varying effect of Pills on different individuals. This also has an influence on the estimated fixed effects:

Fixed Effects:

> fixef(model)
  (Intercept)         Pills         DrugB       DayDay2       DayDay3 DrugB:DayDay2 DrugB:DayDay3 
    6.0148242    -0.5405151     1.2378970    -2.7704121    -1.2936909    -0.3715879    -3.1542060 
> fixef(model2)
  (Intercept)         Pills         DrugB       DayDay2       DayDay3 DrugB:DayDay2 DrugB:DayDay3 
    5.7833549    -0.4712361     1.2989692    -2.7167875    -1.0871647    -0.3765863    -3.2227614 

However, going this route costs your model some degrees of freedom. It's a less parsimonious model, and generally increased complexity in models should only be approached when needed or scientifically justified. You can check the significance of the improvement through the anova() command:

> anova(model, model2)
Data: sim
Models:
model: Pain ~ Pills + Drug * Day + (1 | id)
model2: Pain ~ Pills + Drug * Day + (Pills | id)
       Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
model   9 122.35 134.97 -52.177   104.35                        
model2 11 124.73 140.14 -51.364   102.73 1.625      2     0.4437

Adding in Pills as a random slope wasn't significant by any of the usual measures ($p=0.44$), but don't hang model building on $p$-values. It's just one tool alongside scientific understanding of your problem (in addition to this being mock data).

Summary

It's always difficult to give a definitive answer to whether another person's model is appropriate or not, but what you have seems appropriate from what you've said. I would suggest thinking about random slopes along with random intercepts, if it's appropriate.

Quick comment on:

How does one usually look for significance in the summary of the model (I noticed that p-values are not calculated in lmer)?

Along with @Piotr's comment and links, and the above anova() command can be used to compare nested models, like a model without the drug term.

> model.null <- lmer(Pain ~ Pills + Day + (1 | id), data= sim, REML =F)
> anova(model, model.null)
Data: sim
Models:
model.null: Pain ~ Pills + Day + (1 | id)
model: Pain ~ Pills + Drug * Day + (1 | id)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
model.null  6 124.36 132.77 -56.182   112.36                           
model       9 122.35 134.97 -52.177   104.35 8.0092      3    0.04582 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Compared to a model that doesn't include the Drug term, there would be a significant difference in Pain induced by the effect of Drug ($p=0.046$).

There is also the confint() command, which will profile the likelihood to determine confidence intervals for the individual effects. On your mock data, this produced a convergence error, but presumably on real data, this hopefully will not be an issue.

EDIT

For completeness in one answer, for your third question:

Can lmer models handle unequal groups?

They certainly can. As @Piotr states in their answer, this is one of the many advantages of the mixed model.

Answer 2

2. Try summary(model) and remember, that those are not main effects, but rather treatment contrasts, please refer to great answers here, here and here.
3. One of many advantages of mixed models is that they handle unequal groups well

I'll leave the answer to the first question to someone else, though I do not see any errors in your model.