# glmer in R: Significance estimates are not robust to order of data frame

I'm using a mixed effects model with logistic link function (using lme4 version 1.1-7 in R). However, I noticed that the estimates of significance for fixed effects change depending on the order of the rows in the dataset.

That is, if I run a model on a dataset, I get certain estimate for my fixed effect and it has a certain p-value. I run the model again, and I get the same estimate and p-value. Now, I shuffle the order of rows (the data is not mixed, just the rows are in a different order). Running the model a third time, the p-value is very different.

For the data I have, the estimated p-value for the fixed effect can be between p=0.001 and p=0.08. Obviously, these are crucial differences given conventional significance levels.

I understand that the estimates are just estimated, and there will be differences between values for a number of reasons. However, the magnitude of the differences for my data seem large to me, and I wouldn't expect the order of my dataframe to have this effect (we discovered this problem by chance when a colleague ran the same model but got different results. It turned out they had ordered their data frame.).

Here is the output of my script: (X and Y are binary variables which are contrast-coded and centred, Group and SubGroup are categorical variables)

> # Fit model
> m1 = glmer(X ~Y+(1+Y|Group)+(1+Y|SubGroup),family=binomial(link='logit'),data=d)
> # Shuffle order of rows
> d = d[sample(1:nrow(d)),]
> # Fit model again
> m2 = glmer(X ~Y+(1+Y|Group)+(1+Y|SubGroup),family=binomial(link='logit'),data=d)
> summary(m1)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
['glmerMod']
Family: binomial  ( logit )
Formula: X ~ Y + (1 + Y | Group) + (1 + Y | SubGroup)
Data: d

AIC       BIC    logLik  deviance  df.resid
200692.0  200773.2 -100338.0  200676.0    189910

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.1368 -0.5852 -0.4873 -0.1599  6.2540

Random effects:
Groups       Name        Variance Std.Dev. Corr
SubGroup     (Intercept) 0.2939   0.5421
Y1          0.1847   0.4298   -0.79
Group        (Intercept) 0.2829   0.5319
Y1          0.4640   0.6812   -0.07
Number of obs: 189918, groups:  SubGroup, 15; Group, 12

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.0886     0.1325  -8.214   <2e-16 ***
Y1            0.3772     0.2123   1.777   0.0756 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr)
Y1 0.112
>
> # -----------------
> summary(m2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
['glmerMod']
Family: binomial  ( logit )
Formula: X ~ Y + (1 + Y | Group) + (1 + Y | SubGroup)
Data: d

AIC       BIC    logLik  deviance  df.resid
200692.0  200773.2 -100338.0  200676.0    189910

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.1368 -0.5852 -0.4873 -0.1599  6.2540

Random effects:
Groups       Name        Variance Std.Dev. Corr
SubGroup     (Intercept) 0.2939   0.5422
Y1          0.1846   0.4296   -0.79
Group        (Intercept) 0.2829   0.5318
Y1          0.4641   0.6813   -0.07
Number of obs: 189918, groups:  SubGroup, 15; Group, 12

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.0886     0.1166  -9.334  < 2e-16 ***
Y1            0.3773     0.1130   3.339 0.000841 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr)
Y1 0.074


I'm afraid that I can't attach the data due to privacy reasons.

Both models converge. The difference appears to be in the standard errors, while the differences in coefficient estimates are smaller. The model fit (AIC etc.) are the same, so maybe there are multiple optimal convergences, and the order of the data pushes the optimiser into different ones. However, I get slightly different estimates every time I shuffle the data frame (not just two or three unique estimates). In one case (not shown above), the model did not converge simply because of a shuffling of the rows.

I suspect that the problem lies with the structure of my particular data. It's reasonably large (nearly 200,000 cases), and has nested random effects. I have tried centering the data, using contrast coding and feeding starting values to lmer based on a previous fit. This seems to help somewhat, but I still get reasonably large differences in p-values. I also tried using different ways of calculating p-values, but I got the same problem.

Below, I've tried to replicate this problem with synthesised data. The differences here aren't as big as with my real data, but it gives an idea of the problem.

library(lme4)
set.seed(999)

# make a somewhat complex data frame
x = c(rnorm(10000),rnorm(10000,0.1))
x = sample(x)
y = jitter(x,amount=10)
a = rep(1:20,length.out=length(x))
y[a==1] = jitter(y[a==1],amount=3)
y[a==2] = jitter(x[a==2],amount=1)
y[a>3 & a<6] = rnorm(sum(a>3 & a<6))
# convert to binary variables
y = y >0
x = x >0
# make a data frame
d = data.frame(x1=x,y1=y,a1=a)

# run model

# shuffle order of rows
d = d[sample(nrow(d)),]

# run model again

# show output
summary(m1)
summary(m2)


One solution to this is to run the model multiple times with different row orders, and report the range of p-values. However, this seems inelegant and potentially quite confusing.

The problem does not affect model comparison estimates (using anova), since these are based on differences in model fit. The fixed effect coefficient estimates are also reasonably robust. Therefore, I could just report the effect size, confidence intervals and the p-value from a model comparison with a null model, rather than the p-values from within the main model.

Anyway, has anyone else had this problem? Any advice on how to proceed?

• Two thoughts in order of likeliness: 1. choice of control points in integration (so maybe try raising nAGQ). 2. Starting values depend on some row-sensitive data feature. – conjugateprior Dec 10 '14 at 15:53
• this is surprising, interesting, and slightly worrying. will take a look when I have a chance. Presumably these differences do not affect the profile likelihood confidence intervals? – Ben Bolker Dec 10 '14 at 17:54
• When I run this (with r-devel, 32-bit Linux) and compare summaries with all.equal(coef(summary(m1)),coef(summary(m2))) I get  "Mean relative difference: 0.0004398406". Is this comparable to your result, or are there cross-platform differences? – Ben Bolker Dec 11 '14 at 4:15
• Continuing the conversation on github: github.com/lme4/lme4/issues/262 . There is definitely a difference based on order, but there's something about your data that makes it much, much worse ... – Ben Bolker Dec 11 '14 at 20:23
• Also, I just realised that, for my model with real data, the random slopes and intercepts are exactly correlated (singular fit), but that doesn't explain everything. – Sean Roberts Jan 12 '15 at 14:08