Chapter 5 - Multivariate Multilevel Model for Examining Intraindividual Coupling (Part 2)
Overview
In the previous tutorial we covered how the multivariate multilevel modeling is used to examine dyadic repeated measures data (a version of APIM. In this tutorial, we develop a parallel example that uses teh same model set-up to examine intraindividual coupling. Specifically, we use the same dummy-variable “trick” to model relations of two outcome variables simultaneously. Here, distinguishability is given by the variables inherently.
Outline
- Introduction to the research questions, model, and data.
- Plotting and Restructuring the data.
- The multilevel model.
- Conclusion.
Prelim - Loading libraries used in this script.
library(ggplot2)
library(nlme)
library(psych)
library(plyr)
library(DataCombine)
library(reshape2)1. Introduction to the Research Questions, Model, and Data.
The Research Questions.
For this example, we use the AMIB data. Here, our “dyadic data”, rather than being being about two individuals nested in dyads, are about two variables nested within an individual.
Our interest in this example is to model the intraindividual coupling of positive and negative emotion.
We are going to address: How do the emotions carry-over from day to day? This is being called emotional intertia in the literature.
How does postive emotion influence negative emotions? (protective factor) and, vice versa,
How does negative emotion influence positive emotions? (risk factor)
Basically, this is a cross-lag panel model - except that we are working with diary data, and that means we can prioritize within-person associations.
The Modeling Enterprise.
The basic multilevel model is designed as a model with a univariate outcome. Here again, we “trick” the program into thinking that two (or more) variables are one variable.
The Data.
We use the AMIB data - working from a standard long file to obtain a “double-entry” stacked file with lagged variables, and dummy variables that will be used to “turn on” and “turn off” the double records and invoke parameter estimation for each variable. In the dyadic example, the data were already prepared. Here, we walk through the data preparation process.
Load data and needed libraries.
#set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/AMIBbrief_raw_daily1.csv"
#read in the .csv file using the url() function
daily <- read.csv(file=url(filepath),header=TRUE)
#Little bit of clean-up
var.names.daily <- tolower(colnames(daily))
colnames(daily)<-var.names.daily
#subsetting data
daily1 <- daily[,c("id","day","posaff","negaff")]
# Examine first few rows of the data set
names(daily1)## [1] "id" "day" "posaff" "negaff"2. Plotting & Restructuring the Data.
Plotting the data
Before we begin running our models, it is always a good idea to look at our data.
We look at our two outcomes, posaff and negaff.
#sample descriptives
describe(daily1$negaff)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 1441 2.45 1.04 2.2 2.34 1.04 1 6.9 5.9 0.96 0.77
## se
## X1 0.03describe(daily1$posaff)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 1441 4.12 1.1 4.2 4.15 1.19 1 7 6 -0.25 -0.33
## se
## X1 0.03#histogram
ggplot(data=daily1, aes(x=negaff)) +
geom_histogram(fill="white", color="black") +
labs(x = "Negative Affect")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data=daily1, aes(x=posaff)) +
geom_histogram(fill="white", color="black") +
labs(x = "Positive Affect")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Intraindividual plots
#setting up for changing ids
i <- 106
ggplot(data = subset(daily1, id <=i), aes(x=day, group=id), legend=FALSE) +
#geom_rect(mapping=aes(xmin=day-.5, xmax=day+.5, ymin=0, ymax=10, fill=wrkstrscw), alpha=0.6) +
geom_point(aes(x=day,y = posaff), color="blue", shape=17, size=2) +
geom_line(aes(x=day,y = posaff), color="blue", lty=1, size=1) +
geom_point(aes(x=day,y = negaff), color="red", shape=17, size=2) +
geom_line(aes(x=day,y = negaff), color="red", lty=1, size=1) +
xlab("Day") +
ylab("Affect (Positive = Blue, Negative = Red)") + ylim(1,7) +
scale_x_continuous(breaks=seq(0,8,by=1)) +
facet_wrap( ~ id)
Restructuring the Data Set
Splitting Predictors
First, we split our predictor variables into time-varying and time-invariant components. We calculate the imeans as usual using the plyr functions.
#Calculating imeans
daily.imeans <- ddply(daily1, "id", summarize, imean.posaff=mean(posaff, na.rm=TRUE),
imean.negaff=mean(negaff, na.rm=TRUE))
#Calculating sample-centered versions *Note that this is done in a person-level data file.
daily.imeans$imean.posaff.c <- scale(daily.imeans$imean.posaff,center=TRUE,scale=FALSE)
daily.imeans$imean.negaff.c <- scale(daily.imeans$imean.negaff,center=TRUE,scale=FALSE)
describe(daily.imeans)## vars n mean sd median trimmed mad min max
## id 1 190 318.29 130.44 321.50 318.99 151.23 101.00 532.00
## imean.posaff 2 190 4.10 0.73 4.10 4.10 0.63 2.23 6.04
## imean.negaff 3 190 2.48 0.73 2.41 2.43 0.72 1.11 5.09
## imean.posaff.c 4 190 0.00 0.73 0.00 0.00 0.63 -1.88 1.94
## imean.negaff.c 5 190 0.00 0.73 -0.07 -0.05 0.72 -1.37 2.61
## range skew kurtosis se
## id 431.00 -0.04 -1.09 9.46
## imean.posaff 3.81 0.04 0.08 0.05
## imean.negaff 3.98 0.68 0.45 0.05
## imean.posaff.c 3.81 0.04 0.08 0.05
## imean.negaff.c 3.98 0.68 0.45 0.05#merging "trait" scores back into the *long* data file and calculate "state" scores.
daily1 <- merge(daily1,daily.imeans,by="id")
daily1$posaff.state <- daily1$posaff - daily1$imean.posaff
daily1$negaff.state <- daily1$negaff - daily1$imean.negaffLagged variables
Second, we make some lag variables. These will serve as predictors.
We use the slide() function in the DataCombine library (see http://www.r-bloggers.com/slide-one-function-for-laglead-variables-in-data-frames-including-time-series-cross-sectional-data/). Note that becasue slide() returns a data frame, there is a bit of extra mechanics.
#make a duplicate data set
daily1.new <- daily1
#make lag variables (original)
daily1.new$posaff.lag1 <- slide(daily1, Var = "posaff", GroupVar = "id", slideBy = -1)[,11]##
## Remember to order daily1 by id and the time variable before running.##
## Lagging posaff by 1 time units.daily1.new$negaff.lag1 <- slide(daily1, Var = "negaff", GroupVar = "id", slideBy = -1)[,11]##
## Remember to order daily1 by id and the time variable before running.##
## Lagging negaff by 1 time units.#make lag variables (person-centered)
daily1.new$negaff.state.lag1 <- slide(daily1, Var = "negaff.state", GroupVar = "id", slideBy = -1)[,11]##
## Remember to order daily1 by id and the time variable before running.##
## Lagging negaff.state by 1 time units.daily1.new$posaff.state.lag1 <- slide(daily1, Var = "posaff.state", GroupVar = "id", slideBy = -1)[,11]##
## Remember to order daily1 by id and the time variable before running.##
## Lagging posaff.state by 1 time units.#Checking structure
head(daily1.new,12)## id day posaff negaff imean.posaff imean.negaff imean.posaff.c
## 1 101 0 3.9 3.0 4.5625 1.50000 0.460261
## 2 101 1 3.8 2.3 4.5625 1.50000 0.460261
## 3 101 2 5.1 1.0 4.5625 1.50000 0.460261
## 4 101 3 5.6 1.3 4.5625 1.50000 0.460261
## 5 101 4 4.3 1.1 4.5625 1.50000 0.460261
## 6 101 5 3.9 1.0 4.5625 1.50000 0.460261
## 7 101 6 5.1 1.2 4.5625 1.50000 0.460261
## 8 101 7 4.8 1.1 4.5625 1.50000 0.460261
## 9 102 0 6.3 1.4 5.1750 2.21875 1.072761
## 10 102 1 7.0 1.6 5.1750 2.21875 1.072761
## 11 102 2 6.1 2.6 5.1750 2.21875 1.072761
## 12 102 3 6.2 2.8 5.1750 2.21875 1.072761
## imean.negaff.c posaff.state negaff.state posaff.lag1 negaff.lag1
## 1 -0.9783089 -0.6625 1.50000 NA NA
## 2 -0.9783089 -0.7625 0.80000 3.9 3.0
## 3 -0.9783089 0.5375 -0.50000 3.8 2.3
## 4 -0.9783089 1.0375 -0.20000 5.1 1.0
## 5 -0.9783089 -0.2625 -0.40000 5.6 1.3
## 6 -0.9783089 -0.6625 -0.50000 4.3 1.1
## 7 -0.9783089 0.5375 -0.30000 3.9 1.0
## 8 -0.9783089 0.2375 -0.40000 5.1 1.2
## 9 -0.2595589 1.1250 -0.81875 NA NA
## 10 -0.2595589 1.8250 -0.61875 6.3 1.4
## 11 -0.2595589 0.9250 0.38125 7.0 1.6
## 12 -0.2595589 1.0250 0.58125 6.1 2.6
## negaff.state.lag1 posaff.state.lag1
## 1 NA NA
## 2 1.50000 -0.6625
## 3 0.80000 -0.7625
## 4 -0.50000 0.5375
## 5 -0.20000 1.0375
## 6 -0.40000 -0.2625
## 7 -0.50000 -0.6625
## 8 -0.30000 0.5375
## 9 NA NA
## 10 -0.81875 1.1250
## 11 -0.61875 1.8250
## 12 0.38125 0.9250Let’s plot the updated data - within-person and between-variables all in one.
Can we add between-person differences as well?
#Person-centered
#setting up for changing ids
i <- 106
ggplot(data = subset(daily1.new, id <=i), aes(x=day, group=id), legend=FALSE) +
geom_rect(mapping=aes(xmin=0, xmax=8, ymin=-3, ymax=0, fill=imean.negaff.c), alpha=0.6) +
geom_rect(mapping=aes(xmin=0, xmax=8, ymin=0, ymax=+3, fill=imean.posaff.c), alpha=0.6) +
guides(fill=FALSE) + #turining off the legend
geom_hline(aes(yintercept=0)) +
geom_point(aes(x=day,y = posaff.state), color="blue", shape=17, size=2) +
geom_line(aes(x=day,y = posaff.state), color="blue", lty=1, size=1) +
geom_point(aes(x=day,y = negaff.state), color="red", shape=17, size=2) +
geom_line(aes(x=day,y = negaff.state), color="red", lty=1, size=1) +
geom_point(aes(x=day,y = posaff.state.lag1), color="blue", shape=17, size=2) +
geom_line(aes(x=day,y = posaff.state.lag1), color="blue", lty=2, size=1) +
geom_point(aes(x=day,y = negaff.state.lag1), color="red", shape=17, size=2) +
geom_line(aes(x=day,y = negaff.state.lag1), color="red", lty=2, size=1) +
xlab("Day") +
ylab("Affect (Positive = Blue, Negative = Red)") + #ylim(1,7) +
scale_x_continuous(breaks=seq(0,8,by=1)) +
facet_wrap( ~ id)
The colors and labelling are not great - but the idea is intact. Cool!
NOTICE THAT THERE IS A PROBLEM WITH ID#103. Missing Data is not being treated correctly! We need to make some more adjustments when using slide(). Be careful! (In this case the problem could be corrected by using a full dataset with 190 x 8 = 1520 rows)
Making Double-Entry Data
Ok - now we have all our variables in place and can shape the data into a double-entry set-up. We make use of the meltfunction in the reshape2 package.
#melting data
names(daily1.new)## [1] "id" "day" "posaff"
## [4] "negaff" "imean.posaff" "imean.negaff"
## [7] "imean.posaff.c" "imean.negaff.c" "posaff.state"
## [10] "negaff.state" "posaff.lag1" "negaff.lag1"
## [13] "negaff.state.lag1" "posaff.state.lag1"daily2 <- melt(data=daily1.new,
id.vars=c("id","day","imean.posaff","imean.negaff","imean.posaff.c","imean.negaff.c",
"posaff.state","negaff.state",
"negaff.lag1","posaff.lag1","negaff.state.lag1","posaff.state.lag1"),
na.rm=FALSE, value.name="affect")
#reordering for convenience
daily2 <- daily2[order(daily2$id,daily2$day,daily2$var),]
#look at updated data set
head(daily2, 12)## id day imean.posaff imean.negaff imean.posaff.c imean.negaff.c
## 1 101 0 4.5625 1.5 0.460261 -0.9783089
## 1459 101 0 4.5625 1.5 0.460261 -0.9783089
## 2 101 1 4.5625 1.5 0.460261 -0.9783089
## 1460 101 1 4.5625 1.5 0.460261 -0.9783089
## 3 101 2 4.5625 1.5 0.460261 -0.9783089
## 1461 101 2 4.5625 1.5 0.460261 -0.9783089
## 4 101 3 4.5625 1.5 0.460261 -0.9783089
## 1462 101 3 4.5625 1.5 0.460261 -0.9783089
## 5 101 4 4.5625 1.5 0.460261 -0.9783089
## 1463 101 4 4.5625 1.5 0.460261 -0.9783089
## 6 101 5 4.5625 1.5 0.460261 -0.9783089
## 1464 101 5 4.5625 1.5 0.460261 -0.9783089
## posaff.state negaff.state negaff.lag1 posaff.lag1 negaff.state.lag1
## 1 -0.6625 1.5 NA NA NA
## 1459 -0.6625 1.5 NA NA NA
## 2 -0.7625 0.8 3.0 3.9 1.5
## 1460 -0.7625 0.8 3.0 3.9 1.5
## 3 0.5375 -0.5 2.3 3.8 0.8
## 1461 0.5375 -0.5 2.3 3.8 0.8
## 4 1.0375 -0.2 1.0 5.1 -0.5
## 1462 1.0375 -0.2 1.0 5.1 -0.5
## 5 -0.2625 -0.4 1.3 5.6 -0.2
## 1463 -0.2625 -0.4 1.3 5.6 -0.2
## 6 -0.6625 -0.5 1.1 4.3 -0.4
## 1464 -0.6625 -0.5 1.1 4.3 -0.4
## posaff.state.lag1 variable affect
## 1 NA posaff 3.9
## 1459 NA negaff 3.0
## 2 -0.6625 posaff 3.8
## 1460 -0.6625 negaff 2.3
## 3 -0.7625 posaff 5.1
## 1461 -0.7625 negaff 1.0
## 4 0.5375 posaff 5.6
## 1462 0.5375 negaff 1.3
## 5 1.0375 posaff 4.3
## 1463 1.0375 negaff 1.1
## 6 -0.2625 posaff 3.9
## 1464 -0.2625 negaff 1.0#adding the two dummy indicators
daily2$pos <- ifelse(daily2$variable=="posaff", 1, 0)
daily2$neg <- ifelse(daily2$variable=="negaff", 1, 0)
#adding another indicator
daily2$var <- as.numeric(factor(daily2$variable))
#look at updated data set
head(daily2,12)## id day imean.posaff imean.negaff imean.posaff.c imean.negaff.c
## 1 101 0 4.5625 1.5 0.460261 -0.9783089
## 1459 101 0 4.5625 1.5 0.460261 -0.9783089
## 2 101 1 4.5625 1.5 0.460261 -0.9783089
## 1460 101 1 4.5625 1.5 0.460261 -0.9783089
## 3 101 2 4.5625 1.5 0.460261 -0.9783089
## 1461 101 2 4.5625 1.5 0.460261 -0.9783089
## 4 101 3 4.5625 1.5 0.460261 -0.9783089
## 1462 101 3 4.5625 1.5 0.460261 -0.9783089
## 5 101 4 4.5625 1.5 0.460261 -0.9783089
## 1463 101 4 4.5625 1.5 0.460261 -0.9783089
## 6 101 5 4.5625 1.5 0.460261 -0.9783089
## 1464 101 5 4.5625 1.5 0.460261 -0.9783089
## posaff.state negaff.state negaff.lag1 posaff.lag1 negaff.state.lag1
## 1 -0.6625 1.5 NA NA NA
## 1459 -0.6625 1.5 NA NA NA
## 2 -0.7625 0.8 3.0 3.9 1.5
## 1460 -0.7625 0.8 3.0 3.9 1.5
## 3 0.5375 -0.5 2.3 3.8 0.8
## 1461 0.5375 -0.5 2.3 3.8 0.8
## 4 1.0375 -0.2 1.0 5.1 -0.5
## 1462 1.0375 -0.2 1.0 5.1 -0.5
## 5 -0.2625 -0.4 1.3 5.6 -0.2
## 1463 -0.2625 -0.4 1.3 5.6 -0.2
## 6 -0.6625 -0.5 1.1 4.3 -0.4
## 1464 -0.6625 -0.5 1.1 4.3 -0.4
## posaff.state.lag1 variable affect pos neg var
## 1 NA posaff 3.9 1 0 1
## 1459 NA negaff 3.0 0 1 2
## 2 -0.6625 posaff 3.8 1 0 1
## 1460 -0.6625 negaff 2.3 0 1 2
## 3 -0.7625 posaff 5.1 1 0 1
## 1461 -0.7625 negaff 1.0 0 1 2
## 4 0.5375 posaff 5.6 1 0 1
## 1462 0.5375 negaff 1.3 0 1 2
## 5 1.0375 posaff 4.3 1 0 1
## 1463 1.0375 negaff 1.1 0 1 2
## 6 -0.2625 posaff 3.9 1 0 1
## 1464 -0.2625 negaff 1.0 0 1 2Yes! Now the data are in the double-entry format we need for modeling.
3. The Multilevel Model.
We are now ready to start running our models!
We start with the simple “unconditional means” model.
\[affect_{it} = \beta_{A0i}dummyA_{it} + \beta_{B0i}dummyB_{it} + e_{it}\] \[\beta_{A0i} = \gamma_{A00} + u_{A0i}\] \[\beta_{B0i} = \gamma_{B00} + u_{B0i}\]
with \[e_{it} \sim \mathcal{N}(0,\mathbf{R})\], and \[\mathbf{U_{i}} \sim \mathcal{N}(0,\mathbf{G})\]. where \[\mathbf{R} = \mathbf{I} \left[\begin{array} {r} \sigma^2_{e} \end{array}\right]\] and \[\mathbf{G} = \left[\begin{array} {rr} \sigma^{2}_{uA0} & \sigma_{uA0uB0} \\ \sigma_{uA0uB0} & \sigma^{2}_{uB0} \end{array}\right]\]
model1.fit <- lme(fixed = affect ~ -1 + pos + neg,
random = ~ -1 + pos + neg|id,
data=daily2,
na.action = na.exclude)
summary(model1.fit)## Linear mixed-effects model fit by REML
## Data: daily2
## AIC BIC logLik
## 7856.977 7892.77 -3922.488
##
## Random effects:
## Formula: ~-1 + pos + neg | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## pos 0.6638692 pos
## neg 0.6482248 -0.493
## Residual 0.8478374
##
## Fixed effects: affect ~ -1 + pos + neg
## Value Std.Error DF t-value p-value
## pos 4.105112 0.05333648 2691 76.96630 0
## neg 2.463913 0.05230692 2691 47.10492 0
## Correlation:
## pos
## neg -0.403
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.75033477 -0.59659553 -0.03672567 0.55217304 3.83689841
##
## Number of Observations: 2882
## Number of Groups: 190In this model we only have one time-varying residual - which is implicitly assumed to apply to both positive affect and negative affect. We need to change that. We want the residuals to be structured as
\[\mathbf{R} = \left[\begin{array} {ll} \sigma^2_{eA} & \sigma_{eAeB} \\ \sigma_{eAeB} & \sigma^2_{eB} \end{array}\right]\]
We add two statements to do this.
#adding the proper residual structure
model2.fit <- lme(fixed = affect ~ -1 + pos + neg,
random = ~ -1 + pos + neg | id,
weights=varIdent(form=~1 | var),
corr=corCompSymm(form=~1 | id/day),
data=daily2,
na.action = na.exclude)
summary(model2.fit)## Linear mixed-effects model fit by REML
## Data: daily2
## AIC BIC logLik
## 7428.507 7476.231 -3706.254
##
## Random effects:
## Formula: ~-1 + pos + neg | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## pos 0.6561755 pos
## neg 0.6534759 -0.374
## Residual 0.8806943
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id/day
## Parameter estimate(s):
## Rho
## -0.5363001
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | var
## Parameter estimates:
## 1 2
## 1.0000000 0.9243499
## Fixed effects: affect ~ -1 + pos + neg
## Value Std.Error DF t-value p-value
## pos 4.109555 0.05322494 2691 77.21107 0
## neg 2.463665 0.05228379 2691 47.12102 0
## Correlation:
## pos
## neg -0.404
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.84623955 -0.59782579 -0.04443476 0.54610635 3.94725062
##
## Number of Observations: 2882
## Number of Groups: 190Ok, that accommodates the possibility of different variances, and of possible 3rd variables.
We add the auto-reggression time-varying predictors posaff.state.lag1 and negaff.state.lag1.
Our equation becomes \[affect_{it} = \\ \beta_{A0i}dummyA_{it} + \beta_{A1i}dummyA_{it}posaff_{it-1} + e_{Ait} + \\ \beta_{B0i}dummyB_{it} + \beta_{B1i}dummyB_{it}negaff_{it-1} + e_{Bit}\] \[\beta_{A0i} = \gamma_{A00} + u_{A0i}\] \[\beta_{B0i} = \gamma_{B00} + u_{B0i}\]
with \[\mathbf{R} = \left[\begin{array} {ll} \sigma^2_{eA} & \sigma_{eAeB} \\ \sigma_{eAeB} & \sigma^2_{eB} \end{array}\right]\] and \[\mathbf{G} = \left[\begin{array} {rr} \sigma^{2}_{uA0} & \sigma_{uA0uB0} \\ \sigma_{uA0uB0} & \sigma^{2}_{uB0} \end{array}\right]\]
Note the use of the : for the “interaction” terms.
#adding autoregression predictor
model3.fit <- lme(fixed = affect ~ -1 +
pos + pos:posaff.state.lag1 +
neg + neg:negaff.state.lag1,
random = ~ -1 + pos + neg | id,
weights=varIdent(form = ~1 | var),
corr=corCompSymm(form = ~1 | id/day),
data=daily2,
na.action = na.exclude)
summary(model3.fit)## Linear mixed-effects model fit by REML
## Data: daily2
## AIC BIC logLik
## 6319.05 6377.235 -3149.525
##
## Random effects:
## Formula: ~-1 + pos + neg | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## pos 0.6808219 pos
## neg 0.6594216 -0.313
## Residual 0.8592714
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id/day
## Parameter estimate(s):
## Rho
## -0.5782807
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | var
## Parameter estimates:
## 1 2
## 1.0000000 0.9249552
## Fixed effects: affect ~ -1 + pos + pos:posaff.state.lag1 + neg + neg:negaff.state.lag1
## Value Std.Error DF t-value p-value
## pos 4.081515 0.05567630 2301 73.30794 0
## neg 2.403529 0.05346519 2301 44.95502 0
## pos:posaff.state.lag1 0.149636 0.02554707 2301 5.85728 0
## neg:negaff.state.lag1 0.138608 0.02572431 2301 5.38821 0
## Correlation:
## pos neg ps:..1
## neg -0.362
## pos:posaff.state.lag1 -0.008 -0.004
## neg:negaff.state.lag1 -0.003 -0.012 0.307
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.80177208 -0.56924485 -0.04599583 0.54317950 4.33287762
##
## Number of Observations: 2490
## Number of Groups: 186We add the cross-reggression time-varying predictors negaff.state.lag1 and posaff.state.lag1.
Our equation becomes \[affect_{it} = \\ \beta_{A0i}dummyA_{it} + \beta_{A1i}dummyA_{it}posaff_{it-1} + \beta_{A2i}dummyA_{it}negaff_{it-1}+ e_{Ait} + \\ \beta_{B0i}dummyB_{it} + \beta_{B1i}dummyB_{it}negaff_{it-1} + \beta_{B2i}dummyB_{it}posaff_{it-1} + e_{Bit}\] \[\beta_{A0i} = \gamma_{A00} + u_{A0i}\] \[\beta_{B0i} = \gamma_{B00} + u_{B0i}\]
#adding cross-regression predictors
model4.fit <- lme(fixed = affect ~ -1 +
pos + pos:posaff.state.lag1 + pos:negaff.state.lag1 +
neg + neg:negaff.state.lag1 + neg:posaff.state.lag1,
random = ~ -1 + pos + neg | id,
weights=varIdent(form = ~1 | var),
corr=corCompSymm(form = ~1 | id/day),
data=daily2,
na.action = na.exclude)
summary(model4.fit)## Linear mixed-effects model fit by REML
## Data: daily2
## AIC BIC logLik
## 6322.196 6392.007 -3149.098
##
## Random effects:
## Formula: ~-1 + pos + neg | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## pos 0.6813932 pos
## neg 0.6591309 -0.311
## Residual 0.8563035
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id/day
## Parameter estimate(s):
## Rho
## -0.5781934
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | var
## Parameter estimates:
## 1 2
## 1.000000 0.927447
## Fixed effects: affect ~ -1 + pos + pos:posaff.state.lag1 + pos:negaff.state.lag1 + neg + neg:negaff.state.lag1 + neg:posaff.state.lag1
## Value Std.Error DF t-value p-value
## pos 4.078292 0.05569206 2299 73.22933 0.0000
## neg 2.404015 0.05344993 2299 44.97696 0.0000
## pos:posaff.state.lag1 0.177320 0.03504729 2299 5.05944 0.0000
## pos:negaff.state.lag1 0.107901 0.03814840 2299 2.82846 0.0047
## negaff.state.lag1:neg 0.102557 0.03538434 2299 2.89839 0.0038
## posaff.state.lag1:neg 0.022531 0.03250784 2299 0.69310 0.4883
## Correlation:
## pos neg ps:p..1 ps:n..1 ng..1:
## neg -0.361
## pos:posaff.state.lag1 -0.020 0.011
## pos:negaff.state.lag1 -0.023 0.013 0.532
## negaff.state.lag1:neg 0.013 -0.022 -0.307 -0.578
## posaff.state.lag1:neg 0.012 -0.019 -0.578 -0.307 0.532
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.74583840 -0.56948604 -0.05883869 0.54773833 4.14314073
##
## Number of Observations: 2490
## Number of Groups: 186We add the person-level time-invariant predictors imean.posaff.c and imean.negaff.c.
\[affect_{it} = \beta_{A0i}dummyA_{it} + \beta_{A1i}dummyA_{it}posaff_{it-1} + \beta_{A2i}dummyA_{it}negaff_{it-1}+ e_{Ait} + \\ \beta_{B0i}dummyB_{it} + \beta_{B1i}dummyB_{it}negaff_{it-1} + \beta_{B2i}dummyB_{it}posaff_{it-1} + e_{Bit}\] \[\beta_{A0i} = \gamma_{A00} + \gamma_{A01}imeanposaff_{i} + \gamma_{A02}imeannegaff_{i} + u_{A0i}\] \[\beta_{B0i} = \gamma_{B00} + \gamma_{B01}imeannegaff_{i} + \gamma_{B02}imeanposaff_{i} + u_{B0i}\]
#adding main effects for time-invariant predictors
model5.fit <- lme(fixed = affect ~ -1 +
pos + pos:imean.posaff.c + pos:imean.negaff.c +
pos:posaff.state.lag1 + pos:negaff.state.lag1 +
neg + neg:imean.negaff.c + neg:imean.posaff.c +
neg:negaff.state.lag1 + neg:posaff.state.lag1,
random = ~ -1 + pos + neg | id,
weights=varIdent(form = ~1 | var),
corr=corCompSymm(form = ~1 | id/day),
data=daily2,
na.action = na.exclude)
summary(model5.fit)## Linear mixed-effects model fit by REML
## Data: daily2
## AIC BIC logLik
## 5355.654 5448.71 -2661.827
##
## Random effects:
## Formula: ~-1 + pos + neg | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## pos 1.633020e-05 pos
## neg 2.400985e-05 0
## Residual 8.017601e-01
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id/day
## Parameter estimate(s):
## Rho
## -0.5754094
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | var
## Parameter estimates:
## 1 2
## 1.0000000 0.9279277
## Fixed effects: affect ~ -1 + pos + pos:imean.posaff.c + pos:imean.negaff.c + pos:posaff.state.lag1 + pos:negaff.state.lag1 + neg + neg:imean.negaff.c + neg:imean.posaff.c + neg:negaff.state.lag1 + neg:posaff.state.lag1
## Value Std.Error DF t-value p-value
## pos 4.066356 0.02278661 2295 178.45372 0.0000
## neg 2.436755 0.02114433 2295 115.24390 0.0000
## pos:imean.posaff.c 1.032514 0.03398790 2295 30.37887 0.0000
## pos:imean.negaff.c 0.033348 0.03487640 2295 0.95617 0.3391
## pos:posaff.state.lag1 0.175271 0.03257307 2295 5.38085 0.0000
## pos:negaff.state.lag1 0.109973 0.03544511 2295 3.10263 0.0019
## imean.negaff.c:neg 1.012772 0.03236277 2295 31.29435 0.0000
## imean.posaff.c:neg 0.015787 0.03153832 2295 0.50055 0.6167
## negaff.state.lag1:neg 0.095372 0.03289050 2295 2.89968 0.0038
## posaff.state.lag1:neg 0.020261 0.03022546 2295 0.67032 0.5027
## Correlation:
## pos neg ps:mn.p. ps:mn.n. ps:p..1 ps:n..1
## neg -0.575
## pos:imean.posaff.c -0.005 0.003
## pos:imean.negaff.c 0.045 -0.026 0.398
## pos:posaff.state.lag1 -0.044 0.025 0.022 0.009
## pos:negaff.state.lag1 -0.049 0.028 0.011 0.006 0.533
## imean.negaff.c:neg -0.026 0.045 -0.229 -0.575 -0.005 -0.004
## imean.posaff.c:neg 0.003 -0.005 -0.575 -0.229 -0.013 -0.006
## negaff.state.lag1:neg 0.028 -0.049 -0.006 -0.004 -0.307 -0.575
## posaff.state.lag1:neg 0.025 -0.044 -0.013 -0.005 -0.575 -0.307
## imn.n.: imn.p.: ng..1:
## neg
## pos:imean.posaff.c
## pos:imean.negaff.c
## pos:posaff.state.lag1
## pos:negaff.state.lag1
## imean.negaff.c:neg
## imean.posaff.c:neg 0.398
## negaff.state.lag1:neg 0.006 0.011
## posaff.state.lag1:neg 0.009 0.022 0.533
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.97737894 -0.58357800 -0.05172056 0.54842112 4.11842551
##
## Number of Observations: 2490
## Number of Groups: 186Wait a second! - this model has imean.posaff predicting posaff with \(\gamma_{A01}\) = 1.01, and imean.negaff predicting negaff with \(\gamma_{B01}\) = 1.0. That is a “duh!”. The model doesn’t make sense! Do we actually need the between-person variables at all?
Note the difference here between the standard cross-lag regression model (which is based on between-person associations) and this model (which is based on within-person associations).
We add additional random effects.
We can push one more step though - and add some more random effects (between-person differences). How far we can push I’m not sure, but lets try!
In reality, we build the model up slowly - increasing complexity one random effect at a time. However, here I just jump to as far as could go.
#adding additional random effects
model6.fit <- lme(fixed = affect ~ -1 +
pos + pos:posaff.state.lag1 + pos:negaff.state.lag1 +
neg + neg:negaff.state.lag1 + neg:posaff.state.lag1,
random = ~ -1 +
pos + pos:posaff.state.lag1 + pos:negaff.state.lag1 +
neg + neg:negaff.state.lag1 + neg:posaff.state.lag1 | id,
weights=varIdent(form = ~1 | var),
corr=corCompSymm(form = ~1 | id/day),
data=daily2,
na.action = na.exclude,
control = lmeControl(maxIter=500, opt = "optim"))
summary(model6.fit)## Linear mixed-effects model fit by REML
## Data: daily2
## AIC BIC logLik
## 6332.808 6507.336 -3136.404
##
## Random effects:
## Formula: ~-1 + pos + pos:posaff.state.lag1 + pos:negaff.state.lag1 + neg + neg:negaff.state.lag1 + neg:posaff.state.lag1 | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## pos 0.6824439 pos neg ps:p..1 ps:n..1 ng..1:
## neg 0.6625180 -0.306
## pos:posaff.state.lag1 0.1236745 0.071 -0.092
## pos:negaff.state.lag1 0.2419485 -0.010 -0.178 0.840
## negaff.state.lag1:neg 0.1943738 -0.102 -0.035 -0.771 -0.914
## posaff.state.lag1:neg 0.2016603 -0.487 0.028 -0.477 -0.470 0.726
## Residual 0.8369835
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id/day
## Parameter estimate(s):
## Rho
## -0.5738982
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | var
## Parameter estimates:
## 1 2
## 1.0000000 0.9235012
## Fixed effects: affect ~ -1 + pos + pos:posaff.state.lag1 + pos:negaff.state.lag1 + neg + neg:negaff.state.lag1 + neg:posaff.state.lag1
## Value Std.Error DF t-value p-value
## pos 4.078768 0.05554350 2299 73.43376 0.0000
## neg 2.402182 0.05344377 2299 44.94784 0.0000
## pos:posaff.state.lag1 0.172519 0.03665855 2299 4.70609 0.0000
## pos:negaff.state.lag1 0.116137 0.04375927 2299 2.65400 0.0080
## negaff.state.lag1:neg 0.081494 0.03928389 2299 2.07448 0.0381
## posaff.state.lag1:neg 0.002962 0.03642181 2299 0.08134 0.9352
## Correlation:
## pos neg ps:p..1 ps:n..1 ng..1:
## neg -0.354
## pos:posaff.state.lag1 -0.003 -0.011
## pos:negaff.state.lag1 -0.025 -0.054 0.552
## negaff.state.lag1:neg -0.021 -0.032 -0.359 -0.650
## posaff.state.lag1:neg -0.169 -0.007 -0.558 -0.336 0.559
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.46829673 -0.56399712 -0.05553651 0.54784304 3.90461396
##
## Number of Observations: 2490
## Number of Groups: 186#intervals(model6.fit)
#Error in intervals.lme(model6.fit) :
# cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covarianceWell, we get estimates - but there is a problem with the model. We cannot get confidence intervals. So, we probably would try to find justification for fitting a simpler set of random effects.
4. Conclusion.
This tutorial is meant to illustrate how one can use the multivariate multilevel model to examine intraindividual coupling. We have not given an in-depth treatment, only provided a “seat-of-the-pants” example. When using the model for a real empirical inquiry - more care is needed. Be careful!
Cheers!