1 Introduction

Health-related quality of life is a key outcome in health technology assessments because it is patient-relevant and it is needed to calculate quality-adjusted life years. Quality of life instruments typically measure health problems in multiple domains using ordinal Likert scales. Value sets or valuation functions convert these profiles of ordinal measures into cardinal health-related utilities between 1 (perfect health) and minus infinity, where 0 represents death, and negative values represent health states worse than death. Because 100% quality of life represents perfect health, health state utilities are limited at 1. The lowest possible utility in a local value set further defines a lower limit of health state utilities in a local population. Thus, health state utilities are limited dependent variables. In addition, health state utilities often show gaps between 1 and the next smaller utility in the value set. These gaps occur more frequently in quality of life instruments with few levels in the Likert scales such as the EQ-5D-3L (Mulhern et al. 2018). A last but important particularity of health state utilities is that they can be the consequence of multiple latent classes, or they can exhibit multi-modal marginal densities (Hernández Alava et al. 2014).

Adjusted limited dependent variable mixture models are finite mixtures of normal distributions that account for limits, gaps between 1 and the next smaller utility value and multi-modality (Hernández Alava, Wailoo, and Ara 2012; Hernández Alava et al. 2013, 2014; Hernández Alava and Wailoo 2015; Mukuria et al. 2019). These features can improve empirical fit, parameter identification and predictive accuracy compared to standard regression models. Thus, adjusted limited dependent variable mixture models are particularly useful for mapping studies (Gray, Wailoo, and Hernández Alava 2018; Gray, Hernández Alava, and Wailoo 2018; Dixon, Hollingworth, and Sparrow 2020; Yang et al. 2019; Xu et al. 2020; Fuller et al. 2017; Pennington et al. 2020).

The R ‘aldvmm’ package is an implementation of the adjusted limited dependent variable mixture model proposed by Hernández Alava and Wailoo (2015) using normal component distributions and a multinomial logit model of probabilities of component membership.

The objectives of this vignette are to demonstrate the usage of the ‘aldvmm’ package, show important challenges of fitting adjusted limited dependent variable mixture models and validate the R implementation against the STATA package (Hernández Alava and Wailoo 2015) using publicly available data.

2 Methods

Adjusted limited dependent variable mixture models are finite mixtures of normal distributions in \(K\) components \(c\) with conditional expectations \(E[y|X, c] = X\beta^{c}\) and standard deviations \(\sigma^{c}\). The probabilities of component membership are estimated using a multinomial logit model as \(P[c|X]=exp(X\delta^{c})/\sum_{k=1}^{K}exp(X\delta^{k})\). The model accumulates the density mass of the finite mixture below a minimum value \(\Psi_1\) at the value \(\Psi_1\), and the density mass above a maximum value \(\Psi_{2}\) at 1. If the maximum value \(\Psi_2\) is smaller than 1, the model emulates a value set with a gap between 1 and the next smaller value.

\[\begin{equation} \label{eq:limits} \begin{array}{ll} y_{i}|c =& \begin{cases} \begin{array}{ll} 1 & \text{if } y_{i}|c > \Psi_{2}\\ \Psi_{1} & \text{if } y_{i}|c \leq \Psi_{1}\\ y_{i}|c & \text{if } \Psi_{1} < y_{i}|c \leq \Psi_{2}\\ \end{array} \end{cases} \end{array} \end{equation}\]

In this vignette, we estimate the same models of post-operative EQ-5D-3L utilities as Hernández Alava and Wailoo (2015) and include post-operative Oxford Hip Scores (divided by 10) as the only explanatory variable \(x\).

\[\begin{equation} \label{eq:models} \begin{array}{lrl} \text{Model 1:}& E[y|c, X] &= \beta_{0}^{c} + \beta_{1}^{c}x\\ & P[c|X] &= \text{mlogit}(\delta_{0}^{c})\\ &&\\ \text{Model 2:}& E[y|c, X] &= \beta_{0}^{c} + \beta_{1}^{c}x\\ & P[c|X] &= \text{mlogit}(\delta_{0}^{c} + \delta_{1}^{c}x) \end{array} \end{equation}\]

The aldvmm() function fits an adjusted limited dependent variable mixture model using the likelihood function from Hernández Alava and Wailoo (2015). The function calls optimx::optimr() to minimize the negative log-likelihood using analytical gradients from aldvmm.gr(). The aldvmm() function accepts all optimization methods available in optimx::optimr() except for “nlm”, which requires a different implementation of the likelihood function. The default optimization method is “BFGS”.

