3
$\begingroup$

I am fitting a linear mixed model for the first time. I have a dataset that looks something like this:

Surgeon Handed Illness Speed
1 Left A 8
2 Right B 15
3 Right C 12
4 Left A 10
1 Left B 10
2 Right C 16
3 Right A 14
4 Left B 9
1 Left C 9
2 Right A 18
3 Right B 18
4 Left C 10

I am interested in the influence of patient illness and surgeon's handedness (left vs right) on the speed of the surgery. Each surgeon operates on each disease exactly once.

My key concern: You may notice that surgeon (random effect) implies handedness (fixed effect). My key questions: Does this violate the assumptions of mixed models? Why does the model still report a significant fixed effect for handedness?

When I run this code:

library(reshape2)
library(dplyr)
library(ggplot2)
library(pROC)
library(lme4)

data = read.csv("surgeon.csv")
data$Surgeon = as.factor(data$Surgeon)
model <- lmer(Speed ~ Handed + Illness + (1|Surgeon), data=data)
print(summary(model))

I get this output

Linear mixed model fit by REML ['lmerMod']
Formula: Speed ~ Handed + Illness + (1 | Surgeon)
   Data: data

REML criterion at convergence: 37.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5236 -0.6274  0.1120  0.6498  1.2995 

Random effects:
 Groups   Name        Variance Std.Dev.
 Surgeon  (Intercept) 0.000    0.00    
 Residual             3.458    1.86    
Number of obs: 12, groups:  Surgeon, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)    9.417      1.074   8.771
HandedRight    6.167      1.074   5.744
IllnessB       0.500      1.315   0.380
IllnessC      -0.750      1.315  -0.570

Correlation of Fixed Effects:
            (Intr) HnddRg IllnsB
HandedRight -0.500              
IllnessB    -0.612  0.000       
IllnessC    -0.612  0.000  0.500
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

It surprises me that handedness has a strong effect of $t=5.7$ when there shouldn't be any signal left after the random effect accounts for the surgeon. Where am I going wrong?

Further, I get a singular fit warning. After removing Handedness as variable the warning goes away. How do I handle this? I found a lot about random-random-effect interactions online when people talk about singular fits. But I haven't seen anything yet about a random-fixed-effect interaction.

$\endgroup$
4
  • $\begingroup$ The highest speed of a left handed surgeon is 10, the lowest speed of a right handed surgeon is 12. I'm not surprised there is "Handedness" effect (in this simulated dataset). $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 15:06
  • $\begingroup$ As @ShawnHemelstrand explains (+1), the fit is singular because the variance of the Surgeon random effect is estimated to be zero. Variances and variance components are hard to estimate and in the simulated dataset there are 4 surgeons. The general advice is to treat a grouping variable as random if there are at least 5 groups (here group = surgeon) and even then one might get a singular fit. It doesn't mean that the handedness effect is not recoverable. How many surgeons and illnesses are there in your real dataset; do you get a singular fit warning? $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 15:41
  • $\begingroup$ I'm not sure what you mean by random-fixed-effect interaction; can you clarify? Since each surgeon is either left-handed or right-handed I'm not sure what the interaction would look like, unless the variance (in Speed) is different for different-handed surgeons. You can investigate this graphically first. $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 15:52
  • $\begingroup$ @dipetkov thanks for your replies. In my real dataset I have 20 groups. With interaction effect, I dont mean a statistical interaction. What I mean is that, if we know the surgeon, we know the handedness variable. So after accounting for surgeon, there shouldnt be any effect left for handedness. However, the model summary still shows an fixed effect. $\endgroup$ Commented Oct 16, 2023 at 20:33

2 Answers 2

4
$\begingroup$

Let's consider a simpler mixed effects model without an illness effect: (Illness varies within surgeon and so is irrelevant to the core of this question.)

$$ Y_{ij} = \beta_0 + \beta_1\operatorname{I}_{h(i) = \text{right}} + \left(\epsilon_i + \epsilon_{ij}\right) $$

This model has an intercept, a fixed handedness effect $\beta_1$ (the expected difference in speed between a surgery performed by a right-handed surgeon and a left-handed one), random surgeon effects $\epsilon_i$ and residual effects (or errors) $\epsilon_{ij}$. The $\beta$ coefficients specify $\operatorname{E}(Y)$ and the random effects (in parentheses) specify $\operatorname{Var}(Y)$.

