I ran into (what I think is) an inconsistency when running the same 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 model. 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.