Introductie op multilevel modelleren met gebruik van rstanarm: Een tutorial voor onderwijsonderzoekers

Een tutorial over multilevel Bayes

JoonHo Lee, Nicholas Sim, Feng Ji, and Sophia Rabe-Hesketh, vertaling HarrieJonkman
2022-02-06

Een paar jaar geleden schreven Joonho Lee en collega’s van de Berkley Universiteit onder leiding van Sophia Rabe-Hesketh (van wie ik al eerder veel geleerd heb over multilevel analyse) een mooie tutorial over het gebruik van rstanarm in multilevelanalyse. Omdat deze tutorial in vele opzichten leerzaam is, heb ik deze vertaald. Zo maakte ik mij het werken met rstanarm eigen maar ook hoe je zo’n tutorial in rmarkdown zet. Alle credits gaan natuurlijk naar Joonhoo Lee e.a.. De oorspronkelijke tutorial van 24 april 2018 vind je hier. Later vond ik ook op github de syntax en kon ik diverse eigen fouten weer bijstellen hier.

Introductie

Een veel voorkomend kenmerk van datastructuren in het onderwijs is dat analyse-eenheden (b.v. leerlingen) genest zijn in hogere organisatorische clusters (b.v. scholen). Dit soort structuren leidt tot afhankelijkheid tussen de antwoorden die worden waargenomen voor eenheden binnen dezelfde cluster. Leerlingen in dezelfde school hebben in het algemeen de neiging meer gelijkenissen te vertonen in hun academische en attitudinale kenmerken dan leerlingen die willekeurig uit de bevolking worden gekozen.

Multilevel modellen1 zijn ontworpen om dergelijke afhankelijkheid binnen clusters te modelleren. Multilevelmodellen erkennen het bestaan van gegevensclustering (op twee of meer niveaus) door residuele componenten op elk niveau in de hiërarchie toe te staan. Bijvoorbeeld, een model met twee niveaus dat de groepering van leerlingresultaten binnen scholen toelaat, zou residuen op zowel leerling- als schoolniveau omvatten. De residuele variantie wordt dan verdeeld in een tussen-school component (de variantie van de residuen op schoolniveau) en een binnen-school component (de variantie van de residuen op leerlingniveau).

In deze tutorial laten Joonhoo Lee en collega’s zien hoe je een multilevel lineair model kunt fitten binnen een volledig Bayesiaans raamwerk met rstanarm. Deze handleiding is in de eerste plaats gericht op onderwijsonderzoekers die lme4 in R hebben gebruikt om modellen te fitten op hun gegevens en die mogelijk geïnteresseerd zijn in het leren fitten van Bayesiaanse multilevel modellen. Voor lezers die lme4 nog niet eerder hebben gebruikt, geven we echter een kort overzicht van het gebruik van het pakket voor het fitten van multilevel modellen.

Deze tutorial maakt gebruik van Stan versie 2.17.3 en vereist de volgende R pakketten.

# Pakketten die nodig zijn
library(mlmRev)
library(lme4)
library(rstanarm)
library(ggplot2)

Data voorbeeld

Ze analyseren de Gcsemv dataset (Rasbash et al. 2000) uit het mlmRev pakket in R. De gegevens omvatten de GCSE-examenscores (General Certificate of Secondary Education) van 1.905 leerlingen van 73 scholen in Engeland op een natuurwetenschappelijk vak. De Gcsemv-dataset bestaat uit de volgende 5 variabelen:

# Use example dataset from mlmRev package: GCSE exam score
data(Gcsemv, package = "mlmRev")
summary(Gcsemv)
     school        student     gender      written     
 68137  : 104   77     :  14   F:1128   Min.   : 0.60  
 68411  :  84   83     :  14   M: 777   1st Qu.:37.00  
 68107  :  79   53     :  13            Median :46.00  
 68809  :  73   66     :  13            Mean   :46.37  
 22520  :  65   27     :  12            3rd Qu.:55.00  
 60457  :  54   110    :  12            Max.   :90.00  
 (Other):1446   (Other):1827            NA's   :202    
     course      
 Min.   :  9.25  
 1st Qu.: 62.90  
 Median : 75.90  
 Mean   : 73.39  
 3rd Qu.: 86.10  
 Max.   :100.00  
 NA's   :180     

Twee onderdelen van het examen werden geregistreerd als uitkomstvariabelen: schriftelijk werkstuk (written) en cursuswerkstuk (course). In deze tutorial wordt alleen de totaalscore op het cursuswerkstuk (course) geanalyseerd. Zoals hierboven te zien is, ontbreken er bij sommige waarnemingen waarden voor bepaalde covariaten. Hoewel we de data niet subsetten om alleen volledige gevallen op te nemen om aan te tonen dat rstanarm deze waarnemingen automatisch laat vallen, is het over het algemeen een goed gebruik om dit handmatig te doen indien nodig.

# Maak Male dereferentiecategorie en hernoem de variabele
Gcsemv$female <- relevel(Gcsemv$gender, "M")


# Gebruik alleen de totaalscore op course 
GCSE <- subset(x = Gcsemv, 
               select = c(school, student, female, course))

# Tel unieke scholen en studenten
J <- length(unique(GCSE$school))
N <- nrow(GCSE)

Het pakket rstanarm automatiseert verschillende stappen van datavoorbewerking, waardoor het gebruik ervan sterk lijkt op dat van lme4 en wel op de volgende manier.

# Check structure of data frame
str(GCSE)
'data.frame':   1905 obs. of  4 variables:
 $ school : Factor w/ 73 levels "20920","22520",..: 1 1 1 1 1 1 1 1 1 2 ...
 $ student: Factor w/ 649 levels "1","2","3","4",..: 16 25 27 31 42 62 101 113 146 1 ...
 $ female : Factor w/ 2 levels "M","F": 1 2 2 2 1 2 2 1 1 2 ...
 $ course : num  NA 71.2 76.8 87.9 44.4 NA 89.8 17.5 32.4 84.2 ...

Likelihood inferentie met lmer()

In dit deel bespreken we kort drie lineaire multilevelmodellen die in deze tutorial zullen worden toegepast. We beginnen met een variërend intercept model zonder voorspellers (Model 1), dan gaan we verder met het variërend intercept model met één voorspeller (Model 2), en het variërend intercept en helling model (Model 3).

Model 1: Variërend intercept model zonder voorspellers (Variantiecomponentenmodel)

Beschouw het eenvoudigste multilevel model voor leerlingen \(i=1,...,n\) genest binnen scholen \(j=1,...,J\) en voor wie we examenresultaten als respons hebben. We kunnen een variabel interceptmodel met twee niveaus zonder voorspellers schrijven met de gebruikelijke tweedelige formulering als

\[Y_{ij}=\alpha_{j} + \epsilon_{ij}, \text{ where } \epsilon{ij} \sim N(0,\sigma^2_{y}) \] \[\alpha{_j}=u_{\alpha} + u_{j}, \text { where } \epsilon{ij} \sim N(0,\sigma^2_{y})\] waarin \(yij\) de examenscore is voor de \(ith\) leerling op de \(jth\) school, \(\alpha_{j}\) de variërende intercept voor de jde school, en \(u_{a}\) het algemene gemiddelde voor alle scholen. Als alternatief kan het model in verkorte vorm worden uitgedrukt als

\[y_{ij}=u_{a} + u_{j} = \epsilon_{ij}\] Als we verder aannemen dat de fouten op leerlingniveau \(\epsilon_{ij}\) normaal verdeeld zijn met gemiddelde 0 en variantie \(\sigma^2_{y}\), en dat de variërende intercepten op schoolniveau \(\alpha_j\) normaal verdeeld zijn met gemiddelde \(u_{a}\) en variantie \(\sigma^2_{a}\), dan kan het model worden uitgedrukt als

