Dyadic Multilevel Modeling for Intensive Longitudinal Data

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 their 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 pedagogically 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 dta
  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?

  1. What is the extent of association for the typical male partner and for the typical female partner (i.e., “actor” fixed effects)?
  2. What is the extent of influence of the typical male partner on the typical female partner, and vice versa (i.e., “partner” fixed effects)?
  3. 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 simultaneously 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.

Libraries

Loading libraries used in this script.

library(psych)   #descriptives
library(ggplot2) #plotting
library(nlme)    #mlm analysis
library(dplyr)   #data manipulation
library(tidyr)   #tidy data, reshaping

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

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://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/Bolger_Laurenceau_2013/BL2013_data_dyads_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.

Reshaping the data

To analyze in the multilevel framework, the data need to be reshaped 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.

# Create double-entry data (note that the outcome variables to be combined are commented out of the id.vars list)
BLdyads_reshape <- BLdyads_long %>%
  pivot_longer(cols = c(f_reldis, m_reldis),
               names_to = "variable",
               values_to = "reldis",
               values_drop_na = FALSE)

#create additional variables
BLdyads_reshape <- BLdyads_reshape %>%
        #add the two dummy indicators
  mutate(female = ifelse(variable == "f_reldis", 1, 0),
         male = ifelse(variable == "m_reldis", 1, 0),
         #add another gender indicator
         gender = as.numeric(factor(variable)),
         #create integrated personid variable
         personid = ifelse(female == 1, f_personid, m_personid))

#subset and organize into a new data set for easy viewing and efficiency
#dput(names(BLdyads_reshape)) #to obtain variable list for easy cut&paste
BLdyads_doubleentry <- BLdyads_reshape %>%
  select(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)
## # A tibble: 6 × 16
##   coupleid personid  time time7c gender reldis female  male f_wrkstrs f_wrkstrsc
##      <int>    <int> <int>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>     <int>      <dbl>
## 1        1        1     0  -1.5       1   3.03      1     0         3    0.00952
## 2        1        2     0  -1.5       2   4.46      0     1         3    0.00952
## 3        1        1     1  -1.36      1   4.62      1     0         3    0.00952
## 4        1        2     1  -1.36      2   4.88      0     1         3    0.00952
## 5        1        1     2  -1.21      1   2.85      1     0         3    0.00952
## 6        1        2     2  -1.21      2   4.58      0     1         3    0.00952
## # ℹ 6 more variables: f_wrkstrsc_b <dbl>, f_wrkstrsc_w <dbl>, m_wrkstrs <int>,
## #   m_wrkstrsc <dbl>, m_wrkstrsc_b <dbl>, m_wrkstrsc_w <dbl>

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.

