Introduction: Instructor- or Section-level Effects in Program Evaluation
In the propensity score weight vignette, we
saw how one might evaluate a student support program using propensity
scores to estimate the causal effects of the program. The example used
for illustration involves the writing_center data from the
IRexamples package. Briefly, the data set contains
students’ first enrollment/attempt in ENG 1, a first year college level
writing course, over several years. Students may choose to receive
additional support outside of the classroom at the college’s writing
center by enrolling in a 0.5-unit writing center/conference support
course. Students may use the writing center as a study center where they
have access to computers, dictionaries, thesauri, etc. At the center,
students may also schedule 1:1 conferences with the instructor on duty,
to get assistance and/or feedback on their writing by bringing a writing
prompt or draft of their paper. Students registered for this support
course are required to spend at least 12 hours in the center, and to
have at least 2 formal conferences.
Suppose course success in ENG 1 is the outcome of interest. As described in the previous vignette, a regression model comparing writing center enrollment to non-enrollment with adjustment variables or group comparisons weighted by propensity score weights are common approaches for evaluating the impact of the writing center on course success. However, if we were to carefully consider the data sample and structure, we would recognize the natural clustering of students in sections taught by the same instructor. This clustering may give rise to correlated data within a section or among students taught by the same instructor (i.e., non-independent observations). Inferences obtained by the two aforementioned approaches may be flawed in that estimates for both the writing center impact and the associated standard errors may be inaccurate due to a violation of the independence assumption.
A natural first step that one might consider to account for clustering is to fit a Generalized Linear Model (GLM) with fixed effects for each section or each instructor (i.e., using section and/or instructor flags as adjustment variables). Although intuitive and straightforward, this approach is usually not recommended in practice when the number of clusters are large and the number of observations in each cluster are small because there is insufficient data (small sample size) to estimate the large number of parameters (number of fixed effects plus other adjustment variables). Those trained in econometrics may point to de-meaning the outcome variable and adjustment variables in linear regression / ordinary least squares to avoid estimation of the fixed effects, but the mathematics do not work out as nicely in non-linear models such as GLM’s with non-linear link functions. In the writing center example, we have many sections, but most sections have less than 25 students, as seen below.
library(IRexamples) # contains writing_center data
library(dplyr) # ease with data manipulations
# Load data
data(writing_center)
# Class size frequency
writing_center %>%
group_by(Term, Section_ID) %>%
tally %>%
ungroup %>%
rename(Class_Size=n) %>%
group_by(Class_Size) %>%
tally %>%
ungroup %>%
rename(Num_Sections=n)
## # A tibble: 20 x 2
## Class_Size Num_Sections
## <int> <int>
## 1 10 1
## 2 11 2
## 3 12 2
## 4 13 4
## 5 14 2
## 6 15 8
## 7 16 6
## 8 17 11
## 9 18 20
## 10 19 15
## 11 20 18
## 12 21 27
## 13 22 27
## 14 23 20
## 15 24 25
## 16 25 20
## 17 26 9
## 18 27 5
## 19 30 2
## 20 33 1Popular methods that account for the classroom clustering or
instructor-level effects include Generalized Linear Mixed Models (GLMM)
and Generalized Estimating Equations (GEE; also refered to as
population-averaged or marginal mean models). Both are
extensions of GLM, where the former incorporates random effects or
latent variables to capture the effect of the clusters, and the latter
assumes a correlation structure for students within a section or for
students taught by the same instructor and allows for misspecification
of the correlation structure since inference is based on robust standard
errors. In this vignette, we illustrate the use of Bayesian Multilevel
Models (or Bayesian Hierarchical Models), which in essence will be a
GLMM under a Bayesian framework. We opt for this approach because
estimation and inference under a Bayesian framework for GLMM is more
numerically stable compared to a likelihood-based frequentist framework.
Moreover, the rstanarm package
in R allows for a convenient interface for complex models estimated
using Stan, a robust statistical
computing platform specializing in Hamiltonian Monte Carlo. Additional
details on how to use rstanarm for a Multilevel Model could
be found in this
tutorial.
Example: GLM Not Accounting for Instructor-level Effects
Before we illustrate the use of Bayesian Multilevel Models, let’s
first illustrate how one would fit a GLM on course success using the
writing center flag and some adjustment variables. In particular, we
will fit a logistic regression model (binary outcome with logit link
function). For those not familiar with the writing_center
data, type ?writing_center in the R console after loading
the IRexamples package.
fit_glm <- glm(Main_Course_SuccessFlag ~ 1 + Gender + Ethnicity + FirstGen + Military + FosterYouth + DSPS + FinAid + Units_Attempted_Beg_Of_Term + GPA_Beg_Of_Term + TermUnitsAttempted + K12_Student + First_Time_College_Student + Nonresident_Tuition_Exempt + International + Nonresident + Online + WR_Center
, data=writing_center
, family=binomial(link='logit')
)In the previous code, we call the glm function to model
Main_Course_SuccessFlag as the outcome. On the right hand
side of ~ in the formula specification, 1
represents the intercept, and the remaining terms are predictors in our
model. The WR_Center variable is a flag indicating
enrollment in the writing center, and the other predictors are
adjustment variables in the model. The
family=binomial(link='logit') argument specifies a logistic
regression model. The results are stored in fit_glm, as
seen below.
# Model summary
summary(fit_glm)
##
## Call:
## glm(formula = Main_Course_SuccessFlag ~ 1 + Gender + Ethnicity +
## FirstGen + Military + FosterYouth + DSPS + FinAid + Units_Attempted_Beg_Of_Term +
## GPA_Beg_Of_Term + TermUnitsAttempted + K12_Student + First_Time_College_Student +
## Nonresident_Tuition_Exempt + International + Nonresident +
## Online + WR_Center, family = binomial(link = "logit"), data = writing_center)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2440 0.1614 0.4299 0.6654 2.3091
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.067410 0.448125 -4.613 3.96e-06 ***
## GenderFemale 0.144271 0.339512 0.425 0.670883
## GenderMale -0.160479 0.338189 -0.475 0.635126
## EthnicityAsian 0.524207 0.277177 1.891 0.058593 .
## EthnicityDecline to State -0.054581 0.364994 -0.150 0.881128
## EthnicityHispanic / Latino -0.286383 0.279436 -1.025 0.305428
## EthnicityMixed Ethnicity 0.096240 0.275141 0.350 0.726502
## EthnicityNative American -0.557699 1.589359 -0.351 0.725667
## EthnicityPacific Islander 0.352540 0.999109 0.353 0.724197
## EthnicityWhite, Non-Hispanic 0.237113 0.272097 0.871 0.383521
## FirstGen -0.089381 0.084427 -1.059 0.289745
## Military -0.218572 0.259138 -0.843 0.398971
## FosterYouth 0.160135 0.525722 0.305 0.760671
## DSPS -0.257622 0.136136 -1.892 0.058439 .
## FinAid -0.037308 0.087455 -0.427 0.669671
## Units_Attempted_Beg_Of_Term 0.004585 0.002835 1.617 0.105865
## GPA_Beg_Of_Term 0.328934 0.035889 9.165 < 2e-16 ***
## TermUnitsAttempted 0.230946 0.009936 23.243 < 2e-16 ***
## K12_Student 3.381302 1.014323 3.334 0.000857 ***
## First_Time_College_Student 0.574293 0.115087 4.990 6.04e-07 ***
## Nonresident_Tuition_Exempt 0.002794 0.206427 0.014 0.989200
## International -0.304306 0.176058 -1.728 0.083908 .
## Nonresident 0.609108 0.313673 1.942 0.052155 .
## Online -0.106318 0.131146 -0.811 0.417550
## WR_Center 0.146165 0.095495 1.531 0.125869
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5093.0 on 4726 degrees of freedom
## Residual deviance: 4038.5 on 4702 degrees of freedom
## AIC: 4088.5
##
## Number of Fisher Scoring iterations: 6
# Odds ratio
exp(coef(fit_glm))
## (Intercept) GenderFemale
## 0.1265130 1.1551972
## GenderMale EthnicityAsian
## 0.8517358 1.6891196
## EthnicityDecline to State EthnicityHispanic / Latino
## 0.9468820 0.7509747
## EthnicityMixed Ethnicity EthnicityNative American
## 1.1010228 0.5725249
## EthnicityPacific Islander EthnicityWhite, Non-Hispanic
## 1.4226770 1.2675845
## FirstGen Military
## 0.9144973 0.8036652
## FosterYouth DSPS
## 1.1736688 0.7728871
## FinAid Units_Attempted_Beg_Of_Term
## 0.9633791 1.0045957
## GPA_Beg_Of_Term TermUnitsAttempted
## 1.3894856 1.2597906
## K12_Student First_Time_College_Student
## 29.4090327 1.7758748
## Nonresident_Tuition_Exempt International
## 1.0027982 0.7376349
## Nonresident Online
## 1.8387905 0.8991390
## WR_Center
## 1.1573868
# Confidence intervals
exp(confint(fit_glm))
## 2.5 % 97.5 %
## (Intercept) 0.05344506 0.3109644
## GenderFemale 0.57568249 2.1945993
## GenderMale 0.42544121 1.6134593
## EthnicityAsian 0.96774097 2.8768103
## EthnicityDecline to State 0.46106194 1.9333211
## EthnicityHispanic / Latino 0.42848656 1.2850283
## EthnicityMixed Ethnicity 0.63320071 1.8674079
## EthnicityNative American 0.02623404 18.5547001
## EthnicityPacific Islander 0.22844502 12.9140864
## EthnicityWhite, Non-Hispanic 0.73307162 2.1363931
## FirstGen 0.77530262 1.0795289
## Military 0.48965995 1.3565830
## FosterYouth 0.44268324 3.5375731
## DSPS 0.59319758 1.0118437
## FinAid 0.81142847 1.1433340
## Units_Attempted_Beg_Of_Term 0.99907283 1.0102473
## GPA_Beg_Of_Term 1.29595582 1.4917834
## TermUnitsAttempted 1.23573581 1.2848293
## K12_Student 6.33632156 523.9105661
## First_Time_College_Student 1.41833362 2.2271716
## Nonresident_Tuition_Exempt 0.67405801 1.5161891
## International 0.52422783 1.0460131
## Nonresident 1.02646424 3.5383057
## Online 0.69681674 1.1654993
## WR_Center 0.96084736 1.3972699After adjusting for student demographic and background information, the odds ratio comparing those enrolled in the writing center to those not in the writing center is estimated to be 1.16 (95% CI: 0.96 – 1.40).
Example: Bayesian Multilevel Model Accounting for Instructor-level Effects
In Bayesian inference, besides specifying a statistical model to represent and help answer a question of interest (e.g., logistic regression model), one also specifies a prior distribution for each model parameter (e.g., the regression coefficients, the outcome variance parameter, etc.). A posterior distribution resulting from the likelihood of the statistical model and the prior distribution is used to draw statistical inference. In nearly all non-trivial scenarios, inference is done by obtaining a random sample of the model parameters using the posterior distribution. From this monte carlo sample, one could calculate the posterior mean, median, and quantiles for each parameter. The posterior mean or median could be used as an estimate of each parameter, and the 2.5% and 97.5% quantiles could be used to construct 95% credible intervals for each parameter.
The rstanarm package allows users to fit a Bayesian GLMM
or Bayesian Multilevel Model using the stan_glmer function,
as illustrated below. Due to the computational complexity of sampling
from the posterior distribution, it is not unusual for the computer to
run for 10 minutes before results are returned.
library(rstanarm) # Load package for multilevel models
# Convert to factor instead of numeric
writing_center$Instructor_ID <- factor(writing_center$Instructor_ID)
writing_center$Section_ID <- factor(writing_center$Section_ID)
# Fit
fit_bhm <- stan_glmer(Main_Course_SuccessFlag ~ 1 + Gender + Ethnicity + FirstGen + Military + FosterYouth + DSPS + FinAid + Units_Attempted_Beg_Of_Term + GPA_Beg_Of_Term + TermUnitsAttempted + K12_Student + First_Time_College_Student + Nonresident_Tuition_Exempt + International + Nonresident + Online + WR_Center + (1 | Instructor_ID)
# + (1 | Section_ID) # If we wanted to include section mixed effects
, data=writing_center
, family=binomial(link='logit')
, seed=1000
# , prior=normal(0, 10, autoscale=FALSE)
, init_r=0.1
)
##
## SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.004 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 40 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 156.267 seconds (Warm-up)
## Chain 1: 70.044 seconds (Sampling)
## Chain 1: 226.311 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.002 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 20 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 142.509 seconds (Warm-up)
## Chain 2: 73.586 seconds (Sampling)
## Chain 2: 216.095 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.002 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 20 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 238.512 seconds (Warm-up)
## Chain 3: 71.84 seconds (Sampling)
## Chain 3: 310.352 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.002 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 20 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 243.786 seconds (Warm-up)
## Chain 4: 92.942 seconds (Sampling)
## Chain 4: 336.728 seconds (Total)
## Chain 4:The interface to stan_glmer is quite similar to that of
glm. In the formula argument, we use the exact same
specification as before, but add (1 | Instructor_ID) on the
right hand side of the of ~. This specifies a random
intercept on Instructor_ID, which is a random effect at the
instructor level and accounts for correlation among students taught by
the same instructor. The user could also specify a random effect on the
section instead of the instructor by adding
(1 | Section_ID), or specify random effects for both
instructors and sections with
(1 | Instructor_ID) + (1 | Section_ID). For the purpose of
this vignette and the writing center example, we focus on
instructor-level effects, as specified by
(1 | Instructor_ID).
The seed argument ensures reproducibility of the results
as inference is based on monte carlo samples (data sampled from the
posterior distribution of the model). The default prior
argument is used, but could be modified if the user wishes to impose an
informative prior on the parameters in the model. The
init_r argument helps to start sampling from the posterior
distribution (without specifying this, errrors may arise).
Once a monte carlo sample is obtained from stan_glmer,
we could obtain a summary of the results with the following:
# Summary of all parameters
summary(fit_bhm)
##
## Model Info:
## function: stan_glmer
## family: binomial [logit]
## formula: Main_Course_SuccessFlag ~ 1 + Gender + Ethnicity + FirstGen +
## Military + FosterYouth + DSPS + FinAid + Units_Attempted_Beg_Of_Term +
## GPA_Beg_Of_Term + TermUnitsAttempted + K12_Student + First_Time_College_Student +
## Nonresident_Tuition_Exempt + International + Nonresident +
## Online + WR_Center + (1 | Instructor_ID)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 4727
## groups: Instructor_ID (52)
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) -1.9 0.5 -2.6 -2.0 -1.3
## GenderFemale 0.0 0.3 -0.4 0.0 0.5
## GenderMale -0.3 0.3 -0.7 -0.3 0.2
## EthnicityAsian 0.4 0.3 0.1 0.4 0.8
## EthnicityDecline to State -0.1 0.4 -0.6 -0.1 0.4
## EthnicityHispanic / Latino -0.3 0.3 -0.7 -0.3 0.1
## EthnicityMixed Ethnicity 0.1 0.3 -0.3 0.1 0.4
## EthnicityNative American -0.5 1.8 -2.6 -0.6 1.7
## EthnicityPacific Islander 0.4 1.1 -0.9 0.4 1.8
## EthnicityWhite, Non-Hispanic 0.2 0.3 -0.2 0.2 0.6
## FirstGen -0.1 0.1 -0.2 -0.1 0.0
## Military -0.2 0.3 -0.6 -0.2 0.1
## FosterYouth 0.2 0.6 -0.6 0.1 0.9
## DSPS -0.3 0.1 -0.5 -0.3 -0.1
## FinAid -0.1 0.1 -0.2 -0.1 0.1
## Units_Attempted_Beg_Of_Term 0.0 0.0 0.0 0.0 0.0
## GPA_Beg_Of_Term 0.3 0.0 0.3 0.3 0.4
## TermUnitsAttempted 0.2 0.0 0.2 0.2 0.3
## K12_Student 4.0 1.3 2.6 3.8 5.6
## First_Time_College_Student 0.6 0.1 0.5 0.6 0.8
## Nonresident_Tuition_Exempt 0.0 0.2 -0.3 0.0 0.2
## International -0.3 0.2 -0.6 -0.3 -0.1
## Nonresident 0.7 0.3 0.3 0.7 1.1
## Online -0.3 0.2 -0.5 -0.3 0.0
## WR_Center 0.3 0.1 0.2 0.3 0.4
## b[(Intercept) Instructor_ID:910] -0.2 0.3 -0.6 -0.2 0.1
## b[(Intercept) Instructor_ID:1210] 0.3 0.3 -0.1 0.3 0.7
## b[(Intercept) Instructor_ID:1515] -1.0 0.2 -1.2 -1.0 -0.7
## b[(Intercept) Instructor_ID:1597] 0.5 0.4 0.0 0.5 1.0
## b[(Intercept) Instructor_ID:1612] -0.2 0.2 -0.5 -0.2 0.0
## b[(Intercept) Instructor_ID:1649] 0.0 0.3 -0.3 0.0 0.3
## b[(Intercept) Instructor_ID:1658] -0.6 0.2 -0.9 -0.6 -0.3
## b[(Intercept) Instructor_ID:1660] -0.3 0.3 -0.8 -0.3 0.1
## b[(Intercept) Instructor_ID:1672] 0.1 0.3 -0.3 0.1 0.5
## b[(Intercept) Instructor_ID:1833] -0.3 0.2 -0.6 -0.3 0.0
## b[(Intercept) Instructor_ID:1900] 0.8 0.3 0.4 0.8 1.2
## b[(Intercept) Instructor_ID:1938] 0.4 0.3 0.0 0.4 0.7
## b[(Intercept) Instructor_ID:1953] 0.6 0.3 0.3 0.6 1.0
## b[(Intercept) Instructor_ID:1980] 0.0 0.4 -0.4 0.0 0.5
## b[(Intercept) Instructor_ID:2001] -0.6 0.2 -0.8 -0.6 -0.3
## b[(Intercept) Instructor_ID:2017] 0.6 0.2 0.3 0.6 0.9
## b[(Intercept) Instructor_ID:2019] 0.6 0.3 0.1 0.5 1.0
## b[(Intercept) Instructor_ID:2029] 0.5 0.3 0.1 0.5 1.0
## b[(Intercept) Instructor_ID:2034] -0.2 0.4 -0.7 -0.2 0.3
## b[(Intercept) Instructor_ID:2038] 0.3 0.3 -0.1 0.3 0.6
## b[(Intercept) Instructor_ID:2156] 0.2 0.4 -0.3 0.2 0.7
## b[(Intercept) Instructor_ID:2292] 0.1 0.4 -0.5 0.1 0.6
## b[(Intercept) Instructor_ID:2293] -0.5 0.4 -1.0 -0.5 0.0
## b[(Intercept) Instructor_ID:2321] -0.1 0.2 -0.4 -0.1 0.1
## b[(Intercept) Instructor_ID:2350] -0.6 0.4 -1.1 -0.6 -0.1
## b[(Intercept) Instructor_ID:2351] 0.2 0.3 -0.1 0.2 0.5
## b[(Intercept) Instructor_ID:2357] 0.1 0.3 -0.3 0.1 0.4
## b[(Intercept) Instructor_ID:2358] -0.5 0.4 -1.0 -0.5 0.0
## b[(Intercept) Instructor_ID:3272] 0.8 0.3 0.4 0.8 1.1
## b[(Intercept) Instructor_ID:3316] -0.4 0.3 -0.9 -0.4 0.0
## b[(Intercept) Instructor_ID:3380] 0.3 0.4 -0.3 0.3 0.8
## b[(Intercept) Instructor_ID:3598] 0.6 0.4 0.1 0.6 1.1
## b[(Intercept) Instructor_ID:3984] 0.1 0.3 -0.3 0.1 0.4
## b[(Intercept) Instructor_ID:4550] -0.3 0.3 -0.6 -0.3 0.1
## b[(Intercept) Instructor_ID:4776] -0.1 0.3 -0.5 -0.1 0.3
## b[(Intercept) Instructor_ID:5224] -0.7 0.2 -1.0 -0.7 -0.4
## b[(Intercept) Instructor_ID:6216] -0.6 0.5 -1.2 -0.6 0.0
## b[(Intercept) Instructor_ID:7882] 0.4 0.3 0.0 0.4 0.8
## b[(Intercept) Instructor_ID:8450] 0.4 0.5 -0.2 0.4 1.0
## b[(Intercept) Instructor_ID:8483] -0.1 0.2 -0.4 -0.1 0.1
## b[(Intercept) Instructor_ID:8605] -0.3 0.2 -0.6 -0.3 0.0
## b[(Intercept) Instructor_ID:8722] -0.4 0.4 -0.9 -0.4 0.1
## b[(Intercept) Instructor_ID:8735] -0.1 0.3 -0.5 -0.1 0.2
## b[(Intercept) Instructor_ID:8816] -1.3 0.3 -1.6 -1.2 -0.9
## b[(Intercept) Instructor_ID:8975] -0.6 0.2 -0.8 -0.6 -0.3
## b[(Intercept) Instructor_ID:9014] 0.4 0.3 0.0 0.4 0.8
## b[(Intercept) Instructor_ID:9102] -0.4 0.3 -0.8 -0.4 0.0
## b[(Intercept) Instructor_ID:9412] 0.3 0.2 0.0 0.3 0.6
## b[(Intercept) Instructor_ID:9569] 0.8 0.3 0.4 0.8 1.2
## b[(Intercept) Instructor_ID:9589] 0.4 0.3 0.0 0.3 0.7
## b[(Intercept) Instructor_ID:9927] 0.5 0.4 -0.1 0.5 1.0
## b[(Intercept) Instructor_ID:9935] 0.4 0.2 0.1 0.4 0.7
## Sigma[Instructor_ID:(Intercept),(Intercept)] 0.3 0.1 0.2 0.3 0.5
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.8 0.0 0.8 0.8 0.8
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 1740
## GenderFemale 0.0 1.0 2937
## GenderMale 0.0 1.0 2954
## EthnicityAsian 0.0 1.0 1188
## EthnicityDecline to State 0.0 1.0 1639
## EthnicityHispanic / Latino 0.0 1.0 1222
## EthnicityMixed Ethnicity 0.0 1.0 1170
## EthnicityNative American 0.0 1.0 5255
## EthnicityPacific Islander 0.0 1.0 4188
## EthnicityWhite, Non-Hispanic 0.0 1.0 1149
## FirstGen 0.0 1.0 6328
## Military 0.0 1.0 6401
## FosterYouth 0.0 1.0 6846
## DSPS 0.0 1.0 5501
## FinAid 0.0 1.0 5786
## Units_Attempted_Beg_Of_Term 0.0 1.0 4688
## GPA_Beg_Of_Term 0.0 1.0 3500
## TermUnitsAttempted 0.0 1.0 6135
## K12_Student 0.0 1.0 3266
## First_Time_College_Student 0.0 1.0 3393
## Nonresident_Tuition_Exempt 0.0 1.0 7035
## International 0.0 1.0 6159
## Nonresident 0.0 1.0 6218
## Online 0.0 1.0 3504
## WR_Center 0.0 1.0 4820
## b[(Intercept) Instructor_ID:910] 0.0 1.0 4503
## b[(Intercept) Instructor_ID:1210] 0.0 1.0 4997
## b[(Intercept) Instructor_ID:1515] 0.0 1.0 2714
## b[(Intercept) Instructor_ID:1597] 0.0 1.0 5064
## b[(Intercept) Instructor_ID:1612] 0.0 1.0 2940
## b[(Intercept) Instructor_ID:1649] 0.0 1.0 4047
## b[(Intercept) Instructor_ID:1658] 0.0 1.0 3558
## b[(Intercept) Instructor_ID:1660] 0.0 1.0 5576
## b[(Intercept) Instructor_ID:1672] 0.0 1.0 5086
## b[(Intercept) Instructor_ID:1833] 0.0 1.0 3120
## b[(Intercept) Instructor_ID:1900] 0.0 1.0 4697
## b[(Intercept) Instructor_ID:1938] 0.0 1.0 3602
## b[(Intercept) Instructor_ID:1953] 0.0 1.0 4839
## b[(Intercept) Instructor_ID:1980] 0.0 1.0 5923
## b[(Intercept) Instructor_ID:2001] 0.0 1.0 3518
## b[(Intercept) Instructor_ID:2017] 0.0 1.0 4116
## b[(Intercept) Instructor_ID:2019] 0.0 1.0 5321
## b[(Intercept) Instructor_ID:2029] 0.0 1.0 4585
## b[(Intercept) Instructor_ID:2034] 0.0 1.0 6708
## b[(Intercept) Instructor_ID:2038] 0.0 1.0 4102
## b[(Intercept) Instructor_ID:2156] 0.0 1.0 4735
## b[(Intercept) Instructor_ID:2292] 0.0 1.0 6324
## b[(Intercept) Instructor_ID:2293] 0.0 1.0 5567
## b[(Intercept) Instructor_ID:2321] 0.0 1.0 3456
## b[(Intercept) Instructor_ID:2350] 0.0 1.0 3600
## b[(Intercept) Instructor_ID:2351] 0.0 1.0 3842
## b[(Intercept) Instructor_ID:2357] 0.0 1.0 4300
## b[(Intercept) Instructor_ID:2358] 0.0 1.0 4845
## b[(Intercept) Instructor_ID:3272] 0.0 1.0 4234
## b[(Intercept) Instructor_ID:3316] 0.0 1.0 5221
## b[(Intercept) Instructor_ID:3380] 0.0 1.0 7226
## b[(Intercept) Instructor_ID:3598] 0.0 1.0 5157
## b[(Intercept) Instructor_ID:3984] 0.0 1.0 4942
## b[(Intercept) Instructor_ID:4550] 0.0 1.0 4254
## b[(Intercept) Instructor_ID:4776] 0.0 1.0 4796
## b[(Intercept) Instructor_ID:5224] 0.0 1.0 3771
## b[(Intercept) Instructor_ID:6216] 0.0 1.0 5737
## b[(Intercept) Instructor_ID:7882] 0.0 1.0 4378
## b[(Intercept) Instructor_ID:8450] 0.0 1.0 5625
## b[(Intercept) Instructor_ID:8483] 0.0 1.0 2983
## b[(Intercept) Instructor_ID:8605] 0.0 1.0 3646
## b[(Intercept) Instructor_ID:8722] 0.0 1.0 5696
## b[(Intercept) Instructor_ID:8735] 0.0 1.0 4644
## b[(Intercept) Instructor_ID:8816] 0.0 1.0 3706
## b[(Intercept) Instructor_ID:8975] 0.0 1.0 3220
## b[(Intercept) Instructor_ID:9014] 0.0 1.0 5408
## b[(Intercept) Instructor_ID:9102] 0.0 1.0 5009
## b[(Intercept) Instructor_ID:9412] 0.0 1.0 3965
## b[(Intercept) Instructor_ID:9569] 0.0 1.0 4778
## b[(Intercept) Instructor_ID:9589] 0.0 1.0 4091
## b[(Intercept) Instructor_ID:9927] 0.0 1.0 5488
## b[(Intercept) Instructor_ID:9935] 0.0 1.0 3706
## Sigma[Instructor_ID:(Intercept),(Intercept)] 0.0 1.0 1230
## mean_PPD 0.0 1.0 3934
## log-posterior 0.3 1.0 932
##
## 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).
# Focus on WR_Center coefficient
summary(fit_bhm
, pars = c('WR_Center')
, probs = c(0.025, 0.5, 0.975)
, digits = 2
)
##
## Model Info:
## function: stan_glmer
## family: binomial [logit]
## formula: Main_Course_SuccessFlag ~ 1 + Gender + Ethnicity + FirstGen +
## Military + FosterYouth + DSPS + FinAid + Units_Attempted_Beg_Of_Term +
## GPA_Beg_Of_Term + TermUnitsAttempted + K12_Student + First_Time_College_Student +
## Nonresident_Tuition_Exempt + International + Nonresident +
## Online + WR_Center + (1 | Instructor_ID)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 4727
## groups: Instructor_ID (52)
##
## Estimates:
## mean sd 2.5% 50% 97.5%
## WR_Center 0.30 0.11 0.09 0.30 0.51
##
## MCMC diagnostics
## mcse Rhat n_eff
## WR_Center 0.00 1.00 4820
##
## 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).In the next set of code, we extract the simulated draws and summarize them for our own needs (e.g., inference on the odds ratio comparing writing center to no writing center).
sims <- as.matrix(fit_bhm) # Extract random draws
dim(sims) # Columns are parameters, rows are observations from MCMC
## [1] 4000 78
colnames(sims) # Each column contains draws for a model parameter
## [1] "(Intercept)"
## [2] "GenderFemale"
## [3] "GenderMale"
## [4] "EthnicityAsian"
## [5] "EthnicityDecline to State"
## [6] "EthnicityHispanic / Latino"
## [7] "EthnicityMixed Ethnicity"
## [8] "EthnicityNative American"
## [9] "EthnicityPacific Islander"
## [10] "EthnicityWhite, Non-Hispanic"
## [11] "FirstGen"
## [12] "Military"
## [13] "FosterYouth"
## [14] "DSPS"
## [15] "FinAid"
## [16] "Units_Attempted_Beg_Of_Term"
## [17] "GPA_Beg_Of_Term"
## [18] "TermUnitsAttempted"
## [19] "K12_Student"
## [20] "First_Time_College_Student"
## [21] "Nonresident_Tuition_Exempt"
## [22] "International"
## [23] "Nonresident"
## [24] "Online"
## [25] "WR_Center"
## [26] "b[(Intercept) Instructor_ID:910]"
## [27] "b[(Intercept) Instructor_ID:1210]"
## [28] "b[(Intercept) Instructor_ID:1515]"
## [29] "b[(Intercept) Instructor_ID:1597]"
## [30] "b[(Intercept) Instructor_ID:1612]"
## [31] "b[(Intercept) Instructor_ID:1649]"
## [32] "b[(Intercept) Instructor_ID:1658]"
## [33] "b[(Intercept) Instructor_ID:1660]"
## [34] "b[(Intercept) Instructor_ID:1672]"
## [35] "b[(Intercept) Instructor_ID:1833]"
## [36] "b[(Intercept) Instructor_ID:1900]"
## [37] "b[(Intercept) Instructor_ID:1938]"
## [38] "b[(Intercept) Instructor_ID:1953]"
## [39] "b[(Intercept) Instructor_ID:1980]"
## [40] "b[(Intercept) Instructor_ID:2001]"
## [41] "b[(Intercept) Instructor_ID:2017]"
## [42] "b[(Intercept) Instructor_ID:2019]"
## [43] "b[(Intercept) Instructor_ID:2029]"
## [44] "b[(Intercept) Instructor_ID:2034]"
## [45] "b[(Intercept) Instructor_ID:2038]"
## [46] "b[(Intercept) Instructor_ID:2156]"
## [47] "b[(Intercept) Instructor_ID:2292]"
## [48] "b[(Intercept) Instructor_ID:2293]"
## [49] "b[(Intercept) Instructor_ID:2321]"
## [50] "b[(Intercept) Instructor_ID:2350]"
## [51] "b[(Intercept) Instructor_ID:2351]"
## [52] "b[(Intercept) Instructor_ID:2357]"
## [53] "b[(Intercept) Instructor_ID:2358]"
## [54] "b[(Intercept) Instructor_ID:3272]"
## [55] "b[(Intercept) Instructor_ID:3316]"
## [56] "b[(Intercept) Instructor_ID:3380]"
## [57] "b[(Intercept) Instructor_ID:3598]"
## [58] "b[(Intercept) Instructor_ID:3984]"
## [59] "b[(Intercept) Instructor_ID:4550]"
## [60] "b[(Intercept) Instructor_ID:4776]"
## [61] "b[(Intercept) Instructor_ID:5224]"
## [62] "b[(Intercept) Instructor_ID:6216]"
## [63] "b[(Intercept) Instructor_ID:7882]"
## [64] "b[(Intercept) Instructor_ID:8450]"
## [65] "b[(Intercept) Instructor_ID:8483]"
## [66] "b[(Intercept) Instructor_ID:8605]"
## [67] "b[(Intercept) Instructor_ID:8722]"
## [68] "b[(Intercept) Instructor_ID:8735]"
## [69] "b[(Intercept) Instructor_ID:8816]"
## [70] "b[(Intercept) Instructor_ID:8975]"
## [71] "b[(Intercept) Instructor_ID:9014]"
## [72] "b[(Intercept) Instructor_ID:9102]"
## [73] "b[(Intercept) Instructor_ID:9412]"
## [74] "b[(Intercept) Instructor_ID:9569]"
## [75] "b[(Intercept) Instructor_ID:9589]"
## [76] "b[(Intercept) Instructor_ID:9927]"
## [77] "b[(Intercept) Instructor_ID:9935]"
## [78] "Sigma[Instructor_ID:(Intercept),(Intercept)]"
# Odds ratio of WR_Center
OR_WR_Center <- exp(sims[, 'WR_Center'])
mean(OR_WR_Center) ## Mean of OR_WR_Center
## [1] 1.35359
median(OR_WR_Center) ## Median of OR_WR_Center
## [1] 1.347018
quantile(OR_WR_Center, c(0.025, 0.975)) ## 95% credible interval
## 2.5% 97.5%
## 1.099113 1.657074The above results show the odds ratio estimate and credible intervals for the impact of the writing center accounting for both the adjustment variables and the instructor-level effects. That is, after adjusting for the student background characteristics (those listed in the model) and instructor-level effects, those that enrolled in the writing center are associated with a 1.35-fold increase in the odds for course success in ENG 1 compared to those not enrolled in the writing center (95% credible interval: 1.10 – 1.66.
The following table compares the estimated odds ratios from GLM and Bayesian Multi-level Model:
| Model | Odds Ratio: WR Center vs. Not | Lower Limit | Upper Limit |
|---|---|---|---|
| GLM | 1.16 | 0.96 | 1.40 |
| Bayesian Multi-level Model | 1.35 | 1.10 | 1.66 |
Illustration: Range in Instructor Effects Based on Model
The following shows the different instructor-level intercepts based on the fitted Bayesian GLMM, which further supports the case for accounting for instructor-level effects when evaluating the writing center program.
library(ggplot2)
library(ggthemes)
# Extract simulations for intercept
mu_a_sims <- as.matrix(fit_bhm, pars='(Intercept)')
# Extract simulations for instructor-level intercepts
u_sims <- as.matrix(fit_bhm, regex_pars='\\(Intercept\\) Instructor_ID')
# Combine intercepts
a_sims <- as.numeric(mu_a_sims) + u_sims
# Store mean and quantiles as data frame
a_df <- data.frame(Varying_Intercept=apply(binomial(link='logit')$linkinv(a_sims), 2, mean)
, Q2.5=apply(binomial(link='logit')$linkinv(a_sims), 2, quantile, 0.025)
, Q97.5=apply(binomial(link='logit')$linkinv(a_sims), 2, quantile, 0.975)
) %>%
mutate(Rank=rank(Varying_Intercept))# Plot instructor-level alphas's posterior mean and 95% credible interval
ggplot(data = a_df, aes(x = Rank, y = Varying_Intercept)) +
geom_pointrange(aes(ymin = Q2.5, ymax = Q97.5), position = position_jitter(width = 0.1, height = 0)) +
geom_hline(yintercept = mean(a_df$Varying_Intercept), 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")
#dev.copy(png, filename='instructor_plot.png', width=9, height=5, units='in', res=150)
#dev.off()
References and Further Reading
- Mixed-effects modeling with R by Douglas Bates
- Data Analysis Using Regression and Multilevel/Hierarchical Models by Andrew Gelman and Jennifer Hill
Appendix: R and R Package Versions
This vignette was generated using an R session with the following packages. There may be some discrepancies when the reader replicates the code caused by version mismatch.
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggthemes_4.2.4 ggplot2_3.3.3 knitr_1.33
## [4] rstanarm_2.21.1 Rcpp_1.0.6 dplyr_1.0.6
## [7] IRexamples_0.0.1 RevoUtils_11.0.2 RevoUtilsMath_11.0.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-152 matrixStats_0.59.0 xts_0.12.1
## [4] threejs_0.3.3 rstan_2.21.2 tools_4.0.2
## [7] bslib_0.2.5.1 utf8_1.2.1 R6_2.5.0
## [10] DT_0.18 colorspace_2.0-1 withr_2.4.2
## [13] tidyselect_1.1.1 gridExtra_2.3 prettyunits_1.1.1
## [16] processx_3.5.2 curl_4.3.1 compiler_4.0.2
## [19] cli_2.5.0 shinyjs_2.0.0 labeling_0.4.2
## [22] colourpicker_1.1.0 sass_0.4.0 scales_1.1.1
## [25] dygraphs_1.1.1.6 ggridges_0.5.3 callr_3.7.0
## [28] stringr_1.4.0 digest_0.6.27 StanHeaders_2.21.0-7
## [31] minqa_1.2.4 rmarkdown_2.9 base64enc_0.1-3
## [34] pkgconfig_2.0.3 htmltools_0.5.1.1 lme4_1.1-27
## [37] highr_0.9 fastmap_1.1.0 htmlwidgets_1.5.3
## [40] rlang_0.4.11 rstudioapi_0.13 shiny_1.6.0
## [43] farver_2.1.0 jquerylib_0.1.4 generics_0.1.0
## [46] zoo_1.8-9 jsonlite_1.7.2 crosstalk_1.1.1
## [49] gtools_3.9.2 inline_0.3.19 magrittr_2.0.1
## [52] loo_2.4.1 bayesplot_1.8.1 Matrix_1.3-4
## [55] munsell_0.5.0 fansi_0.5.0 lifecycle_1.0.0
## [58] stringi_1.4.6 yaml_2.2.1 MASS_7.3-54
## [61] pkgbuild_1.2.0 plyr_1.8.6 grid_4.0.2
## [64] parallel_4.0.2 promises_1.2.0.1 crayon_1.4.1
## [67] miniUI_0.1.1.1 lattice_0.20-44 splines_4.0.2
## [70] ps_1.6.0 pillar_1.6.1 igraph_1.2.6
## [73] boot_1.3-28 markdown_1.1 shinystan_2.5.0
## [76] reshape2_1.4.4 codetools_0.2-18 stats4_4.0.2
## [79] rstantools_2.1.1 glue_1.4.2 evaluate_0.14
## [82] V8_3.4.2 RcppParallel_5.1.4 nloptr_1.2.2.2
## [85] vctrs_0.3.8 httpuv_1.6.1 gtable_0.3.0
## [88] purrr_0.3.4 xfun_0.23 mime_0.10
## [91] xtable_1.8-4 later_1.2.0 survival_3.2-11
## [94] rsconnect_0.8.18 tibble_3.1.2 shinythemes_1.2.0
## [97] ellipsis_0.3.2