The model formula in aldvmm() is an object of class “formula” with two parts on the right-hand side of ~. The first part on the left of the | delimiter represents the model of expected values of \(K\) normal distributions. The second part on the right of the | delimiter represents the model of probabilities of component membership from a multinomial logit model.

The ‘aldvmm’ package provides four options for the generation of starting values of the optimization algorithm.

1. “zero”: A vector of zeroes (default).

2. “random”: A vector of standard normal random values.

3. “constant”: Parameter estimates of a constant-only model as starting values for intercepts and standard deviations, and zeroes for all other parameters. The auxiliary models for obtaining starting values are fitted using zero starting values.

4. “sann”: Parameter estimates of a simulated annealing algorithm.

The ‘aldvmm’ package obtains fitted values using the expected value function from Hernández Alava and Wailoo (2015). Covariance matrices and standard errors of parameters are obtained from a numerical approximation of the hessian matrix using numDeriv::hessian(). Standard errors of fitted values in the estimation data \(SE^{fit}_{i}\) and standard errors of predicted values in new data \(SE^{pred}_{i}\) are calculated using the delta method (Dowd, Greene, and Norton 2014; Whitmore 1986). \(G_{i}\) denotes the gradient of a fitted value with respect to changes in parameter estimates, \(\Sigma\) denotes the covariance matrix of parameters, and \(MSE\) denotes the mean squared error of fitted versus observed values in the estimation data. \[\begin{equation} \begin{array}{rl} SE^{fit}_{i} &= \sqrt{G_{i}'\Sigma G_{i}} \end{array} \end{equation}\] \[\begin{equation} \begin{array}{rl} SE^{pred}_{i} &= \sqrt{MSE + G_{i}'\Sigma G_{i}} \end{array} \end{equation}\]

The aldvmm() function returns an object of S3 class “aldvmm” for which methods for generic functions print(), summary(), stats::predict(), stats::coef(), stats::nobs(), stats::vcov(), stats::model.matrix(), stats::formula(), stats::residuals(), stats::update() and sandwich::estfun() are available. Objects of class “aldvmm” can be supplied to sandwich::sandwich(), sandwich::vcovCL(), sandwich::vcovPL(), sandwich::vcovHAC() and sandwich::vcovBS() to estimate robust or clustered standard errors (see example in the appendix). sandwich::estfun() calls aldvmm.sc() to return analytical gradients of the log-likelihood with respect to parameters for each observation, which allows fast computation of robust standard errors. sandwich::vcovBS() allows re-estimating the covariance matrix using bootstrapping with and without clustering, which can be particularly useful in cases with no valid covariance matrix from model fitting.

3 Installation

The latest version of the ‘aldvmm’ package can be installed from cran.

install.packages("aldvmm")

4 Data

We analyze the same publicly available EQ-5D-3L utility data from English patients after hip replacement in 2011 and 2012 (NHS Digital 2013) as Hernández Alava and Wailoo (2015) in their description of the STATA ALDVMM package.

temp <- tempfile()

url <- paste0("https://files.digital.nhs.uk/publicationimport/pub11xxx/",
              "pub11359/final-proms-eng-apr11-mar12-data-pack-csv.zip")

download.file(url, temp)
rm(url)

df <- read.table(unz(description = temp,
                     filename = "Hip Replacement 1112.csv"),
                 sep = ",",
                 header = TRUE)

unlink(temp)
rm(temp)

df <- df[, c("AGEBAND", "SEX", "Q2_EQ5D_INDEX", "HR_Q2_SCORE")]
df <- df[df$AGEBAND != "*" & df$SEX != "*", ]

df$eq5d <- df$Q2_EQ5D_INDEX
df$hr <- df$HR_Q2_SCORE/10

df <- df[stats::complete.cases(df), ]

set.seed(101010101)
df <- df[sample(1:nrow(df), size = nrow(df)*0.3), ]

The data includes 35’166 observations with complete information on patients’ post-operative utilities, Oxford Hip Scores, age and sex. Like Hernández Alava and Wailoo (2015), we draw a 30% sub-sample of 10’549 observations from the population of complete observations of these variables. Although we follow a similar approach in data preparation preparation as Hernández Alava and Wailoo (2015), our random sample is not identical to the data used in their study. Post-operative EQ-5D-3L utilities from English value sets (Dolan 1997) show a bimodal distribution, limits at -0.594 and 1 and a gap between 1 and 0.883 (figure 4.1).

Figure 4.1: Frequency distribution of observed EQ-5D-3L utilities

5 Examples

5.1 Model 1: Default settings

