Adding random effect influences coefficient estimates

Question

I have always been taught that random effects only influence the variance (error), and that fixed effects only influence the mean. But I have found an example where random effects influence also the mean - the coefficient estimate:

require(nlme)
set.seed(128)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, each = n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

# simulate missing data
y[c(1:(n/2), (n*k-n/2):(n*k))] <- NA

m1 <- lm(y ~ x)
summary(m1)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

You can see that the estimated coefficient for x from model m1 is -0.013780, while from model m3 it is 0.0011713 - both significantly different from zero.

Note that when I remove the line simulating missing data, the results are the same (it is full matrix).

Why is that?

PS: please note I am not a professional statistician, so if you are about to respond with a lot of math then please make also some simple summary for dummies :-)

Answer 1

"I have always been taught that random effects only influence the variance (error), and that fixed effects only influence the mean."

As you have discovered, this is only true for balanced, complete (i.e., no missing data) datasets with no continuous predictors. In other words, for the kinds of data/models discussed in classical ANOVA texts. Under these ideal circumstances, the fixed effects and random effects can be estimated independent of one another.

When these conditions do not hold (as they very very often do not in the "real world"), the fixed and random effects are not independent. As an interesting aside, this is why "modern" mixed models are estimated using iterative optimization methods, rather than being exactly solved with a bit of matrix algebra as in the classical mixed ANOVA case: in order to estimate the fixed effects, we have to know the random effects, but in order to estimate the random effects, we have to know the fixed effects! More relevant to the present question, this also means that when data are unbalanced/incomplete and/or there are continuous predictors in the model, then adjusting the random-effects structure of the mixed model can alter the estimates of the fixed part of the model, and vice versa.

Edit 2016-07-05. From the comments: "Could you elaborate or provide a citation on why continuous predictors will influence the estimates of the fixed part of the model?"

The estimates for the fixed part of the model will depend on the estimates for the random part of the model -- that is, the estimated variance components -- if (but not only if) the variance of the predictors differs across clusters. Which will almost certainly be true if any of the predictors are continuous (at least in "real world" data -- in theory it would be possible for this to not be true, e.g. in a constructed dataset).

Answer 2

On first level, I think all you are ignoring shrinkage toward the population values; "the per-subject slopes and intercepts from the mixed-effects model are closer to the population estimates than are the within-subject least squares estimates." [ref. 1]. The following link probably will also be of help (What are the proper descriptives to look at for my mixed-models?), see Mike Lawrence's answer).

Furthermore, I think you are marginally unlucky in your toy example because you have a perfectly balanced design that cause you to have the exact same estimate in the case of no missing values.

Try the following code which has the same process with no missing value now:

 cat <- as.factor(sample(1:5, n*k, replace=T) ) #This should be a bit unbalanced.
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma) 

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits= 7,fixef(m3)) ==  round(digits=7, coef(m1)) #Not this time lad.
 #(Intercept)           x 
 #      FALSE       FALSE 

Where now, because your design is not perfectly balanced you don't have the same coefficient estimates.

Actually if you play along with your missing value pattern in a silly way ( so for instance: y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NA) so your design is still perfectly balanced you'll get the same coefficients again.

 require(nlme)
 set.seed(128)
 n <- 100
 k <- 5
 cat <- as.factor(rep(1:k, each = n))
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
 plot(x, y)

 # simulate missing data in a perfectly balanced way
 y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NA

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits=7,fixef(m3)) ==  round(digits=7, coef(m1)) #Look what happend now...
 #(Intercept)           x 
 #       TRUE        TRUE 

You are marginally misguided by the perfect design of your original experiment. When you inserted the NA's in a non-balanced away you changed the pattern of how much "strength" could the individual subjects borrow from each other.

In short the differences you see are due to shrinkage effects and more specifically because you distorted your original perfectly-balanced design with non-perfectly-balanced missing values.

Ref 1: Douglas Bates lme4:Mixed-effects modeling with R, pages 71-72