APA-ATI Intensive Longitudinal Data: Session M
Dyadic Multilevel Modeling
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
- Introduction to the research questions, model, and data.
- Plotting the Data.
- The Dyadic Multilevel Model.
- Cautions.
- 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)?
- What is the extent of association for the typical male partner and for the typical female partner (i.e., “actor” fixed 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 descriptivesLoading 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.1428571These 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.1428571The 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")
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-peson differences
#females
ggplot(data = subset(BLdyads_doubleentry, time == 0 & female == 1),
aes(x = f_wrkstrsc_b)) +
geom_histogram(fill = "white", color = "black") +
labs(x = "Female Work Stress (dyad-level centered)")
#males
ggplot(data = subset(BLdyads_doubleentry, time == 0 & male == 1), 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 anly, 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
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=corCompSymm(form = ~1 | coupleid/time), # invokes off-diagonal sigma_{e1e2}
data=BLdyads_doubleentry,
control=lmeControl(maxIter=10000, opt="optim"))
summary(model1)## Linear mixed-effects model fit by REML
## Data: BLdyads_doubleentry
## AIC BIC logLik
## 12112.56 12328.14 -6022.281
##
## 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.95867958 female male fml:f__ ml:m__ fml:m__
## male 1.00870418 0.267
## female:f_wrkstrsc_w 0.12802327 0.481 0.083
## male:m_wrkstrsc_w 0.16955643 -0.019 0.180 0.517
## female:m_wrkstrsc_w 0.06872053 0.636 0.494 0.581 -0.067
## male:f_wrkstrsc_w 0.05358497 -0.028 0.147 0.720 0.554 0.413
## Residual 0.99546796
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | coupleid/time
## Parameter estimate(s):
## Rho
## 0.07005602
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | gender
## Parameter estimates:
## 1 2
## 1.0000000 0.8745067
## 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
## Value Std.Error DF t-value p-value
## female 4.641936 0.0991324 4091 46.82562 0.0000
## male 5.094219 0.1035640 4091 49.18909 0.0000
## female:f_wrkstrsc_b 0.763597 0.4395660 4091 1.73716 0.0824
## female:f_wrkstrsc_w 0.161184 0.0253139 4091 6.36742 0.0000
## female:m_wrkstrsc_b 0.471148 0.4197285 4091 1.12251 0.2617
## female:m_wrkstrsc_w -0.007262 0.0227571 4091 -0.31909 0.7497
## m_wrkstrsc_b:male 0.072782 0.4502500 4091 0.16165 0.8716
## m_wrkstrsc_w:male 0.109398 0.0256369 4091 4.26720 0.0000
## f_wrkstrsc_b:male 0.631484 0.4725361 4091 1.33637 0.1815
## f_wrkstrsc_w:male 0.016108 0.0198907 4091 0.80985 0.4181
## Correlation:
## female male fml:f_wrkstrsc_b fml:f_wrkstrsc_w
## male 0.259
## female:f_wrkstrsc_b 0.089 0.023
## female:f_wrkstrsc_w 0.234 0.040 -0.006
## female:m_wrkstrsc_b -0.084 -0.021 0.110 0.002
## female:m_wrkstrsc_w 0.185 0.146 -0.006 0.108
## m_wrkstrsc_b:male -0.021 -0.086 0.029 0.005
## m_wrkstrsc_w:male -0.012 0.116 0.001 0.174
## f_wrkstrsc_b:male 0.022 0.091 0.250 -0.002
## f_wrkstrsc_w:male -0.007 0.038 0.001 0.160
## 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.001
## m_wrkstrsc_b:male 0.248 -0.008
## m_wrkstrsc_w:male 0.000 0.035 -0.004
## f_wrkstrsc_b:male 0.029 -0.003 0.111
## f_wrkstrsc_w:male -0.001 0.040 0.003
## m_wrkstrsc_w: f_wrkstrsc_b:
## 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.113 -0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.3777584357 -0.6550602915 -0.0004825079 0.6439829383 3.4130866726
##
## Number of Observations: 4200
## Number of Groups: 100Note 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
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), #this invokes separate sigma^{2}_{e} for each gender
corr=corAR1(form = ~1 | coupleid/gender/time), #this invokes an AR structure
data=BLdyads_doubleentry,
control=list(maxIter=1000, opt="optim"))
summary(model2)## Linear mixed-effects model fit by REML
## Data: BLdyads_doubleentry
## AIC BIC logLik
## 12121.65 12337.23 -6026.827
##
## 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.95747245 female male fml:f__ ml:m__ fml:m__
## male 1.00844599 0.270
## female:f_wrkstrsc_w 0.13343214 0.451 0.083
## male:m_wrkstrsc_w 0.16989649 -0.018 0.181 0.483
## female:m_wrkstrsc_w 0.06623571 0.694 0.531 0.625 0.130
## male:f_wrkstrsc_w 0.06025820 0.004 0.110 0.788 0.549 0.448
## Residual 0.99515353
##
## 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.8741739
## 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
## Value Std.Error DF t-value p-value
## female 4.641136 0.0990138 4091 46.87362 0.0000
## male 5.094184 0.1035341 4091 49.20298 0.0000
## female:f_wrkstrsc_b 0.737799 0.4396664 4091 1.67809 0.0934
## female:f_wrkstrsc_w 0.161573 0.0255943 4091 6.31283 0.0000
## female:m_wrkstrsc_b 0.481079 0.4196549 4091 1.14637 0.2517
## female:m_wrkstrsc_w -0.007096 0.0226753 4091 -0.31294 0.7543
## m_wrkstrsc_b:male 0.071059 0.4496717 4091 0.15803 0.8744
## m_wrkstrsc_w:male 0.109292 0.0256551 4091 4.26005 0.0000
## f_wrkstrsc_b:male 0.628226 0.4719580 4091 1.33111 0.1832
## f_wrkstrsc_w:male 0.016005 0.0200786 4091 0.79712 0.4254
## Correlation:
## female male fml:f_wrkstrsc_b fml:f_wrkstrsc_w
## male 0.258
## female:f_wrkstrsc_b 0.089 0.022
## female:f_wrkstrsc_w 0.227 0.041 -0.005
## female:m_wrkstrsc_b -0.084 -0.021 0.111 0.003
## female:m_wrkstrsc_w 0.195 0.152 -0.006 0.114
## m_wrkstrsc_b:male -0.020 -0.086 0.028 0.005
## m_wrkstrsc_w:male -0.012 0.117 0.001 0.167
## f_wrkstrsc_b:male 0.022 0.091 0.245 -0.002
## f_wrkstrsc_w:male 0.001 0.032 0.001 0.130
## 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.002
## m_wrkstrsc_b:male 0.244 -0.008
## m_wrkstrsc_w:male 0.001 0.027 -0.003
## f_wrkstrsc_b:male 0.028 -0.003 0.111
## f_wrkstrsc_w:male -0.001 0.044 0.003
## m_wrkstrsc_w: f_wrkstrsc_b:
## 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.123 -0.001
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.3759567188 -0.6561709432 0.0006555097 0.6431144091 3.3927004171
##
## Number of Observations: 4200
## Number of Groups: 100Although 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 possibilites 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.09 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, p = 0.00; females = 0.16, p = 0.00); 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, p = 0.42; females = -0.01, p = 0.75); 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.27, 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. The way in which dyad members’ relationship dissatisfaction changes in response to work stressors is not very similar, as indicated by the -0.02 correlation.
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 consierations when using this model.
The model demonstrated here is for distinguishable dyads - i.e., any pair of people or variables that can be separated by a chosen feature.
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.
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!