|Home | About | Journals | Submit | Contact Us | Français|
This paper presents a continuous-time Bayesian model for analyzing durations of behavior displays in social interactions. Duration data of social interactions are often complex because of repeated behaviors (events) at individual or group (e.g., dyad) level, multiple behaviors (multistates), and several choices of exit from a current event (competing risks). A multilevel, multistate model is proposed to adequately characterize the behavioral processes. The model incorporates dyad-specific and transition-specific random effects to account for heterogeneity among dyads and interdependence among competing risks. The proposed method is applied to child-parent observational data derived from the School Transitions Project to assess the relation of emotional expression in child-parent interaction to risk for early and persisting child conduct problems.
The increasing availability of sequential behavioral data at the microlevel for social interactions has allowed researchers to model behavioral interaction processes between two or more individuals, such as couples, parents and children, and peers (Patterson, 1982; Gardner, 1993; Gottman, 1994; Lindeboom & Kerkhofs, 2000; Dagne et al., 2002, 2003; Howe, Dagne & Brown, 2005). The patterns of social interaction that may be extracted from such sequential behavior data are linked to (1) types of behaviors, (2) recurrence of behaviors, and (3) duration of a behavior act in specified behavioral state before making transition to other defined states over time. To date, little attention has been given to the modeling of patterns of social interaction by considering these features of behavioral data. This paper contributes to this gap of knowledge by developing a Bayesian multilevel, multistate competing risks duration model for behavioral observation data.
In health-related applications, social interaction between therapists and clients, parents and children, husbands and wives, and peers are often construed as critical social processes (Griffin & Gardner, 1989; Snyder et al., 2003; Dagne et a., 2007). Many current theories of mental health focus on social interaction patterns that begin in early childhood as risk factors for mental health problems that emerge many years later in young adulthood (e.g., the theory of antisocial behavior of Patterson, Reid & Dishion, 1992). In such social interactions, the types of behaviors displayed along with their durations by two or more individuals are the main interest of research. The occurrence of particular behaviors (events) of individuals can sometimes be recurrent and complex. That is, there may be repeated events with multi-state origins and multi-state destinations from each origin (competing risks). In this article, we propose multilevel semi-Markov models to analyze such complex duration data by not only incorporating covariates but also taking into account various potential sources of variation across dyads and among transition rates. Typical sources of variations in a complex duration dataset are (1) dyad or individual level heterogeneity and (2) transition-specific random effects. These random effects are allowed to correlate to each other, and follow a multivariate normal distribution.
Thus, we propose hierarchical parametric duration models for analyzing repeated events combined with multiple origin states and competing risks simultaneously. When dealing with many competing risks and allowing for the possibility of intra-dyadic correlation because of repeated events and inter-dependence among competing risks, a Bayesian approach provides greater flexibility in estimation than a frequentist approach (Gelman, Carlin, Stern, & Rubin, 1995). The reason for this that the number of unknown parameters in these types of competing risks duration models could become very large especially under an unrestricted variance-covariance matrix for random effects and a large number of behavior states.
The remainder of the paper is organized as follows: In the next section, we will develop the proposed models, and followed by a description of methods of estimation. Application of the models to real data and discussion of results are given in the Example and Conclusion sections, respectively.
When repeated durations of behavior are observed, such as multiple behavior sojourn times for an individual or a dyad, the durations within a dyad tend to be more alike than those across dyads, thereby generating a correlation at the dyad level. Ignoring this correlation will result in understating the standard errors of estimates for the effects of background covariates (Diggle, Heagerty, Liang, & Zeger, 2002, p. 19). A multilevel model takes into account this correlation directly by using random effects, which are commonly known as frailties. Even after adjusting for the frailties, there may be additional correlations due to unobservable factors (e.g., propensity for exhibiting particular behaviors) that affect transitions from a specific origin state to another specific destination state. To properly account for such correlations, we introduce additional levels of random effects for each origin-destination transition in the modeling process, described below. These random effects are allowed to correlate to each other, and follow a multivariate normal distribution. In addition to the random effects, covariates are also incorporated to determine which covariates might impact the occurrence of a behavior and how long it lasted. These models help us answer questions of why and how a particular behavior may change faster under certain conditions and persist under other conditions.
Let m denote the total number of dyads in a sample, and these dyads move in and out of a finite number of states according to a continuous-time Markov chain. If the chain is in a certain state, it remains there for some time before moving to another state. Transition probabilities or hazard rates are used to characterize this process. Let O represent origin states, and likewise D represents destinations (transitions). Suppose there are Do transitions in which an event in origin state o can end, where o = 1, 2, …, O, and d = 1, 2, …, Do.
The hazard of making a transition to destination d from origin state o after sojourn duration tij for the ith dyad (i = 1, 2, …, m) at the jth episode (j = 1, …, ni) is denoted by hod(tij), and it is given by
The hazard rate measures the instantaneous rate at which a particular behavior display of interest that had lasted time t in origin o was ending. High hazard rate for termination from an origin to one destination implies a low rate for termination to another. Ignoring this dependence can potentially bias inference. Random effects are often used to assess such dependencies (Yeshin, Manton & Vaupel, 1985; Therneau & Grambsch, 2000; Omar et al., 1999). Assuming these random effects operate multiplicatively on the baseline hazard, the hazard function could be re-written as
where ui is the dyadic-specific frailty, which is assumed to have a normal distribution with mean zero and variance τ2, quantifying the dependencies among intra-dyadic duration times and the heterogeneity among dyads; X is a vector of time invariant covariates, such that exp(Xodiβod) is the relative risk associated with covariates; Zodij is the transition design matrix in which zodij is an indicator variable that takes unity for the odi transition and zero otherwise, and bodi are the origin-destination specific random effects which are assumed to be independent of ui and multi-normally distributed with mean zero and covariance matrix Σ of dimension Do. Note that the covariate effects need not be the same for transition from origin o to destination d. hod(tij) is a baseline Weibull hazard function, defined as
where λ > 0 is a scale parameter, gauging the multiplicative effects of duration on the termination of a particular behavior, and the other parameter is γ > 0 which is a shape parameter. If γ < 1, the hazard function is decreasing over time; if γ = 1, there is a constant hazard rate over time, and if γ > 1 then the hazard rate is increasing over time. The hazard rates measure the magnitude of the risk that a dyad will move to another state from current state after staying duration t during child-mother interactions.
The corresponding integrated hazard function is
The survival function becomes
and the density for the duration is thus given by
where t is the sojourn duration of an event before transition to other competing risks.
We use a Bayesian estimation approach to fit model (2) above. Recent advances in Bayesian computation via Markov chain Monte Carlo algorithm (Gelman & Rubin, 1992; Gamerman, 1997) enable us to fit such models in a straightforward manner. The algorithm requires the specification of likelihood function and prior distributions for the unknown parameters. These are discussed below.
For constructing the likelihood function, we assume that there are O and D number of origin states and destination states, respectively. In addition, let δodij = 1 if the transition from the oth state to the dth destination (target behavior) for the ith dyad at the jth episode, and 0 otherwise. ζoi is another indicator variable taking the value 1 if the origin state is o for the ith dyad, otherwise it takes zero. If a transition is observed to a target destination from a given origin state, its contribution to the likelihood function is fod(tij|u, β), but if the transition is to other destinations, its contribution to the likelihood is Sod(tij|u, β). Thus, conditional on the random effects, the likelihood function becomes
where hod(tij|b, u, X, Z, θ) is a Weibull hazard function; , and ηodij = Xiβod + Zodijbodi + ui), ui ~ N(0, τ2), bodi ~ N(0, Σ), and θ = (γ, Σ, τ2, β).
As noted above, our estimation approach is Bayesian and requires the specification of prior distributions for the unknown parameters of the model.
We assume normal prior distributions for the fixed parameters, β, with mean β* and variance Ωβ which are chosen to make the distribution proper but diffuse with large variances. The origin-destination specific random effects, bodi ~ N(0, Σ); and and Σ−1 are distributed as Wishart distributions with hyperparameters ν1, V1 and ν2, V2 where V1 and V2 are positive definite matrices and ν1 and ν2 are degrees of freedom, respectively. The dyad-specific random effects, ui ~ N(0, τ2), where τ2 has an inverse gamma distribution with parameters a1 and c1. The inverse gamma prior distribution is computationally easy to use but it has limitations particularly when the number of units for second level or higher is small (Gelman, 2006). Likewise, γod has a gamma distribution with parameters a2 and c2. These hyperparameters need to be assessed to give noninformative but proper prior distributions. For our example (see Section 4), the assessed values of the hyperparameters are (a1, c1) = (.0001, .0001) for the prior distribution of τ−2, and (a2, c2) = (1, .0001) for γ; β* = 0, and V1 = .001I, V2 = .001I, where I is an identity matrix. We have also tried different values for the hyperparameters to check if the results were sensitive to choices of the hyperparameters. All converged to the same results indicating that the data dominated the inference.
The joint posterior distribution of the parameters of the model in (2) conditional on the data is obtained by combining the likelihood function given in (7) and the above specified prior distributions via Bayes’ theorem, as:
where, ; ; ; ; f(τ2|a2, c2) τ−a1 exp(−b1/2τ2); ; is degrees of freedom and V1 is a scale matrix; f(Σ−1|ν2, V2) |Σ|ν2−q−1)/2 exp(−0.5tr(Σ−1V2), ν2 ≥ q is degrees of freedom and V2 is a scale matrix.
From the joint posterior distribution given in (8) one may obtain marginal posterior distributions of the parameters for inferential purposes. To do this, one can use a freely available software called WinBUGS Version 1.4 (Spiegelhalter, Thomas, & Best, 2000) to get the marginal posterior distributions.
To illustrate the proposed method, we analyze data from the School Transitions Project (STP) which is a study of emotion regulation during parent-child interaction (Snyder et al., 2003). In this study, the interaction of kindergarten children with their mothers was videotaped for one hour on each of two occasions. On each occasion, the mother-child dyad engaged in a series of tasks including playing a novel game, problem solving, engaging in letter and number identification, and planning a fun activity. The videotaped interaction was coded by trained coders using the Specific affect Coding System (SPAFF) (Gottman et al., 1996). The SPAFF codes the onset and offset of both the mother’s and child’s behavior into mutually exclusive and collectively exhaustive categories on a real time.
Three categories for child behavior and two categories for mother behavior were used in these analyses. Child “anger” (designated as A) was defined by the codes for anger, contempt and disgust. Child “sad/fearful” (designated as S) was defined by the codes for sad and fear. Child “neutral” (designated as N) was used to designate periods when behavior does not fit the other codes. Mother “negative” (designated as M) was defined by the codes for anger, contempt, disgust, criticism, domineering, belligerence, threats, whining, stonewalling, defensiveness, fear, and sadness. Any other behavior, which was not negative, was defined mother “non-negative” (designated by R). Thus, a parent can be in one of two behavioral states, M or R, and the child can be in one of three states, A, S or N. The childmother behavioral combinations are formed in which a dyadic state is indicated by the child state abbreviation first followed by the parent state abbreviation (e.g., AM is the dyadic state of child angry-parent negative process). We have focused in this paper on two subsets of the possible transitions that have AM and AR origin (start) states:
Our focus in this paper is to study whether children’s anger exhibits duration dependence in the context of parents’ behavioral states. For this purpose, we use the AM start state. The dependent variable, t, is the sojourn time (duration) in each dyadic state before transition to another state. For example, the sojourn time for a dyad in an AM state (i.e. the dyadic state of child angry-parent negative) is the time it stays in this state before moving to either AR, SM or NM. The latter states are competing transition states (risks). For this duration data set, based on 267 families (dyads), it was observed that the relationship between log(−log(S(t))) and log(t) is approximately linear, indicating that the Weibull hazard rate fits the data reasonably well (see Figures Figures11--2).2). The definition of S(t) was given in equation (5) above.
To account for variation among dyads and interdependence among competing risks, we introduce two kinds of random effects. These random effects are: (1) dyad-specific random effects for accounting repeated events per dyad, and (2) transition-specific random effects for accounting for shared unobservable latent factors across competing risks or transition rates emanating from the same origin behavioral state. More detailed model specifications for the random effects were given in Section 2.
Estimation of the model parameters in (8) was carried out using WinBUGS (Spiegelhalter Thomas, & Best, 2000). Convergence with a three-chain run was achieved after 20,000 iterations as it was judged by the scale reduction factor between .8 and 1.2 (Gelman & Rubin, 1992). Posterior estimates of the parameters were computed based on subsequent 20,000 iterations after convergence and they are given in Tables Tables11--33.
From the outputs given in Table 1, the posterior estimates of the shape parameters (γod) of the Weibull hazard models for transition from AM origin state to the other three competing risks (AR, SM, and NM) are (1.238, 1.311, 1.115), respectively. The estimate of the shape parameter for the hazard rate of AM → AR is greater than unity, suggesting that there is an increased hazard of terminating negative mood and move to non-negative behavior for the mother while the child stays in angry state for longer duration. In the case of the hazard rates of AM → SM and AM → NM, however, interval estimates of their shape parameters contain unity, implying that the duration of AM behavior has constant hazard rates of its being terminated by either SM or NM. That is, the hazard of cooling down for a child from angry mood to sad mood or neutral state when the mother stays in negative mood is a chance occurrence regardless of how long the dyad stays in the current AM state. Such an interpretation, however, requires caution since the presence of random effects in the model may underestimate the hazard rate (Omeri & Johnson, 1993). In terms of the scale parameters (λ = exp(ν)) of the hazard models, the results show that there are different interaction patterns among the competing risks over time. Those patterns may be explained by adding dyad-level covariates in the hazard models (see Section 4.2).
In the middle part of Table 1, the posterior estimate of the standard deviation of the dyadic-specific random effects is given as .222 with a 95% credible interval of (.078, .425). The lower bound is not very close to zero, suggesting that there is a moderate non-zero standard deviation of heterogeneity among families (dyads) in terms of the length of durations of displaying AM behavior. If the standard deviation were zero, there would have been no heterogeneity among families, suggesting that the repeated episodes of behaviors were independent.
In addition to dyadic-level random effects, transition-specific random effects are also included in the hazard models. The estimates of the variances and covariances of these transition-specific random effects are at the bottom of Table 1. The estimated variances for the transitions AM → AR, AM → SM, AM → NM are .122, .378, and .154, respectively. The non-zero variances suggest that all 3 hazard rates of family transitions from AM state to AR, SM and NM are heterogeneous in this sample. The covariance between the hazard rates of transitions AM → SM and AM → NM is −0.850 with a 95% credible interval of (−1.942, −0.089). The interval does not include zero, implying that the correlation is significant. That is, children with a high transition rate from angry to sad are more likely to have low transition rate from angry to neutral behavior when the mothers are exhibiting negative behavior. In other words, when mother stays in negative mode, children who are likely to transition to neutral are much likely to transition to sad.
We have also fitted Weibull hazard models with AR origin state. The results are depicted in Table 2. The 95% credible intervals of the estimates of the shape parameters for the hazard rates of AR → AM, AR → SR, and AR → NR do not include unity, suggesting that, on average, the transitions from AR to the other competing behavioral states are faster as duration increases. When the mother is in a non-negative mood during interaction, it seems that the hazard for the child to move away from angry state is higher.
For this article, covariates hypothesized to be related to emotion expression during mother-child interaction are (1) maternal depression (MDEP), (2) maternal hostility (MHOS), (3) parent rated child overt antisocial behavior (POA) and (4) parent rated child covert antisocial behavior (PCA). These covariates were collected in the fall of the kindergarten year, which was the baseline year. Parent hostility and depression were derived from the Brief Symptom Inventory (Derogatis & Melisaratos, 1983), and measures of child antisocial behavior from the Child Behavior Checklist (Achenbach, 1991).
Our goal here is to characterize parent-child interactions and correlate these with the child’s social antisocial behavior and family characteristics mentioned above. We introduce time-invariate, dyad-level covariates to model (2) to assess the effects of these covariates on the hazard rates. That is, the durations for social behaviors are modeled to determine which background covariates might have impact on the hazard of transitions to AR, SM, or NM from AM state. For example, we are especially interested to estimate the contrast of AM → AR versus AM → NM which represents the difference between transition rates for parent and child representing backing down from an angry confrontation. We hypothesize that whatever difference exists in the overall population, there is mutual irritability. That is, mothers of overt antisocial kids are less likely to “back down,” and kids with hostile mothers are less likely to back down. To test such hypotheses we use the model given in (2) by adding covariates as follows.
where Xodi = (POA, PCA, MDEP, MHOS), βod is 4 × 1 vector of coeffcients for fixed covariates X.
The posterior means and posterior standard deviations of the effects of covariates are displayed in Table 3. For ease of interpretation of the results, we focus on who changes behavioral state, parent or child. So, for example, if the start state is AM (child angry and parent negative) and the end state is AR (child angry and parent non-negative), the transition could be described as “the parent cools down given the child is angry”. A negative (positive) sign means the parent is less (more) likely to cool down given the child is angry if they score high on a predictor compared to if they score low on the predictor.
Before interpreting the estimates in Table 3, it is also important to assess whether the more complex dependent risk model, which incorporates transition-specific random effects, offers a significant improvement over the independent risk model in terms of overall fit. To do this a Bayesian model comparison criterion, DIC (deviance information criterion) (Spiegelhalter et al., 2002), can be used to assess the contribution of the transition-specific random effects. The values of DIC for the dependent risk model and the independent risk model are 1621.69 and 1630.75, respectively. Thus, the former model has lower value of DIC than that of the independent risk model, suggesting a substantial improvement of goodness-of-fit for the dependent risk model.
Turning now attention to the dependent risk model in Table 3, the estimated coeffcient of POA (overt antisocial reported by parent) on the hazard rate of AM → NM is −.863 with a 90% credible interval of (−1.683, −.070). This shows an indication that a child with higher score on overt antisocial stays longer in angry state before cooling down to neutral behavioral state. The effect of PCA (covert antisocial reported by parent) on the hazard rate of AM → NM is 1.630 with a 95% credible interval of (.1591, 3.150). The result suggests that, a typical child with higher score on covert antisocial behavior moves away from angry state to neutral state faster while the mother stays in negative mood. Thus, the sojourn time in angry state is shorter for a child with covert antisocial than a child with overt antisocial. The other credible coeffcient is for the impact of MHOS (maternal hostility) on the hazard rate of AM → AR. The estimated coeffcient is −1.047 with 95% credible interval (−1.86, −.25). It means that a mother with higher score on maternal hostility takes her longer time to terminate negative behavior and move to non-negative behavior when the child displays anger. Thus, kids with hostile mothers are less likely to back down, and mothers of overt antisocial kids are also less likely to back down.
We now turn to an examination of the transition-specific random effects variances and covariances in the bottom part of Table 3. The variance estimates for the random effects of transitions AM → AR, AM → SM, and AM → NM are .122, .378, and .154, respectively. The estimates show that there is evidence of unobserved heterogeneity between families (dyads) in the hazards of terminating the AM behavioral state. In other words, families in the sample show strong variation in how long they stay in negative state (AM) before moving to the competing behavioral states of “child angry & mother non-negative”, “child sad/fearful & mother negative” or “child neutral & mother negative”. The existence of such an unobserved heterogeneity suggests that there may be some important missing factors, such as time-varying characteristics of families, that yet to be accounted for explaining the variation of duration among the dyads. The covariances between the transition-specific random effects of pairs of transitions are not strong. The covariance between the hazards of the transitions from child being angry to sad and neutral is negative and its estimated value is −0.105. The negative covariance indicates that children with a high (low) hazard of moving from angry state to sad state when their mothers show negative emotions during interaction tend to have low (high) hazard rate of moving from angry behavioral state to neutral behavioral state.
In this paper, we have described repeated events (behaviors) and competing risks model with heterogeneity to study emotion regulation of children during child-parent interactions. We have adopted a Bayesian approach and a continuous-time semi-Markov model to describe those social interactions. The paper also extends previous work (e.g., Gardner & Griffin, 1989) in the area of duration modeling in social interaction in three ways. First, transition-specific random effects are incorporated to account for unobserved heterogeneity in hazard rates for moving from one behavioral state to another behavioral state, in addition to dyad-level random effects to account for the within-dyad repeated events of the same behavior. A term for unobserved heterogeneity helps us capture the influence of missing covariates on the transition rates. Second, transition-specific random effects are allowed to be correlated. We have assumed that these random effects follow multivariate normal distribution. Third, for estimating the parameters of the proposed models we have used a multilevel Bayesian method. A Bayesian method is flexible and computationally straightforward to implement, and especially recent advances in Bayesian computational methods (e.g., WinBUGS) have alleviated the diffculty of fitting complex models to a great extent.
In addition, the proposed models in this paper provide a mechanism of assessing the pattern of emotion regulation of children. That is, by considering how long the child or the mother took to make a choice for transition, we can learn about his or her emotion regulation processes. For example, we were able to test the hypothesis of mutual irritability. That is, mothers of overt antisocial kids are less likely to “back down,” and kids with hostile mothers are less likely to back down.
An examination of the transition-specific random effects variances and covariances after controlling time-independent covariates showed that there were unobserved heterogeneities between families (dyads) in the hazard rates for the transitions considered in this paper. This suggests that there may be some unobserved latent factors, such as time-varying characteristics of families, that yet to be accounted for explaining the heterogeneity among the dyads.
The proposed method may be applied to other dyadic processes which involve interpersonal relationships and mutual influences. Examples of such processes are marital couples (Raudenbush, Brennen, & Barnett, 1995), twins (Segal, & Hershberger, 1999), employees and supervisors (Fleenor, McCauley, & Brutus, 1996), doctors and patients (Goldberg, Cohen, & Rubin, 1998), or teachers and students (Dolan et al., 1993).
Our models can also be extended to more complex situations where within-individual between-episode random variations are present. In principle one can incorporate such variations in a straightforward manner. They can be incorporated in a multilevel model fashion by using episode-specific random effects, provided there are relatively high proportion of dyads who experience more than one transition of each behavior. Another potential extension of our method is to incorporate latent class mixtures for identifying subgroups of children from the sample who may show high or low risk for early and persisting child emotional problems.
This work was supported in part by the National Institute of Mental Health (NIMH) and the National Institute on Drug Abuse (NIDA) under grant 5R01-MH40859 and by the NIMH under grant number MH57342.
The authors would like to thank the editor and two anonymous reviewers for their valuable comments.
Getachew A. Dagne, University of South Florida.
James Snyder, Wichita State University.