Group effects

Write \(j = 1, \dots, J\) for the groups. Group \(j\) contributes \(n_j\) observations \(y_j\) with its own design matrix \(D_j\), and has its own coefficient vector \(\beta_j\):

\[ y_j \mid \beta_j \sim N\!\left(D_j \beta_j,\; \sigma^2 I_{n_j}\right). \]

The \(\beta_j\) are the group effects. There are \(J \times p_{\mathrm{re}}\) of them, which is far too many to estimate independently from the data available, and \(\sigma^2\) is the observation variance describing how tightly a group’s observations cluster around its own fitted line.

Population effects

The group effects are not free parameters. They are assumed to be drawn from a common distribution, whose mean is itself a regression on group-level characteristics \(\mathcal{W}_j\):

\[ \beta_j \mid \gamma \sim N\!\left(\mathcal{W}_j \gamma,\; \Psi\right). \]

The vector \(\gamma\) holds the population effects: the coefficients that explain how group effects vary with group characteristics. The matrix \(\Psi\) is the between-group variance, describing how far groups scatter around the population pattern after those characteristics are accounted for. Finally \(\gamma\) receives a prior of its own,

\[ \gamma \sim N(\mu_0, V). \]

Two features of this parameterization matter for what follows. First, it is centered: \(\beta_j\) is the group’s own coefficient vector, not a deviation from an average. Second, both levels have the same shape — a normal likelihood paired with a normal prior — which is what allows a single sampling machine to serve the group level and the population level alike. This is the Lindley–Smith formulation (Lindley and Smith 1972), and Chapter-B01 describes how the engine exploits it.

What is fixed and what has a prior

The group effects \(\beta_j\) and the population effects \(\gamma\) always have priors and are always sampled; they are what a fit is for. The two variance quantities are the ones that may be either estimated or held at a known value, and choosing between them is what separates the chapters of this series:

Fixing both variances is what makes the posterior exactly normal, so that it can be drawn from directly, which is why C02 comes first. C03 then runs an iterative sampler at that same posterior, where the correct answer is already known and any error is visible. Each later row relaxes one more assumption, and the price of each relaxation is a convergence rate that has to be certified rather than assumed.

Variables used

The full dataset carries 32 columns covering pupils, teachers and schools. The model below uses seven of them:

knitr::kable(C01$model_vars[, c("variable", "level", "meaning", "values")])

The level column is worth noting. private_school, title1 and free_reduced_lunch are properties of the school, so they take one value for all pupils in a school and are natural candidates for explaining group effects. The distracted_* ratings vary from pupil to pupil within a school.

Other columns not used here include a treatment indicator, pupil age and sex, two earlier assessment scores, and a set of teacher-level measures:

setdiff(C01$all_columns, C01$model_vars$variable)
#>  [1] "participant_id"       "treat"                "age_months"          
#>  [4] "female"               "kindergarten"         "teacher_id"          
#>  [7] "state"                "esl_observed"         "special_ed_observed" 
#> [10] "new_student"          "distracted_a2"        "score_a1"            
#> [13] "invalid_a1"           "score_a2"             "invalid_a2"          
#> [16] "score_ppvt_age"       "invalid_ppvt"         "t_logins_april"      
#> [19] "t_logins_total"       "t_weeks_used"         "t_words_learned"     
#> [22] "t_financial_struggle" "t_misbehavior"        "t_years_experience"  
#> [25] "score_pct_change"

How the population effects are indexed

Because \(\gamma\) explains group effects, it is not a flat vector. It is grouped by which group effect each coefficient explains:

C01$popef_names
#> $`(Intercept)`
#> [1] "(Intercept)"        "private_school"     "title1"            
#> [4] "free_reduced_lunch"
#> 
#> $distracted_ppvt
#> [1] "(Intercept)"
#> 
#> $distracted_a1
#> [1] "(Intercept)"        "free_reduced_lunch"

Reading these three blocks in turn: the intercept block holds the school-level predictors of a school’s average score. The distracted_ppvt block holds a single entry, its own intercept, which is the average distracted_ppvt slope across schools. The distracted_a1 block holds that same average plus a free_reduced_lunch term — a cross-level interaction, asking whether distraction matters more in poorer schools.

This indexing is why the summary tables in later chapters label rows component::parameter, and why lining them up against a flat lmer coefficient vector takes some care.

The priors

Prior_Setup_GLMM() calibrates priors from the data rather than asking for hyperparameters directly:

ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01)
pf <- pfamily_list(ps)

For \(\gamma\) it produces the normal prior of Section 1, with pop.pwt setting its weight. At 0.01 that is about one percent of the weight of the data, so it regularizes without materially moving the answer:

C01$pop_prior
#>         component          parameter prior_mean    prior_sd
#> 1     (Intercept)        (Intercept)   88.37564  78.4502504
#> 2     (Intercept)     private_school    0.00000  96.8526878
#> 3     (Intercept)             title1    0.00000 100.2776090
#> 4     (Intercept) free_reduced_lunch    0.00000   1.5705244
#> 5 distracted_ppvt        (Intercept)    0.00000  17.5064432
#> 6   distracted_a1        (Intercept)    0.00000  41.4190827
#> 7   distracted_a1 free_reduced_lunch    0.00000   0.5841062

The prior means are zero apart from the intercept, and the standard deviations are large relative to the coefficients they cover, which is what makes this prior weak.

pfamily_list(ps) then wraps each component of \(\Psi\). A dNormal component declares \(\tau^2_k\) known, so no sampling step is needed for it; dIndependent_Normal_Gamma estimates it instead. The same choice arises one level down for \(\sigma^2\). Which are fixed is precisely the table in Section 1.

Where a quantity is fixed, it is worth knowing what it is fixed at. The natural assumption is the REML estimate, and that assumption is wrong:

tab <- C01$fixed_vs_reml
tab$ratio <- round(tab$calibrated / tab$reml, 3)
tab
#>         component calibrated      reml ratio
#> 1     (Intercept)  159.13613 191.78503 0.830
#> 2 distracted_ppvt    6.73747  26.03416 0.259
#> 3   distracted_a1   11.47851  41.05084 0.280
#> 4        Residual  388.80365 366.40576 1.061

The between-school variances are calibrated to values well below the REML ones, while the residual runs the other way, because Prior_Setup_GLMM() uses a pooled marginal estimate rather than the classical residual:

c(source = C01$dispersion_source,
  classical_reml = round(C01$dispersion_classical, 2))
#>                      source              classical_reml 
#> "pooled S_marg (A12 3.3.4)"                    "366.41"

This does not affect whether a sampler is correct: any fixed \(\Psi\) and \(\sigma^2\) still define a legitimate posterior, and reaching that posterior is what the later chapters verify. It does mean the draws are not interchangeable with an lmer plug-in analysis, which is why the reference comparisons there are judged against lmer’s own standard errors rather than against exact agreement. Chapter-B05 covers the calibration itself.

With the model specified and the priors calibrated, the first question is what happens when both variances are simply held at these values. That is Chapter-C02.