RealVAMS : An R Package for Fitting a Multivariate Value-added Model ( VAM )

We present RealVAMS, an R package for fitting a generalized linear mixed model to multimembership data with partially crossed and partially nested random effects. RealVAMS utilizes a multivariate generalized linear mixed model with pseudo-likelihood approximation for fitting normally distributed continuous response(s) jointly with a binary outcome. In an educational context, the model is referred to as a multidimensional value-added model, which extends previous theory to estimate the relationships between potential teacher contributions toward different student outcomes and to allow the consideration of a binary, real-world outcome such as graduation. The simultaneous joint modeling of continuous and binary outcomes was not available prior to RealVAMS due to computational difficulties. In this paper, we discuss the multidimensional model, describe RealVAMS, and demonstrate the use of this package and its modeling options with an educational data set.


Introduction
Originally developed in Broatch and Lohr (2012), the RealVAMS model is a multivariate generalized linear mixed model (GLMM) that extends previous theory by allowing the simultaneous joint modeling of continuous and binary outcomes within a value-added framework.Although the model may be used in a variety of contexts (e.g., sports, medicine, etc.), it is particularly useful in education for studying relationships between teachers and students.For example, the RealVAMS model can estimate potential teacher contributions toward a variety of student outcomes, including quantitative scores from different subjects and different assessments, as well as categorical success outcomes such as graduation status.In addition, the RealVAMS model estimates the relationships between the predicted teacher effects for these various measures of student success.
The RealVAMS package employs the RealVAMS model to analyze multimembership data.Multimembership data (Browne et al., 2001) poses unique modeling and estimation issues.Bates et al. (2015) describes the random effects that result from this type of data structure as "partially crossed" and "partially nested" because subjects are not perfectly nested under one higher level effect and may be nested under multiple, different higher level effects.The newly released lme4 (Bates et al., 2015) adequately analyzes this type of effect, but cannot model a continuous outcome and a binary outcome at the same time.Specifically, glmer cannot utilize both a normal and binary link in a single function call.
Modeling both types of responses is essential for capturing a broader picture of the latent effects.If a relationship between the latent effects on the outcome variables exists, this simultaneous estimation improves prediction of outcomes.In an application of the RealVAMS model to the context of sports, Broatch and Karl (2018) found that the joint model benefits from significantly improved median log-loss and absolute residuals of cross-validation predictions if the estimated "team-level" random effects are correlated.However if the team-level random effects are not correlated, the joint model does not show improved predictive power.These "team effects" are analogous to the "teacher effects" in an education setting such that if different educational outcomes are somewhat correlated, the joint model will likely lead to improved predictions of academic achievement and more accurate estimates of teacher effectiveness.However, if the effects are not correlated, the joint model will likely not be an improvement over separate models that can be run in lme4.
Similar to other value-added models (VAMs) used within an educational context, the RealVAMS model aims to estimate the effects of educational factors such as teachers, schools, and districts on student learning while controlling for prior student achievement and other covariates, if available.These statistical models typically measure the correlation between different educational factors and student achievement, and do not directly measure causation (American Statistical Association, 2014).Consequently, "teacher effect" estimates represent unexplained classroom-level heterogeneity (Lockwood et al., 2007) and should be "viewed within the context of quality improvement" (American Statistical Association, 2014, p. 2) rather than used for high-stakes evaluation purposes.Although most VAMs are based on standardized test scores alone (American Statistical Association, 2014), the RealVAMS model enables users to explore the potential contributions of educational factors and programs to other types of student outcomes, providing a more nuanced characterization of student success to help enhance the quality and fidelity of educational programs and instruction.
The combination of different types of responses allows measurement of educational effectiveness using additional indicators of student success, not just standard test scores.RealVAMS also estimates The R Journal Vol.10/1, July 2018 ISSN 2073-4859 the relationships among the multidimensional random effects on all continuous and binary outcomes.Therefore, in addition to obtaining multidimensional real-world estimates of value added by teachers or programs, education researchers can simultaneously explore the relationships among the effects on typical continuous measures of student achievement and real-world outcomes, which are frequently binary.In this paper, we present the RealVAMS model and discuss the technical challenges overcome by RealVAMS when using multimembership data.We also demonstrate the use of RealVAMS (Karl et al., 2015) with an educational data set.

RealVAMS model
This section discusses the unique features of the multidimensional RealVAMS model and how the model extends current theory.

