Preliminaries

Loading libraries used in this script

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

The Research Questions

For this example, we use the 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:

1. How do the emotions carry-over from day to day? This is a construct called emotional inertia in the literature.

2. How does positive emotion influence negative emotions? (protective factor) and, vice versa,

3. How does negative emotion influence positive emotions? (risk factor)

Basically, this is a cross-lag panel model (or APIM-type model) - except that we are working with intensive longitudinal data, which provides opportunity to prioritize examination of within-person (rather than between-person) associations.

The Multivariate Modeling Enterprise

The basic multilevel model is designed as a model with a univariate outcome. So, we “trick” the model/program into thinking that two (or more) variables are one variable. So clever!

The AMIB Data

We use The AMIB Data, a multiple time-scale data set that has been shared for teaching purposes. The data include person-level dispositions, daily diary assessments, and ecological momentary assessments that were obtained after everyday social interactions (event-contingent sampling).

The data are shared with intent to facilitate learning about analysis of intensive longitudinal data (from experience sampling, daily diary, or EMA designs). The data are posted only for teaching purposes and should not be used for research.

Use requires citation and acknowledgement of both this website and the following paper:

We are 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.

#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/AMIB/AMIB_daily.csv"
#read in the .csv file using the url() function
daily <- read.csv(file=url(filepath), header=TRUE)

#subsetting data to variables of interest
daily1 <- daily %>%
  select(id, day, posaff, negaff)

#examine first few rows of the data set
names(daily1)

## [1] "id"     "day"    "posaff" "negaff"

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
psych::describe(daily1$negaff)

##    vars    n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 1441 2.45 1.04    2.2    2.34 1.04   1 6.9   5.9 0.96     0.77 0.03

psych::describe(daily1$posaff)

##    vars    n mean  sd median trimmed  mad min max range  skew kurtosis   se
## X1    1 1441 4.12 1.1    4.2    4.15 1.19   1   7     6 -0.25    -0.33 0.03

#histogram
daily1 %>%
  ggplot(aes(x=negaff)) +
  geom_histogram(fill="red", color="black") +
  labs(x = "Negative Affect")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