BLdyads_doubleentry %>%
  ggplot(aes(x = reldis)) +
  geom_histogram(fill = "white", color = "black") + 
  labs(x = "Relationship Dissatisfaction") +
  #creating a separate plot for each gender
  facet_grid(. ~ 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.

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

It looks like there is quite a lot of day-to-day variability!

Finally, we examine a histogram of the dyad-level (between-dyad) time-invariant variables f_wrkstrsc_b and m_wrkstrsc_b.

#between-person differences
#female identifying
BLdyads_doubleentry %>%
  filter(time == 0 & female == 1) %>%
  ggplot(aes(x = f_wrkstrsc_b)) +
  geom_histogram(fill = "white", color = "black") + 
  labs(x = "Female Work Stress (dyad-level centered)")

#males identifying
BLdyads_doubleentry %>%
  filter(time == 0 & male == 1) %>%
  ggplot(aes(x = f_wrkstrsc_b)) +
  geom_histogram(fill = "white", color = "black") + 
  labs(x = "Male Work Stress (dyad-level centered)")

3. The Dyadic Multilevel Model

We are now ready to start running the dyadic model!

We’ll construct a model looking at the within-person and between-person associations of relationship dissatisfaction (reldis) with female and male work stressors (f_wrkstrs and m_wrkstrs - which have been centered and separated into within-person and between-person components (f_wrkstrsc_w,m_wrkstrsc_w and f_wrkdstrsc_b,m_wrkdstrsc_b respectively).

We use the two dummy variables (male and female) to turn on and off the parameters. The parameters invoked with \(female\) are associated with the female member of the dyad’s relationship dissatisfaction (within reldis), and parameters invoked with \(male\) are associated with the male member of the dyad’s relationship dissatisfaction (within reldis).

Level 1 of the multilevel model is as follows:

\[RelDis_{it} = \beta_{0F}Female_{it} + \beta_{1F}Female_{it}*FemaleWorkstressWithin_{it} + \beta_{2F}Female_{it}*FemaleWorkstressBetween_{it} + \\ \beta_{3F}Female_{it}*MaleWorkstressWithin_{it} + \beta_{4F}Female_{it}time7c_{it} + \\ \beta_{0M}Male_{it} + \beta_{1M}Male_{it}*MaleWorkstressWithin_{it} + \beta_{2M}Male_{it}*MaleWorkstressBetween_{it} +\\ \beta_{3M}Male_{it}*FemaleWorkstressWithin_{it} + \beta_{4M}Male_{it}time7c_{it} + \\ e_{Fi} +e_{Mi}\]

And Level 2 of the multilevel model: \[\beta_{0F} = \gamma_{00F} + u_{0Fi} \\ \beta_{0M} = \gamma_{00M} + u_{0Mi} \\ \beta_{1F} = \gamma_{10F} + u_{1Fi} \\ \beta_{1M} = \gamma_{10M} + u_{1Mi} \]

Noting that our random effect matrices also expand, \[\mathbf{R} = \left[\begin{array} & \sigma^2_{eF} & \sigma_{eFeM} \\ \sigma_{eFeM} & \sigma^2_{eM} \end{array}\right]\]

where \(\sigma_{eFeM}\) is the residual covariance between male and female relationship dissatisfaction.
and \[\left[\begin{array} {rrrr} \sigma^{2}_{u0F} & \sigma_{u0Fu0M} & \sigma_{u0Fu1F} & \sigma_{u0Fu1M} \\ \sigma_{u0Mu0F} & \sigma^{2}_{u0M} & \sigma_{u0Mu1F} & \sigma_{u0Mu1M} \\ \sigma_{u1Fu0} & \sigma_{u1Fu0M} & \sigma^{2}_{u1F} & \sigma_{u1Fu1M} \\ \sigma_{u1Mu0} & \sigma_{u1Mu0M} & \sigma_{u1Mu1F} & \sigma^{2}_{u1M} \end{array}\right]\]

where the matrix is blocks of between-dyad associations, some among males only, some among females only, and some across genders. These are sample-level, between-dyad relations, and should be interpreted appropriately.

Fitting the full model with both between-dyad and within-dyad actor and partner effects

model1 <- lme(fixed = reldis ~ -1 +
                female + #intercept
                female:f_wrkstrsc_b + female:f_wrkstrsc_w + #actor-effects
                female:m_wrkstrsc_b + female:m_wrkstrsc_w + #partner-effects
                male + #intercept
                male:m_wrkstrsc_b + male:m_wrkstrsc_w + #actor-effects
                male:f_wrkstrsc_b + male:f_wrkstrsc_w + #partner-effects 
                male:time7c + female:time7c,            #fixed effects of time for each gender
              random = ~ -1 + 
                 female + male + #intercepts
                 female:f_wrkstrsc_w + male:m_wrkstrsc_w + #actor-effects
                 female:m_wrkstrsc_w + male:f_wrkstrsc_w | coupleid, #partner-effects
              weights=varIdent(form = ~1 | gender),    # invokes separate sigma^{2}_{e} for each gender
              corr=corSymm(form = ~1 | coupleid/time), # invokes off-diagonal sigma_{e1e2} symmetric structure for level 1 residuals
              data=BLdyads_doubleentry,
              control=lmeControl(maxIter=10000, opt="optim")) 

summary(model1)
## Linear mixed-effects model fit by REML
##   Data: BLdyads_doubleentry 
##        AIC      BIC    logLik
##   12126.74 12354.98 -6027.369
## 
## Random effects:
##  Formula: ~-1 + female + male + female:f_wrkstrsc_w + male:m_wrkstrsc_w +      female:m_wrkstrsc_w + male:f_wrkstrsc_w | coupleid
##  Structure: General positive-definite, Log-Cholesky parametrization
##                     StdDev     Corr                                
## female              0.95915554 female male   fml:f__ ml:m__ fml:m__
## male                1.00937344  0.265                              
## female:f_wrkstrsc_w 0.12822848  0.484  0.083                       
## male:m_wrkstrsc_w   0.17026001 -0.015  0.180  0.531                
## female:m_wrkstrsc_w 0.06787238  0.634  0.505  0.557  -0.071        
## male:f_wrkstrsc_w   0.05179949 -0.042  0.149  0.703   0.566  0.387 
## Residual            0.99565514                                     
## 
## Correlation Structure: General
##  Formula: ~1 | coupleid/time 
##  Parameter estimate(s):
##  Correlation: 
##   1    
## 2 0.071
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | gender 
##  Parameter estimates:
##         1         2 
## 1.0000000 0.8743963 
## Fixed effects:  reldis ~ -1 + female + female:f_wrkstrsc_b + female:f_wrkstrsc_w +      female:m_wrkstrsc_b + female:m_wrkstrsc_w + male + male:m_wrkstrsc_b +      male:m_wrkstrsc_w + male:f_wrkstrsc_b + male:f_wrkstrsc_w +      male:time7c + female:time7c 
##                         Value Std.Error   DF  t-value p-value
## female               4.640480 0.0991961 4089 46.78088  0.0000
## male                 5.095043 0.1036419 4089 49.16007  0.0000
## female:f_wrkstrsc_b  0.768676 0.4396487 4089  1.74839  0.0805
## female:f_wrkstrsc_w  0.161576 0.0253286 4089  6.37918  0.0000
## female:m_wrkstrsc_b  0.470131 0.4198417 4089  1.11978  0.2629
## female:m_wrkstrsc_w -0.007121 0.0227347 4089 -0.31321  0.7541
## m_wrkstrsc_b:male    0.072740 0.4504006 4089  0.16150  0.8717
## m_wrkstrsc_w:male    0.109341 0.0256842 4089  4.25711  0.0000
## f_wrkstrsc_b:male    0.630237 0.4727048 4089  1.33326  0.1825
## f_wrkstrsc_w:male    0.015914 0.0198440 4089  0.80195  0.4226
## male:time7c          0.011914 0.0222100 4089  0.53641  0.5917
## female:time7c       -0.022295 0.0252508 4089 -0.88294  0.3773
##  Correlation: 
##                     female male   fml:f_wrkstrsc_b fml:f_wrkstrsc_w
## male                 0.257                                         
## female:f_wrkstrsc_b  0.089  0.023                                  
## female:f_wrkstrsc_w  0.236  0.040 -0.006                           
## female:m_wrkstrsc_b -0.084 -0.021  0.110            0.002          
## female:m_wrkstrsc_w  0.182  0.147 -0.006            0.103          
## m_wrkstrsc_b:male   -0.020 -0.086  0.029            0.005          
## m_wrkstrsc_w:male   -0.010  0.117  0.001            0.180          
## f_wrkstrsc_b:male    0.022  0.091  0.248           -0.002          
## f_wrkstrsc_w:male   -0.010  0.037  0.001            0.155          
## male:time7c          0.001  0.015  0.000            0.000          
## female:time7c        0.018  0.001 -0.002           -0.014          
##                     fml:m_wrkstrsc_b fml:m_wrkstrsc_w m_wrkstrsc_b:
## male                                                               
## female:f_wrkstrsc_b                                                
## female:f_wrkstrsc_w                                                
## female:m_wrkstrsc_b                                                
## female:m_wrkstrsc_w  0.000                                         
## m_wrkstrsc_b:male    0.247           -0.008                        
## m_wrkstrsc_w:male    0.000            0.035           -0.004       
## f_wrkstrsc_b:male    0.028           -0.003            0.111       
## f_wrkstrsc_w:male   -0.001            0.037            0.003       
## male:time7c          0.000            0.000            0.001       
## female:time7c       -0.001           -0.005            0.003       
##                     m_wrkstrsc_w: f_wrkstrsc_b: f_wrkstrsc_w: ml:tm7
## male                                                                
## female:f_wrkstrsc_b                                                 
## female:f_wrkstrsc_w                                                 
## female:m_wrkstrsc_b                                                 
## female:m_wrkstrsc_w                                                 
## m_wrkstrsc_b:male                                                   
## m_wrkstrsc_w:male                                                   
## f_wrkstrsc_b:male   -0.002                                          
## f_wrkstrsc_w:male    0.112        -0.002                            
## male:time7c         -0.004         0.000        -0.014              
## female:time7c       -0.001         0.000        -0.001         0.071
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -3.364402e+00 -6.528970e-01 -4.230411e-06  6.447156e-01  3.385325e+00 
## 
## Number of Observations: 4200
## Number of Groups: 100

Note that the residual structure is not quite the same as in Bolger & Laurenceau (2013). It is not clear exactly how to get lme() structure to match. There is not a straightforward way, but it seems like it is possible. We just haven’t found the exact combination of structures to use yet.

But, here is another close variant.

model2 <- lme(fixed = reldis ~ -1 +
                female + #intercept
                female:f_wrkstrsc_b + female:f_wrkstrsc_w + #actor-effects
                female:m_wrkstrsc_b + female:m_wrkstrsc_w + #partner-effects
                male + #intercept
                male:m_wrkstrsc_b + male:m_wrkstrsc_w + #actor-effects
                male:f_wrkstrsc_b + male:f_wrkstrsc_w + #partner-effects
                male:time7c + female:time7c,            #fixed effects of time for each gender
              random = ~ -1 +
                female + male + #intercepts
                female:f_wrkstrsc_w + male:m_wrkstrsc_w | coupleid, #actor-effects
                female:m_wrkstrsc_w + male:f_wrkstrsc_w | coupleid, #partner-effects
              weights=varIdent(form = ~1 | gender), #invokes separate sigma^{2}_{e} for each gender
              corr=corAR1(form = ~1 | coupleid/gender/time), #invokes an AR(1) structure for level 1 residuals
              data=BLdyads_doubleentry,
              control=list(maxIter=1000, opt="optim")) 
## Warning in female:m_wrkstrsc_w: numerical expression has 4200 elements: only
## the first used
## Warning in male:f_wrkstrsc_w: numerical expression has 4200 elements: only the
## first used
summary(model2)
## Linear mixed-effects model fit by REML
##   Data: BLdyads_doubleentry 
##   Subset: female:m_wrkstrsc_w + male:f_wrkstrsc_w | coupleid 
##        AIC      BIC    logLik
##   12122.27 12280.77 -6036.136
## 
## Random effects:
##  Formula: ~-1 + female + male + female:f_wrkstrsc_w + male:m_wrkstrsc_w | coupleid
##  Structure: General positive-definite, Log-Cholesky parametrization
##                     StdDev    Corr                
## female              0.9595876 female male   fml:__
## male                1.0091299  0.271              
## female:f_wrkstrsc_w 0.1283527  0.474  0.074       
## male:m_wrkstrsc_w   0.1675532 -0.012  0.182  0.509
## Residual            0.9980118                     
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | coupleid/gender/time 
##  Parameter estimate(s):
## Phi 
##   0 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | gender 
##  Parameter estimates:
##         1         2 
## 1.0000000 0.8740942 
## Fixed effects:  reldis ~ -1 + female + female:f_wrkstrsc_b + female:f_wrkstrsc_w +      female:m_wrkstrsc_b + female:m_wrkstrsc_w + male + male:m_wrkstrsc_b +      male:m_wrkstrsc_w + male:f_wrkstrsc_b + male:f_wrkstrsc_w +      male:time7c + female:time7c 
##                         Value Std.Error   DF  t-value p-value
## female               4.642207 0.0992790 4089 46.75919  0.0000
## male                 5.097476 0.1036463 4089 49.18148  0.0000
## female:f_wrkstrsc_b  0.831764 0.4466742 4089  1.86213  0.0627
## female:f_wrkstrsc_w  0.160809 0.0253658 4089  6.33962  0.0000
## female:m_wrkstrsc_b  0.448591 0.4267527 4089  1.05117  0.2932
## female:m_wrkstrsc_w -0.005117 0.0217190 4089 -0.23558  0.8138
## m_wrkstrsc_b:male    0.042081 0.4553011 4089  0.09242  0.9264
## m_wrkstrsc_w:male    0.109203 0.0255261 4089  4.27808  0.0000
## f_wrkstrsc_b:male    0.709805 0.4774740 4089  1.48658  0.1372
## f_wrkstrsc_w:male    0.013253 0.0191494 4089  0.69211  0.4889
## male:time7c          0.011450 0.0222361 4089  0.51495  0.6066
## female:time7c       -0.024637 0.0252996 4089 -0.97379  0.3302
##  Correlation: 
##                     female male   fml:f_wrkstrsc_b fml:f_wrkstrsc_w
## male                 0.260                                         
## female:f_wrkstrsc_b  0.090  0.025                                  
## female:f_wrkstrsc_w  0.232  0.037 -0.004                           
## female:m_wrkstrsc_b -0.085 -0.023  0.108           -0.004          
## female:m_wrkstrsc_w -0.001  0.000 -0.003            0.019          
## m_wrkstrsc_b:male   -0.023 -0.087  0.029            0.000          
## m_wrkstrsc_w:male   -0.008  0.116  0.001            0.169          
## f_wrkstrsc_b:male    0.024  0.092  0.266            0.000          
## f_wrkstrsc_w:male    0.000 -0.001  0.000            0.000          
## male:time7c          0.000  0.015  0.000            0.001          
## female:time7c        0.018  0.000 -0.002           -0.014          
##                     fml:m_wrkstrsc_b fml:m_wrkstrsc_w m_wrkstrsc_b:
## male                                                               
## female:f_wrkstrsc_b                                                
## female:f_wrkstrsc_w                                                
## female:m_wrkstrsc_b                                                
## female:m_wrkstrsc_w  0.010                                         
## m_wrkstrsc_b:male    0.266            0.001                        
## m_wrkstrsc_w:male    0.002           -0.001           -0.004       
## f_wrkstrsc_b:male    0.029            0.000            0.109       
## f_wrkstrsc_w:male   -0.002            0.005            0.004       
## male:time7c         -0.001            0.000            0.001       
## female:time7c       -0.004           -0.005            0.000       
##                     m_wrkstrsc_w: f_wrkstrsc_b: f_wrkstrsc_w: ml:tm7
## male                                                                
## female:f_wrkstrsc_b                                                 
## female:f_wrkstrsc_w                                                 
## female:m_wrkstrsc_b                                                 
## female:m_wrkstrsc_w                                                 
## m_wrkstrsc_b:male                                                   
## m_wrkstrsc_w:male                                                   
## f_wrkstrsc_b:male   -0.002                                          
## f_wrkstrsc_w:male    0.013        -0.001                            
## male:time7c         -0.004         0.000        -0.015              
## female:time7c       -0.001         0.000         0.000        -0.001
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.362712205 -0.648833344  0.000713101  0.657788257  3.322554053 
## 
## Number of Observations: 4200
## Number of Groups: 100

Although the structure is not exactly matched to the example in the book, the ones we are fitting are among the structures used in the literature. Across variants, the interpretations do not change in meaningful ways, suggesting that looking across the range of possibilities is fine.

Brief Interpretation

The results for the dyadic model (model1) are described briefly below.

Fixed effects: On an average day, males and females relationship dissatisfaction is 5.1 and 4.64 respectively, on a 0 - 10 scale.

On days when there is an additional work stressor than usual, relationship dissatisfaction is expected to increase for both genders (males = 0.11, females = 0.16); these are sometimes referred to as the “actor effects”.

On days when a person’s partner has an additional work stressor than usual, relationship dissatisfaction is not expected to change for either gender (males = 0.02, females = -0.01); these are sometimes referred to as the “partner effects”.

Random effects: The standard deviation around males and females average relationship dissatisfaction is 1.01 and 0.96, respectively. Additionally, males and females relationship dissatisfaction reports are correlated 0.26, indicating that members of the same dyad often have relatively similar reports of relationship dissatisfaction. The standard deviation around males and females reactivity to work stressors (i.e., slope of work stressors) is 0.17 and 0.13, respectively.

A more thorough interpretation and formal write-up of the results can be found on pages 165 - 171 of Bolger and Laurenceau (2013) Chapter 8: Design and Analysis of Intensive Longitudinal Longitudinal Study of Distingiushable Dyads.

4. Cautions

While we only provide a brief description of the multivariate multilevel model for repeated measures dyadic data, we want to a highlight a few considerations when using this model.

  1. The model demonstrated here is for distinguishable dyads - i.e., any pair of people or variables that can be separated by a chosen feature.

  2. The model demonstrated here assumed a linear association between work stressors and relationship dissatisfaction. Additional parameters would need to be added to the model to test non-linear relationships.

  3. The model (model1) demonstrated here did not allow for an autoregressive error structure - i.e., correlations between error terms from day-to-day within dyad-member were ignored (but see model2 for adding an autoregressive structure).

5. Conclusion

In sum, this tutorial follows and builds upon the example provided in Chapter 8: Design and Analysis of Intensive Longitudinal Longitudinal Study of Distingiushable Dyads of Bolger and Laurenceau (2013). We provide brief explanation of (1) the underlying model, (2) the code to run these models in R, and (3) the interpretation of the results. More detailed information about these (and related) analyses can be found in Bolger and Laurenceau (2013), Laurenceau and Bolger (2005), and Kenny, Kashy, and Cook (2006).

Thanks for dyading!

Citations

Bolger, N., & Laurenceau, J.-P. (2013). Chapter 8: Design and Analysis of Intensive Longitudinal Longitudinal Study of Distingiushable Dyads. In Intensive longitudinal methods: An introduction to diary and experience sampling research (pp. xv, 256). The Guilford Press.

Kenny, D. A., Kashy, D. A., & Cook, W. L. (2006). Dyadic data analysis. The Guilford Press.

Laurenceau, J.-P., & Bolger, N. (2005). Using Diary Methods to Study Marital and Family Processes. Journal of Family Psychology, 19(1), 86–97. https://doi.org/10.1037/0893-3200.19.1.86

Pinheiro, J. C., Bates, D., & R Core Team. (2025). nlme: Linear and Nonlinear Mixed Effects Models (Version 3.1-167). https://doi.org/10.1007/b98882

R Core Team. (2024). R: A Language and Environment for Statistical Computing. Foundation for Statistical Computing. https://www.R-project.org/

Revelle, W. (2024). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University. https://CRAN.R-project.org/package=psych

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag. https://ggplot2.tidyverse.org/

Wickham, H., François, R., Henry, L., Müller, K., & Vaughan, D. (2023). dplyr: A Grammar of Data Manipulation (Version 1.1.4). https://CRAN.R-project.org/package=dplyr

Wickham, H., Vaughan, D., & Girlich, M. (2024). tidyr: Tidy Messy Data (Version 1.3.1). https://CRAN.R-project.org/package=tidyr