The random variables $\epsilon_i$ are independent and identically distributed with zero mean and variance $\theta^2$ and the residuals are iid $\text{Normal}(0, \sigma^2)$. The $\epsilon_i$s and $\epsilon_{ij}$s are also mutually independent.

The definition of the random effects implies that in expectation half the left-handed surgeons have positive $\epsilon_i$ and the other half — negative $\epsilon_i$. Same for the right-handed surgeons. (The $\epsilon_i$s do not depend on handedness.) Since the average surgeon effect is 0 for both left-handed and right-handed surgeons, we can estimate the handedness effect even though handedness is a group-level variable (ie. each surgeon is either left-handed or right-handed): we calculate a (weighted) average of the observed speed of right handed surgeons and subtract from it a (weighted) average of the observed speed of left handed surgeons. Here the design is balanced — there are three observations from each surgeon and two surgeons in each group, so there is no need to weight the observations and we take the difference between the two simple averages. (Check that this calculation gives 15.5 - 9.33 = 6.17, the same coefficient estimate reported in the summary table.)

This is the theory but then there is the practice of estimating the parameters: the coefficients $\beta=(\beta_0,\beta_1)$ and the variances $\theta^2$ and $\sigma^2$. For the sample dataset of 4 surgeons with 3 observations per surgeon, the lme4::lmer fitting procedure estimates the surgeon variance $\theta$ as 0. It also prints a "singular fit" warning. This means that the surgeon variance $\theta^2$ is much smaller than the residual variance $\sigma^2$ and difficult to estimate with only 4 surgeons.

A simulation may help to illustrate the theory. I use the lme4::simulate.merMod function to simulate $Y_{ij}$ for the 4 surgeon × 3 observations dataset. I specify $\theta = \sigma = 1$ and repeat the simulation three times.

I show only the estimate random effects; keep in mind that the first and fourth surgeons are left handed and the second and third surgeons are right handed.

2 out of 3 simulations the estimated $\widehat{\epsilon}_i$ are non-zero and once we get all 0s with a singular fit message. Note that when the surgeon effects are non-zero, one left-handed surgeon has a positive effect and the other one negative, with the same absolute value. So simulation matches theory: the average random effect for left-handed surgeons is 0, and so for right-handed surgeons. The population-level fixed handedness effect is estimable with this mixed effects model.

library("lme4")

data <- data.frame(
  Surgeon = c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L),
  Handed = c("Left", "Right", "Right", "Left", "Left", "Right", "Right", "Left", "Left", "Right", "Right", "Left"),
  Illness = c("A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C"),
  Speed = c(8L, 15L, 12L, 10L, 10L, 16L, 14L, 9L, 9L, 18L, 18L, 10L),
  stringsAsFactors = FALSE
)

simulate_surgery_speed <- function(seed = NULL) {
  data$Speed <- simulate(
    ~ Handed + (1 | Surgeon),
    newdata = data,
    newparams = list(
      beta = c(`(Intercept)` = 0, HandedRight = 2),
      theta = c(`Surgeon.(Intercept)` = 1),
      sigma = 1
    ),
    seed = seed
  )$sim_1
  data
}

#>   Surgeon Handed
#>         1   Left
#>         2  Right
#>         3  Right
#>         4   Left
ranef(
  lmer(Speed ~ Handed + (1 | Surgeon), data = simulate_surgery_speed(seed = 2))
)
#> $Surgeon
#>   (Intercept)
#> 1   0.5055828
#> 2  -1.0915321
#> 3   1.0915321
#> 4  -0.5055828
ranef(
  lmer(Speed ~ Handed + (1 | Surgeon), data = simulate_surgery_speed(seed = 5))
)
#> $Surgeon
#>   (Intercept)
#> 1   -0.109308
#> 2    1.025426
#> 3   -1.025426
#> 4    0.109308
ranef(
  lmer(Speed ~ Handed + (1 | Surgeon), data = simulate_surgery_speed(seed = 3))
)
#> boundary (singular) fit: see help('isSingular')
#> $Surgeon
#>   (Intercept)
#> 1           0
#> 2           0
#> 3           0
#> 4           0