daily1 %>%
  ggplot(aes(x=posaff)) +
  geom_histogram(fill="royalblue2", 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

#plot
daily1 %>%
  filter(id <= i) %>%
  ggplot(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, linewidth=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, linewidth=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

#Calculating iMeans
daily.imeans <- daily1 %>%
  dplyr::group_by(id) %>%
  dplyr::summarise(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 <- daily.imeans %>%
  mutate(imean.posaff.c = scale(imean.posaff, center=TRUE, scale=FALSE),
         imean.negaff.c = scale(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.22   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 <- daily1 %>%
  left_join(daily.imeans, by="id") %>%
  mutate(posaff.state = posaff - imean.posaff,
         negaff.state = negaff - imean.negaff)

#make a duplicate data set
daily1.new <- daily1

#make lag variables (original)
daily1.new <- daily1.new %>%
  dplyr::group_by(id) %>%
  mutate(posaff.lag1 = lag(posaff, 1),
         negaff.lag1 = lag(negaff, 1),
         #make lag variables (person-centered)
         negaff.state.lag1 = lag(negaff.state, 1),
         posaff.state.lag1 = lag(posaff.state, 1))

#Checking structure
head(daily1.new, 12)

## # A tibble: 12 × 14
## # Groups:   id [2]
##       id   day posaff negaff imean.posaff imean.negaff imean.posaff.c[,1]
##    <int> <int>  <dbl>  <dbl>        <dbl>        <dbl>              <dbl>
##  1   101     0    3.9    3           4.56         1.5               0.460
##  2   101     1    3.8    2.3         4.56         1.5               0.460
##  3   101     2    5.1    1           4.56         1.5               0.460
##  4   101     3    5.6    1.3         4.56         1.5               0.460
##  5   101     4    4.3    1.1         4.56         1.5               0.460
##  6   101     5    3.9    1           4.56         1.5               0.460
##  7   101     6    5.1    1.2         4.56         1.5               0.460
##  8   101     7    4.8    1.1         4.56         1.5               0.460
##  9   102     0    6.3    1.4         5.18         2.22              1.07 
## 10   102     1    7      1.6         5.18         2.22              1.07 
## 11   102     2    6.1    2.6         5.18         2.22              1.07 
## 12   102     3    6.2    2.8         5.18         2.22              1.07 
## # ℹ 7 more variables: imean.negaff.c <dbl[,1]>, posaff.state <dbl>,
## #   negaff.state <dbl>, posaff.lag1 <dbl>, negaff.lag1 <dbl>,
## #   negaff.state.lag1 <dbl>, posaff.state.lag1 <dbl>

#Person-centered 
#setting up for changing ids
i <- 106

#plot
daily1.new %>%
  filter(id <= i) %>%
  ggplot(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, linewidth=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, linewidth=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, linewidth=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, linewidth=1) +
  xlab("Day") + 
  ylab("Affect (Positive = Blue, Negative = Red)") + #ylim(1,7) +
  scale_x_continuous(breaks=seq(0,8,by=1)) + 
  facet_wrap( ~ id)

#reshaping data
dput(names(daily1.new))

## c("id", "day", "posaff", "negaff", "imean.posaff", "imean.negaff", 
## "imean.posaff.c", "imean.negaff.c", "posaff.state", "negaff.state", 
## "posaff.lag1", "negaff.lag1", "negaff.state.lag1", "posaff.state.lag1"
## )

daily2 <- daily1.new %>%
  pivot_longer(cols = c("posaff", "negaff"),
               names_to = "variable",
               values_to = "affect",
               values_drop_na = FALSE)

#look at updated data set
head(daily2, 12)

## # A tibble: 12 × 14
## # Groups:   id [1]
##       id   day imean.posaff imean.negaff imean.posaff.c[,1] imean.negaff.c[,1]
##    <int> <int>        <dbl>        <dbl>              <dbl>              <dbl>
##  1   101     0         4.56          1.5              0.460             -0.978
##  2   101     0         4.56          1.5              0.460             -0.978
##  3   101     1         4.56          1.5              0.460             -0.978
##  4   101     1         4.56          1.5              0.460             -0.978
##  5   101     2         4.56          1.5              0.460             -0.978
##  6   101     2         4.56          1.5              0.460             -0.978
##  7   101     3         4.56          1.5              0.460             -0.978
##  8   101     3         4.56          1.5              0.460             -0.978
##  9   101     4         4.56          1.5              0.460             -0.978
## 10   101     4         4.56          1.5              0.460             -0.978
## 11   101     5         4.56          1.5              0.460             -0.978
## 12   101     5         4.56          1.5              0.460             -0.978
## # ℹ 8 more variables: posaff.state <dbl>, negaff.state <dbl>,
## #   posaff.lag1 <dbl>, negaff.lag1 <dbl>, negaff.state.lag1 <dbl>,
## #   posaff.state.lag1 <dbl>, variable <chr>, affect <dbl>

#adding two dummy indicators for the two variables of interest
daily2 <- daily2 %>%
  mutate(pos = ifelse(variable=="posaff", 1, 0),
         neg = ifelse(variable=="negaff", 1, 0),
         #adding another variable indicator 
         var = as.numeric(factor(variable, 
                                 levels = c("posaff", "negaff"),
                                 labels = c(1, 2))))

#look at updated data set
head(daily2, 12)

## # A tibble: 12 × 17
## # Groups:   id [1]
##       id   day imean.posaff imean.negaff imean.posaff.c[,1] imean.negaff.c[,1]
##    <int> <int>        <dbl>        <dbl>              <dbl>              <dbl>
##  1   101     0         4.56          1.5              0.460             -0.978
##  2   101     0         4.56          1.5              0.460             -0.978
##  3   101     1         4.56          1.5              0.460             -0.978
##  4   101     1         4.56          1.5              0.460             -0.978
##  5   101     2         4.56          1.5              0.460             -0.978
##  6   101     2         4.56          1.5              0.460             -0.978
##  7   101     3         4.56          1.5              0.460             -0.978
##  8   101     3         4.56          1.5              0.460             -0.978
##  9   101     4         4.56          1.5              0.460             -0.978
## 10   101     4         4.56          1.5              0.460             -0.978
## 11   101     5         4.56          1.5              0.460             -0.978
## 12   101     5         4.56          1.5              0.460             -0.978
## # ℹ 11 more variables: posaff.state <dbl>, negaff.state <dbl>,
## #   posaff.lag1 <dbl>, negaff.lag1 <dbl>, negaff.state.lag1 <dbl>,
## #   posaff.state.lag1 <dbl>, variable <chr>, affect <dbl>, pos <dbl>,
## #   neg <dbl>, var <dbl>

Unconditional means model

We start with the simple “unconditional means” model.

\[affect_{it} = \beta_{A0i}dummyA_{it} + \beta_{B0i}dummyB_{it} + e_{it}\] (Notice that there is no intercept in the model). \[\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]\]

Coding the model. We put a “-1” where we would usually put a “1” for the intercept in order to remove the default.

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: 190

Fixing Residual Structure

In the previous 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), #to get different for each var
                  corr=corSymm(form=~1 | id/day), #to allow them to correlate
                  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.6561721 pos   
## neg      0.6534680 -0.374
## Residual 0.8806941       
## 
## Correlation Structure: General
##  Formula: ~1 | id/day 
##  Parameter estimate(s):
##  Correlation: 
##   1     
## 2 -0.536
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | var 
##  Parameter estimates:
##         1         2 
## 1.0000000 0.9243557 
## Fixed effects:  affect ~ -1 + pos + neg 
##        Value  Std.Error   DF  t-value p-value
## pos 4.109555 0.05322472 2691 77.21140       0
## neg 2.463665 0.05228332 2691 47.12144       0
##  Correlation: 
##     pos   
## neg -0.404
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.84620289 -0.59782678 -0.04443837  0.54610645  3.94722848 
## 
## Number of Observations: 2882
## Number of Groups: 190

Ok, that accommodates the possibility of different variances, and for unknown 3rd variable “common cause” (via the correlated residuals).

Auto-Regressions

We add the auto-regression 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}\] and \[\beta_{A1i} = \gamma_{A10}\] \[\beta_{B1i} = \gamma_{B10}\] (Notice that these autoregression do not yet have random effects on them; no \(u_{A1i}\) or \(u_{B1i}\))

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. We do this becasue we need control of what main effects and interactions are being added. We are overriding the defaults with out “trick”.

#adding auto-regression 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=corSymm(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.6808490 pos   
## neg      0.6594594 -0.313
## Residual 0.8592680       
## 
## Correlation Structure: General
##  Formula: ~1 | id/day 
##  Parameter estimate(s):
##  Correlation: 
##   1     
## 2 -0.578
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | var 
##  Parameter estimates:
##         1         2 
## 1.0000000 0.9249473 
## Fixed effects:  affect ~ -1 + pos + pos:posaff.state.lag1 + neg + neg:negaff.state.lag1 
##                          Value  Std.Error   DF  t-value p-value
## pos                   4.081514 0.05567805 2301 73.30562       0
## neg                   2.403529 0.05346759 2301 44.95301       0
## pos:posaff.state.lag1 0.149636 0.02554709 2301  5.85727       0
## neg:negaff.state.lag1 0.138607 0.02572411 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.80189130 -0.56926009 -0.04598771  0.54317530  4.33285365 
## 
## Number of Observations: 2490
## Number of Groups: 186

Cross Regressions

We add the cross-regression 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}\] and \[\beta_{A1i} = \gamma_{A10}\] \[\beta_{B1i} = \gamma_{B10}\] and \[\beta_{A2i} = \gamma_{A20}\] \[\beta_{B2i} = \gamma_{B20}\]

