Overview

Repeated measures data, often obtained from experience sampling or daily diary studies, require analytic methods that accommodate the inherent nesting in the data. One special case is repeated measures data obtained from dyads (e.g., couples, parent/child). Dyadic data are structured such that repeated measures are nested within persons, and persons are nested within dyads. The dyads (higher-level units) are assumed to be independent, while the members of the dyad and the repeated measures are non-independent. Multivariate multilevel modeling is one suitable technique to effectively handle this type of data.

In this tutorial, we follow the example from Bolger and Laurenceau (2013) Chapter 8: Design and Analysis of Intensive Longitudinal Longitudinal Study of Distingiushable Dyads. Specifically, we will use the simulated dyadic process data set (p. 150). The data were simulated to represent 100 dual-career heterosexual couples where each partner provided diary reports twice daily over the course of 21 consecutive days. The first report is an end-of-workday report that includes (a) the number of stressors that occurred at work, and (b) work dissatisfaction. The data are complete and clean - with Bolger and Laurenceau’s goal being to “create a pedagologically uncomplicated data set.”

Note that we are using distinguishable dyads in this example. Distinguishability is typically determined conceptually based on a stable differentiating characteristic (e.g., gender, age, role).

Also, “dyad units” can be persons or variables, where in the latter case the framework is used to accommodate study of two variables (e.g., emotion and behavior) nested within persons. Procedures for preparing bivariate (and multivariate) data and dyadic (and multiperson) data are the same.

Outline

  1. Introduction to the research questions, model, and data.
  2. Plotting the Data.
  3. The Dyadic Multilevel Model.
  4. Cautions.
  5. Conclusion.

1. Introduction to the Research Questions, Model, and Data.

The Research Questions.

We are going to address:

  • Is the number of daily work stressors associated with end-of-day relationship satisfaction?
    • What is the extent of association for the typical male partner and for the typical female partner (i.e., “actor” fixed effects)?
    • What is the extent of influence of the typical male partner on the typical female partner, and vice versa (i.e., “partner” fixed effects)?
    • Is there heterogeneity in the strength of the association across male partners in couples and across female partners in couples (random effects)?

Notice that, so far, the questions are stated separately for the two types of dyad members (males and females). The dyadic longitudinal model provides for another type of question - questions related to non-independence.

  • Are dyad members’ relationship satisfactions on any given day related, after accounting for other explanatory variables?

Here, these relations manifest as correlations/covariances between intercepts, slopes, and residuals.

The Modeling Enterprise.

The multilevel model is designed as a model with a univariate outcome. The ability to model multiple outcomes simultaneoulsy used to be a distinguishing feature of structural equation models (SEM). However, researches discovered that the multilevel model can be adapted for examination of multivariate outcomes quite easily. One simply has to “trick” the program into thinking that two (or more) variables are one variable.

The Data.

The data are organized as follows:

  • There should be N (number of individuals)*measurement occasions rows per dyad in the data set. In this case, we should have a data set with 4,200 rows.

  • Columns:
    • Couple ID
    • Person ID (e.g., 1 = partner 1, 2 = partner 2)
    • Time (e.g., day within the daily diary study)
    • Centered version of time (optional)
    • Gender (or whatever feature is distinguishing the dyad; in this case, 0 = male and 1 = female)
    • An indicator variable for each partner (e.g., husband, wife) - dichotomous (0/1)
    • Outcome variable (in this case, “reldis”)
    • Predictor variable (in this case, “wrkstrs”)
    • Centered version of the predictor variable (“wrkstrsc”)
    • Trait component of the predictor variable (“wrkstrsb”)
    • State component of the predictor variable (“wrkstrsw”)

    • State component of the predictor variable for each gender (which we will need to create)

    • for a total of 14 columns in this data set

Preliminaries

Loading Libraries

Loading libraries used in this script.

library(dplyr) #for data manipulation
library(reshape2) #for data manipulation
library(ggplot2) #for plotting
library(nlme) #for mlm analysis
library(psych) #for descriptives

