# 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
> # Shuffle order of rows
> d = d[sample(1:nrow(d)),]
> # Fit model again
> 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. Commented Dec 10, 2014 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? Commented Dec 10, 2014 at 17:54
• 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 ... Commented Dec 11, 2014 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. Commented Jan 12, 2015 at 14:08
• Did this problem ever get sorted out, out of interest? Have you tried using other packages with different methods of computing GLMMs? Commented Jun 9, 2015 at 11:43

Following up on this question: The effect noticed is still present in more recent R and lme4 versions. The testing below helped me draw two conclusions:

1. lme4::glmer will produce the same answer with the same seed and the same dataset
2. Re-ordering the dataset introduces numerical precision effects that are larger in some cases than expected.

One problem in the original poster's question was that the fit of the simulated dataset was singular. This singular fit likely caused additional numerical instability. The examples below do not have that issue.

require(lme4)

R.version\$version.string
#> [1] "R version 4.1.2 (2021-11-01)"
utils::packageVersion("lme4")
#> [1] '1.1.31'

set.seed(1893)
cbpp_shuffle <- cbpp[sample(1:nrow(cbpp), replace = FALSE), ]
set.seed(284938)
gm0 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial)
set.seed(284938)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial)
set.seed(284938)
gm2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp_shuffle, family = binomial)

# setting the random seed with the same data order, yields the coefficients and p-value
all(coef(gm0)$$herd - coef(gm1)$$herd == 0)
#> [1] TRUE
all(summary(gm0)$$coefficients[,4] - summary(gm1)$$coefficients[,4] == 0)
#> [1] TRUE

# re-ordering the dataset causes differences in the range of numerical precision
#   on coefficients and p-value
all(abs(coef(gm1)$$herd - coef(gm2)$$herd) < 1E-9)
#> [1] TRUE
all(summary(gm1)$$coefficients[,4] - summary(gm2)$$coefficients[,4] < 1E-8)
#> [1] TRUE

# Creating a non-singular example like the original poster's
set.seed(999)

d <- data.frame(x = c(rnorm(10000, 0, 1), rnorm(10000, 0.1, 1)),
a = rep(1:20, each = 1000))
d$$y <- 6*d$$x + rnorm(20000, 0, 1)
d$$y <- d$$y + rep(rnorm(20, 0, 2), each = 1000)
d$$y <- d$$y > 4

d_shuffle <- d[sample(1:nrow(d), replace = FALSE),]

set.seed(10339)
m0 = glmer(y ~ x + (1|a), data = d, family = binomial(link='logit'))
set.seed(10339)
m1 = glmer(y ~ x + (1|a), data = d, family = binomial(link='logit'))
set.seed(10339)
m2 = glmer(y ~ x + (1|a), data = d_shuffle, family = binomial(link='logit'))

# setting the random seed with the same data order, yields the same coefficients and p-value
all(coef(m0)$$a - coef(m1)$$a == 0)
#> [1] TRUE
all(summary(m0)$$coefficients[,4] - summary(m1)$$coefficients[,4] == 0)
#> [1] TRUE

# re-ordering the dataset causes differences in the range of numerical precision
#   on coefficients and p-value
all(abs(coef(m1)$$a - coef(m2)$$a) < 1E-6)
#> [1] TRUE
all(summary(m1)$$coefficients[,4] - summary(m2)$$coefficients[,4] < 1E-4)
#> [1] TRUE


Created on 2023-12-08 with reprex v2.0.2