[R-sig-ME] Why lme4 doesn't throw an error for illegitimate an model

[Phillip Alday]

Potentially something that is uninteresting or nonsensical or not
physically real ... the math for a Klein bottle is well defined, but I
doubt you'll find one in the three spatial dimensions we experience. ;)


Here it is actually capturing some aspects of the different sectors:

> summary(coef(lmer(math ~ (sector | sch.id), data =
hsb))$sch.id[,"sector"])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-2.0087 -0.1117  0.4981  0.4908  1.0710  2.8777

> summary(lmer(math ~ sector + (1|sch.id), data=hsb))
Linear mixed model fit by REML ['lmerMod']
Formula: math ~ sector + (1 | sch.id)
   Data: hsb

REML criterion at convergence: 47080.1

Scaled residuals:
    Min      1Q  Median      3Q     Max
-3.0130 -0.7523  0.0253  0.7602  2.7472

Random effects:
 Groups   Name        Variance Std.Dev.
 sch.id   (Intercept)  6.677   2.584
 Residual             39.151   6.257
Number of obs: 7185, groups:  sch.id, 160

Fixed effects:
            Estimate Std. Error t value
(Intercept)  11.3930     0.2928  38.907
sector        2.8049     0.4391   6.388


Basically, it's the (shrunken) sector by-school adjustment
from-the-zeroed-out-fixed/population-effect . Since you've omitted
sector from the fixed effects, then the adjustments correspond to
(shrunken) estimates (well, predictions, see previous fine print) of
what the total population + sch.id-level effect would be.  Because the
sch.id-level effect is in reality ideally/approximately(*) zero, these
work out to be approximations to the population level effects.

(*) "ideally/approximately" because the variability between schools may
differ between sectors. In other words, the public schools may have
adjustments distributed as N(public_mean, public_sd) and the private
schools may have adjustments distributed as N(private_mean, private_sd),
where the two SDs aren't equal. The means are handled by the fixed
effects, the SDs by the random effects. This is what you get from this
model:

> summary(lmer(math ~ sector + (1+sector|sch.id), data=hsb))
Linear mixed model fit by REML ['lmerMod']
Formula: math ~ sector + (1 + sector | sch.id)
   Data: hsb

REML criterion at convergence: 47080.1

Scaled residuals:
     Min       1Q   Median       3Q      Max
-3.01392 -0.75219  0.02518  0.76045  2.74806

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 sch.id   (Intercept)  6.7346  2.5951
          sector       0.5322  0.7295   -0.17
 Residual             39.1513  6.2571
Number of obs: 7185, groups:  sch.id, 160

Fixed effects:
            Estimate Std. Error t value
(Intercept)  11.3930     0.2939  38.762
sector        2.8048     0.4387   6.394

> summary(coef(lmer(math ~ sector + (sector | sch.id), data =
hsb))$sch.id[,"sector"])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  2.531   2.741   2.806   2.805   2.860   3.104

I'll end with a question which I won't answer but it will help you to
think about why fitting these models might be useful: how is this
related to heteroskedacity?

Phillip

On 7/10/20 11:59 pm, Simon Harmel wrote:
> Also Phillip, what does `sector` in the output of `VarCorr(mn)` below
> denote, now that you say this model is mathematically defined?
> 
> mn <- lmer(math ~ ses +  (sector | sch.id), data = hsb)
> 
>> VarCorr(mn)
>  Groups   Name        Std.Dev. Corr  
>  sch.id <http://sch.id>   (Intercept) 2.0256        
>           sector      1.3717   -0.071
>  Residual             6.0858 
>