We first fit model 1 with the default “BFGS” optimization method and “zero” initial values. The values 0.883 and -0.594 in the argument ‘psi’ represent the maximum and minimum values smaller than 1 in the English value set (Dolan 1997). As the data shows a bi-modal distribution (figure 4.1), we estimate a mixture of 2 normal distributions (‘ncmp’ = 2). aldvmm() returns an object of class “aldvmm”.

library("aldvmm")

fit <- aldvmm::aldvmm(eq5d ~ hr | 1,
                      data = df,
                      psi = c(0.883, -0.594),
                      ncmp = 2)

summary(fit)

pred <- predict(fit,
                se.fit = TRUE,
                type = "fit")

We obtain a summary table of regression results using the generic function summary(). The model converges at a log-likelihood of -706 and an Akaike information criterion value of -1’399 (table 5.1). The function aldvmm() returns the negative log-likelihood and Akaike information criterion, and smaller values suggest better fit.

The coefficients of the intercepts and covariates for the expected values \(E[y|c, X]\) of the normal distributions can be interpreted as marginal effects on component means. ‘lnsigma’ denotes the natural logarithm of the estimated standard deviation \(\sigma^{c}\). The coefficients of covariates in the multinomial logit model of probabilities of component membership are log-transformed relative probabilities. Our model only includes two components, and the multinomial logit model collapses to a binomial logit model. The intercept of -0.7281 means that the average probability of an observation in the data to belong to component 1 is \(\text{exp}(-0.7281)\) or 0.4828255 times the probability to belong to component 2.

Table 5.1: Regression results from model 1 with “BFGS” optimization method and “zero” starting values

		Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
E[y|X, c]
Comp1	(Intercept)	-0.4312	0.0224	-19.2471	0	-0.4751	-0.3873
	hr	0.3136	0.0065	47.9444	0	0.3008	0.3264
	lnsigma	-1.2479	0.0215	-57.9753	0	-1.2901	-1.2057
Comp2	(Intercept)	0.2358	0.0069	34.1725	0	0.2223	0.2493
	hr	0.1459	0.0019	76.359	0	0.1422	0.1497
	lnsigma	-2.4623	0.0178	-138.1184	0	-2.4972	-2.4273
P[c|X]
Comp1	(Intercept)	-0.7281	0.0607	-12.0038	0	-0.847	-0.6092
N = 10549	ll = -706	AIC = -1399	BIC = -1348

We obtain expected values of observations in the estimation data using the generic function predict(). Standard errors of fitted (estimation data) or predicted (new data) values are calculated using the delta method. Expected values exhibit a smoother distribution than observed values and do not show a gap between 1 and 0.883, because they are weighted averages of component distributions and 1 (figure 5.1).

Figure 5.1: Expected values from base case model

Expected values of observations in the estimation data can also be used to calculate average incremental effects or average treatment effects on the treated. Average treatment effects on the treated compare predictions for treated individuals to predictions for the same individuals without the effect of the treatment indicator. Standard errors of average treatment effects on the treated can be calculated using the delta method (see example code in the appendix).

A visual inspection of mean residuals over deciles of expected values (see example code of the modified Hosmer-Lemeshow test in the appendix) shows that model 1 fits the data poorly (figure 5.2). Although the test is overpowered with large data, and confidence bands should not be interpreted directly, the absolute over-prediction among individuals with low expected values is quite large.

Figure 5.2: Mean residuals over deciles of expected values, BFGS with zero starting values

We can use the ‘sandwich’ package (Zeileis 2006) to calculate robust or clustered standard errors from objects of class “aldvmm” (see example code in the appendix).

5.2 Model 1: Comparison of optimization methods

Hernández Alava and Wailoo (2015) suggested that the likelihood function of the adjusted limited dependent variable mixture model with the English EQ-5D-3L data might have multiple local optima, and that the estimation is sensitive to initial values. We thus fit model 1 with all optimization algorithms and methods for generating initial values available in aldvmm() to assess the sensitivity of model fits to optimization settings and to find the maximum likelihood estimates.

init.method <- c("zero", "random", "constant", "sann")

optim.method <- c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "nlminb", "Rcgmin", 
                  "Rvmmin","hjn")

fit1_all <- list()

for (i in init.method) {
  for (j in optim.method) {
    set.seed(101010101) # Seed for random starting values
    fit1_all[[i]][[j]] <- aldvmm::aldvmm(eq5d ~ hr | 1,
                                         data         = df,
                                         psi          = c(0.883, -0.594),
                                         ncmp         = 2,
                                         init.method  = i,
                                         optim.method = j)
  }
}