#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.6813933 pos   
## neg      0.6591403 -0.311
## Residual 0.8563014       
## 
## Correlation Structure: Compound symmetry
##  Formula: ~1 | id/day 
##  Parameter estimate(s):
##        Rho 
## -0.5781921 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | var 
##  Parameter estimates:
##         1         2 
## 1.0000000 0.9274466 
## 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.05569204 2299 73.22936  0.0000
## neg                   2.404015 0.05345052 2299 44.97646  0.0000
## pos:posaff.state.lag1 0.177320 0.03504720 2299  5.05945  0.0000
## pos:negaff.state.lag1 0.107901 0.03814830 2299  2.82847  0.0047
## negaff.state.lag1:neg 0.102558 0.03538424 2299  2.89840  0.0038
## posaff.state.lag1:neg 0.022531 0.03250774 2299  0.69311  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.74585606 -0.56948320 -0.05883887  0.54773629  4.14314587 
## 
## Number of Observations: 2490
## Number of Groups: 186

Person-Level Predictors

We 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}\] and \[\beta_{A1i} = \gamma_{A10}\] \[\beta_{B1i} = \gamma_{B10}\] and \[\beta_{A2i} = \gamma_{A20}\] \[\beta_{B2i} = \gamma_{B20}\]

#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=corSymm(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.136189e-05 pos 
## neg      2.478979e-05 0   
## Residual 8.017601e-01     
## 
## Correlation Structure: General
##  Formula: ~1 | id/day 
##  Parameter estimate(s):
##  Correlation: 
##   1     
## 2 -0.575
## 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 imn.n.:
## 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   0.398 
## negaff.state.lag1:neg  0.028 -0.049 -0.006   -0.004   -0.307  -0.575   0.006 
## posaff.state.lag1:neg  0.025 -0.044 -0.013   -0.005   -0.575  -0.307   0.009 
##                       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                  
## negaff.state.lag1:neg  0.011        
## posaff.state.lag1:neg  0.022   0.533
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.97737894 -0.58357800 -0.05172056  0.54842112  4.11842550 
## 
## Number of Observations: 2490
## Number of Groups: 186

Wait 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).

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! Can we include all the random effects for

and \[\beta_{A1i} = \gamma_{A10} + u_{A1i}\] \[\beta_{B1i} = \gamma_{B10} + u_{B1i}\] and \[\beta_{A2i} = \gamma_{A20} + u_{A2i}\] \[\beta_{B2i} = \gamma_{B20} + u_{B2i}\]

In reality, we build the model up slowly - increasing complexity one random effect at a time. However, here I just jump as far as we could go (and then regret it :-).

#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")) #added in attempt to obtain convergence
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.6625179 -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.8369836                                     
## 
## 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.05344376 2299 44.94785  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.46829629 -0.56399714 -0.05553655  0.54784299  3.90461447 
## 
## 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-covariance

Well, we get estimates - but there is a problem with the model. We cannot get confidence intervals. So, we could try to find justification for fitting a simpler set of random effects; or we can try to do all of this in a Bayesian framework. Look for that in a future tutorial.