7
$\begingroup$

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?

$\endgroup$
  • 5
    $\begingroup$ 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. $\endgroup$ – Michael M Oct 1 '13 at 16:04
7
$\begingroup$

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.

Updates:

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)
[1] 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)
[1] 0

Updates 2:

Previously I used ML in lmer because I found REML variance estimate for the random intercept seems off the true value in this extreme case. I know that it may not make sense to include both fixed and random group-specific intercepts in a single model, but it can be an interesting example. Note that the only different between this model and the random intercept model is a larger design matrix for the fixed effects.

The lmer example below uses REML, and the model is the same as Model fm2 above. The estimated random effects are all close to zero, and their standard deviation is very close to zero. I know the standard deviation of the estimated random effects is not a perfect estimate of the variance of the random effects, but the two should correspond somehow. But the REML variance estimate for the random effects is 33, which is very off zero.

> fm3 <- lmer(Reaction ~ Days + Subject + (1 | Subject), sleepstudy, REML = T, verbose = T)
> 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)
[1] 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
$\endgroup$
  • $\begingroup$ +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). $\endgroup$ – amoeba Apr 18 '17 at 18:12
  • 1
    $\begingroup$ The behaviour in fm3 is really weird. All random effect estimates are zero to machine precision, but the reported variance is high. I agree with you that this does not make sense. Can it be that the [wrong] st. deviation reported by VarCorr for Subject somehow corresponds to the st. deviation of the fixed effect for Subject? If so, this could be a bug in the VarCorr that gets confused by Subject being both random and fixed. $\endgroup$ – amoeba Apr 19 '17 at 19:40
  • $\begingroup$ Hmm, I checked the standard deviation of fixed effects for Subject, and it seems to be equal to 37.95, which is a bit different from what VarCorr reports (even though in the same ballpark). $\endgroup$ – amoeba Apr 19 '17 at 19:44
  • $\begingroup$ @amoeba I suspect that. The variance component is coded as a parameter in the likelihood function and is updated at each iteration. I guess that the REML version of lmer does not converge to the true value. The error message from Stata may provide some clue. $\endgroup$ – Randel Apr 19 '17 at 19:50

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.