1
$\begingroup$

I'm new to hierarchical models and am learning to use the lme4 package. My understanding is that the fixed effects generated from the lmer() function are suppose to match the coefficients from lm(). However, I can't seem to get the two sets of coefficients to match. Below is an example using the "mpg" data set from ggplot2. If anyone could help me identify the problem with my code or the source of the issue, I'd greatly appreciate it. Thank you!

lm.mpg = lm(hwy ~ displ + year, data = mpg)
summary(lm.mpg)
lm.mpg$coefficients
hlm.mpg = lmer(hwy ~ displ + year + (1 | class),
               data = mpg,
               control=lmerControl(optCtrl=list(maxfun=100000)))
fixef(hlm.mpg)
$\endgroup$
  • 6
    $\begingroup$ These aren't the same models. The second one has a random intercept for each class, so there is no reason that the fixed effects should be the same. $\endgroup$ – The Laconic Sep 8 '18 at 21:47
  • $\begingroup$ @TheLaconic, my understanding was that the fixed effects for other variables from a model with a random intercept for each class would still match the coefs for the linear model w/o the random effect... $\endgroup$ – mgsberger Sep 8 '18 at 23:27
  • 1
    $\begingroup$ @mgsberger this is not correct. In lm() the coefficients are estimated using ordinary least squares; these assume that the measurements within each class are uncorrelated. In a linear mixed model fitted by lmer() the coefficients are estimated using generalized least squares that account for the correlations in the measurements within each class. $\endgroup$ – Dimitris Rizopoulos Sep 9 '18 at 11:23
  • $\begingroup$ Thank you, @DimitrisRizopoulos! This is very helpful. Is there anyway to convert the estimation method to OLS? I didn't see anything in the reference materials. Also, I have seen lmer() and ln() output the same coefficients. Does this occur when the OLS classical assumptions are correct? $\endgroup$ – mgsberger Sep 9 '18 at 17:16
  • $\begingroup$ AFAIK it is not possible to set the estimation method to OLS for mixed models because this is not supported by standard theory of mixed models. It can happen that the OLS estimates are very close to the GLS estimates when for example the correlations in the repeated measurements are very low. In settings with correlations, the quality of OLS is always inferior to the one of GLS. In particular, if you have no missing data of missing data that missing completely at random, the OLS are unbiased but lose in efficiency. If you have missing at random missing data, the OLS are expected to be biased. $\endgroup$ – Dimitris Rizopoulos Sep 9 '18 at 19:09
1
$\begingroup$

As stated in comments, the two models and estimators are different :

  1. lm.mpg = lm(hwy ~ displ + year, data = mpg) provides the ML estimates for $\alpha$, $\beta_{displ}$, $\beta_{year}$ and $\sigma^2$ for: $$ hwy(k) = \alpha + displ(k) \cdot \beta_{displ} + year(k) \cdot \beta_{year} + \epsilon_k , $$ with $\epsilon_k \sim N(0, \sigma^2)$ ($\sigma$ being unknown).

  2. lmer(hwy ~ displ + year + (1 | class),data = mpg) provides the REML for $\alpha_{displ}$, $\alpha_{year}$, $\sigma$ and $\sigma_{class}$ for $$ hwy(k) = displ(k) \cdot \beta_{displ} + year(k) \cdot \beta_{year} + \epsilon_k + \alpha_{class}(class(k)) $$ with $\alpha_{class}(class) \sim N(0,\sigma^2_{class})$ and $\epsilon_k \sim N(0, \sigma^2)$ ($\sigma$ and $\sigma_{class}$ being unknown).

and thus are not expected to give the same result.

$\endgroup$
  • $\begingroup$ Thank you for the detailed response, @peuhp. However, this was a little over my head. :) I've taken a bootcamp machine learning course and have self studied since that point, so my expertise is not that deep. I recognize the coefficients, y-intercepts, and the normal error terms, but I'm not clear on the difference between alpha and beta. Also, is αclass(class(k)) just referencing the multilevel nature of class in the model? Thank you in advance for your patience! $\endgroup$ – mgsberger Sep 11 '18 at 2:09
  • $\begingroup$ @mgsberger (I corrected the notation in my answer that was quite messy, sorry for that). So $\alpha$ is the y-intercept and the $\beta$ are two slopes (they resp. multiply the displacement and year coefficient associated to kth outcome hwy(k)). In the second model, you don't have only one y-intercept but one y-intercept per class, so for a given outcome $hwy(k)$ the model adds the intercept $\alpha_{class}(class(k))$ where $class(k)$ is the class for the $k$ observations. Making a picture of the cumulative effects of the different components can help. Hope it helps! $\endgroup$ – peuhp Sep 11 '18 at 8:19

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.