\[y_{ij}∼ N(\alpha_{j},\sigma^2_{y})\] \[a_{j}∼ N(u_{a},\sigma^2_{a})\]

Dit model kan dan worden aangepast met lmer(). We specificeren een intercept (de voorspeller “1”) en laten deze variëren met de niveau-2-identifier (school). We specificeren ook de REML = FALSE optie om maximum likelihood (ML) schattingen te krijgen in plaats van de standaard restricted maximum likelihood (REML) schattingen.

M1 <- lmer(formula = course ~ 1 + (1 | school), 
           data = GCSE, 
           REML = FALSE)
summary(M1)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: course ~ 1 + (1 | school)
   Data: GCSE

     AIC      BIC   logLik deviance df.resid 
 14111.4  14127.7  -7052.7  14105.4     1722 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.9693 -0.5101  0.1116  0.6741  2.7613 

Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept)  75.24    8.674  
 Residual             190.77   13.812  
Number of obs: 1725, groups:  school, 73

Fixed effects:
            Estimate Std. Error t value
(Intercept)    73.72       1.11    66.4

Onder Fixed effects zien we dat het intercept \(u_{a}\), gemiddeld over de populatie van scholen, wordt geschat op 73,72. Onder Random-effects zien we dat de standaardafwijking tussen de scholen \(\sigma_{a}\) wordt geschat op \(8.67\) en de standaardafwijking binnen de scholen \(\sigma_{y}\) op \(13.81\).

Model 2: Variërend intercept model met een enkele voorspeller

Het variërende interceptmodel5 met een indicatorvariabele voor het vrouwzijn xij kan worden geschreven als

\[Y_{ij}∼N(a_{j}+\beta x{ij}, \sigma^2_{y}),\]
\[a_{j}∼N(u_{a}, \sigma^2_{a})\].

De vergelijking van de gemiddelde regressielijn voor alle scholen is \(u_{ij}=u_{α}+βxij\). De regressielijnen voor specifieke scholen zullen evenwijdig zijn aan de gemiddelde regressielijn (met dezelfde helling β), maar verschillen wat betreft het intercept \(a_{j}\). Dit model kan geschat worden door vrouwelijk toe te voegen aan de formule in de lmer() functie, waardoor alleen het intercept per school varieert en de “helling” voor vrouwelijkheid constant blijft voor alle scholen

M2 <- lmer(formula = course ~ 1 + female + (1 | school), 
           data = GCSE, 
           REML = FALSE)
summary(M2)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: course ~ 1 + female + (1 | school)
   Data: GCSE

     AIC      BIC   logLik deviance df.resid 
 14017.4  14039.2  -7004.7  14009.4     1721 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.7809 -0.5401  0.1259  0.6795  2.6753 

Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept)  76.65    8.755  
 Residual             179.96   13.415  
Number of obs: 1725, groups:  school, 73

Fixed effects:
            Estimate Std. Error t value
(Intercept)   69.730      1.185   58.87
femaleF        6.739      0.678    9.94

Correlation of Fixed Effects:
        (Intr)
femaleF -0.338

De gemiddelde regressielijn over de scholen wordt dus geschat als \(\hat{\mu}_{ij}=69.73+ 6.74 x_{ij}\), waarbij \(\sigma_\_alpha\) en \(\sigma_y\) worden geschat als respectievelijk \(8.76\) en \(13.41\). Als we deze schattingen van \(\mu_alpha\), \(\beta\), \(\sigma^2_{y}\), en \(\sigma^2_{alpha}\) als de ware parameterwaarden beschouwen, kunnen we de **Best Linear Unbiased Predictions (BLUPs) voor de fouten op schoolniveau \(\hat{u}_j = \hat{\alpha}_{j} - \hat{\mu}_{\alpha}\) verkrijgen.

De BLUPs zijn equivalent met de zogenaamde Empirical Bayes (EB)-voorspelling, die het gemiddelde is van de posterieure verdeling van \(u_{j}\) gegeven alle geschatte parameters, alsmede de willekeurige variabelen \(y_{ij}\) en \(x_{ij}\) voor het cluster. Deze voorspellingen worden “Bayes” genoemd omdat ze gebruik maken van de vooraf gespecificeerde prioriteitsverdeling1 \(u_j \sim N(\mu_\alpha, \sigma^2_\alpha)\), en bij uitbreiding \(u_j \sim N(0, \sigma^2_\alpha)\), en “Empirisch” genoemd omdat de parameters van deze prior, \(\mu_\alpha\) en \(\sigma^2_{\alpha}\), naast \(\beta\) en \(\sigma^2_{y}\), geschat worden uit de data.

In vergelijking met de ML-benadering (Maximum Likelihood - maximale waarschijnlijkheid), waarbij waarden voor \(u_j\) worden voorspeld door alleen de geschatte parameters en gegevens van cluster \(j\) te gebruiken, houdt de EB-benadering bovendien rekening met de voorafgaande verdeling van \(u_{j}\), en levert zij voorspelde waarden op die dichter bij \(0\) liggen (een verschijnsel dat wordt beschreven als shrinkage of partial pooling). Om te zien waarom dit verschijnsel shrinkage wordt genoemd, drukken we de uit EB voorspelling verkregen schattingen voor \(u_j\) gewoonlijk uit als \(\hat{u}_j^{text{EB}} = \hat{R}_j\hat{u}_j^{\text{ML}}\) waarbij \(\hat{u}_j^{\text{ML}}\) de ML schattingen zijn, en \(\hat{R}_j = \frac{\sigma_alpha^2}{n_j}\) de zogenaamde Shrinkage factor is.

head(ranef(M2)$school)
      (Intercept)
20920 -10.1702110
22520 -17.0578149
22710   7.8007260
22738   0.4871012
22908  -8.1940346
23208   4.4304453

Deze waarden schatten hoeveel het intercept naar boven of beneden verschoven is in bepaalde scholen. Bijvoorbeeld, in de eerste school in de dataset is het geschatte intercept ongeveer 10.17 lager dan gemiddeld, zodat de schoolspecifieke regressielijn \((69.73 - 10.17) + 6.74 x_{ij}\) is.

Gelman&Hill (2006) karakteriseren multilevel modellering als partial pooling (ook wel shrinkage genoemd), wat een compromis is tussen twee uitersten: complete pooling, waarbij de clustering helemaal niet in het model wordt meegenomen, en no pooling, waarbij voor elke school aparte intercepten worden geschat als coëfficiënten van dummy-variabelen. De geschatte schoolspecifieke regressielijnen in het bovenstaande model zijn gebaseerd op partial pooling schattingen. Om dit aan te tonen, schatten we eerst de intercept en de helling voor elke school op drie manieren:

# Complete-pooling regression
pooled <- lm(formula = course ~ female,
             data = GCSE)
a_pooled <- coef(pooled)[1]   # complete-pooling intercept
b_pooled <- coef(pooled)[2]   # complete-pooling slope

# No-pooling regression
nopooled <- lm(formula = course ~ 0 + school + female,
               data = GCSE)
a_nopooled <- coef(nopooled)[1:J]   # 73 no-pooling intercepts              
b_nopooled <- coef(nopooled)[J+1]

# Partial pooling (multilevel) regression
a_part_pooled <- coef(M2)$school[, 1]
b_part_pooled <- coef(M2)$school[, 2]

Vervolgens plotten we de gegevens en schoolspecifieke regressielijnen voor een selectie van acht scholen met behulp van de volgende commando’s.2:

# (0) Assen plaatsen & scholen kiezen
y <- GCSE$course
x <- as.numeric(GCSE$female) - 1 + runif(N, -.05, .05)
schid <- GCSE$school
sel.sch <- c("65385",
             "68207",
             "60729",
             "67051",
             "50631",
             "60427",
             "64321",
             "68137")

