I found that glmer in lme4 package sometimes gives slightly different results, e.g. AIC/BIC and/or fixed/random factors, when model is a bit complex and the order of data is shuffled. Therefore, I'm picking up the best results in terms of AIC or BIC. Does this make sense?
The second question is where is the divergence coming from? In theory, the order of data should not matter as long as the i.i.d assumption is made. The different results suggest that glmer uses (perhaps iterative) algorithms that depend on initial values.
I would be grateful if you could provide any feedback or suggestions.
For your information, here is the sample code I tested, assuming that multiple practitioners make binary diagnoses based on four types of features in medial images.
# data generation
nPracts <- 40 # number of practitioners
nImages <- 100 # number of medial images
var1 <- rep(seq(-20,20,10), 20) # 4 image features
var2 <- c(rep(-20,20), rep(-10,20), rep(0,20), rep(10,20), rep(20,20))
var3 <- rep(seq(-20,25,5), 10)
var4 <- c(rep(10,20), rep(0,20), rep(-20,20), rep(20,20), rep(10,20))
fix1 <- 0.05 # 4 fixed factors
fix2 <- 0.5
fix3 <- 0.1
fix4 <- 0.3
rand1 <- rnorm(nPracts, mean=0, sd=0.1) # 4 random factors (practitioner-dependent)
rand2 <- rnorm(nPracts, mean=0, sd=0.05)
rand3 <- rnorm(nPracts, mean=0, sd=0.01)
rand4 <- rnorm(nPracts, mean=0, sd=0.03)
data <- data.frame()
for (j in 1:nPracts){
for (i in 1:nImages){
z <- var1[i]*(fix1+rand1[j])+var2[i]*(fix2+rand2[j])+var3[i]*(fix3+rand3[j])+var4[i]*(fix4+rand4[j])
pr <- 1/(1+exp(-z))
y <- rbinom(1, 1, pr) # binary diagnosis
data<- rbind(data, c(j, var1[i], var2[i], var3[i], var4[i], y))
}
}
colnames(data) <- c('practitioner', 'var1', 'var2', 'var3', 'var4', 'diagnosis')
data$diagnosis<- factor(data$diagnosis)
data$practitioner <- factor(data$practitioner)
# refit model 5 times to see the divergence
for (k in 1:5){
set.seed(123+k)
data2 <- data[sample(nrow(data)),] # change data order
fit <- glmer(diagnosis~ var1+var2+var3+var4+(0+var1|practitioner)+(0+var2|practitioner)+(0+var3|practitioner)+(0+var4|practitioner),
data=data2, family=binomial(link="logit"))
print(fit)
}
These are two samples results from the last block. They are slightly different.
AIC BIC logLik deviance df.resid
511.9535 568.5999 -246.9767 493.9535 3991
Random effects:
Groups Name Std.Dev.
practitioner var1 0.099567
practitioner.1 var2 0.080622
practitioner.2 var3 0.002188
practitioner.3 var4 0.046851
Number of obs: 4000, groups: practitioner, 40
Fixed Effects:
(Intercept) var1 var2
0.39222 0.05796 0.56366
var3 var4
0.10046 0.31254
AIC BIC logLik deviance df.resid
512.0151 568.6616 -247.0076 494.0151 3991
Random effects:
Groups Name Std.Dev.
practitioner var1 0.100646
practitioner.1 var2 0.080578
practitioner.2 var3 0.006207
practitioner.3 var4 0.044535
Number of obs: 4000, groups: practitioner, 40
Fixed Effects:
(Intercept) var1 var2
0.26154 0.05492 0.54916
var3 var4
0.10028 0.30373