|Home | About | Journals | Submit | Contact Us | Français|
Mixed-effects logistic regression models are described for analysis of longitudinal ordinal outcomes, where observations are observed clustered within subjects. Random effects are included in the model to account for the correlation of the clustered observations. Typically, the error variance and the variance of the random effects are considered to be homogeneous. These variance terms characterize the within-subjects (i.e., error variance) and between-subjects (i.e., random-effects variance) variation in the data. In this article, we describe how covariates can influence these variances, and also extend the standard logistic mixed model by adding a subject-level random effect to the within-subject variance specification. This permits subjects to have influence on the mean, or location, and variability, or (square of the) scale, of their responses. Additionally, we allow the random effects to be correlated. We illustrate application of these models for ordinal data using Ecological Momentary Assessment (EMA) data, or intensive longitudinal data, from an adolescent smoking study. These mixed-effects ordinal location scale models have useful applications in mental health research where outcomes are often ordinal and there is interest in subject heterogeneity, both between- and within-subjects.
The ordinal logistic regression model, described as the proportional odds model by McCullagh (1980), is a popular model for analyzing ordinal outcomes. For multilevel data, where observations are nested within clusters (e.g., classes, schools, clinics) or are repeatedly assessed across time, mixed-effects regression models are often used to account for the dependency inherent in the data (Hedeker and Gibbons, 2006), and several authors have developed mixed-effects models for ordinal outcome data. Ezzet and Whitehead (1991) and Agresti and Lang (1993) describe random-intercepts proportional odds models. Hedeker and Gibbons (1994) describe both an ordinal logistic and probit model with multiple random effects. Tutz and Hennevogl (1996) propose similar mixed models that additionally allow the model thresholds to be considered as random effects. Besides use of the logit link function, other authors have developed ordinal mixed models utilizing the probit link (Harville and Mee, 1984; Jansen, 1990; Mehta et al., 2004; Saei et al., 1996) and complementary log-log link (Ten Have, 1996).
Most of the models for ordinal outcomes referenced above include the proportional odds assumption (or its equivalent under the probit or complementary log-log link function) for model covariates. For an ordinal response with C categories, this assumption states that the effect of an explanatory variable is the same across the C − 1 cumulative logits of the model, or proportional across the cumulative odds. In previous papers (Hedeker and Mermelstein, 1998, 2000), we have described an extension to the mixed proportional odds model to allow for non-proportional odds for a subset of the explanatory variables. A similar extension is described in Saei and McGilchrist (1998), who allow non-proportional time effects in panel studies. These developments follow the extension due to Peterson and Harrell (1990) of the fixed-effects proportional odds model. In this model, explanatory variables are allowed to have varying effects on the C − 1 cumulative logits. Thus, for a particular explanatory variable, C − 1 regression coefficients are estimated. These additional parameters reflect different location effects of the explanatory variables. This extended model has recently been applied succesfully in several articles (Wakefield et al., 2001; Xie et al., 2001; Freels et al., 2002; Fielding et al., 2003), and a similar Bayesian multilevel model is described in Ishwaran (2000). Fielding et al. (2003) additionally allow the random-effect parameters to have non-proportional effects.
A somewhat different extension of the proportional odds model is described by Tosteson and Begg (1988). Here, in the context of receiver operating characteristic (ROC) analysis, the scale of the effects of explanatory variables is allowed to vary. In other words, the underlying variance of the logistic distribution can vary as a function of model covariates. McCullagh and Nelder (1989) refer to this extended model for ordinal data as a generalized “rational” model. Toledano and Gatsonis (1996) use this extension in describing generalized estimating equations (GEE) analysis of correlated ROC data, while Ishwaran and Gatsonis (2000) build upon this approach using Bayesian methods.
For cross-sectional data, Cox (1995) brought together these extensions of the proportional odds model into what he termed location-scale cumulative odds models. Hedeker et al. (2006) built upon this approach within a mixed model framework for longitudinal ordinal data. The inclusion of scale parameters within the mixed model is particularly advantageous because it allows modeling of both the within-subjects (WS) and between-subjects (BS) variances. In this regard, Hedeker et al. (2008) described a mixed model for variance modeling of continuous longitudinal data that also included a random subject effect to the WS variance model. Here, we extend this to longitudinal ordinal data. Specifically, our model will include a log-linear structure for both the WS and BS variance, allowing covariates to influence both sources of variation. Also, as in Hedeker et al. (2008), a random subject effect is included in the WS variance specification to allow the WS variance to vary at the subject level, above and beyond the influence of covariates on this variance. Data from an adolescent EMA smoking study are used to illustrate the mixed ordinal location-scale model.
The article is organized as follows. Section 2 describes data from Ecological Momentary Assessment (EMA) procedures and lists some pertinent mental health and smoking studies that have employed this approach to data collection. Section 3 presents details of the EMA study that we will use to illustrate our proposed mixed ordinal location-scale model. Section 4 presents the model in detail. Estimation aspects are described in Section 5. Application of our model to the smoking EMA data are presented in Section 6. Finally, in Section 7, we discuss and summarize features of the model and our application.
Modern data collection procedures, such as ecological momentary assessments (EMA, (Stone and Shiffman, 1994; Smyth and Stone, 2003)), experience sampling (de Vries, 1992; Scollon et al., 2003; Feldman Barrett and Barrett, 2001), and diary methods (Bolger et al., 2003), have been developed to record the momentary events and experiences of subjects in daily life (Bolger et al., 2003). These procedures yield relatively large numbers of subjects and observations per subject, and data from such designs are sometimes referred to as intensive longitudinal data (Walls and Schafer, 2006). Such designs are in keeping with the “bursts of measurement” approach described by Nesselroade and McCollam (2000), who called for such an approach in order to assess intra-individual variability. As noted by Nesselroade and McCollam (2000), such bursts of measurement increase the research burden in several ways; however, they are necessary for studying intra-individual variation and to explain why subjects differ in variability rather than solely in mean level (Bolger et al., 2003). In this article we describe data from a natural history study of adolescent smoking, using EMA, where interest was on determinants of the variation in the adolescents’ moods.
In mental health research, EMA methods have been used in studying pediatric affective disorders (Axelson et al., 2003), eating disorders (Boseck et al., 2007; le Grange et al., 2002), drug abuse (Epstein et al., 2009), schizophrenia (Granholm et al., 2008; Kimhy et al., 2006), borderline personality disorder (Trull et al., 2008), stress and anxiety (de Vries et al., 2001; Yoshiuchi et al., 2008), and sexual abuse (Simonich et al., 2004). Similarly, in smoking research, some EMA studies include those studying relapse in people who are quitting smoking (Shiffman, 2005), relapse among adolescent smokers (Gwaltney et al., 2008), examining the urge to smoke (O’Connell et al., 1998), and our own EMA studies on adolescents (Mermelstein et al., 2002, 2007).
Data from EMA studies are inherently multilevel with, for example, (level-1) observations nested within (level-2) subjects. Thus, linear mixed models (LMMs, aka multilevel or hierarchical linear models) are increasingly used for EMA data analysis (Walls and Schafer, 2006; Schwartz and Stone, 2007). A basic characteristic of these models is the inclusion of random subject effects into regression models in order to account for the influence of subjects on their repeated observations. The variance of these random effects indicate the degree of variation that exists in the population of subjects, or the between-subjects variance. Analogously, the error variance characterizes how much variation exists within a subject, or the within-subjects variance. These variances are usually treated as being homogeneous across subject groups or levels of covariates.
In EMA studies, it is common to have up to thirty or forty observations per subject, and this allows greater modeling opportunities than what conventional LMMs allow. In particular, one very promising extended approach is the modeling of both between-subject (BS) and within-subject (WS) variances as a function of covariates, in addition to their effect on overall mean levels. For example, if a smoker’s mood is the outcome, then one can consider the effect of covariates on their mood level (e.g., how happy/sad are they on average), as well as on their variation in mood (e.g., how labile/erratic is their mood). Or, one can examine mood changes when a person smokes in terms of the mean (does mood improve?) and variance (does mood stabilize?), and what variables might be related to those smoking-related changes of mood level and variation. Thus, by allowing WS variance to be a function of covariates, we can more directly examine the hypothesis that smoking helps to regulate mood (Russell et al., 1974). Indeed, in a recent review, Parrott (2006) stated that mood vacillation and its relationship to nicotine dependency is an important topic for future research.
Data from a natural history study of adolescent smoking motivated the application of the mixed effects ordinal location scale model. Students included in this study were either in 9th or 10th grade at baseline, 55.1% female, and self-reported on a screening questionnaire 6–8 weeks prior to baseline that they had smoked at least one cigarette in their lifetime. The majority (57.6%) had smoked at least one cigarette in the past month at baseline. Written parental consent and student assent were required for participation. A total of 461 students completed the baseline measurement wave. The study utilized a multi-method approach to assess adolescents in terms of self-report questionnaires, a week-long time/event sampling method via palmtop computers (EMA), and in-depth interviews.
Here, we focus on the EMA data. Adolescents carried the hand-held computers with them at all times during a seven consecutive day data collection period. They were trained to both respond to random prompts from the computers and to event record (initiate a data collection interview) smoking episodes. We refer to the former as random prompts and the latter as smoking events. Questions included ones about place, activity, companionship, mood, and other subjective items. The hand-held computers date and time-stamped each entry. For the analyses reported, we treated the responses obtained from the random prompts. In all, there were 14,105 random prompts obtained from the 461 students with an approximate average of 30 prompts per student (range = 7 to 71). We also used information from the self-initiated smoking events as covariates (this is described in Section 6).
The mood outcome we considered was a subject’s assessment of their negative mood before the prompt signal. Specifically, subjects were asked to rate their pre-prompt mood response to the item Sad: I felt sad. This was rated on a 10-point Likert scale with “10” representing very high levels of the attribute (i.e., sadness). Over all prompts, and ignoring the clustering of the data, the category frequencies are listed in Table 1. As can be seen, the distribution is highly skewed with the first category (i.e., the lowest level of sadness) as the modal response.
Of interest is the degree of heterogeneity in this mood measure in terms of both WS and BS variation. Furthermore, it is of interest to examine whether certain covariates can explain some of the variation in these two sources of heterogeneity, over and above the influence of these covariates on the mean response. It is also reasonable to allow random subject effects for both the mean response (to allow for subjects with different average levels of mood) and for a subject’s WS variance (to allow for different levels of mood consistency). These considerations led to the application of the mixed location scale model.
Let i = 1, …, N denote the subjects, j = 1, …, ni the occasions, c = 1, 2, …, C denote the response categories, and Yij the ordinal response associated with subject i and occasion j. Ordinal regression models are often specified in terms of the cumulative comparisons of the ordinal outcome. For this, define the cumulative probabilities for the C categories of the ordinal outcome Y as , where pijk denote the individual category probabilities. The mixed-effects logistic regression model for the cumulative probabilities (Hedeker and Gibbons, 1994) is given in terms of the cumulative logits λijc (c = 1, …, C − 1) as:
where xij is the p × 1 covariate vector and β is the p × 1 vector of regression parameters. The regressors can either be at the subject level, vary across occasions, or be interactions of subject-level and occasion-level variables. The random subject effect υi indicates the influence of individual i on his/her repeated assessments. The population distribution of these random effects is usually assumed to be a normal distribution with zero mean and variance . The model also includes C-1 strictly increasing model thresholds γc (i.e., γ1 < γ2 <·< γC−1).
As written above, a positive coefficient for a regressor indicates that as values of the regressor increase so do the odds that the response is greater than or equal to c. This agrees with the formulation in McCullagh (1980), though it is not the only way to write the ordinal model (e.g., see Fielding et al. (2003)). Also, note that the relationship between the regressors x and the cumulative logits does not depend on c (hence the regression coefficients β do not carry the c subscript). McCullagh (1980) calls this assumption of identical odds ratios across the C − 1 cut-offs the proportional odds assumption.
Ordinal regression models are often motivated and described using the “threshold concept” (Bock, 1975). This is also termed a latent variable model for ordinal variables (Long, 1997). For this, it is assumed that a continuous latent variable y underlies the observed ordinal response Y. Typically, the continuous latent variable y is assumed to follow either a normal or logistic distribution, leading to ordinal probit or logistic regression, respectively. Under the threshold concept, the observed ordinal outcome Yij = c if γc−1 ≤ yij < γc for the latent variable (with γ0 = −∞ and γC = ∞). In this article we will consider the logistic formulation.
In the above model, represents the between-subjects (BS) variance. To extend the model to also include within-subject (WS) variance, we note that ordinal regression models including scaling effects are common in ROC analysis (Tosteson and Begg, 1988). Here we will add these terms within the mixed model, namely,
To allow covariates to influence the BS andWS variances, we can utilize a log-linear representation, as has been described in the context of heteroscedastic (fixed-effects) regression models (Harvey, 1976; Aitkin, 1987), namely,
The variances are subscripted by i and j to indicate that their values change depending on the values of the covariates ui and wij (and their coefficients). The number of parameters associated with these variances does not vary with i or j. ui will usually include a (first) column of ones for the reference BS variance (α0). Thus, the BS variance equals exp α0 when the subject-level covariates ui equal 0, and is increased or decreased as a function of these covariates and their coefficients α. Specifically, for a particular covariate u*, if α*> 0, then the BS variance increases as α* increases (and vice versa if α* < 0). Note that the exponential function ensures a positive multiplicative factor for any finite value of α, and so the resulting variance is guaranteed to be positive. The WS variance is modeled in the same way, although, for identification, an intercept cannot be included as one of the w variables. However, both time-varying and subject-varying covariates can influence the WS variance. For this reason, the covariate vector is indicated as wij for the WS variance. Thus, this model allows both subject-varying and time-varying covariates to influence the WS variance, but only subject-varying variables to influence the BS variance. The coefficients in α and τ indicate the degree of influence on the BS and WS variances, respectively, and the ordinary random intercept model is obtained as a special case if α = τ = 0 for all covariates in ui and wij (i.e., excluding the reference variance α0).
We can further allow the WS variance to vary across subjects, above and beyond the contribution of covariates, namely,
where the random subject (scale) effects ωi are distributed in the population of subjects with mean 0 and variance . The idea for this is akin to the inclusion of the random (location) effects in equation (1), namely, covariates do not account for all of the reasons that subjects differ from each other. The υi parameters in (1) indicate how subjects differ in terms of their means and the ωi parameters in (4) indicate how subjects differ in variation, beyond the effect of covariates. Notice that taking logs in (4) yields , which indicates that if the distribution of ωi is specified as normal, then the random effects serve as log normal subject-specific perturbations of the WS variance. In other words, the WS variances follow a log normal distribution at the individual level. The skewed, nonnegative nature of the log normal distribution makes it a reasonable choice for representing variances, and it has been used in many diverse research areas for this purpose (Fowler and Whitlock, 1999; Leonard, 1975; Renò and Rizza, 2003; Shenk and Burnhamb, 1998; Vasseur, 1999).
In this model, υi is a random effect which influences an individual’s mean, or location, and ωi is a random effect which influences an individual’s variance, or (square of the) scale. Thus, we dub the model with both types of random effects as a mixed-effects location scale model. These two random effects are correlated with covariance parameter συω. This covariance parameter indicates the degree to which the random location and scale effects are associated.
Visually, Figure 1 presents the pertinent details of the model. In this figure, the average logit across all subjects is depicted with the solid horizontal line, and the lines of two subjects are presented as dotted lines. Data points for these two subjects are also included in the plot. In a given dataset, there will be as many dotted lines as there are subjects, but for simplicity here we only plot two subjects. Also, for simplicity, consider all covariates as subject-varying only. The effect of covariates x on the logit is represented by β; this effect influences the solid line by either raising or lowering it as a function of the covariates. Each dotted line is indicative of a person’s random location effect υi, which indicates how a subject deviates from the average logit response. In Figure 1, one subject is above and another subject is below the average line. The heterogeneity in these dotted lines is indicative of BS variance: if the dotted lines are close together then there is not much subject heterogeneity, conversely if the dotted lines are spread out then more heterogeneity is indicated. The effect of covariates u on this heterogeneity is represented in the model as α. The degree of variation of a person’s data points around their line is the WS variance. If the points are quite close to a subject’s line, then that subject has low WS variance (i.e., the lower subject in Figure 1). Conversely, if a subject’s data points vary greatly around that person’s line then more WS variation is indicated (i.e., the upper subject in Figure 1). Covariates w influence this WS variation through the coefficients τ. Finally, the model posits that covariates do not explain all of WS variance by including the random scale effect ωi.
Parameter estimation can be solved using maximum likelihood (ML) or variants of ML. Such solutions are iterative ones that can be numerically quite intensive. Here, the approach is sketched; further details can be found in Hedeker and Gibbons (1994) and Fahrmeir and Tutz (2001). Let Y i denote the vector of responses from subject i, and let the vector θi represent the two random effects (i.e., υi and ωi). The probability of any response pattern Y i (of size ni), conditional on the random effects θ, is equal to the product of the probabilities of the level-1 responses:
and Ψ(·) represents the logistic cumulative distribution function (cdf). The assumption that a subject’s responses are independent given the random effects (and therefore can be multiplied to yield the conditional probability of the response vector) is known as the conditional independence assumption. The marginal density of Y i in the population is expressed as the following integral of the conditional likelihood (·)
where f(θ) represents the distribution of the random effects, namely a bivariate normal density. Whereas (5) represents the conditional probability, (7) indicates the unconditional probability for the response vector of subject i. The marginal log-likelihood from the sample of N subjects is then obtained as . Maximizing this log-likelihood yields ML estimates, which are sometimes referred to as maximum marginal likelihood estimates because they are obtained by maximizing the marginal likelihood.
In order to solve the likelihood solution, integration over the random-effects distribution must be performed. As a result, estimation is much more complicated than in models for continuous normally-distributed outcomes where the solution can be expressed in closed form. Various approximations for evaluating the integral over the random-effects distribution have been proposed in the literature; many of these are reviewed in Rodríguez and Goldman (1995). Perhaps the most frequently used methods are based on first- or second-order Taylor expansions. Marginal quasi-likelihood (MQL) involves expansion around the fixed part of the model, whereas penalized or predictive quasi-likelihood (PQL) additionally includes the random part in its expansion (Goldstein and Rasbash, 1996). Unfortunately, these procedures yield estimates of the regression coefficients and random effects variances that are biased towards zero in certain situations, especially for the first-order expansions (Breslow and Lin, 1995). To remedy this, Raudenbush et al. (2000) proposed an approach that uses a combination of a fully multivariate Taylor expansion and a Laplace approximation. This method yields accurate results and is computationally fast. Also, as opposed to the MQL and PQL approximations, the deviance obtained from this approximation can be used for likelihood-ratio tests.
Numerical integration can also be used to perform the integration over the random-effects distribution (Bock and Lieberman, 1970). Specifically, if the assumed distribution is normal, Gauss-Hermite quadrature can approximate the above integral to any practical degree of accuracy. Additionally, like the Laplace approximation, the numerical quadrature approach yields a deviance that can be readily used for likelihood-ratio tests. The integration is approximated by a summation on a specified number of quadrature points for each dimension of the integration. An issue with the quadrature approach is that it can involve summation over a large number of points, especially as the number of random-effects is increased. To address this, methods of adaptive quadrature have been developed that use a few number of points per dimension that are adapted to the location and dispersion of the distribution to be integrated (Rabe-Hesketh et al., 2002).
An approach that is often used in econometrics and transportation research uses simulation methods to integrate over the random-effects distribution (Stern, 1997; Train, 2003). When used in conjunction with maximum likelihood estimation, it is termed “maximum simulated likelihood” or “simulated maximum likelihood.” The idea behind this approach is to draw a number of values from the random-effects distribution, calculate the likelihood for each of these draws, and average over the draws to obtain a solution. Thus, this method maximizes a simulated sample likelihood instead of an exact likelihood, but can be considerably faster than quadrature methods, especially as the number of random effects increases (Haan and Uhlendorff, 2006).
Bayesian approaches, such as the use of Gibbs sampling (Geman and Geman, 1984) and related methods (Tanner, 1996), can also be used to integrate over the random effects distribution. For multilevel data, this is described in detail by Draper (2008) and Gelman and Hill (2006).
Several computer programs provide maximum likelihood estimates of the mixed-effects proportional odds model, including SAS PROC NLMIXED, Stata (StataCorp, 1999), HLM (Raudenbush et al., 2004), MLwiN (Rasbash et al., 2004), LIMDEP (Greene, 2002), GLLAMM (Rabe-Hesketh et al., 2001), Mplus (Muthén and Muthén, 2001), and MIXOR (Hedeker and Gibbons, 1996). Some of these programs, however, only allow random-intercepts models, and few allow the inclusion of the scaling parameters described here. Also, these programs differ in how the integration over the random effects is performed. For the analyses presented in this article, SAS PROC NLMIXED, which utilizes adaptive quadrature for integration of the random effects, was used and is detailed in the appendix. In terms of software for maximum simulated likelihood, LIMDEP (Greene, 2002) has included this estimation approach for several types of outcome variables, including nominal and ordinal, and Haan and Uhlendorff (2006) describe a Stata routine. For estimation under the Bayesian approach, the freeware BUGS software program (Spiegelhalter et al., 1995) is popular; Marshall and Spiegelhalter (2001) provides an example of multilevel modeling using BUGS, including some syntax and discussion of the program.
To begin, we estimated a random-intercepts proportional odds model. Namely, a model without any covariates for the variances (i.e., no ui or wij variables), but with covariates in xij. Subject-level covariates included Smoker (an indicator of whether the student is a current smoker, coded no=0 or yes=1; this was determined based on whether or not the subject provided at least one smoking event during the week-long data collection period), Male (coded 0=female or 1=male), and PropSmk which equaled the proportion of occasions (both random prompts and smoking events) that were smoking events. This variable is an indication of the level of a subject’s smoking. It should be noted that the mood outcome, Sad, is only from the random prompts and not from the smoking events. In terms of prompt-level covariates, we considered whether the subject was alone or not (coded 0=not alone or 1=alone) at the time of the random prompt. For this, we created both a between-subjects and within-subjects version (AloneBS and AloneWS) as described in Neuhaus and Kalbfleisch (1998), namely,
Notice that AloneBS, the first term on the right-hand side, equals the proportion of random prompts in which a subject was alone, and AloneWS, the latter term on the right-hand side, is the prompt-specific deviation relative to this proportion (i.e., it equals either 0—AloneBS or 1—AloneBS if the subject was not alone or was alone, respectively, for the given random prompt).
Table 2 lists the estimates of this model (for space, the threshold estimates are not listed). As can be seen, sadness is significantly less for males, and increased for smokers, loners (i.e., subjects with higher levels of AloneBS), and when subjects are alone. Although smoking level has a negative effect on sadness, this effect is not significant. In terms of the BS variance, this is estimated to be . This can be expressed as an intraclass correlation (ICC) by noting that the variance for the underlying standard logistic distribution equals π2/3 (Agresti, 2002), namely ICC = 2.625/(2.625 + π2/3) = .44. Thus, the data are correlated within subjects to a moderate degree.
Next, we added these covariates into the loglinear models of the BS and WS variances (albeit, the WS variable AloneWS was not included in the BS variance model). For this, including the the random location effect (i.e., υi), we estimated models with and without the additional random scale effect (i.e., ωi). The results for these two models are listed in Table 3.
In terms of the effects on location (i.e., the βs), inspection of Table 3 reveals that the same basic conclusions are found as in the proportional odds mixed model, although the estimates are appreciably smaller and p-values larger. In terms of the variance modeling, first, inspection of the BS variance effects indicates that males are considerably less heterogeneous as a group than females. This highly significant negative effect on sadness heterogeneity is observed for both models in Table 3. For the WS variance effects, difference between the two models of Table 3 emerge. The model without the random scale effect suggests that all of the covariates are significantly related to WS variation in sadness: smokers are more varied in their sadness responses, while males, loners, being alone, and smoking level are indicative of significantly less varied (i.e., more consistent) mood responses. However, the model that includes the random subject scale effect yields only a marginally significant result for smoking level (p < .09) and a non-significant result for the BS alone effect (p < .14). Additionally, the remaining variables have less significant effects in the model that includes random subject scale effects. As can be seen, a primary reason for this is that the standard errors for the WS variance effects are considerably larger in the model with random scale effects relative to the model without them. The model AIC value, in —2 log L metric, equals 44814 (without random scale effect) and 43322 (with random scale effect), supporting selection of the latter model.
Because of the focus on smoking in this study, it is particularly interesting to note that while being a smoker significantly increases the WS variance, the effect of smoking level (Propsmk) on WS variance is negative (though marginally significant in the random scale model). This suggests that while, overall, smokers are more varied in their responses, this variation diminishes to some degree as smoking level increases. In other words, it is the lower level smokers that are the most varied in their reported levels of sadness. This resonates with the hypothesis that smoking helps to regulate mood (Russell et al., 1974), and also the review of Parrott (2006), who indicated the need for research on mood vacillation and its relationship to nicotine dependency. To get a better sense of this, based on the WS estimates of these two variables, one can calculate the value of PropSmk that would be required to equalize the positive overall Smoker effect on WS variance. Namely, for the model without random scale effects, it is .325/.909 = .358, while in the model with random scale effects it is .371/1.116 = .332. In this dataset of 461 subjects, 234 subjects were classified as smokers in terms of the variable Smoker. Of these 234 subjects, the median value of PropSmk equaled .081, while the 90% percentile equaled .3, and the 95% percentile equaled .367. Thus, based on the model estimates and the percentiles of PropSmk, most smokers elicited more varied responses than non-smokers.
Finally, as can be seen from Table 3, the random scale variance and covariance parameters are both highly significant. Thus, there is clear evidence that the WS variance varies across individuals, above and beyond the contribution of the covariates in the WS variance model. In other words, subjects differ in terms of their variation in sadness. The covariance parameter συω is estimated to be negative (and is highly significant). This suggests that subjects who are relatively high in terms of their mean sadness level are less varied across prompts in their sad responses.
This article has illustrated how mixed models can be used to model differences in variances, and not just means, across subject covariates. As such, these models can help to identify predictors of both within-subjects and between-subjects variation, and to test hypotheses about these variances. Additionally, by including a random subject effect on the WS variance, this model can examine the degree to which subjects are heterogeneous in terms of their variation on the outcome variable. Our example with sadness clearly shows that subjects are quite heterogeneous in terms of their mood variation, as one might expect.
More applications of this class of models clearly exist. For example, many questions of both normal development and the development of psychopathology address the issue of variability or stability in emotional responses to various situations and/or contexts. Often, a concern is with the range of responses an individual gives to a variety of stimuli or situations, and not just with the overall mean level of responsivity. These models also allow us to examine hypotheses about cross-situational consistency of responses as well.
In this paper, we have only considered the case of a single random subject effect for location. This could be generalized to allow multiple location random effects. For example, it is typical in longitudinal studies in which time is a factor, to consider a random subject intercept as well as random time trend parameters. However, for EMA data, there is not necessarily a notion that a person has some kind of systematic trend over the random prompts. In any case, the model could clearly be extended to allow multiple random location effects, and SAS PROC NLMIXED could still be used to estimate such a model.
Modern data collection procedures, such as EMA and/or real-time data captures, usually provide a fair amount of both WS and BS data, and so give rise to the opportunity for modeling of both WS and BS variances as a function of covariates. Clearly, these data from this EMA study, and the questions of the investigators, motivated the development of the model presented in this article. One might wonder about how much WS and BS data are necessary for estimation and variance modeling purposes. For random coefficient models, Longford (1993) noted the difficulty with providing general guidelines about the degree of complexity, for the variation part of a model, that a given dataset could support. This would also seem to be true here. Nonetheless, carrying out some simulations with relatively small sample sizes (e.g., 20 subjects with 5 observations each) gives the general impression that the primary issue is that the estimation algorithm does not often converge, but instead has estimation difficulties of one sort or another, in small sample situations.
As this is a relatively new modeling technique, certain limitations and cautions should be mentioned. First, our model assumes that the random location effects are normally distributed and that the random scale effects are lognormally distributed. It is unclear how robust this model is to violations of these assumptions. To some extent, this can be examined empirically using the approach of Liu and Yu (2008) for estimating models with non-normal random effects. Because of the focus on the variance of the ordinal outcome, the number of ordinal categories and relative frequencies within the categories might be an aspect that needs further examination. Finally, attention should be paid to outliers and influential observations, as these might have undue effects on estimation of the model parameters, especially the variance parameters.
Below is a sample of syntax necessary to run the mixed ordinal location scale model described in this article. In this syntax, uppercase letters are used for SAS specific syntax and lowercase letters are used for user defined entities. In terms of the variables used in this syntax, y denotes the outcome, x1 denotes a prompt- or time-varying covariate, x2 denotes a subject-level or time-invariant covariate, and id is a subject identifier. For simplicity, here we consider an ordinal outcome y with only three categories; the two thresholds are named gamma1 and gamma2, the two cumulative logits are clogit1 and clogit2, and the two cumulative probabilities are cprob1 and cprob2. The random location effect is named u1 and the random scale effect is named u2. The model for the mean response is summarized by mean, with the regression coefficients (β) named beta1 and beta2. The model for the BS variance is given by bsvar, with alpha0 indicating the reference BS variance (i.e., the between-subjects variance when the covariate x2 equals 0), in ln units, and alpha1 characterizing how this variance varies with x2. Similarly, for the model of the within-subjects (WS) variance, wsvar is modeled with coefficients tau1 and tau2 specified for the two WS variance influences x1 and x2, respectively.
PROC NLMIXED GCONV=1e-12; PARMS beta1=−.5 beta2=.3 gamma1=−1 gamma2=0 alpha0=1 alpha1=0 tau1=0 tau2=0 scalevar=.05 cov=0; mean = beta1*x1 + beta2*x2 + u1; bsvar = EXP(alpha0 + alpha1*x2); wsvar = EXP(tau1*x1 + tau2*x2 + u2); clogit1 = (gamma1 − mean)/SQRT(wsvar); clogit2 = (gamma2 − mean)/SQRT(wsvar); cprob1 = 1/(1 + EXP(−clogit1)); cprob2 = 1/(1 + EXP(−clogit2)); IF (y=1) THEN p = cprob1; ELSE IF (y=2) THEN p = cprob2 − cprob1; ELSE IF (y=3) THEN p = 1 − cprob2; logl = LOG(p); MODEL y ~ GENERAL(logl); RANDOM u1 u2 ~ NORMAL([0,0], [bsvar, cov, scalevar]) SUBJECT=id; RUN;
Users must provide starting values for all parameters on the PARMS statement. To do so, it is beneficial to run the model in stages using estimates from a prior stage as starting values and setting the additional parameters to zero or some small value. For example, one can start by estimating a random-intercepts ordinal model with fixed effects (β), BS variance (alpha0), and threshold parameters (gamma1, gamma2). Estimates of these parameters can then be specified as starting values in a model that adds in the WS variance parameters τ, and then the BS variance parameters α (or vice versa). Finally, the full model with the additional parameters (or scalevar) and συω (or cov) can be estimated. In practice, this approach works well with PROC NLMIXED, which sometimes has difficulties in converging to a solution for complex models. Also, in our experience, it seems that specifying a small starting value for the second random effect variance ( or scalevar) helps model convergence. Furthermore, for complex models, it is sometimes the case that the default convergence criteria is not strict enough. In the above syntax, the convergence criteria is specified as GCONV=1e-12 on the PROC NLMIXED statement. The results in this article did change a bit as this stricter criteria was applied, relative to the default specification, however the results did not change beyond this level. It would seem that this level is a reasonable choice, however it probably should be examined on a case-by-case basis.
*This work was supported by National Cancer Institute grant 5PO1 CA98262. The authors thank Siu Chi Wong for assisting with data analysis.
Donald Hedeker, University of Illinois at Chicago, School of Public Health (MC 923), 1603 W. Taylor Street, Chicago, IL 60612-4336.
Hakan Demirtas, University of Illinois at Chicago, School of Public Health (MC 923), 1603 W. Taylor Street, Chicago, IL 60612-4336.
Robin J. Mermelstein, University of Illinois at Chicago, Institute for Health Research and Policy (MC 275), 1747 West Roosevelt, Chicago, Illinois 60608.