# (1) Subset 8 van de scholen; genereer data frame
df <- data.frame(y, x, schid)
df8 <- subset(df, schid %in% sel.sch)

# (2) Aangeven van schattingen van volledige-pooling, geen-pooling, gedeeltelijke pooling 
df8$a_pooled <- a_pooled 
df8$b_pooled <- b_pooled
df8$a_nopooled <- a_nopooled[df8$schid]
df8$b_nopooled <- b_nopooled
df8$a_part_pooled <- a_part_pooled[df8$schid]
df8$b_part_pooled <- b_part_pooled[df8$schid]

# (3) Plot van hoe regressie fit voor de 8 scholen
ggplot(data = df8, 
       aes(x = x, y = y)) + 
  facet_wrap(facets = ~ schid, 
             ncol = 4) + 
  theme_bw() +
  geom_jitter(position = position_jitter(width = .05, 
                                         height = 0)) +
  geom_abline(aes(intercept = a_pooled, 
                  slope = b_pooled), 
              linetype = "solid", 
              color = "blue", 
              size = 0.5) +
  geom_abline(aes(intercept = a_nopooled, 
                  slope = b_nopooled), 
              linetype = "longdash", 
              color = "red", 
              size = 0.5) + 
  geom_abline(aes(intercept = a_part_pooled, 
                  slope = b_part_pooled), 
              linetype = "dotted", 
              color = "purple", 
              size = 0.7) + 
  scale_x_continuous(breaks = c(0, 1), 
                     labels = c("male", "female")) + 
  labs(title = "Schatting van Volledige-pooling, Geen-pooling en Gedeeltelijke pooling",
       x = "", 
       y = "Totale score op cursuswerk paper")+theme_bw( base_family = "serif")
De blauw-gestreepte, rood-gestreepte en paars-gestippelde lijnen geven respectievelijk de volledige pooling, no-pooling en gedeeltelijke pooling schattingen weer. We zien dat de geschatte schoolspecifieke regressielijn van de gedeeltelijke pooling-schattingen tussen de volledige pooling- en de geenpooling-regressielijn ligt. Er is meer pooling (paarse stippellijn dichter bij blauwe ononderbroken lijn) op scholen met een kleine steekproefomvang.

Figure 1: De blauw-gestreepte, rood-gestreepte en paars-gestippelde lijnen geven respectievelijk de volledige pooling, no-pooling en gedeeltelijke pooling schattingen weer. We zien dat de geschatte schoolspecifieke regressielijn van de gedeeltelijke pooling-schattingen tussen de volledige pooling- en de geenpooling-regressielijn ligt. Er is meer pooling (paarse stippellijn dichter bij blauwe ononderbroken lijn) op scholen met een kleine steekproefomvang.

Model 3: Variërend intercept en slope model met een enkele voorspeller

We breiden nu het variërende interceptmodel met één voorspeller uit om zowel het intercept als de helling willekeurig te laten variëren tussen scholen met behulp van het volgende model \[y_{ij} = \alpha_j + \beta_j x_{ij} +\epsilon_{ij},\] \[\alpha_j = \mu_\alpha + u_j,\] \[\beta_j = \mu_\beta + v_j,\] or in a reduced form as \[y_{ij} = \mu_\alpha + \mu_\beta x_{ij} + u_j + v_j x_{ij} + \epsilon_{ij}\] where \(\epsilon_{ij} \sim N(0, \sigma_{y}^{2})\) and \(\left( \begin{matrix} u_j \\ v_j \end{matrix} \right) \sim N\left( \left( \begin{matrix} 0 \\ 0 \end{matrix} \right) ,\left( \begin{matrix} { \sigma }_{ \alpha }^{ 2 } & \rho { \sigma }_{ \alpha }{ \sigma }_{ \beta } \\ \rho { \sigma }_{ \alpha }{ \sigma }_{ \beta } & { \sigma }_{ \beta }^{ 2 } \end{matrix} \right) \right)\).]:

\[y_{ij}\sim N(\alpha_{j}+\beta_{j}x_{ij} , \sigma_y ^2 ),\] \[\left( \begin{matrix} \alpha _{ j } \\ \beta _{ j } \end{matrix} \right) \sim N\left( \left( \begin{matrix} { \mu }_{ \alpha } \\ { \mu }_{ \beta } \end{matrix} \right) , \left( \begin{matrix} { \sigma }_{ \alpha }^{ 2 } & \rho { \sigma }_{ \alpha }{ \sigma }_{ \beta } \\ \rho { \sigma }_{ \alpha }{ \sigma }_{ \beta } & { \sigma }_{ \beta }^{ 2 } \end{matrix} \right) \right).\]

Merk op dat we nu variatie hebben in de \(\alpha_{j}\)’s en de \(\beta_{j}\)’s, en ook een correlatie parameter \(\rho\) tussen \(\alpha_{j}\) en \(\beta_{j}\). Dit model kan gefit worden met gebruik van lmer() en wel als volgt:

M3 <- lmer(formula = course ~ 1 + female + (1 + female | school), 
           data = GCSE, 
           REML = FALSE)
summary(M3) 
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: course ~ 1 + female + (1 + female | school)
   Data: GCSE

     AIC      BIC   logLik deviance df.resid 
 13983.4  14016.2  -6985.7  13971.4     1719 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.6886 -0.5222  0.1261  0.6529  2.6729 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 school   (Intercept) 102.93   10.146        
          femaleF      47.94    6.924   -0.52
 Residual             169.79   13.030        
Number of obs: 1725, groups:  school, 73

Fixed effects:
            Estimate Std. Error t value
(Intercept)   69.425      1.352  51.338
femaleF        7.128      1.131   6.302

Correlation of Fixed Effects:
        (Intr)
femaleF -0.574

In dit model wordt de residuele standaardafwijking binnen de school geschat als \(\hat{\sigma}_{y}=\) 13.03. De geschatte standaardafwijkingen van de schoolintercepten en de schoolhellingen zijn respectievelijk \(\hat{\sigma}_{\alpha}= 10.15\) en \(\hat{\sigma}_{\beta}= 6.92\). De geschatte correlatie tussen variërende intercepts en hellingen is \(\hat{\rho} = -0.52\). We kunnen een soortgelijke code als die in paragraaf 2.2 gebruiken om de gegevens en de schoolspecifieke regressielijnen voor een selectie van acht scholen te plotten.

Bayesiaanse inferentie voor Model 1

We kunnen snel en eenvoudig veel multilevel modellen fitten met behulp van de lmer() functie in R. Zoals eerder vermeld, zijn functies zoals lmer() gebaseerd op een combinatie van maximale waarschijnlijkheid (ML) schatting van de model parameters en empirische Bayes (EB) voorspellingen van de variërende intercepts en/of hellingen. In sommige gevallen, wanneer het aantal groepen klein is of wanneer het model veel variërende coëfficiënten of niet-nested componenten bevat, kan de ML-benadering echter minder goed werken, ten dele omdat er misschien niet genoeg informatie is om de variantieparameters nauwkeurig te schatten. In dergelijke gevallen levert de Restricted Maximum Likelihood (REML) Estimation redelijker gevolgtrekkingen op.

Een volledig Bayesiaanse benadering levert ook in deze gevallen redelijke gevolgtrekkingen op, met als bijkomend voordeel dat bij de voorspelling van de variërende intercepten en hellingen rekening wordt gehouden met alle onzekerheid in de parameterschattingen en met de daarmee samenhangende onzekerheid. Dit is een van de vele redenen waarom men geïnteresseerd zou moeten zijn in volledig Bayesiaanse schattingen. Andere redenen worden besproken in paragraaf 3.4.3. Men kan beginnen met het snel fitten van vele specificaties bij het bouwen van een model met behulp van de lmer() functie, en dan gebruik maken van de flexibiliteit van een volledig Bayesiaanse benadering met rstanarm om simulaties te verkrijgen die de onzekerheid over coëfficiënten, voorspellingen en andere grootheden die van belang zijn samenvatten.