$\endgroup$
3
  • $\begingroup$ thank you so much for your answer. This is very helpful. I noticed that when I fit a lm(Speed ~ Handed, data=data) the fixed effect is the same. Is that always the case? My understanding now is that the random effect merely partitions the residuals from lm further into random effect and a smaller lmer residual while removing dependencies from the residual structure "making them more iid" again. Because of that, the presence of a random effect does not influence the fixed effects fitted. Correct? $\endgroup$ Commented Oct 18, 2023 at 12:52
  • $\begingroup$ I'm glad the answer is helpful, this topic is not easy. It's not always the case that lm(Speed ~ Handed) and lmer(Speed ~ Handed + (1 | Surgeon)) give the same estimate for handedness. It would be the same if either holds: (a) the design is balanced; or (b) the Surgeon random effects are all estimated to be 0. $\endgroup$
    – dipetkov
    Commented Oct 18, 2023 at 13:48
  • 1
    $\begingroup$ I find simulations really useful to help me understand, and the simulate() function makes this really easy. (How didn't know about it before?!) Vary one setting at a time, run several simulations and look at the fixed and random effects, compare lm with lmer, etc. For example, you can add one more row for Surgeon 1 to make the design imbalanced. $\endgroup$
    – dipetkov
    Commented Oct 18, 2023 at 13:48
3
$\begingroup$

First, you may want to consider something like a Poisson distribution if your outcome variable is right skewed, as it is impossible to have negative values for speed. In that case, you may need to switch to glmer using the family=poisson argument or something analogous for your data. However, if your outcome variable is essentially normally distributed and has good variance properties, this may not be an issue.

Second, it appears plausible that this would occur, as handedness perfectly predicts which surgeon is involved. In other words, surgeons here only have one unique value, either left or right, and this does not change across conditions. Thus there is no variance to estimate when surgeon is entered as a random effect with handedness fixed. You can imagine the same would happen if you entered surgeon ID as a random effect and their name as a fixed effect...their name will never change across surgeons, so with both entered into the regression this means you can't estimate the differences.

You can see the exact same thing happens if I simulate this below:

#### Load Library ####
library(lmerTest)
library(tidyverse)

#### Assign IDs ####
id <- factor(
  rep(c("Jerry","Larry","Moe","Curly","Doug"),
    each=20)
  )

#### Assign Handedness ####
hand <- ifelse(id %in% c("Jerry","Larry","Moe"),
               "Left", "Right")
hand.numeric <- ifelse(hand == "Left",0,1)

#### Simulate Normally Distributed DV ####
speed <- hand.numeric + rnorm(100)

#### Store as Tibble ####
tib <- tibble(id,hand,speed)
tib

#### Fit ####
fit <- lmer(speed ~ hand + (1|id),
            data = tib)
summary(fit)

The variance/SD is zero and the fit is singular:

Linear mixed model fit by REML. t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: speed ~ hand + (1 | id)
   Data: tib

REML criterion at convergence: 290.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.90454 -0.65359 -0.01887  0.70280  2.26338 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.000    0.000   
 Residual             1.044    1.022   
Number of obs: 100, groups:  id, 5

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)
(Intercept) -0.04899    0.13191 98.00000  -0.371    0.711
handRight    1.05726    0.20857 98.00000   5.069 1.89e-06
               
(Intercept)    
handRight   ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr)
handRight -0.632
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Edit

It seems my example was a bit confusing, so I will just plot the data I have here. If you just dummy code handedness and plot it against speed, you will get something like this:

tib %>% 
  ggplot(aes(x=hand.numeric,
             y=speed))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()+
  labs(x="Handedness",
       y="Speed",
       title = "Handedness Associations with Speed")

enter image description here

If I now try to group the regression by each person (ID), the regression lines disappear:

tib %>% 
  ggplot(aes(x=hand.numeric,
             y=speed,
             group=id))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()+
  labs(x="Handedness",
       y="Speed",
       title = "Handedness Associations with Speed")

enter image description here

If we facet the plot instead, we can see why:

tib %>% 
  ggplot(aes(x=hand.numeric,
             y=speed))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()+
  labs(x="Handedness",
       y="Speed",
       title = "Handedness Associations with Speed")+
  facet_wrap(~id)

enter image description here

Handedness can only belong to one person, and so grouping the data (random intercepts) makes it so there is no variation in data to estimate.

Edit 2

Dipetkov noted that not including illness probably taints the simulation I used for this example. Not including actual random effect variance also probably complicates matters. However, we can still add in the predictor illness and random effect variance and still end up with problems.

#### Load Library ####
set.seed(123)
library(lmerTest)
library(tidyverse)
library(faux)

#### Assign IDs ####
id <- factor(
  rep(c("Jerry","Larry","Moe","Curly","Doug"),
      each=20)
)

#### Assign Handedness ####
hand <- ifelse(id %in% c("Jerry","Larry","Moe"),
               "Left", "Right")
