inertia_coordination_V03
2019-10-17
As for all rties analyses, the first step for an Inertia-Coordination analysis is to follow the instructions in “overview_data_prep” to visualize and prepare the data. We include only the minimal required steps here:
The Inertia-Coordination model represents the within-person pattern of inertia, which is defined as the extent to which a person’s state can be predicted from his or her own state at a prior time point (a.k.a. auto-regression), and the between-person pattern of coordination, which is defined as the extent to which one partner’s state can be predicted from their partner’s state either concurrently or time-lagged (Reed, Randall, Post, & Butler, 2013). Specifically, the time-series state variable is predicted by: 1) separate intercepts for each partner, 2) a person’s own state variable at a prior time point, which gives two “inertia” estimates, one for each partner, and 3) the person’s partner’s state variable at the same prior time point, which gives two “coordination” estimates, again one for each partner. This model is identical to the “Stability-Influence” model (Thorson, West & Mendes, 2017), but we prefer the terms “inertia” and “coordination” rather than “stability” and “influence” because they have fewer connotations that may or may not be appropriate for a given research design. The lm model used by the rties functions is:
lm(obs_deTrend ~ -1 + dist0 + dist1 + dist0:obs_deTrend_Lag + dist0:p_obs_deTrend_Lag + dist1:obs_deTrend_Lag + dist1:p_obs_deTrend_Lag, na.action=na.exclude, data=datai)
where “obs_deTrend” is the observed state variable with individual linear trends removed (e.g., it is the residuals from each person’s state variable predicted from time). The “-1” , “dist0” and “dist1” work together to implement a two-intercept model, whereby the overall intercept is omitted and instead separate intercepts are estimated for the level-0 and level-1 of the distinguisher variable provided by the user (for a discussion of this approach see: Kenny, Kashy, & Cook, 2006). The terms “dist0:obs_deTrend_Lag” and “dist1:obs_deTrend_Lag” estimate the inertia parameters for the partners scored 0 and 1 respectively on the distinguishing variable (e.g., “obs_deTrend_Lag” is the person’s own de-trended observed state variable lagged by how ever many steps the user specifies). Similarly, “dist0:p_obs_deTrend_Lag” and “dist1:p_obs_deTrend_Lag” estimate the coordination parameters for the partners scored 0 and 1 respectively on the distinguishing variable (e.g., “p_obs_deTrend_Lag” is the person’s partner’s de-trended observed state variable lagged by how ever many steps the user specifies). The model is estimated separately for each dyad (e.g., “datai” is the data from couple “i”).
If we consider the parameters of the model in isolation from each other, positive inertia estimates imply slower fluctuations of the state variable, while negative inertia estimates imply that the observed variable is oscillating back and forth between each lag (see Figures in overview_data_prep). For the between-person coordination parameters, a positive estimate implies an in-phase pattern, such that when one partner is high on the observed state variable, so is their partner at the specified lag, while a negative parameter implies an anti-phase pattern, such that when one partner is high the other partner is low at the specified lag (Randall, Post, Reed, & Butler, 2013; Reed et al., 2013; Wilson et al., 2018). It is important to realize, however, that the parameters do not act in isolation - they work together as a set to determine potentially complex trajectories for both partners over time. As with any non-linear model, it is impossible to inspect the inertia and coordination parameters and draw any conclusions about the dynamic pattern implied without actually plotting it using the full set of parameter estimates. The rties package provides a variety of functions for doing so that are described below.
Sample Size Considerations
There are two sample size considerations for any of the models implemented in rties. The first pertains to the number of observations per dyad that are required, while the second is the number of dyads required. The first consideration comes into play when we estimate the dynamics one dyad at a time. Greater complexity requires finer-grained measurement of time and hence more observations per dyad. One advantage of the Inertia-Coordination model is that it is fairly simple and hence requires relatively few observations per dyad. The exact number will depend on how much variance there is over time, both within-people and between-partners, but it is likely to provide good results with as few as 10 observations per dyad (someone should do a simulation study of this!).
The second sample size consideration comes into play when we use latent profiles based on the estimated dynamic parameters for each dyad to either predict the system variable, or be predicted by it (these steps are described in detail in sections below). In both cases, the system variable can be either a dyadic variable (e.g., both partners have the same score, as in relationship length) or an individual variable (e.g., partners can have different scores, as in age). In the case of predicting a dyadic system variable, a regular regression model is appropriate for normally distributed system variables, or a generalized linear model for non-normal ones (any of the families available for glm can be used). In this case, the shared system variable is predicted by the categorical latent profile membership and you can use your favorite rule of thumb along the lines of “n observations for a one-predictor regression model” to choose your sample size. Alternately, you could conduct a power analysis for the linear (e.g., regular regression) or generalized linear model you choose.
The situation is more complicated when the system variable is assessed at the individual level, or when it is the predictor of the latent profiles. In the former case, the system variable is predicted using a cross-sectional random-intercept dyadic model with the latent profile membership, the distinguisher variable and their interaction as fixed effects. For these models, it should be possible to use the R package simr to estimate the necessary sample size to obtain desired power. In the latter case, profile membership is predicted using either a binomial (for 2 latent profiles) or multinomial (for more than 2 latent profiles) model, with either just the system variable as a predictor (for dyadic system variables), or each partner’s system variable and their interaction as predictors (for individual system variables). For these models, it should be possible to use G-Power, or some other software for logistic and multinomial regression, to assess the needed sample size for a given power level.
Choosing the Lag Length
One complexity in using the Inertia-Coordination model is that the results are highly dependent upon the chosen lag length (for discussion of a similar issue when choosing the measurement interval see Boker & Nesselroade, 2002). This dependence on the lag makes interpretation problematic unless one has a strong theory about the temporal processes at work. At a minimum, it is a good idea to inspect the auto-correlation (relevant to inertia) and cross-correlation (relevant to coordination) plots for the observed state variables. If the state variables are oscillating at all, the auto- and cross- correlations will also oscillate depending on lag. For example, the correlations can vary from negative, to zero, to positive depending on the lag. To facilitate making an informed decision about lag length, rties provides functions to produce auto- and cross- correlation plots. Here we plot a chosen subset of dyads to show the strong impact that is sometimes found due to lag, especially for cross-correlation.
temp <- subset(data1, couple %in% c(2,5,27,31,47,60,103,110))
autoCorPlots(basedata=temp, dyadId="couple", personId="person", obs_name="dial", time_name="time")
In addition to visual inspection of the auto- and cross- correlations, choosing the lag relies on a combination of theory, prior research, and how quickly you expect the phenomenon of interest to be changing (see Thorson, West & Mendes, 2017, for additional discussion). For example, if the process is fairly stable, auto-correlation will be high across a long lag time and will dominate any results, because when estimated with the Inertia-Coordination model, the stronger the autocorrelation is, the weaker the cross-correlation will be. This results from the behavior of multiple regression models, where the coordination parameters can only account for independent variance in the outcome after accounting for inertia. Thus part of the decision depends upon how much of the observed state variables behavior you want to be explained by within-person auto-correlation processes versus how much you prefer to prioritize between-person cross-correlations. Another consideration is how many observation time points there were. For each lag you lose one observation point (e.g., if the lag is two steps, then the parameter estimates will be based on the total number of observations per dyad minus 2). One strength of rties is that it makes it very easy to alter the lag length and observe the impact on the results, which is helpful for developing an intuitive understanding of inertia and coordination.
Another approach provided for the lag is the “absMaxCC” option. If this is chosen for setting the lag, then the absolute maximum cross-correlation is found for each dyad and the lag at which that occurs is used to lag their data. This approach gives priority to between-partner interdependencies in the data, which may be justified when the research focus is on that component of the system rather than the within-person autocorrelation.
Assessing Model Fit
Having chosen a lag length (here we chose a lag of 5; see “overview_data_prep”" for our reasoning), we complete the data preparation with the “dataPrep” function.
The next step, which is often neglected in the literature, is to assess how well different variants of the Inertia-Coordination model fit the observed temporal data. Our ultimate goal is to either predict outcomes of interest from the dynamics assessed by the model, or to test whether other variables predict those dynamics. Either way, the results are only meaningful if the Inertia-Coordination model does, in fact, realistically represent the dynamics of the system. We therefore provide a function, “indivInertCoordCompare”“, that fits three versions of the model (inertia-only, coordination-only, inertia and coordination) to each dyad’s data and returns a list that includes: 1) the adjusted R^2 for each dyad for each of the 3 models (“R2inert”, “R2coord” and “R2inertCoord” e.g., how well each model predicts the observed temporal trajectories of the data), and 2) differences between the R2 for each model, where “R2dif_I_C” is the R2 for the inertia-only model minus the coordination-only model, “R2dif_IC_I” (inertia-coordination minus inertia), and “R2dif_IC_C” (inertia-coordination minus coordination-only). The function takes the name of the processed dataframe, here “data2”, and the results can be accessed with the “summary”" function:
In this example, we see that the full model accounts on average for about 37% of the variance in the observed state variables, while the inertia-only model accounts for almost as much (36%) and the coordination-only model accounts for only 10% on average and adds only about 1% to the inertia-only model when included in the full model. These results suggest that coordination is not important for describing the dynamics of the system, over and above inertia. It is still possible, however, that the small amount of variance explained by coordination may predict system variables, such as relationship quality or health, better than the larger amount of variance explained by inertia.
In addition to the model comparison results, we provide a function, “indivInertCoordPlots”“, to plot the observed state variables for each dyad, superimposed with their model predicted trajectories. The function takes as arguments the name of the processed dataframe, which of the 3 variants of the inertia-coordination model to plot the results for (”inert“,”coord“,”inertCoord“), optional strings providing labels for the level-0 and level-1 of the distinguishing variable, as well as the name of the observed state variable (default is”observed“) and minimum/maximum values for the y-axis in terms of quantiles (default is the observed minimum and maximum of the state variable). These plots provide a visual tool for assessing model fit at the individual dyad level. The figures can be accessed from the returned list called “plots” and they are also automatically saved as a .pdf file in the working directory (this process takes awhile). The following figures show examples of poor and good fit respectively:
Finally, we provide a function that plots the residuals from fitting the model to each dyad. As with any regression model, a good fit is indicated by normally distributed residuals with a mean of zero and suggests that the model assumptions have been met. In our use case, it provides evidence that the state variables are normally distributed enough to be appropriate for this model. The figures can be accessed from the returned list called “plots” and they are also automatically saved as a .pdf file in the working directory. Here we show a few examples of such plots.
Generating Dynamic Parameter Estimates
The next step is to use the function “indivInertCoord” to generate the dynamic parameter estimates for each dyad. The function takes the name of the processed dataframe (“data2” in this example) and the name of the version of the Inertia-Coordination model that you want the parameter estimates for (“inert”, “coord”, or “inertCoord”). It returns a dataframe containing the parameter estimates, called “params”, for use in the latent profile analysis (see below). In the returned dataframe, “int0” and “int1” are the estimated intercepts for the partner scored 0 amd 1 on the distinguisher respectively (e.g., their estimated state variable at time zero). The variables “inert0” and “inert1” are the inertia estimates for each partner and “coord0” and “coord1” are their coordination estimates.
Latent Profile Analysis (LPA)
The next step is to use the parameter estimates generated by “indivInertCoord” as input to a latent profile analysis. Because the models in rties all represent non-linear dynamics, the behavior of the dyadic system cannot be intuited by examining individual parameter estimates from the model. Rather, they operate together as a set to produce potentially complex temporal trajectories of both partner’s state variables (emotional experience in this example). Thus when it comes to using the dynamic parameter estimates to predict, or be predicted by, system variables such as conversation quality, we wish the parameters to be operating as a set, not as isolated individual predictors or outcomes such as would be the case if we used them as predictors in a regression model. Instead, we use them as input to a latent profile analysis in order to estimate qualitatively distinct groups of dyads based on the dynamics assessed by the Inertia-Coordination model. In other words, we wish to group the dyads based on their similarities and differences in terms of patterns of inertia and coordination.
The latent profile analysis makes use of the mclust package (Scrucca L., Fop M., Murphy T. B. and Raftery A. E., 2016). Specifically, the rties function “inspectProfiles” provides a wrapper to some of the functionality provided by mclust. LPA is a form of Gaussian Mixture Model, which assumes an underlying multivariate normal distribution and groups cases based on a set of observed variables (in our case the Inertia-Coordination estimates for each dyad). In much of the literature on LPA the focus is on trying to estimate the “true” number of underlying profiles by using fit statistics such as the BIC, but simulation studies tend to show that this process doesn’t work very well unless you have a huge sample and large effect size (see for example Tein, Coxe & Cham, 2013). Thus, in the context of rties we completely ignore fit statistics and suggest choosing the number of profiles based on: 1) interpretability in terms of meaningful dynamics, 2) the number of cases assigned to each profile, and 3) how well separated the profiles are. The “inspectProfiles” function provides results relevant to each of the these criteria, which we discuss below. We suggest considering results for 2, 3 and 4 profiles. More than that is unlikely to be interpretable, unless you have a theory that suggests 5 or more distinct dynamic patterns. Note that mclust automatically considers a set of different models, given the number of profiles, that differ based on how the variances and covariances of the input variables are handled (e.g., fixed, varying, zero). The results that are returned are for the best-fitting of those models based on the BIC. As an example, the following syntax runs the “inspectProfiles” function for 3 profiles and returns the profile membership scores in a dataframe for subsequent analyses. It also produces a frequency table and a number of plots, which we discuss below.
#> [[1]]
#> [[1]][[1]]
#>
#> [[1]][[2]]
#>
#> [[1]][[3]]
#>
#>
#> [[2]]
#> [[2]][[1]]
#>
#> [[2]][[2]]
#>
#> [[2]][[3]]
#>
#>
#> [[3]]
#> [[3]][[1]]
#>
#> [[3]][[2]]
#>
#> [[3]][[3]]
In choosing the number of profiles, we first consider the number of dyads included in each profile and prefer solutions where there are at least about 10% of the dyads in the smallest profile. The rationale here is that we do not want a solution that is driven by a very small portion of our data, since that is unlikely to be a robust finding. In this example, the first part of the output (before the plots) is a frequency table which shows that there are 15 dyads in Profile 1, 5 in Profile 2 and 59 in Profile 3. Thus this is not ideal, since there are only about 8% of the dyads in the second profile, but it is not terrible either.
Next we consider the first plot above showing the projection of the data onto the clustering solution. We prefer solutions where the profiles are clearly separated. In this example, we see adequate separation, but there is more overlap of the green with the other two than is ideal. The second plot shows the means for each parameter for each profile. We can see the profiles are fairly similar for “coord0” and differ the most from each other on “inert0”.
Finally, and most importantly, we consider how meaningful the dynamics of each profile are, which can be established by inspecting the model predicted trajectories of the state variable for each profile. If plotted directly, the Inertia-Coordination model predicts simple exponential growth or decline, which by itself would not be a good representation of most interpersonal temporal data. For example, the time-series emotional data in our example clearly do not show simple exponential growth/decline, but rather fluctuate up and down over time. The model only becomes a good representation of our data when random noise is added at each temporal step. In other words, local exponential growth/decline becomes a fluctuating trajectory when the system is randomly perturbed at each time step. In order to visualize the model predicted dynamics, therefore, we provide 3 exemplar plots for each profile, based on the Inertia-Coordination parameter estimates for that profile in combination with random noise at each time step (see the plots above). We prefer solutions where it can be seen with the naked-eye that there are distinct dynamics characterizing the profiles. In this example, the dynamics are not as clearly separated as would be ideal, but we can see that Profile 1 is characterized by moderate amplitude anti-phase oscillations, Profile 2 is almost flat and anti-phase, and Profile 3 shows a moderate amplitude pattern that hints at an in-phase trajectory.
Due to the importance of the plots for deciding the number of profiles, we also provide a stand-alone function, “inertCoordPlotTraj” that produces the same type of plots as “inspectProfiles”, but with more optional arguments for controlling the presentation. The following syntax provides an example. The three most useful arguments are “minMax,” which allows you to scale the y-axis by setting the minimum and maximum values based on quantiles (in this example, we set them to the 5th and 95th percentiles of the observed state variable), “time_length” which sets the x-axis to however many time steps you would like, and “numPlots” which sets the number of exemplar plots to produce.
Once you have decided on the number of profiles you would like to use, you should re-run “inspectProfiles” with that number set. It will return a dataframe (called “lpaData”" in this example) that includes the profile membership scores to be used as a predictor or outcome of any system variables of interest. The function “makeFullData” then combines the profile membership data with your original dataframe containing the system variables you would like to explore as predictors and outcomes of the profile membership.
Finally, if the number of profiles you have chosen provides a good representation of the data, you should be able to at least fantasize that you can see the prototypical trajectories reflected in the raw data. Of course, the raw data has a lot of noise which makes it hard to see patterns, but a good LPA solution should produce groups of dyads that appear to at least have some similarity to each other and to the prototypical trajectories. In order to assess this, we provide the function “plotDataByProfile” that produces plots of the raw data grouped by profile membership. Here we show a few examples from each profile. Profiles 1 and 3, which were the ones with more dyads in them, show some suggestion of the anti- versus in- phase trajectories predicted by the model.
Predicting the System Variable From the Profiles
The next step in the analysis is to use each dyad’s profile membership to predict the system variable using the “sysVarOut” function. The system variable can be either dyadic (sysVarType = “dyadic”), where both partners have the same score (e.g., relationship length) or individual (sysVarType = “indiv”), where the partners can have different scores (e.g., age). It takes as arguments the name of the dataframe containing the profile membership scores combined with your original dataframe (created by the “makeFullData” function), the name of the column in the dataframe containing the variable you would like to use as the system variable, and whether the system variable is “dyadic” or “individual”. For dyadic system variables, the only predictor is profile membership and the model is a regular regression model since all variables are at the level of the dyad. If the system variable is individual then the model is a random-intercept dyadic model and 3 models are estimated: 1) the main effect of profile membership (“profile”), 2) main effects of profile membership and the distinguishing variable (“profilePlusDist”), and 3) the interaction of profile membership and the distinguishing variable (“profileByDist”). If the system variable is not normally distributed, any of the generalized linear models supported by glm (for dyadic system variables) or glmmPQL (for individual system variables) are available by specifying the “family” distribution (use ?sysVarOut for more information).
For normally distributed system variables, the function returns a list including the lm or lme objects containing the full results for each model (called “models”). These results can be inspected using the usual “summary” function and the models can be compared using the “anova” function. Similarly, for non-normal system variables, the function returns a list of the glm or glmmPQL objects containing the full results for the models. By default, the function also displays a boxplot of the system variable for each profile, histograms of the residuals and plots of the predicted values against observed values for each model, but these can be turned off by setting “printPlots=F”. In this example, we find a significant effect of sex, but no effect of profile (see summary and anova output). The boxplot helps to interpret this result and shows that the system variable (engage) is predicted to be higher for women than men, especially in Profiles 1 & 3.
Predicting Profile Membership From the System Variable
Finally, you can turn the direction of prediction around and use the system variable to predict couples’ profile membership. The function “sysVarIn” accomplishes this. It takes as arguments the name of the dataframe containing the profile membership scores combined with your original dataframe (created by the “makeFullData” function), the name of the column in the dataframe containing the variable you would like to use as the system variable, whether the system variable is “dyadic” or “individual”, the number of profiles, and optional names to use in the plots for the levels of the distinguisher and the system variable. If there are 2 profiles, then binomial regression models are used. If there are more than 2 profiles then multinomial regression is used. For dyadic system variables, a couple’s shared score is the only predictor of their profile membership (the model is called “sysVarMain”). For individual system variables, two models are tested, one with the main effects of both partner’s system variable (“sysVarMain”) and one with the main effects and their interaction (“sysVarInteract”). In both cases an intercept-only model is included as a comparison point (called “base”).
The function returns a list of the full model results and by default produces a boxplot of the system variable for each profile and plots of profile membership against the system variable(s), but these can be turned off by setting printPlots=F. Due to these being either logistic or multinomial models, interpretation is somewhat more complicated than for the sysVarOut results. In this example, it is a multinomial model since we used 3 profiles. To aid intrepretation, we exponentiate the coefficients from the sysVarInteract model to get odds ratios and then calculate p-values for them. The results here do not show any evidence of an interaction effect.
sysIn <- sysVarIn(fullData=fullData, sysVar_name="relstress", sysVarType="indiv", n_profiles=3, dist0name="men", dist1name="women", printPlots=F)
#> # weights: 6 (2 variable)
#> initial value 86.790371
#> final value 55.943729
#> converged
#> # weights: 12 (6 variable)
#> initial value 83.494534
#> iter 10 value 50.698250
#> final value 50.698038
#> converged
#> # weights: 15 (8 variable)
#> initial value 83.494534
#> iter 10 value 48.466257
#> final value 48.452427
#> converged
#> Model names are base, sysVarMain and sysVarInteract
summary(sysIn$models$sysVarInteract)
#> Call:
#> nnet::multinom(formula = profileN ~ sysVar0 * sysVar1, data = basedata,
#> na.action = na.exclude)
#>
#> Coefficients:
#> (Intercept) sysVar0 sysVar1 sysVar0:sysVar1
#> 1 -0.7085458 -0.8507039 -1.9305469 2.089534
#> 2 1.5716090 0.6324667 0.2452257 -1.759158
#>
#> Std. Errors:
#> (Intercept) sysVar0 sysVar1 sysVar0:sysVar1
#> 1 1.1210380 1.552565 2.211627 1.722723
#> 2 0.6742623 1.041372 1.363606 1.662565
#>
#> Residual Deviance: 96.90485
#> AIC: 112.9049
exp(coef(sysIn$models$sysVarInteract))
#> (Intercept) sysVar0 sysVar1 sysVar0:sysVar1
#> 1 0.4923597 0.4271142 0.1450688 8.0811463
#> 2 4.8143883 1.8822478 1.2779097 0.1721898
z <- summary(sysIn$models$sysVarInteract)$coefficients/summary(sysIn$models$sysVarInteract)$standard.errors
p <- (1 - pnorm(abs(z), 0, 1)) * 2
round(p, 3)
#> (Intercept) sysVar0 sysVar1 sysVar0:sysVar1
#> 1 0.527 0.584 0.383 0.225
#> 2 0.020 0.544 0.857 0.290