In dit gedeelte laten we zien hoe model 1 kan worden ingepast en geëvalueerd met behulp van het pakket rstanarm. Als alternatief, maar dat komt niet in deze tutorial aan de orde, kan men ook een handgeschreven programma maken in Stan en het uitvoeren met behulp van het rstan pakket.

Gebruik van het rstanarm pakket

Veel relatief eenvoudige modellen kunnen worden aangepast met behulp van het rstanarm pakket zonder enige code te schrijven in de Stan taal. Het rstanarm pakket is een “wrapper” voor het rstan pakket waarmee de meest gebruikte regressiemodellen kunnen worden geschat met behulp van Markov Chain Monte Carlo (MCMC) en toch kunnen worden gespecificeerd met de gebruikelijke R modelleersyntaxis. Onderwijs onderzoekers kunnen Bayesiaanse schatting gebruiken voor multilevel modellen met slechts minimale veranderingen in hun bestaande code met lmer().

Bijvoorbeeld, Model 1 met standaard prior verdelingen voor \(\mu_{\alpha}\), \(\sigma_{\alpha}\), en \(\sigma_{y}\) kan worden gespecificeerd met het rstanarm pakket door stan_ toe te voegen aan de lmer aanroep:

M1_stanlmer <- stan_lmer(formula = course ~ 1 + (1 | school), 
                         data = GCSE,
                         seed = 349)

Deze stan_lmer() functie is qua syntax gelijk aan lmer(), maar in plaats van een maximum likelihood schatting uit te voeren, wordt een Bayesiaanse schatting uitgevoerd via MCMC. Omdat elke stap in de MCMC schatting random trekkingen uit de parameter ruimte inhoudt, voegen we een seed optie toe om ervoor te zorgen dat stan_lmer elke keer dat de code wordt uitgevoerd, dezelfde resultaten geeft.

Prior distributies

Model 1 is een variabel interceptiemodel met normaal verdeelde leerlingresiduen en intercepten op schoolniveau: \(y_{ij} \sim N(\alpha_{j}, \sigma_{y}^{2}),\) en \(\alpha_{j},\sim N(\mu_{alpha}, \sigma_{alpha}^{2})\). De normale verdeling voor de \(\alpha{j}\)’s kan worden beschouwd als een prioriteitsverdeling voor deze variërende intercepten. De parameters van deze prior verdeling, \(\mu_{\alpha}\) en \(\sigma_{\alpha}\), worden geschat uit de gegevens bij gebruik van maximum likelihood schatting. Bij volledige Bayesiaanse inferentie hebben alle hyperparameters (\(\mu_{\alpha}\) en \(\sigma_{\alpha}\)), samen met de andere niet-gemodelleerde parameters (in dit geval, \(\sigma_{y}\)) ook een priorverdeling nodig.

Hier gebruiken we de standaard prior verdelingen voor de hyperparameters in stan_lmer door geen prior opties op te geven in stan_lmer() functie. De standaard priors zijn bedoeld als zwak informatief in de zin dat ze gematigde regularisatie bieden [Regularisatie kan worden beschouwd als een techniek om ervoor te zorgen dat schattingen binnen een acceptabel bereik van waarden worden begrensd] en helpen bij het stabiliseren van de berekening. Opgemerkt moet worden dat de auteurs van rstanarm suggereren om niet te vertrouwen op rstanarm om de standaard prior voor een model te specificeren, maar eerder om de priors expliciet te specificeren, zelfs als ze inderdaad de huidige standaard zijn, aangezien updates van het pakket andere defaults mee kunnen krijgen.

Ten eerste wordt, alvorens rekening te houden met de schaal van de variabelen, \(\mu_{alpha}\) een normale priorverdeling gegeven met gemiddelde 0 en standaardafwijking 10. Dat wil zeggen, \(mu_{alpha} \sim N(0, 10^2)\). De standaardafwijking van deze prioriteitsverdeling, 10, is vijf keer zo groot als de standaardafwijking van de respons indien deze gestandaardiseerd zou zijn. Dit zou een dichte benadering moeten zijn van een niet-informatieve prior over het door de waarschijnlijkheid ondersteunde bereik, die in gevolgtrekkingen zou moeten geven die vergelijkbaar zijn met die verkregen met maximale waarschijnlijkheidsmethoden indien even zwakke priors worden gebruikt voor de andere parameters.

Ten tweede wordt de (ongeschaalde) prior voor \(\sigma_{y}\) ingesteld op een exponentiële verdeling met de ‘rate’-parameter op 1.

Ten derde, om een prior voor de varianties en covarianties van de variërende (of “willekeurige”) effecten te specificeren, zal rstanarm deze matrix ontbinden in een correlatiematrix van de variërende effecten en een functie van hun varianties. Omdat er in dit voorbeeld slechts één variërend effect is, reduceert de standaard (ongeschaalde) prior voor \(\sigma_{\alpha}\) die rstanarm gebruikt tot een exponentiële verdeling met de rate parameter op 1.

Ook moet worden opgemerkt dat rstanarm de priors zal schalen tenzij de autoscale = FALSE optie wordt gebruikt. Na het fitten van een model met stan_lmer, kunnen we de gebruikte priors controleren door de prior_summary() functie op te roepen.

# Een samenvatting krijgen van de priors die gebruikt worden 
prior_summary(object = M1_stanlmer)
Priors for model 'M1_stanlmer' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 73, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 73, scale = 41)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.061)

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details
# Hiermee krijgen we de SD van de uitkomst 
sd(GCSE$course, na.rm = TRUE)
[1] 16.32096