hand.numeric <- ifelse(hand == "Left",0,1)
illness <- rbinom(100,1,.5)

#### Simulate Normally Distributed DV ####
speed <- hand.numeric + (20*illness) + rnorm(100)

#### Store as Tibble ####
tib <- tibble(id,hand,illness)
sigma <- rnorm(100,sd=10)
b0 <- 30
b1 <- (.70*hand.numeric)
b2 <- (.60*illness)

#### Add Random Effects ####
full.tib <- tib %>% 
  add_ranef("id", u0s = 1) %>% 
  mutate(speed = b0 + b1 + b2 + u0s + sigma) %>% 
  select(-u0s)
full.tib

#### Fit ####
fit <- lmer(speed ~ factor(hand) + factor(illness) + (1|id),
            data = full.tib)
summary(fit)

We can see that with the variance sufficiently low, this can cause issues with estimation.

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: speed ~ factor(hand) + factor(illness) + (1 | id)
   Data: full.tib

REML criterion at convergence: 721.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4935 -0.6775 -0.1736  0.5496  3.3797 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept)  0.00    0.000   
 Residual             88.72    9.419   
Number of obs: 100, groups:  id, 5

Fixed effects:
                  Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)        30.6163     1.5204 97.0000  20.137   <2e-16 ***
factor(hand)Right   0.2310     1.9238 97.0000   0.120    0.905    
factor(illness)1    0.1834     1.8883 97.0000   0.097    0.923    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) fct()R
fctr(hnd)Rg -0.525       
fctr(llns)1 -0.600  0.033
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

However, adjusting the variance helps alleviate this issue, as now it is allowed to converge with enough random effect variance to work with (this is done by simply changing u0s here to $10$):

#### Add Random Effects ####
full.tib <- tib %>% 
  add_ranef("id", u0s = 10) %>% 
  mutate(speed = b0 + b1 + b2 + u0s + sigma) %>% 
  select(-u0s)
full.tib

#### Fit ####
fit <- lmer(speed ~ factor(hand) + factor(illness) + (1|id),
            data = full.tib)
summary(fit)

Shown below:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: speed ~ factor(hand) + factor(illness) + (1 | id)
   Data: full.tib

REML criterion at convergence: 730

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4961 -0.6742 -0.1720  0.5576  3.2304 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 36.49    6.041   
 Residual             90.69    9.523   
Number of obs: 100, groups:  id, 5

Fixed effects:
                  Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)       26.83283    3.81879  3.40691   7.027  0.00386 **
factor(hand)Right 11.67157    5.84748  3.00047   1.996  0.13985   
factor(illness)1   0.05786    1.97115 95.36978   0.029  0.97664   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) fct()R
fctr(hnd)Rg -0.615       
fctr(llns)1 -0.249  0.011

So the structure of your data can also play a very large role here.

$\endgroup$
9
  • $\begingroup$ Thanks for your reply. This is only a dummy datasets, which I set up specifically for the purpose to showcase the relationship between Handedness and Surgeon, so I am not worrying about Normal vs Poisson too much. I do not understand your point "handedness perfectly predicts which surgeon is involved". Handedness does not imply a surgeon. Multiple surgeons are left and right handed. "Right" implies Surgeon "2" or "3". $\endgroup$ Commented Oct 16, 2023 at 12:13
  • $\begingroup$ I think I worded that a bit poorly, but I have edited my answer to better explain what I mean visually. $\endgroup$ Commented Oct 16, 2023 at 13:31
  • $\begingroup$ @ShawnHemelstrand One way in which your simulation doesn't reflect the OP's setup is that it doesn't have anything like "Illness" which is a within-surgeon variable. With more doctors & illnesses I suspect that the mixed model may be able to estimate the surgeon random effect just fine. $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 15:02
  • $\begingroup$ Thinking more about it, the illustration doesn't need illnesses but to simulate actual clustered data. With rnorm(100) you simulate measurement errors but no random group effects. So it's reassuring that the estimated group effect variance is 0. PS: simstudy is supposed to make simulation a bit easier but I've discovered it has quite a learning curve. $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 15:16
  • 2
    $\begingroup$ After writing paragraphs, I figured out a one-sentence summary: The "disappearing" regression line and the one-sided (pun intended) facets illustrate the fact that within-surgeon differences / basic estimators cannot recover the group-level handedness effect. We have to use between-surgeon differences instead. This turned out to be a very interesting thread! $\endgroup$
    – dipetkov
    Commented Oct 17, 2023 at 7:20

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.