The maximum likelihood varies considerably across optimization methods and initial values which confirms the sensitivity of the model to changes in these settings (table 5.2). The most frequent log-likelihood is 706.32, but the Hooke and Jeeves Pattern Search Optimization (“hjn”) with “zero” initial values converges at a log-likelihood 4’295.62.

The optimization methods “Nelder-Mead”, “CG”, “L-BFGS-B”, “Rcgmmin” and “hjn” are particularly sensitive to changes of starting values. Only the method “nlminb” converges at the same log-likelihood of 706.32 regardless of initial values. However, even the “nlminb” algorithm converges at different parameters with starting values from a constant-only model than with the other starting values despite the same maximum log-likelihood.

Table 5.2: Log-likelihood by optimization method

	Nelder-Mead	BFGS	CG	L-BFGS-B	nlminb	Rcgmin	Rvmmin	hjn
zero	-245.7	706.3	-634.8	-4513.3	706.3	-634.8	706.3	4295.6
random	213.9	-616.8	696.9	-2830.2	706.3	706.3	-609.6	-627.8
constant	706.3	706.3	-634.8	-634.8	706.3	-634.8	706.3	1703.3
sann	706.3	706.3	701.4	706.3	706.3	706.3	706.3	706.3

The computation times of optimization routines were relatively short and did not vary considerably across methods expect for the “hjn” algorithm and simulated annealing “sann” starting values which requires more computation time (table 5.3).

Table 5.3: Estimation time [minutes] by optimization method

	Nelder-Mead	BFGS	CG	L-BFGS-B	nlminb	Rcgmin	Rvmmin	hjn
zero	0.43	0.11	0.18	0.06	0.11	1.22	0.10	2.48
random	0.39	0.17	2.24	0.07	0.11	1.36	0.28	0.41
constant	0.43	0.15	0.23	0.07	0.15	1.22	0.12	2.68
sann	2.54	2.36	4.50	2.47	2.38	3.57	2.36	4.23

Parameter estimates differ considerably across selected optimization algorithms (table 5.4). The default “BFGS” - “zero” approach converges at the same parameters as the “nlminb” - “zero” method. The solution of the “hjn” method is rather extreme with no effect of the Oxford Hip Score, a very small standard deviation in component 2 and a very low probability of membership of component 2. Also, the Modified Hosmer-Lemeshow test did not show better fit of the “hjn” solution compared to the other algorithms (figure 7.2 in the appendix). The solution of the “hjn” method suggests that a single-component model should be explored as well.

Table 5.4: Regression results of model 1 with zero starting values in BFGS, Nelder-Mead, nlminb and hjn algorithms

		BFGS	Nelder-Mead	nlminb	hjn
E[y|X, c]
Comp1	(Intercept)	-0.4312	-0.0426	-0.4307	0.5263
	hr	0.3136	0.2194	0.3135	0.0881
	lnsigma	-1.2479	-1.8595	-1.2481	-1.2898
Comp2	(Intercept)	0.2358	4.1853	0.2358	0.691
	hr	0.1459	-1.0081	0.1459	0
	lnsigma	-2.4623	-0.0249	-2.4624	-11.5422
P[c|X]
Comp1	(Intercept)	-0.7281	3.5636	-0.728	2.1902
ll		-706.32	245.69	-706.32	-4295.62

As adjusted limited dependent variable mixture models are frequently used for tasks that rely on predictions, we also compare expected values from the “BFGS”, “Nelder-Mead” and “hjn” methods to expected values from the “nlminb” method. Expected values from the “Nelder-Mead” method are smaller, and values from the “hjn” method are larger than “nlminb” values among observations with the lowest 80% of expected values. Predictions from “BFGS” match predictions from the “nlminb” model fit closely (figure 5.3).

Figure 5.3: Deviation of expected values from BFGS, Nelder-Mead and hjn versus nlminb with zero starting values

Based on the comparison of log-likelihoods, parameter estimates and predicted values, we deem the “nlminb” optimization method with “zero” starting values the most robust approach for the used data and model.

init <- c(0,    0,   0,   0,    0,    0,    0.7283)
lo   <- c(-Inf, -Inf, -3,  -Inf, -Inf, -3, -Inf)
hi   <- c(Inf,  Inf,  Inf, Inf,  Inf,  Inf,  Inf)

fit1_cstr <- aldvmm::aldvmm(eq5d ~ hr | 1,
                            data = df,
                            psi = c(0.883, -0.594),
                            ncmp = 2,
                            init.est = init,
                            init.lo = lo,
                            init.hi = hi)