Loading Data

These data are a modified version of the data that accompanies Bolger and Laurenceau (2013) - which we made in order to (a) demonstrate the data manipulation procedures and (b) accommodate examination of both “actor” and “partner” effects. (The book only does actor effects).

# Set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/BLdyads_long.csv"

# Read in the .csv file using the url() function
BLdyads_long <- read.csv(file=url(filepath),header=TRUE)
head(BLdyads_long)
##   coupleid time     time7c f_personid f_reldis f_wrkstrs f_wrkstrsc
## 1        1    0 -1.5000000          1     3.03         3  0.0095238
## 2        1    1 -1.3571429          1     4.62         3  0.0095238
## 3        1    2 -1.2142857          1     2.85         3  0.0095238
## 4        1    3 -1.0714286          1     6.40         4  1.0095238
## 5        1    4 -0.9285714          1     2.54         1 -1.9904762
## 6        1    5 -0.7857143          1     5.16         2 -0.9904762
##   f_wrkstrsc_b f_wrkstrsc_w m_personid m_reldis m_wrkstrs m_wrkstrsc
## 1   -0.3238095    0.3333333          2     4.46         3  0.0095238
## 2   -0.3238095    0.3333333          2     4.88         3  0.0095238
## 3   -0.3238095    0.3333333          2     4.58         3  0.0095238
## 4   -0.3238095    1.3333333          2     4.49         1 -1.9904762
## 5   -0.3238095   -1.6666667          2     5.04         3  0.0095238
## 6   -0.3238095   -0.6666667          2     4.87         3  0.0095238
##   m_wrkstrsc_b m_wrkstrsc_w
## 1   -0.1333333    0.1428571
## 2   -0.1333333    0.1428571
## 3   -0.1333333    0.1428571
## 4   -0.1333333   -1.8571429
## 5   -0.1333333    0.1428571
## 6   -0.1333333    0.1428571

These data are long-format data with multiple rows per dyad and both the females’ and males’ scores on the same row.

Reformatting the data

To analyze in the multilevel framework, the data need to be melted so that the outcome variable for females (f_reldis) and males (m_reldis) are combined into a single variable along with two dummy indicators (female and male) that indicate which dyad member’s score is in each row of data. The data thus become “double-entry” or “stacked”: two records for each observation. The data file is twice as long, and we structure the model to “turn on” and “turn off” the double records to invoke parameter estimation for each variable.

Melting into double-entry data.

# Create double-entry data (note that the outcome variables to be combined are commented out of the id.vars list)
BLdyads_melt <- melt(data = BLdyads_long,
                     id.vars = c("coupleid","time","time7c",
                                 "f_personid",#"f_reldis",
                                 "f_wrkstrs","f_wrkstrsc","f_wrkstrsc_b","f_wrkstrsc_w",
                                 "m_personid",#"m_reldis",
                                 "m_wrkstrs","m_wrkstrsc","m_wrkstrsc_b","m_wrkstrsc_w"),
                     na.rm = FALSE, 
                     value.name = "reldis") #naming new variable

# Reorder rows for convenience 
BLdyads_melt <- BLdyads_melt[order(BLdyads_melt$coupleid, BLdyads_melt$time, BLdyads_melt$var), ]

# Add the two dummy indicators 
BLdyads_melt$female <- ifelse(BLdyads_melt$variable == "f_reldis", 1, 0)
BLdyads_melt$male <- ifelse(BLdyads_melt$variable == "m_reldis", 1, 0)

# Add another gender indicator
BLdyads_melt$gender <- as.numeric(factor(BLdyads_melt$variable))

# Create integrated personid variable
BLdyads_melt$personid <- ifelse(BLdyads_melt$female == 1, BLdyads_melt$f_personid, BLdyads_melt$m_personid)

