3
$\begingroup$

I have two models - one is including a categorical covariate as a fixed effect, the other includes it as a random effect:

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

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

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

As you can see, both models m2 and m3 produce exactly the same coefficient estimate for x (including SE). Also the residual standard error is the same. The same result is produced when I simulate some missing data:

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

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

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

So can we say in general that adding effect as random will have the same impact on the other coefficients and the overall inference as adding it as fixed? If not, can you please provide a simple example (or change the provided one) when this fails?

$\endgroup$
1
  • $\begingroup$ @user11852, now you can retract it (edited). I admit that the core principles of the answer might be the same, but not at first sight at least :-) Thanks for your help. $\endgroup$ – Innate Imunity is The Way Jul 24 '13 at 16:42
1
$\begingroup$

The ideas behind the answer are the same as in this answer. Again you are a victim of your design. Try the following code:

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

 m2 <- lm(y ~ cat + x)
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
 round(digits=7,fixef(m3)[2]) == round(digits=7, coef(m2)[6]) 
 #Not in this case but results are really close because of your design matrix
 #If you did the same comparison with your original structure you would get TRUE.
 #    x 
 #FALSE 

 #y[c(1:(n/2), (n*k-n/2):(n*k))] <- NA
 # simulate missing data in an obviously unbalanced way
 ym <- y;  ym[1:66] <- NA
 m2m <- lm(y ~ cat + x) 
 m3m <- lme(ym ~ x, random = ~ 1|cat, na.action = na.omit)

 round(digits=7,fixef(m3m)[2]) == round(digits=7, coef(m2m)[6]) 
 #The results are further apart because your design changed more significantly.

As I said to you here, you are seem to ignore shrinkage in general and is one of the main reasons of you getting different estimates. Also, it needs to be emphasized that the main difference between a fixed and a random effect is that the first one is an unknown constant to estimate while a random effect is a random variable; it does not make sense to "estimate a value" for it.

$\endgroup$
2
  • 1
    $\begingroup$ I would argue that it does make sense to "estimate a value for it" in some cases. $\endgroup$ – Dason Aug 23 '13 at 21:51
  • $\begingroup$ @Dason: Obviously it makes sense to estimate the realizations of the random effects. The whole idea is though that these values are picked at random from a zero-meaned distribution of possible batch/grouping/clusters effects. Choosing to "estimate a value" though as if the levels you see are some fixed factors of some sort is misleading. $\endgroup$ – usεr11852 Aug 23 '13 at 22:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.