summary(fit1_cstr)

fit <- aldvmm::aldvmm(eq5d ~ hr,
                      data = df,
                      psi = c(0.883, -0.594),
                      ncmp = 1,
                      init.method = "zero",
                      optim.method = "nlminb")

summary(fit)

init <- c(-.40293118, .30502755, .22614716, .14801581, -.70755741, 0, 
          -1.2632051, -2.4541401)

fit2 <- aldvmm::aldvmm(eq5d ~ hr | hr,
                       data = df,
                       psi = c(0.883, -0.594),
                       ncmp = 2,
                       init.est = init,
                       optim.method = "nlminb")

summary(fit2)

6 Comparison to STATA results

To validate the R implementation of adjusted limited dependent variable mixture models, we estimate the four models presented in Hernández Alava and Wailoo (2015) as reference cases in R and STATA. The STATA and R code for model estimation is included in the appendix.

1. Model 1 with default options

2. Model 1 with parameter constraints

3. Model 1 with initial values from constant-only model

4. Model 2 with user-defined initial values

The parameter estimates and standard errors obtained in R are very similar to the results from STATA (table 6.1 and table 6.2). R did not obtain any standard errors in reference model 2 while STATA returns standard errors for the first component and the probability of belonging to component 1. Although reference models 1 and 3 yield different parameter estimates, they converge at the same log-likelihood which further supports the hypothesis of multiple local optima of the likelihood. The relative ordering of models is consistent across platforms.

Fitted values show very similar marginal distributions on both platforms (figure 6.1). R does not return fitted values in reference case 2. The summary statistics of differences in fitted values between R and STATA suggest that individual predictions are quite similar across platforms as well (table 6.3).

Standard errors of fitted values differ visibly between platforms (figure 6.2 and table 6.4). The difference is particularly pronounced in reference case 1, but the standard errors from STATA seem quite extreme compared to all other reference cases. R does not return standard errors of fitted values in reference case 2.

Figure 6.1: Fitted values in R and STATA

Figure 6.2: Standard errors of fitted values in R and STATA

The comparison of the R and STATA packages shows that the R implementation sometimes behaves differently than the STATA package, but the results are not indicative of technical errors in the R implementation.

7 Discussion

Adjusted limited dependent variable mixture models are powerful tools for regression analysis of health state utilities. Unlike standard regression models, adjusted limited dependent variable mixture models account for limits, gaps and multi-modal distributions.

The comparison of different optimization methods with EQ-5D-3L utility data from English patients after hip replacement in 2011 and 2012 (NHS Digital 2013) showed that the likelihood function can be challenging to maximize and can converge at local optima or extreme solutions. Parameter estimates varied considerably across optimization methods and even across methods with the same maximum log-likelihood. However, fitted values were very similar across the four reference cases which suggests that the model is more robust for the identification of incremental and average marginal effects than for parameter identification.

The ‘aldvmm’ package offers a variety of optimization algorithms and methods for generating initial values which is an important strength in such challenging situations. It is essential to assess different optimization algorithms and methods for initial values before interpreting the parameter estimates or predictions of adjusted limited dependent variable mixture models.

The analysis of the EQ-5D-3L utility data also suggests that simpler models with fewer components should be considered when multi-component models are difficult to fit. Even single-component adjusted limited dependent variable mixture models can improve fit compared to traditional regression techniques because they account for limits and gaps.

Although coefficients can be interpreted as marginal effects within each component, they cannot be interpreted in terms of overall expected values. Thus, average marginal effects and average treatment effects need to be calculated from predictions using the generic function predict(). Standard errors of marginal effects or average treatment effects can be calculated using the delta method.

When no valid covariance matrix of parameters can be obtained or the user is in doubt of the validity of standard errors, the function sandwich::vcovBS() can be used to obtain bootstrapped standard errors. Due to the large variability of results across model fits, a large number of iterations might be needed which increases computation time. Other functions from the ‘sandwich’ package can also be used to estimate robust or clustered standard errors.

In situations with repeated utility measures, the ‘aldvmm’ package only allows fixed effects estimations with group and time fixed effects which can be an important limitation in the analysis of clinical data. However, time fixed effects can be an appropriate modeling strategy in the presence of general time trends and dynamic selection in the population, e.g. because health state utilities decrease over time and treated individuals survive longer and thus are over-represented in later measurements.

Possible extensions of the ‘aldvmm’ package could include additional component distributions, a mixed model implementation for repeated measures and an implementation of functions for calculating average marginal effects and their standard errors.