# Subset and organize into a new data set for easy viewing and efficiency
#dput(names(BLdyads_melt)) #to obtain variable list for easy cut&paste
BLdyads_doubleentry <- BLdyads_melt[ ,c("coupleid", "personid","time", "time7c",
                                "gender", "reldis","female", "male",
                                "f_wrkstrs", "f_wrkstrsc", "f_wrkstrsc_b", "f_wrkstrsc_w",
                                "m_wrkstrs", "m_wrkstrsc", "m_wrkstrsc_b", "m_wrkstrsc_w")]

head(BLdyads_doubleentry)
##      coupleid personid time    time7c gender reldis female male f_wrkstrs
## 1           1        1    0 -1.500000      1   3.03      1    0         3
## 2101        1        2    0 -1.500000      2   4.46      0    1         3
## 2           1        1    1 -1.357143      1   4.62      1    0         3
## 2102        1        2    1 -1.357143      2   4.88      0    1         3
## 3           1        1    2 -1.214286      1   2.85      1    0         3
## 2103        1        2    2 -1.214286      2   4.58      0    1         3
##      f_wrkstrsc f_wrkstrsc_b f_wrkstrsc_w m_wrkstrs m_wrkstrsc
## 1     0.0095238   -0.3238095    0.3333333         3  0.0095238
## 2101  0.0095238   -0.3238095    0.3333333         3  0.0095238
## 2     0.0095238   -0.3238095    0.3333333         3  0.0095238
## 2102  0.0095238   -0.3238095    0.3333333         3  0.0095238
## 3     0.0095238   -0.3238095    0.3333333         3  0.0095238
## 2103  0.0095238   -0.3238095    0.3333333         3  0.0095238
##      m_wrkstrsc_b m_wrkstrsc_w
## 1      -0.1333333    0.1428571
## 2101   -0.1333333    0.1428571
## 2      -0.1333333    0.1428571
## 2102   -0.1333333    0.1428571
## 3      -0.1333333    0.1428571
## 2103   -0.1333333    0.1428571

The data set BLdyads_doubleentry is now ready for analysis.

2. Plotting the Data.

Before we begin running the models, it is always a good idea to look at the data.

We start with examining the distribution of our outcome variable end-of-day relationship dissatisfaction, reldis. Let’s look at the histogram by gender.

ggplot(data = BLdyads_doubleentry, aes(x = reldis)) +
  geom_histogram(fill = "white", color = "black") + 
  labs(x = "Relationship Dissatisfaction") +
  facet_grid(. ~ gender) # creating a separate plot for each gender

The outcome variable for each gender looks approximately normally distributed, which is good news for when we run our models.

Next, let’s plot a few dyads’ reports of relationship dissatisfaction through the course of the study. Since the the predictor variable has already been split into time-varying (state) and time-invariant (trait) components, we use the time-varying predictor wrkstrsc_w as the “background” context variable.

ggplot(data = subset(BLdyads_doubleentry, coupleid <= 9), aes(x = time, group = personid), legend = FALSE) +
  geom_rect(mapping = aes(xmin = time-.5, xmax = time+.5, ymin = 0, ymax = 5, fill = f_wrkstrsc_w), alpha = 0.6) + # creating rectangles in the background of the plot colored by female work stressors
  geom_rect(mapping = aes(xmin = time-.5, xmax = time+.5, ymin = 5, ymax = 10, fill = m_wrkstrsc_w), alpha = 0.6) + # creating rectangles in the background of the plot colored by male work stressors
  geom_point(aes(x = time, y = reldis, color = factor(gender)), shape = 17, size = 3) + # creating a different colored point for each gender
  geom_line(aes(x = time, y = reldis, color = factor(gender)), lty = 1, size=1) + # creating a different colored line for each gender
  xlab("Time") + 
  ylab("Relationship Dissatisfaction") + ylim(0, 10) +
  scale_x_continuous(breaks=seq(0, 20, by = 5)) + 
  facet_wrap( ~ coupleid) +# creating a separate plot for each dyad
  theme(legend.position = "none")