Normal multimembership model
For simplicity, we begin by introducing the RealVAMS model for the case where all responses for subject i are continuous, then extend this model to include a binary outcome.To begin, the RealVAMS model presented in a traditional mixed model framework is The vector of responses for subject i consists of t potentially different continuous responses that do not require time ordering or scaling.The [t × (p + 1)] matrix X i gives the p covariates for subject i.In an educational context, these covariates may include time-invariant covariates such as gender and ethnicity, as well as time-varying covariates such as participation in free lunch programs.The latent effect j = 1, . . ., m for responses k = 1, . . ., t is represented by the vector γ.For latent effect j, γ j = (γ j1 , . . ., γ jt ) is the t-vector of effects, where γ jk represents the latent effect of j on response k.Thus, the vector γ = [γ 1 , . . ., γ m ] is the concatenation of the individual m vectors.The (t × tm) matrix S i is the multimembership design matrix.In an educational context, this matrix indicates which teacher instructs student i for each of the t responses since the continuous response is often measured in multiple years.In the model, S can be expanded to allow for fractional membership.The error term ε i is assumed to follow a normal distribution with mean 0 and variance R i , where R i is unrestricted and allows the subject responses to be correlated over the t responses.Additionally, all γ j and ε i are assumed to be uncorrelated.The model also assumes that γ is normally distributed with mean 0 and variance G = blockdiag(G 1 , . . ., G m ) where and all G j are initially assumed to be equal.
The multimembership structure precludes a factorization of the design and covariance matrices across subjects, which makes the non-sparse estimation routines inefficient and non-scalable for both linear and nonlinear mixed models.Thus, RealVAMS uses sparse matrix routines provided by Matrix (Bates and Maechler, 2015).Broatch and Lohr (2012) use the pseudo-likelihood approach to obtain approximations to the maximum likelihood estimators; they then adopt the penalized quasi-likelihood approach in SAS PROC GLIMMIX (SAS Institute Inc., 2008) to approximate maximum likelihood estimates (Breslow and Clayton, 1993;Wolfinger and O'Connell, 1993).However, the SAS estimation in presented in Broatch and Lohr (2012) is limited to less than 30 random effects, which is not feasible or efficient for a large scale model.

Normal-binary multimembership model
The model executed by RealVAMS simultaneously analyzes normal-binary responses with correlated random effects.To illustrate, model (1) can be extended to a GLMM that also includes a binary response.In this case, the continuous response model ( 1) is adapted by defining the binary response as an unobservable latent trait ỹ.The binary response is defined as r ij = 1 if the latent variable ỹij > 0. For example, if college entry was the binary response of interest, then the response would equal 1 if the student entered college and 0 otherwise.A response r ij = 1 is therefore equivalent to the subject's underlying latent trait exceeding some threshold, ỹij > 0. Thus, ỹi = X i β + S i γ + εi . (2) The R Journal Vol.10/1, July 2018 ISSN 2073-4859 where γ ∼ N(0, G) and ε ∼ N(0, R).To maintain the identifiability in the parameter space, we take R i to be a correlation matrix fixed at 1.This requirement is only necessary for the variance of the latent process for the binary response; the variance components for the continuous responses are unrestricted.The other terms in the model are defined as in model ( 1).The GLMM contains the linear mixed model inside the inverse link function: where g(•) is a multivariate probit link function for a binary response, following the recommendation of McCulloch (1994) and Rabe-Hesketh and Skrondal (2001).
The computational difficulty in estimating the random effects and correlations between those effects increases with the inclusion of different distributional outcomes.Simultaneously modeling a binary response requires the inclusion of a nonlinear link function in the integrand of the expression for the log-likelihood.No closed-form solutions are available for these integrals.Due to the multimembership structure, the dimension of the intractable integral is equal to the number of random effects, with at least one random effect for each effect/response combination.Thus, RealVAMS uses the pseudo-likelihood (PL) method of Wolfinger and O'Connell (1993) to linearize the likelihood function.We use an EM algorithm to maximize the likelihood function, iteratively calculate the pseudo-response, and then re-run EM to maximize.PL is the estimation routine used to handle the probit link, and EM is used within each iteration for maximization instead of Newton-Raphson.This approach differs from the Bayesian estimation techniques used in Zhang et al. (2016).Specifically, we use a doubly-iterative pseudolikelihood routine: the outer pseudolikelihood procedure is given in Wolfinger and O'Connell (1993) and the inner loop follows Karl et al. (2013).For convenience, the comments within the RealVAMS source code identify specific equation numbers in Karl et al. (2013) and Wolfinger and O'Connell (1993) .
The random effects in the RealVAMS model are correlated such that the off-diagonal elements of G j are not assumed to be 0, and there is a potential for high correlation among the random effects.The covariance matrix of the random teacher effects, G, may be near (or on) the boundary of the parameter space.The Newton-Raphson (NR) methods are known to fail to converge in these situations (Demidenko, 2004).The Cholesky factorization is not an option for the RealVAMS model due to the structure imposed on G. RealVAMS presents an efficient alternative: an EM algorithm to maximize the pseudo-likelihood function.Karl et al. (2013) document how the EM algorithm leads to a more stable maxima in the presence of near-singular covariance matrices.These highly efficient estimation routines have already been demonstrated to be effective on this type of multimembership model.Consequently, the resulting RealVAMS package is capable of handling models with tens of thousands of random effects.

Application of the RealVAMS package
Although the RealVAMS model may be used in a variety of contexts in which multimembership data occur (Karl, 2012), we choose to present it within an educational context.In this context, the model is generally referred to as a multidimensional value-added model (VAM).VAMs are popular accountability methods for estimating the potential effects of educational factors such as teachers and schools on student achievement.VAMs typically estimate the potential contributions of these various factors on test scores.Through unexplained classroom-level heterogeneity, existing VAMs measure potential teacher contributions toward the concepts, skills, and abilities directly measured by student exams (Koretz, 2008;McCaffrey et al., 2003) and provide a limited view of how classrooms may more broadly contribute to student success.The addition of a binary response in the RealVAMS model helps to address this issue.
The remainder of this section discusses the three main aspects to utilizing the RealVAMS package: data requirements, model specification, and output.To illustrate its features, we present an analysis of an educational data set.

Data requirements
The RealVAMS model extends previous theory to estimate the relationships between potential teacher contributions toward different student outcomes, often exams, and to allow the consideration of a real-world binary outcome such as graduation.Therefore, RealVAMS requires that data include at least one year of a continuous response (e.g., standardized test scores) and one binary outcome coded as 0 or 1 (e.g., college entry) to create the response vector.These required responses must be specified as separate data frames.At minimum, the continuous response data frame score.datamust contain four different columns: "y" containing numeric student scores, "student" containing unique The R Journal Vol.10/1, July 2018 ISSN 2073-4859 student identification numbers, "teacher" containing unique teacher identification numbers linking students to their respective teachers, and "year" containing the numeric year or numeric time period indicator.Optionally, additional covariates may exist for modeling purposes and can be included in the score.data.Similarly, the binary response data frame outcome.datamust contain: "r" containing the binary response coded as a 0 or 1, "student," "teacher," and "year" containing the numeric year or numeric time period indicator.
To demonstrate the use of the RealVAMS package, we apply the model to examine the effect of teacher effectiveness on college entry and standardized test scores.Student and classroom data were obtained for a single cohort of 11th grade students in a large public school district.Available outcome data included standardized test scores and whether or not each student entered college after graduation.For simplicity, students missing either test scores or college enrollment status were excluded from the analysis.The final data set included 803 students whose records were linked to 71 different teachers.Just as in traditional models, student background and baseline performance scores were included in the model specification.Student information regarding their gender (SEX), ethnicity (ETH), gifted classification (GIF), special needs classification (SPE), and/or English-Language learner classification (ESL) were included.The data set also contained students' math (PLM) and reading (PLR) sub-scale scores from the PLAN assessment, a national assessment administered to 10th grade students.The PLM and PLR scores were included as a baseline measure of student achievement.
The 11th grade state test scores (y) and college entry status (r) served as the continuous and binary outcome responses, respectively.The original data were separated into two different data frames for input to the RealVAMS function: "score" data for the state test scores and "outcome" data for college entry status.The two data frames are linked by teacher ID and student ID.To further illustrate the package requirements for multi-year data, the RealVAMS help file includes a simple simulated example.The example score data frame example.score.dataincludes 1875 observations on 625 students over 3 years, with 25 different teachers in each year.The example score data frame includes the four required variables: • student represents the unique students, • teacher represents the unique teacher identification, • year is a numeric vector containing the years 1, 2 and 3, and • y represents the student score.
The data frame also contains an optional continuous covariate cont_var.Correspondingly, the example binary outcome data frame example.outcome.dataincludes 625 observations and the required two variables: r is composed of 0s and 1s representing a simulated binary outcome measured on students, and student represents the unique student identification that corresponds to the student identification used in outcome.data.

Model specification
The RealVAMS function accommodates different value-added modeling specifications for analyzing multimembership data, including the RealVAMS models (1) and ( 2).The following syntax illustrates how to analyze the joint continuous and binary outcomes using the RealVAMS function.
> res=RealVAMS(score.data,outcome.data,persistence = "CP" score.fixed.effects= formula(~as.factor(SEX)+as.factor(ETH)+as.factor(SPE)+as.factor(ESL)+as.factor(GIF)+PLM+PLR),outcome.fixed.effects=formula(~as.factor(SEX)+as.factor(ETH)+as.factor(SPE)+as.factor(ESL)+as.factor(GIF)+PLM+PLR),school.effects= FALSE, REML = TRUE, max.iter.EM = 10, outcome.family= binomial(link = "probit"), tol1 = 1e-07, max.PQL.it= 30, pconv = .Machine$double.eps*1e9, var.parm.hessian= TRUE, verbose = TRUE) Within the RealVAMS function, the persistence option designates the type of persistence coefficients that are modeled for teacher score effects.Users may choose between complete persistence "CP," as illustrated here, or variable persistence "VP" of teacher score effects.The CP model assumes that teacher effects persist undiminished into the future, whereas the VP model assumes that teacher effects in future years are multiples of the effect in the current year (Lockwood et al., 2007).The multipliers α in the VP model are called persistence parameters and are estimated in the model.In contrast, the CP model fixes the persistence parameters at 1 (Lockwood et al., 2007).The teacher binary outcome effects are modeled with complete persistence, regardless of the selection for the continuous score effects.
The RealVAMS function allows users to account for school effects using the school.effectsoption (not used in the example above).If school.effects=TRUE,correlated random school-level effects are fit in the continuous and binary response models with zero-persistence (e.g., a student's score in each year is associated with the current school attended, and their binary outcome is associated with the last school the student attended).To use the school effects option, the school ID should be included as a column labeled "schoolID" in both the score.dataand the outcome.datadata frames.
By default, the RealVAMS function fits the pseudo-response using restricted maximum likelihood (REML) estimation.If REML=FALSE, maximum likelihood estimation is used instead.In addition, users may specify the link function used for the binary outcome response by using the outcome.familyoption.
Fixed effects are included in the score.dataand outcome.datadata frames and specified in the model using score.fixed.effectsand outcome.fixed.effects.Each of these options utilize the functionality of R's formula class such that categorical variables should be wrapped in an as.factor statement.For example, the RealVAMS syntax includes fixed effects for the factors (SEX, ETH, SPE, ESL, and GIF), continuous covariates (PLM, PLR), and separate intercepts for the continuous score outcome and the binary outcome.To fit a no-intercept model, +0 may be be specified in the formula.To fit an intercept-only model, formula(∼ 1) may be used.In addition, an interaction between any two covariates may be specified using the * operator: as.factor(SEX) * PLM, or equivalently, as.factor(SEX) + PLM + as.factor(SEX):PLM.
Other options allow users to specify the efficiency of the estimation process.The model is linearized using a pseudo-likelihood approach (Wolfinger and O'Connell, 1993), and the resulting multiple membership linear mixed model is estimated via an EM algorithm (Karl et al., 2013).The argument max.iter.EM controls the maximum number of EM iterations during each pseudo-likelihood iteration, with tol1 representing the convergence tolerance for the EM algorithm during each interior pseudo-likelihood iteration.By default, convergence is declared for each interior iteration when (l k − l k−1 ) /l k < tol1, where l k is the log-likelihood at iteration k, max.PQL.itspecifies the maximum number of outer pseudo-likelihood iterations, and pconv represents the convergence criterion for these outer iterations, similar to the PCOV option in SAS PROC GLIMMIX.By default, the Hessian is used to generate covariance matrices, providing standard errors for the estimated parameters.However, setting the var.parm.hessian to FALSE reduces the run time.When verbose is TRUE, model information is printed at each iteration.

Output
There are many output options available in RealVAMS.All parameters of the RealVAMS models (1) and ( 2) are available for output, including covariance and persistence estimates where appropriate.

Fixed effects
As with a standard mixed model, the function RealVAMS returns the fixed effect estimates β for the given formulae.The fixed effect estimates and standard errors can be obtained by using summary or parameters.These effects are reported in a manner similar to a standard VAM or lme4.The RealVAMS model estimates teacher effects on multiple outcomes which results in a separate fixed effect for each of the two response variables: the continuous standardized state test scores and the binary indicator for college entry.For example, when all other variables are held constant, Hispanic American students are estimated to score on average 9.849 points lower on the state standardized test than white students (95% confidence interval [−17.274, −2.423]) and are 16.7% less likely than white students to enter college (95% confidence interval [−29%, −2%]).Note that the point estimate of −0.4333 presented in the table has been transformed from the probit scale.
These separate estimates for the fixed effects allow one to isolate the marginal effects of a covariate on each of the responses.For some covariates, the fixed effect estimates that are statistically significant for one response are also statistically significant for the other response.However, this need not be the case and may help to illuminate useful information about how a covariate's estimated conditional relationship with one response may differ from its estimated conditional relationship with another response.For example, when all other variables are held constant, students who are classified as gifted are estimated to score on average 6.78 points higher on the state standardized test than those who are not (95% confidence interval [1.93, 11.627]), but there is no statistically distinguishable difference in their probability for entering college (95% confidence interval [−17%, 3%]).In other words, holding everything else constant, a student classified as gifted is likely to perform better on the exam, but is not more likely to attend college than another student who is not classified as gifted.

Random effects
Although other value-added methods may consider the teacher effects as fixed effects, the goal of the RealVAMS model is to estimate the "teacher effects" using random effects γ and the relationship between these effects Ĝj .The presented example assumes that each teacher teaches only one year.If there is more than one year of data, the function appends "(year z)" to each teacher identification number to keep these effects separate where z is the year in which the teacher taught.Similarly, the program appends _outcome to indicate the random effect estimate for the binary outcome.The estimated random "teacher" effects (EBLUP) and their standard errors can be obtained using teach.effects:> res$teach.effects[1:4,] The R Journal Vol.These estimates show that teacher 1 in year 1 had an estimated random effect on students' standardized test scores of −19.909 (standard error = 8.697), and an estimated effect of −0.278 on the probit scale (standard error = 0.188) on the probability of students' college entry.These random effects are estimated relative to other teachers in the sample so that the teacher effects are centered on the distributional mean of zero.In this example, we assume complete persistence with persistence = "CP", but changing this argument to persistence = "VP" will give the estimated persistence parameter α in the res$persistence_parameters function output.
Figure 1a and 1b display the caterpillar plots for the estimated teacher effects for the continuous and the binary responses with 90% confidence intervals given by > plot(res) The plot function also displays the conditional residuals and predicted probability of a positive response by response level (omitted here).

Covariance matrices
A key feature of RealVAMS is that it estimates the relationship between the teacher effects for the two types of responses by calculating the covariance matrix for teacher j, Ĝj .In general, Ĝj is a blockdiagonal matrix with blocks estimated for each year.The blockdiagonal components of G j for the one year of data presented in the example are estimated as: > res$teach.covTeachers: Ĝj = 264.753.3765 3.3765 0.0691 .
In the example, the off-diagonal or [1, 2] component of each of the block matrices gives the covariance between the continuous response in [1, 1] and the binary response in [2,2].The matrix above indicates that higher teacher effect estimates on the standardized assessment tend to associate with relatively higher teacher effect estimates on the probability of college entry (r G = 0.79).The offdiagonal, ĝ12 = 3.3765, is significantly different from zero (p = 0.007), showing that there is a significant positive relationship between the teacher effects for the probability of college entry and performance on the state standardized test.The high correlation justifies the need for the multidimensional model as opposed to two separate models.
The R Journal Vol.10/1, July 2018 ISSN 2073-4859 The RealVAMS function can also provide an estimate of the intra-student relationship between the continuous response and the binary response.Let Ri denote the estimated error covariance matrix for student i for the continuous and binary responses.
The estimated relationship between the multidimensional teacher effects can also be compared to other variance component estimates, allowing investigation of the relative contributions of each random effect.For example, the estimated student covariance matrix, R i , has larger variance component estimates than G j .Fitting separate models for the responses assumes that they are independent and the correlation in both G and R is 0. One can compare the RealVAMS model to standard separate models in lme4 or SAS PROC GLIMMIX since the RealVAMS function includes an independent.responsesoption.Setting independent.responses=TRUEfixes the correlations in both R and G to 0 and should lead to results that are identical to those from running separate models, allowing for comparison of the two estimation procedures.

Conclusion
RealVAMS is an R package for fitting a multivariate value-added model with a binary outcome and a normally distributed, continuous outcome.This package overcomes computational difficulties that arise when simultaneously modeling joint binary and continuous outcomes, and may be used in a variety of contexts.In particular, RealVAMS is useful in education to estimate the relationship between the potential teacher contributions and to allow the mixture of binary and continuous outcomes in the same model.By using models that provide a broader picture of student success, these results may help enhance the quality of educational programs.

Table 1 :
student year SEX ETH SPE ESL GIF PLM PLR resp.typeteacher First six observations of score data.Note the required columns student, year, resp.type, and y.

Table 2 :
First six observations of binary outcome data.Note the required columns student, year, resp.type, and r.
To return all parameter estimates including the fixed effects, use

Table 3 :
Fixed effects estimates from a RealVAMS model for state test scores (continuous outcome) and probability of college entry (binary outcome).Note: Fixed effects estimates for the probability of college entry are reported on the probit scale.