Can a variable be included in a mixed model as a fixed effect and as a random effect at the same time?

I am building a mixed model, where the same 4 groups are tested repeatedly 5 weeks on the measure M. I want to test the effect of group on the measure M. Can I build a model as following (let's say with R):

lmer(M ~ Time + Group + (1|Group), data = mydata)

where group is both a fixed and random effect at the same time? Imposing two effects of a group seems contradictory conceptually. But what if I want to test the group effect?

• No. The fixed group effects already take care of the fact that objects in the same group might be more similar with respect to M than objects between groups. However, if you measure M several times for the same objects, you might consider a random object intercept to correct p-values etc. for possible intra-object correlations. – Michael M Oct 1 '13 at 16:04

First to make sure Group is defined as a factor. Your current model is $$M_{ij}=\alpha_i+\beta_1\mathrm{Time}_{ij}+u_i+e_{ij},$$ where $i$ denotes Group and $j$ the time points. If you run your model and test with ranef(), you will find that it is hard to distinguish $\alpha_i$ and $u_i$ in the estimation, thus $u_i$ almost equal 0.

Two possible alternative models are:

1. Random effects model, $$M_{ij}=\beta_0+\beta_1\mathrm{Time}_{ij}+u_i+e_{ij},$$ where $\beta_0$ is the average intercept, $u_i$ is the random individual deviance from the average intercept for each group and is assumed to follow a distribution.
2. Fixed effects model (in the context of econometrics), $$M_{ij}=\alpha_i+\beta_1\mathrm{Time}_{ij}+e_{ij},$$ where $\alpha_i$ is the fixed individual intercept.

I use the sleepstudy data set as an illustration. If the grouping variable (a factor) is included as a covariate (Model fm2 below), both the random effects and its variance tend to be zero. The intuitive explanation is that, $\alpha_i$ and $u_i$ basically model the same quantity (the group specific intercept), although one is assumed fixed and one is random. The majority of the variability is first absorbed by the fixed intercepts ($\alpha_i$), so the random intercepts $u_i$ tend to be all zero. The code and results are listed below.

> fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy, REML = F)
> re1 = as.matrix(ranef(fm1)$Subject) > summary(as.vector(re1)) Min. 1st Qu. Median Mean 3rd Qu. Max. -77.570 -7.460 5.701 0.000 15.850 71.920 > VarCorr(fm1) Groups Name Std.Dev. Subject (Intercept) 36.012 Residual 30.895 > sd(re1)  35.76336 > fm2 <- lmer(Reaction ~ Days + Subject + (1 | Subject), sleepstudy, REML = F) > re2 = as.matrix(ranef(fm2)$Subject)
> summary(as.vector(re2))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0       0       0       0       0       0
> VarCorr(fm2)
Groups   Name        Std.Dev.
Subject  (Intercept)  0.00
Residual             29.31
> sd(re2)
 0

> re3 = as.matrix(ranef(fm3)$Subject) > summary(as.vector(re3)) Min. 1st Qu. Median Mean 3rd Qu. Max. -3.998e-12 -9.220e-13 -6.549e-13 -7.717e-13 -4.567e-13 6.893e-13 > VarCorr(fm3) Groups Name Std.Dev. Subject (Intercept) 33.050 Residual 30.991 > sd(re3)  9.391908e-13 I also tested in Stata and it becomes more interesting. The mixed command uses an EM algorithm but it cannot converge and thus gives a very large estimate. In my understanding, REML and ML should not differ so much in this case. There may be some numerical issues. Given that the estimates rely on iterations, I will think more about this when I have more time. . mixed reaction days i.subject || subject:, reml Performing EM optimization: Performing gradient-based optimization: could not calculate numerical derivatives -- discontinuous region with missing values encountered could not calculate numerical derivatives -- discontinuous region with missing values encountered Computing standard errors: standard-error calculation failed Mixed-effects REML regression Number of obs = 180 Group variable: subject Number of groups = 18 Obs per group: min = 10 avg = 10.0 max = 10 Wald chi2(18) = 169.64 Log restricted-likelihood = -805.65036 Prob > chi2 = 0.0000 ... ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ subject: Identity | var(_cons) | 105231.5 . . . -----------------------------+------------------------------------------------ var(Residual) | 960.4566 . . . ------------------------------------------------------------------------------ LR test vs. linear model: chibar2(01) = 2.3e-13 Prob >= chibar2 = 1.0000 Warning: convergence not achieved; estimates are based on iterated EM • +1. Why did you write "thus$u_i$almost equal 0" in the first paragraph? Why "almost"? Can it be that this model, when estimated e.g. with lmer, will yield non-zero variance for$u_i\$? I thought it's impossible, but somebody just reported this happening in another question (stats.stackexchange.com/questions/273278, see comments). – amoeba Apr 18 '17 at 18:12