Zoals hierboven te zien is, worden de schalen van de priors voor \(\mu_{\alpha}\) en \(\sigma_y\) op respectievelijk \(163,21\) en \(16,32\) gezet na herschaling. Aangezien de standaard prior voor het intercept normaal is met een schaalparameter van \(10\), is de herschaalde prior ook normaal maar met een schaalparameter van \(\text{scale} \times \text{SD}(y) = 10 \times 16.321= 163.21\). Aangezien de standaard prior voor \(\sigma_y\) exponentieel is met een snelheidsparameter van \(1\) (of gelijkwaardig, de schaalparameter \(\text{scale} = \frac{1}{text{rate} = 1\)), is de herschaalde prior eveneens exponentieel met een schaalparameter van \(\text{scale} \maal \text{SD}(y) = 1 maal 16,321= 16,32\).

Directe output van stan_lmer

Posterior medianen en posterior mediaan absolute deviaties

We kunnen een snelle samenvatting van de fit van Model 1 weergeven door de print methode op de volgende manier te gebruiken:

print(M1_stanlmer, digits = 2)
stan_lmer
 family:       gaussian [identity]
 formula:      course ~ 1 + (1 | school)
 observations: 1725
------
            Median MAD_SD
(Intercept) 73.68   1.14 

Auxiliary parameter(s):
      Median MAD_SD
sigma 13.82   0.24 

Error terms:
 Groups   Name        Std.Dev.
 school   (Intercept)  8.87   
 Residual             13.82   
Num. levels: school 73 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Hier is de puntschatting van \(\mu_{\alpha}\) uit stan_lmer \(73.75\) en dit komt overeen met de mediaan van de posterior trekkingen. Dit is vergelijkbaar met de ML schatting uit lmer. De puntschatting voor \(`sigma_{_alpha}\) van stan_lmer is \(8.87\), die groter is dan de ML schatting (\(8.67\)). Dit verschil kan deels komen doordat de ML benadering in lmer() geen rekening houdt met de onzekerheid in \(\mu_{\alpha}\) bij het schatten van \(\sigma_{\alpha}\). De REML benadering (\(8.75\)) in lmer() houdt, zoals eerder vermeld, wel rekening met deze onzekerheid.

Bij gebruik van stan_lmer worden standaardfouten verkregen door de mediaan absolute afwijking (MAD) van elke trekking ten opzichte van de mediaan van die trekkingen te beschouwen. Het is bekend dat ML de neiging heeft om onzekerheden te onderschatten, omdat het gebaseerd is op puntschattingen van hyperparameters. Full Bayes daarentegen propageert de onzekerheid in de hyperparameters over alle niveaus van het model en levert adequatere onzekerheidsschattingen op. Zie ook Brown e.a. (2006) voor verdere discussie.

Posterior gemiddelden, posterior standaard deviaties, 95% credible interval en Monte Carlo fouten

summary(M1_stanlmer, 
        pars = c("(Intercept)", "sigma", "Sigma[school:(Intercept),(Intercept)]"),
        probs = c(0.025, 0.975),
        digits = 2)

Model Info:
 function:     stan_lmer
 family:       gaussian [identity]
 formula:      course ~ 1 + (1 | school)
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 1725
 groups:       school (73)

Estimates:
                                        mean   sd     2.5%   97.5%
(Intercept)                            73.67   1.12  71.49  75.92 
sigma                                  13.82   0.24  13.36  14.30 
Sigma[school:(Intercept),(Intercept)]  78.67  15.55  53.46 114.10 

MCMC diagnostics
                                      mcse Rhat n_eff
(Intercept)                           0.05 1.01  572 
sigma                                 0.00 1.00 4579 
Sigma[school:(Intercept),(Intercept)] 0.61 1.00  641 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Het is de moeite waard om op te merken dat bij gebruik van de summary methode, de schatting voor de standaardafwijking \(sigma_y\) het gemiddelde is van de posterior trekkingen van de parameter. Dit in tegenstelling tot de mediaan van de posterior trekkingen die we krijgen bij gebruik van de print methode. Een voordeel van het gebruik van de mediaan is dat de schatting voor \(\sigma_y^2\) gewoon het kwadraat is van de schatting voor \(\sigma_y\) als het aantal steekproeven oneven is. Dit is niet het geval bij gebruik van het gemiddelde. In dit geval, en meer algemeen wanneer we andere functies van de parameters moeten evalueren, moeten we de posterior trekkingen rechtstreeks benaderen. Dit wordt beschreven in het volgende deel.

Onder Diagnostics, verwijzen we de lezer naar Paragraaf 5 voor meer informatie over Rhat en n_eff. De waarden onder mcse vertegenwoordigen schattingen voor de Monte Carlo standaardfouten, die de willekeurigheid vertegenwoordigen die geassocieerd is met elke MCMC schattingsrun. Dat wil zeggen, met dezelfde dataset, herhaaldelijk gebruik van een MCMC benadering om een parameter te schatten levert schattingen op met een standaardafwijking gelijk aan de Monte Carlo standaardfout.

Andere output van stan_lmer

Zoals gezegd, kunnen gebruikers er de voorkeur aan geven om direct met de posterior trekkingen te werken om schattingen van meer complexe parameters te verkrijgen. Om dit te doen, moeten gebruikers ze handmatig benaderen vanuit het stan_lmer object. We laten zien hoe dit moet in de context van het maken van vergelijkingen tussen individuele scholen.

Toegang tot de simulaties en samenvattende resultaten

Gebaseerd op de standaard instellingen, genereert stan_lmer 4 MCMC-ketens van 2.000 iteraties elk. De helft van deze iteraties in elke keten wordt gebruikt als warming-up/burn-in (om de keten te laten convergeren naar de posterior verdeling), en daarom gebruiken we slechts 1.000 steekproeven per keten. Deze door MCMC gegenereerde steekproeven worden geacht getrokken te zijn uit de posterior verdelingen van de parameters in het model. Wij kunnen deze steekproeven gebruiken voor voorspellingen, om de onzekerheid samen te vatten en ‘çredible intervals’ (geloofwaardige intervallen) te schatten voor elke functie van de parameters.

Om toegang te krijgen tot de posterior trekkingen voor alle parameters, passen we de methode as.matrix() toe op het stanreg object M1_stanlmer. Dit geeft een \(S\) bij \(P\) matrix, waarbij \(S\) de grootte is van de posterior steekproef (of gelijkwaardig, het aantal MCMC iteraties na warm-up) en \(P\) het aantal parameters/kwantiteiten. Door deze matrix te manipuleren kunnen we een matrix genereren voor de variërende intercepts \(\alpha_{j}\) en vectoren met de trekkingen voor de within standaardafwijking en de between variantie. Merk op dat om de juiste kolommen voor de parameter van belang te selecteren, het nuttig is om de kolomnamen van de matrix sims te onderzoeken.

Een meer directe benadering voor het verkrijgen van de posterior trekkingen voor specifieke parameters is gebruik te maken van de ingebouwde functionaliteit van de as.matrix methode voor stanreg objecten. Wanneer de as.matrix methode wordt toegepast op een stanreg object, kan de gebruiker ofwel een optionele karaktervector van parameternamen specificeren, of een optionele karaktervector van reguliere expressies3 om de posterior trekkingen van alleen de parameters waarin ze geïnteresseerd zijn te extraheren. Bijvoorbeeld, omdat de parameter die het totale gemiddelde representeert is gelabeld met (Intercept), kunnen we de posterior trekkingen van alleen deze parameter extraheren door de optie pars = "(Intercept)" op te nemen. En omdat de parameters die de 73 schoolfouten representeren allemaal de string b[(Intercept) school: bevatten, kunnen we alle parameters die deze string bevatten extraheren door de optie regex_pars = "b[(Intercept) school:” te gebruiken.

# Extraheren van de posterior trekkingen voor alle parameters
sims <- as.matrix(M1_stanlmer)
dim(sims)
[1] 4000   76
para_name <- colnames(sims)
para_name
 [1] "(Intercept)"                          
 [2] "b[(Intercept) school:20920]"          
 [3] "b[(Intercept) school:22520]"          
 [4] "b[(Intercept) school:22710]"          
 [5] "b[(Intercept) school:22738]"          
 [6] "b[(Intercept) school:22908]"          
 [7] "b[(Intercept) school:23208]"          
 [8] "b[(Intercept) school:25241]"          
 [9] "b[(Intercept) school:30474]"          
[10] "b[(Intercept) school:35270]"          
[11] "b[(Intercept) school:37224]"          
[12] "b[(Intercept) school:47627]"          
[13] "b[(Intercept) school:50627]"          
[14] "b[(Intercept) school:50631]"          
[15] "b[(Intercept) school:60421]"          
[16] "b[(Intercept) school:60427]"          
[17] "b[(Intercept) school:60437]"          
[18] "b[(Intercept) school:60439]"          
[19] "b[(Intercept) school:60441]"          
[20] "b[(Intercept) school:60455]"          
[21] "b[(Intercept) school:60457]"          
[22] "b[(Intercept) school:60501]"          
[23] "b[(Intercept) school:60729]"          
[24] "b[(Intercept) school:60741]"          
[25] "b[(Intercept) school:63619]"          
[26] "b[(Intercept) school:63833]"          
[27] "b[(Intercept) school:64251]"          
[28] "b[(Intercept) school:64321]"          
[29] "b[(Intercept) school:64327]"          
[30] "b[(Intercept) school:64343]"          
[31] "b[(Intercept) school:64359]"          
[32] "b[(Intercept) school:64428]"          
[33] "b[(Intercept) school:65385]"          
[34] "b[(Intercept) school:66365]"          
[35] "b[(Intercept) school:67051]"          
[36] "b[(Intercept) school:67105]"          
[37] "b[(Intercept) school:67311]"          
[38] "b[(Intercept) school:68107]"          
[39] "b[(Intercept) school:68111]"          
[40] "b[(Intercept) school:68121]"          
[41] "b[(Intercept) school:68125]"          
[42] "b[(Intercept) school:68133]"          
[43] "b[(Intercept) school:68137]"          
[44] "b[(Intercept) school:68201]"          
[45] "b[(Intercept) school:68207]"          
[46] "b[(Intercept) school:68217]"          
[47] "b[(Intercept) school:68227]"          
[48] "b[(Intercept) school:68233]"          
[49] "b[(Intercept) school:68237]"          
[50] "b[(Intercept) school:68241]"          
[51] "b[(Intercept) school:68255]"          
[52] "b[(Intercept) school:68271]"          
[53] "b[(Intercept) school:68303]"          
[54] "b[(Intercept) school:68321]"          
[55] "b[(Intercept) school:68329]"          
[56] "b[(Intercept) school:68405]"          
[57] "b[(Intercept) school:68411]"          
[58] "b[(Intercept) school:68417]"          
[59] "b[(Intercept) school:68531]"          
[60] "b[(Intercept) school:68611]"          
[61] "b[(Intercept) school:68629]"          
[62] "b[(Intercept) school:68711]"          
[63] "b[(Intercept) school:68723]"          
[64] "b[(Intercept) school:68805]"          
[65] "b[(Intercept) school:68809]"          
[66] "b[(Intercept) school:71927]"          
[67] "b[(Intercept) school:74330]"          
[68] "b[(Intercept) school:74862]"          
[69] "b[(Intercept) school:74874]"          
[70] "b[(Intercept) school:76531]"          
[71] "b[(Intercept) school:76631]"          
[72] "b[(Intercept) school:77207]"          
[73] "b[(Intercept) school:84707]"          
[74] "b[(Intercept) school:84772]"          
[75] "sigma"                                
[76] "Sigma[school:(Intercept),(Intercept)]"
# Verkrijgen van school-niveau varyiërende intercept a_j
# trekking voor het algemeen gemiddelde
mu_a_sims <- as.matrix(M1_stanlmer, 
                       pars = "(Intercept)")
# trekkingen voor 73 scholen van de school-niveua fouten 
u_sims <- as.matrix(M1_stanlmer, 
                    regex_pars = "b\\[\\(Intercept\\) school\\:")
# trekkingen van alle 73 school variërende intercepten               
a_sims <- as.numeric(mu_a_sims) + u_sims          

# Verkrijgen van sigma_y en sigma_alpha^2
# trekkingen van sigma_y
s_y_sims <- as.matrix(M1_stanlmer, 
                       pars = "sigma")
# trekkingen van sigma_alpha^2
s__alpha_sims <- as.matrix(M1_stanlmer, 
                       pars = "Sigma[school:(Intercept),(Intercept)]")

Verkrijgen van gemiddelden, standaard deviaties, medianen en 95% geloofwaardigheids intervallen

In a_sims hebben we 4.000 posterior trekkingen (van alle 4 ketens) voor de variërende intercepten \(\alpha_{j}\) van de 73 scholen opgeslagen. De eerste kolom van de matrix van 4.000 bij 73 is bijvoorbeeld een vector van 4.000 posterior simulatietrekkingen voor de variërende intercept van de eerste school (School 20920). Een kwantitatieve manier om de posterior kansverdeling van deze 4.000 schattingen voor \(1,1pha_{1}\) samen te vatten is het onderzoeken van hun quantielen.

# Computeer gemiddelde, SD, mediaan en 95% geloofwaardigheids interval van de varyiërende intercepten

# Posterior gemiddelde en SD van elke alpha
a_mean <- apply(X = a_sims,     # posterior gemiddelde
                MARGIN = 2,
                FUN = mean)
a_sd <- apply(X = a_sims,       # posterior SD
              MARGIN = 2,
              FUN = sd)

# Posterior mediaan en 95% geloofwaardigheids interval
a_quant <- apply(X = a_sims, 
                 MARGIN = 2, 
                 FUN = quantile, 
                 probs = c(0.025, 0.50, 0.975))
a_quant <- data.frame(t(a_quant))
names(a_quant) <- c("Q2.5", "Q50", "Q97.5")

# Combineer samenvattende statistieken van posterior simulatie trekkingen
a_df <- data.frame(a_mean, a_sd, a_quant)
round(head(a_df), 2)
                            a_mean a_sd  Q2.5   Q50 Q97.5
b[(Intercept) school:20920]  63.63 4.40 54.97 63.57 72.17
b[(Intercept) school:22520]  57.15 1.79 53.72 57.15 60.60
b[(Intercept) school:22710]  81.64 3.15 75.47 81.64 87.78
b[(Intercept) school:22738]  73.03 4.27 64.48 73.11 81.50
b[(Intercept) school:22908]  66.51 5.28 56.47 66.47 76.70
b[(Intercept) school:23208]  79.23 2.81 73.60 79.22 85.02

Wij kunnen een rupsplot maken om de volledige Bayes-schattingen voor de schoolafhankelijke intercepts in rangorde te tonen, samen met hun 95% credible intervallen.

# Sorteer dataframe die een geschatte alfa gemiddelde en sd voor elke school omvatten
a_df <- a_df[order(a_df$a_mean), ]
a_df$a_rank <- c(1 : dim(a_df)[1])  # een vector van de schoolranking 

# Plot school-niveau alfas posterior gemiddelde en 95% credible interval
ggplot(data = a_df, 
       aes(x = a_rank, 
           y = a_mean)) +
  geom_pointrange(aes(ymin = Q2.5, 
                      ymax = Q97.5),
                  position = position_jitter(width = 0.1, 
                                             height = 0)) + 
  geom_hline(yintercept = mean(a_df$a_mean), 
             size = 0.5, 
             col = "red") + 
  scale_x_continuous("Rank", 
                     breaks = seq(from = 0, 
                                  to = 80, 
                                  by = 5)) + 
  scale_y_continuous(expression(paste("varying intercept, ", alpha[j]))) + 
  theme_bw( base_family = "serif")

Dezelfde aanpak kan natuurlijk worden gevolgd om 95% geloofwaardige intervallen te genereren voor \(\sigma_y\) en \(\sigma_\alpha\).

Vergelijkingen maken tussen individuele scholen

Het hebben van steekproeven van alle parameters en variërende intercepten uit hun gezamenlijke posterior verdeling maakt het gemakkelijk om inferenties te trekken over functies van deze parameters.

In onderwijsonderzoek en in de onderwijspraktijk is het vaak interessant om de scholen in de data met elkaar te vergelijken. Relevante vragen zijn bijvoorbeeld (1) wat is het verschil tussen de gemiddelden van school A en school B, (2) presteert school A beter dan school B en (3) wat is de rangorde van deze scholen binnen de steekproef. Wanneer niet-Bayesiaanse methoden worden gebruikt, kunnen wij proberen dergelijke vergelijkingen te maken op basis van empirische Bayes- (of Best Linear Unbiased-) voorspellingen van de variërende intercepten. Maar het zal in het algemeen onmogelijk zijn om de onzekerheid uit te drukken voor niet-lineaire functies zoals rangschikkingen. Zie ook Goldstein en Spiegelhalter (2006) voor verdere discussie.

Hier zullen we twee scholen vergelijken als voorbeeld: Scholen 60501 (de 21ste school) en 68271 (de 51ste school). We hebben al 4.000 posterior simulatietrekkingen voor beide scholen. Om conclusies te trekken over het verschil tussen de gemiddelde scores van de twee scholen, kunnen we eenvoudigweg het verschil nemen tussen de twee vectoren van trekkingen \(\alpha_{51} - \alpha_{21}\).

# Het verschil tussen de twee schoolgemiddelden (school #21 en #51)
school_diff <- a_sims[, 21] - a_sims[, 51]

Wij kunnen de posterior verdeling van het verschil als volgt onderzoeken met beschrijvende statistieken en een histogram:

# Onderzoek verschillen van twee distributies 
mean <- mean(school_diff)
sd <- sd(school_diff)
quantile <- quantile(school_diff, probs = c(0.025, 0.50, 0.975))
quantile <- data.frame(t(quantile))
names(quantile) <- c("Q2.5", "Q50", "Q97.5")
diff_df <- data.frame(mean, sd, quantile)
round(diff_df, 2)
  mean   sd  Q2.5  Q50 Q97.5
1 5.12 4.48 -3.53 5.12 14.02
# Histogram van de verschillen 
ggplot(data = data.frame(school_diff), 
       aes(x = school_diff)) + 
  geom_histogram(color = "black", 
                 fill = "gray", 
                 binwidth = 0.75) + 
  scale_x_continuous("Score verschil tussen twee scholen: #21, #51",
                     breaks = seq(from = -20, 
                                  to = 20, 
                                  by = 10)) + 
  geom_vline(xintercept = c(mean(school_diff),
                            quantile(school_diff, 
                                     probs = c(0.025, 0.975))),
             colour = "red", 
             linetype = "longdash") + 
  geom_text(aes(5.11, 20, label = "mean = 5.11"), 
            color = "red", 
            size = 4) + 
  geom_text(aes(9, 50, label = "SD = 4.46"), 
            color = "blue", 
            size = 4) + 
  theme_bw( base_family = "serif") 

prop.table(table(a_sims[, 21] > a_sims[, 51]))

  FALSE    TRUE 
0.12425 0.87575 

Het verwachte verschil komt uit op 5,11 met een standaardafwijking van 4,46 en een grote bandbreedte van onzekerheid. Het 95% geloofwaardigheidsinterval is [-3.64, 13.66], dus we zijn er 95% zeker van dat de ware waarde van het verschil tussen de twee scholen binnen het bereik ligt, gegeven de gegevens.

We kunnen ook het deel van de tijd bepalen dat School 60501 een hoger gemiddelde heeft dan School 68271:

prop.table(table(a_sims[, 21] > a_sims[, 51]))

  FALSE    TRUE 
0.12425 0.87575 

Dit betekent dat de posterior waarschijnlijkheid dat School 60501 beter is dan School 68271 87.6% is. Elk paar scholen binnen de steekproef van scholen kan op deze manier vergeleken worden.

Bayesiaanse inferentie voor Model 2 en 3

Model 2: Een voorspeller op studentenniveau toevoegen

Onderzoekers zouden de variërende interceptmodellen kunnen uitbreiden met waargenomen verklarende variabelen op het niveau van de leerling \(x_{ij}\), in dit voorbeeld een indicatorvariabele voor vrouw. Een eenvoudig variërend interceptmodel met één voorspeller op het niveau van de leerling kan worden geschreven als \(y_{ij} \N(\alpha_{j} + \beta x_{ij}, \sigma_{y}^{2})\) en \(\alpha_{j} \N(\mu_{alpha}, \sigma_{alpha}^{2})\). We gebruiken niet-informatieve prioriteitsverdelingen voor de hyperparameters (\(\mu_{\alpha}\) en \(\sigma_{\alpha}\)) zoals gespecificeerd in het variërende interceptmodel zonder voorspellers. Bovendien krijgt de regressiecoëfficiënt \(\beta\) een normale prioriteitsverdeling met gemiddelde 0 en standaardafwijking 100. Dit betekent, ruwweg, dat we verwachten dat deze coëfficiënt in het bereik \((-100, 100)\) ligt, en als de ML schatting in dit bereik ligt, geeft de prior verdeling zeer weinig informatie voor de inferentie.

Het bovenstaande model kan als volgt worden gefit met de stan_lmer() functie in het rstanarm pakket:

M2_stanlmer <- stan_lmer(formula = course ~ female + (1 | school), 
                         data = GCSE, 
                         prior = normal(location = 0, 
                                        scale = 100,
                                        autoscale = FALSE),
                         prior_intercept = normal(location = 0, 
                                                  scale = 100, 
                                                  autoscale = FALSE),
                         seed = 349)
prior_summary(object = M2_stanlmer)
Priors for model 'M2_stanlmer' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 100)

Coefficients
 ~ normal(location = 0, scale = 100)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.061)

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details
M2_stanlmer
stan_lmer
 family:       gaussian [identity]
 formula:      course ~ female + (1 | school)
 observations: 1725
------
            Median MAD_SD
(Intercept) 69.7    1.3  
femaleF      6.7    0.7  

Auxiliary parameter(s):
      Median MAD_SD
sigma 13.4    0.2  

Error terms:
 Groups   Name        Std.Dev.
 school   (Intercept)  9      
 Residual             13      
Num. levels: school 73 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Merk op dat in plaats van de standaard priors in stan_lmer, \(\mu_{\alpha}\) en \(\beta\) normale prior verdelingen krijgen met gemiddelde 0 en standaardafwijking 100 door de argumenten prior en prior_intercept op te geven als normal(location = 0, scale = 100, autoscale = FALSE). Om te voorkomen dat stan_lmer de prior schaalt, moeten we ervoor zorgen dat het argument autoscale = FALSE wordt toegevoegd.

De puntschattingen van \(_mu_{alpha}\), \(\beta\), en \(\sigma_{y}\) zijn bijna identiek aan de ML-schattingen van de lmer() fit. Echter, deels omdat ML de onzekerheid over \(\mu_{alpha}\) negeert bij het schatten van \(\sigma_{alpha}\), is de Bayesiaanse schatting voor \(\sigma_{alpha}\) (\(9,0\)) groter dan de ML-schatting (\(8,8\)), net als bij model 1.

Model 3: Variërende slopes over scholen toevoegen

We also use stan_lmer to fit Model 3 using the command below. Note that here, we use the default priors which are mostly similar to what was done in Model 1. Additionally, we are also required to specify a prior for the covariance matrix \(\Sigma\) for \(\alpha_j\) and \(\beta_j\) in this Model. stan_lmer decomposes this covariance matrix (up to a factor of \(\sigma_y\)) into (i) a correlation matrix \(R\) and (ii) a matrix of variances \(V\), and assigns them separate priors as shown below.

\[ \begin{aligned} \Sigma &= \left(\begin{matrix} \sigma_\alpha^2 & \rho\sigma_\alpha \sigma_\beta \\ \rho\sigma_\alpha\sigma_\beta&\sigma_\beta^2 \end{matrix} \right)\\ &= \sigma_y^2\left(\begin{matrix} \sigma_\alpha^2/\sigma_y^2 & \rho\sigma_\alpha \sigma_\beta/\sigma_y^2 \\ \rho\sigma_\alpha\sigma_\beta/\sigma_y^2 & \sigma_\beta^2/\sigma_y^2 \end{matrix} \right)\\ &= \sigma_y^2\left(\begin{matrix} \sigma_\alpha/\sigma_y & 0 \\ 0&\sigma_\beta/\sigma_y \end{matrix} \right) \left(\begin{matrix} 1 & \rho\\ \rho&1 \end{matrix} \right) \left(\begin{matrix} \sigma_\alpha/\sigma_y & 0 \\ 0&\sigma_\beta/\sigma_y \end{matrix} \right)\\ &= \sigma_y^2VRV. \end{aligned} \]

De correlatiematrix \(R\) is een 2 bij 2 matrix met 1-en op de diagonaal en \(rho\)’s op de off-diagonaal. stan_lmer kent er een LKJ^[Voor meer details over de LKJ verdeling, zie hier en hier prior aan toe, met regularisatieparameter 1 (Lewandowski et all., 2009). Dit komt overeen met het toekennen van een uniforme prior voor \(rho\). Hoe groter de regularisatieparameter is dan 1, hoe meer de verdeling voor \(\rho\) de waarde 0 aanneemt.

De matrix van (geschaalde) varianties \(V\) kan eerst worden samengevat in een vector van (geschaalde) varianties, en vervolgens ontleed in drie delen, \(J\), \(\tau^2\) en \(\pi\) zoals hieronder getoond. \[ \left(\begin{matrix} \sigma_\alpha^2/\sigma_y^2 \\ \sigma_\beta^2/\sigma_y^2 \end{matrix} \right) = 2\left(\frac{\sigma_\alpha^2/\sigma_y^2 + \sigma_\beta^2/\sigma_y^2}{2}\right)\left(\begin{matrix} \frac{\sigma_\alpha^2/\sigma_y^2}{\sigma_\alpha^2/\sigma_y^2 + \sigma_\beta^2/\sigma_y^2} \\ \frac{\sigma_\beta^2/\sigma_y^2}{\sigma_\alpha^2/\sigma_y^2 + \sigma_\beta^2/\sigma_y^2} \end{matrix} \right)= J\tau^2 \pi. \]
In deze formulering is \(J\) het aantal variërende effecten in het model (hier \(J=2\)), kan \(Jtau^2\) worden beschouwd als een gemiddelde (geschaalde) variantie over de variërende effecten \(Jalpha_j\) en \(Jbeta_j\), en is \(Jpi\) een niet-negatieve vector die sommeert tot 1 (een zogenaamde simplex/probabiliteitsvector). Een symmetrische Dirichlet4 verdeling met concentratieparameter ingesteld op 1 wordt dan gebruikt als de prior voor \(\pi\). Standaard impliceert dit een gezamenlijk uniforme prior over alle simplexvectoren van dezelfde grootte. Een schaalinvariante Gamma-voorrang met vorm- en schaalparameters beide op 1 wordt dan toegekend voor \(\tau\). Dit komt overeen met het toekennen van de exponentiële verdeling met de snelheidsparameter op 1 die consistent is met de prioriteit toegekend aan \(\sigma_y\) als prioriteit.

M3_stanlmer <- stan_lmer(formula = course ~ female + (1 + female | school), 
                         data = GCSE,
                         seed = 349)
prior_summary(object = M3_stanlmer)
Priors for model 'M3_stanlmer' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 73, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 73, scale = 41)

Coefficients
  Specified prior:
    ~ normal(location = 0, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 0, scale = 83)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.061)

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details
M3_stanlmer
stan_lmer
 family:       gaussian [identity]
 formula:      course ~ female + (1 + female | school)
 observations: 1725
------
            Median MAD_SD
(Intercept) 69.4    1.3  
femaleF      7.1    1.1  

Auxiliary parameter(s):
      Median MAD_SD
sigma 13.0    0.2  

Error terms:
 Groups   Name        Std.Dev. Corr 
 school   (Intercept) 10.3          
          femaleF      7.2     -0.49
 Residual             13.0          
Num. levels: school 73 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Hier zien we dat de puntschattingen voor \(\mu_{\alpha}\) en \(\sigma_{y}\) identiek zijn aan de ML-schattingen uit lmer() fit. De puntschatting voor $\(_bèta\) is iets anders in dit model (7.14 vergeleken met 7.13). Verder is, net als bij de vorige twee modellen, de Bayesiaanse schatting voor \(\sigma_{\alpha}\) (10.3) groter dan de ML schatting (10.15). Daarnaast zijn de Bayesiaanse schattingen voor \(\sigma_{\beta}\) (7.2) en \(\rho\) (-0.49) groter dan de overeenkomstige ML schattingen (respectievelijk 6.92 en -0.52).

Evalueren van model convergentie

Standaard zullen alle rstanarm modelleerfuncties 4 willekeurig geïnitialiseerde Markov-ketens laten lopen, elk gedurende 2000 iteraties (inclusief een opwarmperiode van 1000 iteraties). Alle ketens moeten naar de doelverdeling convergeren om geldige conclusies te kunnen trekken. De diagnostica die wij gebruiken om te beoordelen of de ketens naar de posterior verdeling zijn geconvergeerd, zijn de statistieken \(\hat{R}\) en \(N_{text{eff}}\) [@gelman1992inference]. Aan elke parameter zijn de statistiek \(\hat{R}\) en \(N_{text{eff}}\) verbonden. Zoals eerder gezien, worden deze statistieken automatisch gegenereerd bij gebruik van de summary methode.

De \(\hat{R}\) is in wezen de verhouding tussen-keten variantie en binnen-keten variantie analoog aan ANOVA. De statistiek \(\hat{R}\) zou kleiner dan 1,1 moeten zijn als de ketens zijn geconvergeerd. Om een plot te zien van de \(\hat{R}\) waarden over de parameters kunnen we de plot methode gebruiken voor stanreg object M1_stanlmer:

plot(M1_stanlmer, "rhat")
plot(M1_stanlmer, "ess")

De statistiek \(N_{text{eff}}\) geeft het effectieve aantal trekkingen van de simulatie weer. Als de trekkingen onafhankelijk zouden zijn, zou \(N_{text{eff}}\) het aantal opgeslagen trekkingen zijn, 4.000 (4 ketens $ maal$ (2.000 iteraties - 1.000 iteraties voor warmup)), maar \(N_{text{eff}}\) is meestal kleiner dan 4.000 omdat Markov chain simulaties de neiging hebben autocorrelatie te vertonen. De \(N_{text{eff}}\) statistiek zou typisch minstens 100 moeten zijn voor alle parameters. We kunnen een histogram plotten van de verhouding tussen de effectieve steekproefgrootte en de totale steekproefgrootte:

plot(M1_stanlmer, "ess")

Literatuur

Browne, William J, David Draper, and others. 2006. “A Comparison of Bayesian and Likelihood-Based Methods for Fitting Multilevel Models.” Bayesian Analysis 1 (3): 473–514.

Gelman, Andrew, and Jennifer Hill. 2006. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.

Gelman, Andrew, and Donald B Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Statistical Science. JSTOR, 457–72.

Goldstein, Harvey, and David J Spiegelhalter. 1996. “League Tables and Their Limitations: Statistical Issues in Comparisons of Institutional Performance.” Journal of the Royal Statistical Society. Series A (Statistics in Society), 385–443.

Lewandowski, Daniel, Dorota Kurowicka, and Harry Joe. 2009. “Generating Random Correlation Matrices Based on Vines and Extended Onion Method.” Journal of Multivariate Analysis 100 (9). Elsevier: 1989–2001.

Rasbash, Jon, William Browne, Harvey Goldstein, Min Yang, Ian Plewis, Michael Healy, Geoff Woodhouse, David Draper, Ian Langford, and Toby Lewis. 2000. “A User’s Guide to Mlwin.” University of London, Institute of Education, Centre for Multilevel Modelling.


  1. We werken meer over prioriteitsverdelingen uit in hoofdstuk 3. 2↩︎

  2. Gebruikers die niet vertrouwd zijn met de syntaxis van ggplot2 verwijzen we naar hier↩︎

  3. Voor meer informatie over reguliere expressies, zie hier↩︎

  4. De Dirichlet-verdeling is een multivariate veralgemening van de bètaverdeling met één concentratieparameter, die kan worden geïnterpreteerd als voorafgaande tellingen van een multinomiale willekeurige variabele (de simplexvector in onze context), zie voor details hier.↩︎