Why does the log-likelihood change when a variable is linearly transformed in a hierarchical model?

Question

I ran into (what I think is) an inconsistency when running a random-intercept model (using the lmer function in the lme4 package in R).

Here is what I do: I first run a model with a set of covariates; then I run the same model re-scaling (linearly transforming) one of the regressors. To my knowledge, this should change only the coefficient of the variable that is linearly transformed. And indeed, this is what happens when I run this "experiment" with a simple linear regression model and with a logistic model.

This code replicates the "normal" behaviour:

# Create three random independent variables
x1 <- rnorm(20)
x2 <- rnorm(20)
x3 <- as.factor(sample(0:2, 20, replace = TRUE))
# Their random coefficients
coef1 <- runif(1, -1, 1)
coef2 <- runif(1, -1, 1)
# Create a continuous dependent variable and a binomial one
y1 <- coef1 * x1 + coef2 * x2 + runif(20)
y2 <- y1
y2[which(y1 > quantile(y1, 0.5))] <- 1
y2[which(y1 <= quantile(y1, 0.5))] <- 0
# Finally, a linear transformation of x1
x1.trans <- x1*3

So, let us run an OLS model:

lm <- lm(y1 ~ x1 + x2 + x3)
summary(lm)
# OLS model with one variable linearly transformed
lm.bis <- lm(y1 ~ x1.trans + x2 + x3)
summary(lm.bis)

The coefficients of x1 and x1.trans are different, but the R-square of the two models is the same:

summary(lm)$r.sq == summary(lm.bis)$r.sq

The same with a logistic model:

logm <- glm(y2 ~ x1 + x2, family="binomial")
summary(logm)
logm.bis <- glm(y2 ~ x1.trans + x2, family="binomial")
summary(logm.bis)

Even in this case, the log-likelihood of the two models is the same:

logLik(logm) == logLik(logm.bis)

So far, so good. However, when I do the same with a hierarchical model, the log-likelihood (and consequently the AIC and BIC) of the two models are different, although the coefficient of the transformed variable remains significant with the same z value and the other coefficients are the same.

# Multilevel model
mm <- lmer(y1 ~ x1 + x2 + (1 | x3))
summary(mm)
mm.bis <- lmer(y1 ~ x1.trans + x2 + (1 | x3))
summary(mm.bis)
logLik(mm) == logLik(mm.bis) ### FALSE! ###

Why? Also the "REML criterion at convergence" is obviously different. I don't understand this result. This is probably due to my moderate knowledge of the math of hierarchical models. I'd be very happy if some of you could show me what's the trick here.

Since we then use AIC and BIC to compare models, I am puzzled by the fact that a simple transformation that shouldn't change anything makes one model better (or worse) than another.

Answer 1

A short answer for now; hopefully someone else will provide more about "why".

When using REML, likelihood values are not comparable between models with different fixed effects. Thus AIC/BIC/LRT tests are not either. I believe lme4 correctly modifies these values to the ML values when you do these kind of comparisons, but it's worth checking.

To have models that do have comparable likelihood values, fit with ML.