Overview

This tutorial extends examination ofintraindividual covariation with the multilevel modeling framework. Specifically, this tutorial demonstrates how the multilevel model is used to articulate and test hypotheses about dynamic characteristics from experience sampling data.

Outline

In this session we cover …

A. Preliminary preparation and description of two variables and the dynamic characteristic of interest
B. Setting up the Multilevel Model C. Using the Mutlilevel Model to answering a series of research questions

Q1: Do individuals differ in their “baseline” levels of outcome variable?
Q2: Does the average outcome vary by a specific between-person variable(s)?
Q3: Does the outcome vary by a specific within-person variable?
Q4: Does the within-person association differ across persons?
Q5: Are the between-person differences in the within-person association moderated by a between-person predictor?

The main illustration uses a binary predictor, with the final model also given for a continuous predictor.
The overall set-up of the models follows Bolger & Laurenceau (2013) Chapters 4 and 5.

Prelim - Loading libraries used in this script.

library(psych)
library(ggplot2)
library(plyr)
library(nlme)
library(effects)

Prelim - Reading in Repeated Measures Data

Note that we are working from a long file. For your own data, there may be some steps to get to this point.

#Setting the working directory
#setwd("~/Desktop/Fall_2017")  #Person 1 Computer
#setwd("~/Desktop/Fall_2017")  #Person 2 Computer

#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

Everything we do today uses a long file.

A. Preliminary preparation and description of two variables and the dynamic characteristic of interest

We consider two variables, negaff (negaff) and a new variable “stress” (derived from pss).

names(daily)
##  [1] "id"      "day"     "date"    "slphrs"  "weath"   "lteq"    "pss"    
##  [8] "se"      "swls"    "evalday" "posaff"  "negaff"  "temp"    "hum"    
## [15] "wind"    "bar"     "prec"
#sample descriptives
describe(daily$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.03
describe(daily$pss)
##    vars    n mean   sd median trimmed  mad min max range  skew kurtosis
## X1    1 1445 2.61 0.68   2.75    2.64 0.74   0   4     4 -0.35     0.13
##      se
## X1 0.02
#histograms
ggplot(data=daily, aes(x=negaff)) +
  geom_histogram(fill="white", color="black") + 
  labs(x = "Negative Affect")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 17 rows containing non-finite values (stat_bin).

ggplot(data=daily, aes(x=pss)) +
  geom_histogram(fill="white", color="black") + 
  labs(x = "Perceived Stress Scale (low = more stressed") + 
  geom_vline(xintercept=2.75, size=2, color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 13 rows containing non-finite values (stat_bin).

For today’s illustration, we want to have a binary predictor variable stress.
So, we create one from the pss variable using an “arbitrary” cut-point (the median). Days with pss < 2.75 are coded as stress days, and days with pss >= 2.75 are no-stress days (the pss is scored such that low sores are indicative of more stress).
Note that as a general principle, “degrading” the level of measurement (from interval to categorical) is a loss of information. We do it for didactic purposes.

#making a new binary variable
daily$stress <- ifelse(daily$pss < 2.75, 1, 0)

Check that ifelse() dealt with missing data the way we want.

table(daily$stress,daily$pss, useNA = "always")
##       
##          0 0.25 0.75   1 1.25 1.5 1.666666667 1.75   2 2.25 2.333333333
##   0      0    0    0   0    0   0           0    0   0    0           0
##   1      2    4    9  17   23  57           4   90 124  168           2
##   <NA>   0    0    0   0    0   0           0    0   0    0           0
##       
##        2.5 2.666666667 2.75   3 3.25 3.333333333 3.5 3.75   4 <NA>
##   0      0           0  206 194  138           2 121   41  37    0
##   1    205           1    0   0    0           0   0    0   0    0
##   <NA>   0           0    0   0    0           0   0    0   0   13
describe(daily$stress)
##    vars    n mean  sd median trimmed mad min max range skew kurtosis   se
## X1    1 1445 0.49 0.5      0    0.49   0   0   1     1 0.05       -2 0.01
#historgram of the factor variable
ggplot(data=daily, aes(x=factor(stress))) +
  geom_bar(stat="count",fill="white", color="black") 

That looks like we want.

Setting up the substantive inquiry. Defining the dynamic characteristic of interest

OK - Lets briefly consider the substantive framework of our inquiry. We define a person-level dynamic characteristic, stress reactivty as the extent to which a person’s level of daily negative affect differs between stress days (=1) and no-stress days (=0), i.e., stress reactivity is the intraindividual covariation of daily negative affect and daily stress.

Formally,
\[NegAffect_{it} = \beta_{0i} + \beta_{1i}stressday_{it} + e_{it}\]

where \(\beta_{1i}\) is individual i’s level of stress reactivity.
Note that we are giving an explicit name to the coefficient. We are “Naming the Betas” (see Ram & Grimm, 2015).

Let’s look at the stress reactivity phenomena. We make a bivariate longitudinal time-series plots …
For Person A

#ggplot version .. see also http://ggplot.yhathq.com/docs/index.html
ggplot(data = subset(daily, id ==101), aes(x=day), legend=FALSE) +
  geom_rect(mapping=aes(xmin=day-.5, xmax=day+.5, ymin=1, ymax=7, fill=stress), alpha=0.8) +
  geom_point(aes(x=day,y = negaff), shape=17, size=3,colour="red") +
  geom_line(aes(x=day,y = negaff), lty=1, size=1,colour="red") +
  xlab("Day") + 
  ylab("Negative Affect") + ylim(1,7) +
  scale_x_continuous(breaks=seq(0,7,by=1)) 

For Person B

ggplot(data = subset(daily, id ==123), aes(x=day), legend=FALSE) +
  geom_rect(mapping=aes(xmin=day-.5, xmax=day+.5, ymin=1, ymax=7, fill=stress), alpha=0.8) +
  geom_point(aes(x=day,y = negaff), shape=17, size=3,colour="red") +
  geom_line(aes(x=day,y = negaff), lty=1, size=1,colour="red") +
  xlab("Day") + 
  ylab("Negative Affect") + ylim(1,7) +
  scale_x_continuous(breaks=seq(0,7,by=1)) 

From the time-series plots we get a sense that Negative Affect is higher on stress days compared to no-stress days, but to varying degree.

Lets look more explicitly at the within-person relation of negaff and stress.

#Scatter Plots by id
i <- 101
ggplot(data=daily[which(daily$id ==i),], aes(x=stress, y=negaff)) + 
  geom_point(color="red") +
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="red") +
  xlab("Stress Day") + scale_x_continuous(breaks=seq(0,1,by=1)) +
  ylab("Negative Affect") + ylim(1,7) + 
  ggtitle(paste("ID =",i)) 

i <- 123
ggplot(data=daily[which(daily$id ==i),], aes(x=stress, y=negaff)) + 
  geom_point(color="red") +
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="red") +
  xlab("Stress Day") + scale_x_continuous(breaks=seq(0,1,by=1)) +
  ylab("Negative Affect") + ylim(1,7) + 
  ggtitle(paste("ID =",i)) 

Quite informative. Let’s look at all the indvidual-level regressions at once.

ggplot(data=daily, aes(x=stress, y=negaff, group=id, color="gray"), legend=FALSE) + 
  geom_point(color="gray") + guides(color=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=FALSE, lty=1, size=.5, color="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="red") +
  xlab("Stress Day") + scale_x_continuous(breaks=seq(0,1,by=1)) +
  ylab("Negative Affect")  + ylim(1,7) +
  ggtitle("Stress Reactivity (Mix of Within & Between)")
## Warning: Removed 20 rows containing non-finite values (stat_smooth).

## Warning: Removed 20 rows containing non-finite values (stat_smooth).
## Warning: Removed 20 rows containing missing values (geom_point).

Also very informative, BUT, we have to be careful because the plot (and the regressions within it) have mixed together between-person and within-person associations.

B. Setting up the Multilevel Model

Now that we have defined the dynamic characteristic of interest - stress reactivity - we can operationalize it using the multilevel model.

The nlme package contains tools for fitting mixed effects models (MLM, HLM). We will also make use of a supplementary package, effects that provides tools for calculating and plotting model-based predictions.

The lme() function fits the MLMs The fixed argument takes the fixed model The random argument takes the random model The data argument specifies the data sources The na.action argument specifies what to do with missing data

We set up the basic model. The requires choosing which variable will the “outcome” and which variable will be the “predictor”. We choose negaff as the outcome variable, and assume it is normal, not because it is, but because that is how it has been done in almost all the prior literature on stress reactivity.

1. First we fit the Unconditional Means model to see how much between-person and how much within-person variance there is in our outcome variable.

\[y_{it} = \beta_{0i} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i}\]

where \[e_{it} \tilde{} N(0,\sigma^{2}_{e})\], and \[u_{i} \tilde{} N(0,\sigma^{2}_{u0})\].

For clarity, lets write out the full variance covariance matrix of the within-person residuals (spanning across the t = 8 days of observation), \[\left[\begin{array} {rrrrrrrr} \sigma^{2}_{e} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \sigma^{2}_{e} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \sigma^{2}_{e} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \sigma^{2}_{e} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \sigma^{2}_{e} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \sigma^{2}_{e} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \sigma^{2}_{e} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma^{2}_{e} \end{array}\right]\], This is the homoscedasticity of errors assumption.

And, the full variance covariance matrix of the between-person residuals, \[\left[\begin{array} {r} \sigma^{2}_{u0} \end{array}\right]\]

Ok - lets fit the model using the lme() function.

#Use the constant '1' to designate the intercepts
um.fit <- lme(fixed= negaff ~ 1, 
              random= ~ 1|id, 
              data=daily,
              na.action=na.exclude)
summary(um.fit)
## Linear mixed-effects model fit by REML
##  Data: daily 
##       AIC      BIC   logLik
##   3839.46 3855.277 -1916.73
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)  Residual
## StdDev:   0.6534748 0.8140798
## 
## Fixed effects: negaff ~ 1 
##                Value Std.Error   DF  t-value p-value
## (Intercept) 2.463675 0.0522867 1251 47.11858       0
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8739086 -0.6122847 -0.1608485  0.4657840  3.9394042 
## 
## Number of Observations: 1441
## Number of Groups: 190
#Extract the random effects
VarCorr(um.fit)
## id = pdLogChol(1) 
##             Variance  StdDev   
## (Intercept) 0.4270294 0.6534748
## Residual    0.6627260 0.8140798

Compute the ICC … \[ICC_{between} = \frac{\sigma^{2}_{u0}}{\sigma^{2}_{u0} + \sigma^{2}_{e}}\]

RandomEffects <- as.numeric(VarCorr(um.fit)[,1])
ICC_between <- RandomEffects[1]/(RandomEffects[1]+RandomEffects[2]) 
ICC_between
## [1] 0.391858

So, between-person variance = 39.19%, and within-person variance = 100 - 39.19 = 60.81%

2. We need to “split”" our predictor variable into “trait” and “state” components.

We calculate the imeans using the plyr package.

daily.imeans <- ddply(daily, "id", summarize, imean.stress=mean(stress, na.rm=TRUE),
                                              imean.negaff=mean(negaff, na.rm=TRUE))

Note that we really only need to do all of this for the predictor variable. However, we arealso doing it for the outcome in order to have those scores available for some later checking and plotting.

In the stress literature, the proportion of days that are stress days is an indicator of Exposure, a person-level variable. Exposure is a construct, just like Stress Reactivity. Together, they are components of the Daily Stress Process model.

For a specific analysis it may be useful to rename the “trait” variable for conceptual clarity, with something like …
daily.imeans$exposure <- daily.imeans$imean.stress.
However, for the purposes of our didactic presentation here we stick with the imean.xxx labeling.

We calculate a sample-centered versions for later use. *Note that this is done in a person-level data file.

daily.imeans$imean.stress.c <- scale(daily.imeans$imean.stress,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.stress      2 190   0.49   0.31   0.50    0.49   0.37   0.00   1.00
## imean.negaff      3 190   2.48   0.73   2.41    2.43   0.72   1.11   5.09
## imean.stress.c    4 190   0.00   0.31   0.01    0.00   0.37  -0.49   0.51
## 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.stress     1.00  0.03    -1.19 0.02
## imean.negaff     3.98  0.68     0.45 0.05
## imean.stress.c   1.00  0.03    -1.19 0.02
## imean.negaff.c   3.98  0.68     0.45 0.05

We then merge these “trait” scores back into the long data file and calculate the “state” scores.

daily <- merge(daily,daily.imeans,by="id")
daily$negaff.state <- daily$negaff - daily$imean.negaff

Note that we have not calculated deviation scores for the binary stress variable. Rather than defining the “equilibrium” state as the person-specific mean, we define the “equilibrium” as the no-stress, stress = 0, state.
This all depends on the conceptualization, so be aware of the assummption - and its later implications for how we interpret the intercept terms.

Let’s plot the state and trait scores to see if it looks right.

#Trait Scores
ggplot(data = daily[which(daily$id <=105),], aes(x = day, y = imean.negaff, group = id, color=factor(id))) +
  geom_line(data=daily[which(daily$id <= 105 & daily$imean.negaff !="NA"),]) +
  xlab("Day") + 
  ylab("Trait Negative Affect") + ylim(1,7) +
  scale_x_continuous(breaks=seq(0,7,by=1)) 

#State Scores (note change of the ylim)
ggplot(data = daily[which(daily$id <=105),], aes(x = day, y = negaff.state, group = id, color=factor(id))) +
  geom_line(data=daily[which(daily$id <= 105 & daily$negaff.state !="NA"),]) +
  xlab("Day") + 
  ylab("State Negative Affect") + ylim(-3,3) +
  scale_x_continuous(breaks=seq(0,7,by=1)) 

OK - all looking good. Basically, the plots provide a visualiation of the variance decomposition that we did for calculating the ICC.

C. Using the Mutlilevel Model to answering a series of research questions.

We use a series of questions to work up to the full model that is used in most investigations of intraindividual covariation.
Note that we try to match the Bolger and Laurenceau (2013) model configuration as close as possible. The ordering of the questions follows Schwartz and Stone (2007).

Q1. Do people differ in their “baseline” levels of outcome?

e.g., Are there between-person differences in the intercept?

There are a couple of ways to operationalize this question.

We can do it as above, using the Unconditional Means Model

\[y_{it} = \beta_{0i} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i}\]

where \[e_{it} \tilde{} N(0,\sigma^{2}_{e})\], and \[u_{i} \tilde{} N(0,\sigma^{2}_{u0})\].

#Unconditional Means Model
um.fit <- lme(fixed= negaff ~ 1, 
              random= ~ 1|id, 
              data=daily,
              na.action=na.exclude)
summary(um.fit)
## Linear mixed-effects model fit by REML
##  Data: daily 
##       AIC      BIC   logLik
##   3839.46 3855.277 -1916.73
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)  Residual
## StdDev:   0.6534748 0.8140798
## 
## Fixed effects: negaff ~ 1 
##                Value Std.Error   DF  t-value p-value
## (Intercept) 2.463675 0.0522867 1251 47.11858       0
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8739086 -0.6122847 -0.1608485  0.4657840  3.9394042 
## 
## Number of Observations: 1441
## Number of Groups: 190

Note that this is how we specified the model last time. It assumes that there are no time-dependencies in the data - the i.i.d. assumption - and that the “baseline” level is the average. Note that the study design influences the choice here.
Trends that manifest over the short-term may not be meaningful over the long-term.

We might also operationalize the question using a Linear Growth Model with Lag 1 Auto-Regression.
Bolger & Laurenceau (2013) use an expanded base model that takes care of some time dependencies. In particular they examine whether the outcome has any linear time-related trends, and any (intrinsic) auto-regression time-dependencies (carryover from day-to-day). Specifically,

\[y_{it} = \beta_{0i} + \beta_{1i}time_{it} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i}\] \[\beta_{1i} = \gamma_{10} + u_{1i}\] and where the residual variance-covarianc structures are now …

for the within-person residuals (spanning across the t = 8 days of observation), \[\left[\begin{array} {rrrrrrrr} \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 & \phi^6 & \phi^7 \\ \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 & \phi^6 \\ \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 \\ \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 \\ \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 \\ \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 \\ \phi^6 & \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi \\ \phi^7 & \phi^6 & \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} \end{array}\right]\], where \(\phi\) is the lag-1 autoregression.

And, the full variance covariance matrix of the between-person residuals, \[\left[\begin{array} {rr} \sigma^{2}_{u0} & \sigma_{u0u1} \\ \sigma_{u1u0} & \sigma^{2}_{u1} \end{array}\right]\]

Note that we now have two variances and a covariance.

What is the time variable in our sample data? day

describe(daily$day) 
##    vars    n mean  sd median trimmed  mad min max range skew kurtosis   se
## X1    1 1458 3.48 2.3      3    3.47 2.97   0   7     7 0.01    -1.24 0.06

We add day as a predictor into our model, and also the autoregression on the within-person residuals. Note that here we are not directly interested in the effects of time or the autoregression - we are simply using them as “covariates” to get rid of some systematic time-dependencies.

#Linear Growth Model with Lag 1 Auto-Regression 
lg.fit <- lme(fixed= negaff ~ 1 + day, 
              random= ~ 1 + day|id, 
              correlation = corAR1(),
              data=daily,
              na.action=na.exclude)
summary(lg.fit)
## Linear mixed-effects model fit by REML
##  Data: daily 
##        AIC      BIC    logLik
##   3752.929 3789.831 -1869.465
## 
## Random effects:
##  Formula: ~1 + day | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 0.6220962 (Intr)
## day         0.0376319 -0.186
## Residual    0.8285062       
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | id 
##  Parameter estimate(s):
##       Phi 
## 0.2444986 
## Fixed effects: negaff ~ 1 + day 
##                  Value  Std.Error   DF  t-value p-value
## (Intercept)  2.6875518 0.06475834 1250 41.50125       0
## day         -0.0639393 0.01120603 1250 -5.70579       0
##  Correlation: 
##     (Intr)
## day -0.594
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.7497480 -0.6188429 -0.1531224  0.4611812  3.8039166 
## 
## Number of Observations: 1441
## Number of Groups: 190
intervals(lg.fit)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                   lower        est.       upper
## (Intercept)  2.56050479  2.68755183  2.81459886
## day         -0.08592399 -0.06393929 -0.04195459
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                             lower       est.     upper
## sd((Intercept))       0.502201738  0.6220962 0.7706140
## sd(day)               0.005232782  0.0376319 0.2706323
## cor((Intercept),day) -0.729281053 -0.1859620 0.5011850
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.1628294 0.2444986 0.3228402
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.7813905 0.8285062 0.8784629

Interpretation of output …

As we did earlier for the unconditional means model, we can obtain and plot the individual predicted scores from the linear growth model.

#getting predicted scores
daily$pred.negaff.lg.fit <- predict(lg.fit)
#plotting predicted scores
ggplot(data = daily, aes(x = day, y = pred.negaff.lg.fit, group = id)) +
  geom_line(color="gray40") +
  geom_point(data=daily[which(daily$day == 0),]) + #to obtain predicted intercept scores
  xlab("Day") + 
  ylab("Predicted Negative Affect") + #ylim(1,7) +
  scale_x_continuous(breaks=seq(0,7,by=1)) 
## Warning: Removed 5 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_point).

This is informative. But, we must remember that our interest is really in stress, not in time.

Our question was … Q1. Do people differ in their “baseline” levels of outcome? e.g., Are there between-person differences in the intercept?

This model has explicitly removed time-dependencies in the data - thus the “baseline” level is the intercept - which we emphasized in the plot with a black dot.

The statistical information for answering Q1 is in the output. Does the 95% confidence interval for the standard deviation, \sigma_{u0} include zero? It does not. We have evidence to reject the null hypothesis (i.e., we reject the hypothesis that all intercepts are the same).

What does the intercept represent? The expected level of negaff when everything else = 0. So here, “baseline” is a day on which everything in the model is zero. Stress is not in the model - so “baseline” is a time = 0 day with average stress (whatever that is).
Conclusion - much care has to be taken in the interpretation of model parameters - and the assumptions embedded in the model.

Ok - lets move on …

Q2. Does the average outcome vary by a specific between-person variable(s)?

e.g., Adding in a predictor of between-person differences in intercept

For our example this is the question … Is there a significant relation between individuals’ expected level of Negative Affect and their level of stress Exposure? This is a between-person association question.

We have already created the centered person-level exposure variable, imean.stress.c and put it in our long data set.

Our model becomes … \[y_{it} = \beta_{0i} + \beta_{1i}time_{it} + e_{it}\] \[\beta_{0i} = \gamma_{00} + \gamma_{01}Z_{i} + u_{0i}\] \[\beta_{1i} = \gamma_{10} + u_{1i}\] where

for the within-person residuals (spanning across the t = 8 days of observation), \[\left[\begin{array} {rrrrrrrr} \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 & \phi^6 & \phi^7 \\ \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 & \phi^6 \\ \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 \\ \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 \\ \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 \\ \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 \\ \phi^6 & \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi \\ \phi^7 & \phi^6 & \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} \end{array}\right]\],

And, the between-person residuals structure is \[\left[\begin{array} {rr} \sigma^{2}_{u0} & \sigma_{u0u1} \\ \sigma_{u1u0} & \sigma^{2}_{u1} \end{array}\right]\]

Let’s fit the model to the data.

lg.fit.q2 <- lme(fixed= negaff ~ 1 + day + imean.stress.c, 
                 random= ~ 1 + day|id, 
                 correlation = corAR1(),
                 data=daily,
                 na.action=na.exclude)
summary(lg.fit.q2)
## Linear mixed-effects model fit by REML
##  Data: daily 
##        AIC      BIC    logLik
##   3681.348 3723.516 -1832.674
## 
## Random effects:
##  Formula: ~1 + day | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept) 0.41444533 (Intr)
## day         0.04348276 0.109 
## Residual    0.82437720       
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | id 
##  Parameter estimate(s):
##       Phi 
## 0.2329255 
## Fixed effects: negaff ~ 1 + day + imean.stress.c 
##                     Value  Std.Error   DF  t-value p-value
## (Intercept)     2.6859101 0.05483946 1250 48.97769       0
## day            -0.0642326 0.01121055 1250 -5.72966       0
## imean.stress.c  1.3317519 0.13627198  188  9.77275       0
##  Correlation: 
##                (Intr) day   
## day            -0.634       
## imean.stress.c  0.002  0.000
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3508764 -0.6119010 -0.1276897  0.4665249  3.8250373 
## 
## Number of Observations: 1441
## Number of Groups: 190
intervals(lg.fit.q2)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                      lower       est.       upper
## (Intercept)     2.57832259  2.6859101  2.79349768
## day            -0.08622617 -0.0642326 -0.04223904
## imean.stress.c  1.06293328  1.3317519  1.60057056
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                            lower       est.     upper
## sd((Intercept))       0.32319451 0.41444533 0.5314599
## sd(day)               0.01714445 0.04348276 0.1102835
## cor((Intercept),day) -0.39023865 0.10931595 0.5591465
## 
##  Correlation structure:
##         lower      est.    upper
## Phi 0.1595173 0.2329255 0.303772
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.7834373 0.8243772 0.8674565

Interpretation of output

Is there a significant relation between individuals’ expected level of Negative Affect and their level of stress Exposure?
\(\gamma_{01} = 1.33\), with 95%CI = [1.06, 1.60]

Q3. Does the outcome vary by a specific within-person variable?

Is change in the outcome related to changes in predictor? e.g., Adding in a “slope”" but no between-person differences in slope

For our example this is the question … Is there, on average, a significant within-person association between stress day (a time-varying predictor) and Negative Affect? This is a within-person association question.

We have already created our “state” variable stress in the long data set. (Recall the binary coding of 0,1)

Our model becomes … \[y_{it} = \beta_{0i} + \beta_{1i}time_{it} + \beta_{2i}stress_{it} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i}\] \[\beta_{1i} = \gamma_{10} + u_{1i}\] \[\beta_{2i} = \gamma_{20}\] where the within-person residual structure is as above, and, the between-person residuals structure is still \[\left[\begin{array} {rr} \sigma^{2}_{u0} & \sigma_{u0u1} \\ \sigma_{u1u0} & \sigma^{2}_{u1} \end{array}\right]\].

Let’s fit the model to the data.

lg.fit.q3 <- lme(fixed= negaff ~ 1 + day + stress, 
                 random= ~ 1 + day|id, 
                 correlation = corAR1(),
                 data=daily,
                 na.action=na.exclude)
summary(lg.fit.q3)
## Linear mixed-effects model fit by REML
##  Data: daily 
##        AIC      BIC    logLik
##   3493.065 3535.217 -1738.533
## 
## Random effects:
##  Formula: ~1 + day | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 0.5512037 (Intr)
## day         0.0619794 -0.326
## Residual    0.7365948       
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | id 
##  Parameter estimate(s):
##       Phi 
## 0.1619872 
## Fixed effects: negaff ~ 1 + day + stress 
##                  Value  Std.Error   DF  t-value p-value
## (Intercept)  2.2942572 0.06058907 1246 37.86586       0
## day         -0.0628972 0.01035319 1246 -6.07515       0
## stress       0.7904272 0.04618193 1246 17.11551       0
##  Correlation: 
##        (Intr) day   
## day    -0.573       
## stress -0.377  0.006
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.5753647 -0.6184302 -0.1157149  0.5001331  3.6965697 
## 
## Number of Observations: 1438
## Number of Groups: 190
intervals(lg.fit.q3)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                   lower        est.      upper
## (Intercept)  2.17538934  2.29425720  2.4131251
## day         -0.08320883 -0.06289722 -0.0425856
## stress       0.69982427  0.79042720  0.8810301
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                            lower       est.     upper
## sd((Intercept))       0.44285860  0.5512037 0.6860555
## sd(day)               0.03294649  0.0619794 0.1165965
## cor((Intercept),day) -0.66039004 -0.3257310 0.1168782
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.0795809 0.1619872 0.2421945
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.6976842 0.7365948 0.7776754

Interpretation of output

Is there, on average, a significant within-person association between stress day (a time-varying predictor) and Negative Affect? \(\gamma_{20} = 0.79\), with 95%CI = [0.70, 0.88]

Q4. Does the within-person association differ across persons?

e.g., Allowing for between-person differences in slope

For our example this is the question … Are there between-person differences in stress reactivity? This is a question about between-person variance.

Our model becomes … \[y_{it} = \beta_{0i} + \beta_{1i}time_{it} + \beta_{2i}stress_{it} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i}\] \[\beta_{1i} = \gamma_{10} + u_{1i}\] \[\beta_{2i} = \gamma_{20} + \boldsymbol{u_{2i}}\] where the within-person residual structure is as above, but the between-person residuals structure is now expanded \[\left[\begin{array} {rrr} \sigma^{2}_{u0} & \sigma_{u0u1} & \sigma_{u0u2} \\ \sigma_{u1u0} & \sigma^{2}_{u1} & \sigma_{u1u2}\\ \sigma_{u2u0} & \sigma_{u2u1} & \sigma^{2}_{u2} \end{array}\right]\].

Lets fit the model to the data.

lg.fit.q4 <- lme(fixed= negaff ~ 1 + day + stress, 
                 random= ~ 1 + day + stress|id, 
                 correlation = corAR1(),
                 data=daily,
                 na.action=na.exclude)
summary(lg.fit.q4)
## Linear mixed-effects model fit by REML
##  Data: daily 
##        AIC      BIC    logLik
##   3426.434 3484.392 -1702.217
## 
## Random effects:
##  Formula: ~1 + day + stress | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr         
## (Intercept) 0.51107550 (Intr) day   
## day         0.06837691 -0.643       
## stress      0.58491639 -0.079  0.265
## Residual    0.68709871              
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | id 
##  Parameter estimate(s):
##      Phi 
## 0.142763 
## Fixed effects: negaff ~ 1 + day + stress 
##                  Value  Std.Error   DF  t-value p-value
## (Intercept)  2.2849287 0.05704976 1246 40.05150       0
## day         -0.0670127 0.01015746 1246 -6.59739       0
## stress       0.8147939 0.06285115 1246 12.96387       0
##  Correlation: 
##        (Intr) day   
## day    -0.678       
## stress -0.299  0.104
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.4051201 -0.6029285 -0.1307576  0.4849364  4.3291812 
## 
## Number of Observations: 1438
## Number of Groups: 190
intervals(lg.fit.q4)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                   lower        est.       upper
## (Intercept)  2.17300446  2.28492866  2.39685286
## day         -0.08694029 -0.06701269 -0.04708508
## stress       0.69148818  0.81479394  0.93809971
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                               lower        est.      upper
## sd((Intercept))          0.39526228  0.51107550  0.6608224
## sd(day)                  0.04183239  0.06837691  0.1117651
## sd(stress)               0.47299167  0.58491639  0.7233260
## cor((Intercept),day)    -0.83206691 -0.64305558 -0.3202410
## cor((Intercept),stress) -0.33616060 -0.07898019  0.1891617
## cor(day,stress)         -0.13936921  0.26499152  0.5936093
## 
##  Correlation structure:
##          lower     est.     upper
## Phi 0.05640137 0.142763 0.2270041
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.6493411 0.6870987 0.7270518

Interpretation of output

Q4. Does the within-person association differ across persons? \(\sigma_{u2} = 0.58\), with 95% CI = [0.47, 0.72]

Lets use the parameters to get our “prediction equation”.

\[negaff_{it} = \beta_{0i} + \beta_{1i}time_{it} + \beta_{2i}stress_{it} + e_{it}\] \[\beta_{0i} = 2.2849287 + u_{0i}\] \[\beta_{1i} = -0.0670127 + u_{1i}\] \[\beta_{2i} = 0.8147939 + u_{2i}\]

We can use those equations to calculate the expected scores for the protypical person.

Lets make a plot of this result.
We make use of the effect() function in the effects package.

We are not interested in the day effects (in this case day is a covariate). We are interested in the stress effect.

#term="..." designates the specific effect we are interested in
#mod= designates the model object
#xlevels = is the list of values that we want to evaluate the equation at.
ef4 <- effect(term="stress", mod=lg.fit.q4, xlevels=list(stress=c(0,1)))
summary(ef4)
## 
##  stress effect
## stress
##        0        1 
## 2.051270 2.866064 
## 
##  Lower 95 Percent Confidence Limits
## stress
##        0        1 
## 1.968796 2.741435 
## 
##  Upper 95 Percent Confidence Limits
## stress
##        0        1 
## 2.133743 2.990692

Everything we need is in there!

We can convert to a data frame and plot it!

preddata4 <- as.data.frame(ef4)
#plotting the effect evaluation (with standard error ribbon)
ggplot(data=preddata4, aes(x=stress, y=fit), legend=FALSE) + 
  geom_point() +
  geom_line() +
  geom_ribbon(aes(ymin=lower, ymax=upper), alpha=.3) +
  geom_errorbar(aes(x=stress, ymin=lower, ymax=upper), width=.05) +
  xlab("Stress Day") + scale_x_continuous(breaks=seq(0,1,by=1)) +
  ylab("Predicted Negative Affect") + ylim(1,7) +
  ggtitle("Prediction from Model 4")

Q5. Are the between-person differences in the within-person association moderated by a between-person predictor?

e.g., Adding in a between-person variable to explain differences in slope

Our model becomes … \[y_{it} = \beta_{0i} + \beta_{1i}time_{it} + \beta_{2i}stress_{it} + e_{it}\] \[\beta_{0i} = \gamma_{00} + \gamma_{01}Z_{i} + u_{0i}\] \[\beta_{1i} = \gamma_{10} + u_{1i}\] \[\beta_{2i} = \gamma_{20} + \gamma_{21}Z_{i} + u_{2i}\] where the within-person residual structure is as above, and the between-person residuals structure is the expanded \[\left[\begin{array} {rrr} \sigma^{2}_{u0} & \sigma_{u0u1} & \sigma_{u0u2} \\ \sigma_{u1u0} & \sigma^{2}_{u1} & \sigma_{u1u2}\\ \sigma_{u2u0} & \sigma_{u2u1} & \sigma^{2}_{u2} \end{array}\right]\].

Note that \(Z_{i}\) has not been included as a predictor of \(\beta_{1i}\).

For our example this is the question … Is the within-person association between Stress Day (a time-varying predictor) and Negative Affect (i.e., stress reactivity) moderated by Exposure? AND Is there a between-person association of Exposure with Negative Affect?
(We do not have an interest in whether Exposure moderates the within-person across day linear trends.)

We combine the above models together.

lg.fit.q5 <- lme(fixed= negaff ~ 1 + imean.stress.c + day + stress + imean.stress.c*stress, 
                 random= ~ 1 + day + stress|id, 
                 correlation = corAR1(),
                 data=daily,
                 na.action=na.exclude)
summary(lg.fit.q5)
## Linear mixed-effects model fit by REML
##  Data: daily 
##        AIC      BIC    logLik
##   3412.315 3480.793 -1693.158
## 
## Random effects:
##  Formula: ~1 + day + stress | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr         
## (Intercept) 0.44659666 (Intr) day   
## day         0.06843386 -0.619       
## stress      0.57886272 -0.061  0.289
## Residual    0.68976191              
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | id 
##  Parameter estimate(s):
##       Phi 
## 0.1512948 
## Fixed effects: negaff ~ 1 + imean.stress.c + day + stress + imean.stress.c *      stress 
##                            Value  Std.Error   DF  t-value p-value
## (Intercept)            2.3525507 0.05744176 1245 40.95541   0.000
## imean.stress.c         0.6776827 0.15054913  188  4.50141   0.000
## day                   -0.0671249 0.01020935 1245 -6.57485   0.000
## stress                 0.7275144 0.06590980 1245 11.03803   0.000
## imean.stress.c:stress -0.1135894 0.24799955 1245 -0.45802   0.647
##  Correlation: 
##                       (Intr) imn.s. day    stress
## imean.stress.c         0.322                     
## day                   -0.642  0.003              
## stress                -0.316 -0.229  0.107       
## imean.stress.c:stress -0.228 -0.433 -0.011 -0.112
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3939391 -0.6078796 -0.1371780  0.4884681  4.3747875 
## 
## Number of Observations: 1438
## Number of Groups: 190
intervals(lg.fit.q5)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                             lower        est.       upper
## (Intercept)            2.23985735  2.35255069  2.46524402
## imean.stress.c         0.38070009  0.67768274  0.97466539
## day                   -0.08715431 -0.06712489 -0.04709547
## stress                 0.59820789  0.72751443  0.85682096
## imean.stress.c:stress -0.60013258 -0.11358939  0.37295379
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                               lower        est.      upper
## sd((Intercept))          0.33110931  0.44659666  0.6023648
## sd(day)                  0.04169526  0.06843386  0.1123196
## sd(stress)               0.47155243  0.57886272  0.7105934
## cor((Intercept),day)    -0.82837665 -0.61876096 -0.2571451
## cor((Intercept),stress) -0.24421228 -0.06123264  0.1259575
## cor(day,stress)         -0.07548979  0.28911843  0.5855333
## 
##  Correlation structure:
##          lower      est.    upper
## Phi 0.06550236 0.1512948 0.234867
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.6515856 0.6897619 0.7301749

Interpretation of output.

The moderation is not significant. However, if it were, we would want to plot the interaction.
We can again use the effect function.
We are interested in the stress*day effects (which includes the main effects as well).

#xlevels = is the list of values that we want to evaluate at.
#obtaining the standard deviation
describe(daily.imeans$imean.stress.c)
##    vars   n mean   sd median trimmed  mad   min  max range skew kurtosis
## X1    1 190    0 0.31   0.01       0 0.37 -0.49 0.51     1 0.03    -1.19
##      se
## X1 0.02
#calculating the effects
ef5 <- effect(term="imean.stress.c*stress", mod=lg.fit.q5, 
              xlevels=list(imean.stress.c=c(-0.31, +0.31), stress=c(0,1)))
summary(ef5)
## 
##  imean.stress.c*stress effect
##               stress
## imean.stress.c        0        1
##          -0.31 1.908419 2.671146
##          0.31  2.328582 3.020884
## 
##  Lower 95 Percent Confidence Limits
##               stress
## imean.stress.c        0        1
##          -0.31 1.812675 2.454044
##          0.31  2.178475 2.862714
## 
##  Upper 95 Percent Confidence Limits
##               stress
## imean.stress.c        0        1
##          -0.31 2.004163 2.888248
##          0.31  2.478689 3.179054

Everything we need is in there!

We can convert to a data frame and plot it!

#convert to dataframe
preddata5 <- as.data.frame(ef5)
#plotting the effect evaluation (with standard error ribbon)
ggplot(data=preddata5, aes(x=stress, y=fit, group=imean.stress.c), legend=FALSE) + 
  geom_point() +
  geom_line() +
  geom_ribbon(aes(ymin=lower, ymax=upper), alpha=.3) +
  geom_errorbar(aes(x=stress, ymin=lower, ymax=upper), width=.05) +
  xlab("Stress Day") + scale_x_continuous(breaks=seq(0,1,by=1)) +
  ylab("Predicted Negative Affect") + ylim(1,7) +
  ggtitle("Prediction from Model 5")

Yeah - that is cool! We will need to annotate the plot and/or add a legend, but all the information is in there. We can now expand the model in all kinds of interesting ways - adding in additional person-level predictors and/or additional time-varying predictors.

Continuous predictor.

All of the above was done with a categorical time-varying predictor, the stress variable we created.

Basically, everything is the same when we use a continuous time-varying predictor.

We do it in abbreviated from here.

But, we did not do the detrending the way that Bolger and Laurenceau do. To do that explicitly, Lets take the final model from last time and add the time variables.

Our model follows the general form \[y_{it} = \beta_{0i} + \beta_{1i}time_{it} + \beta_{2i}x_{it} + e_{it}\] \[\beta_{0i} = \gamma_{00} + \gamma_{01}Z_{i} + u_{0i}\] \[\beta_{1i} = \gamma_{10} + \gamma_{11}Z_{i} + u_{1i}\] (Note that this time we have also included \(Z_{i}\) as a moderater of the time trend. Adding another hypothesis in for fun.) \[\beta_{2i} = \gamma_{20} + \gamma_{21}Z_{i} + u_{2i}\] with the within-person residual structure with an auto-regressive structure (spanning across the t = 8 days of observation), \[\left[\begin{array} {rrrrrrrr} \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 & \phi^6 & \phi^7 \\ \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 & \phi^6 \\ \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 & \phi^5 \\ \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 & \phi^4 \\ \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 & \phi^3 \\ \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi & \phi^2 \\ \phi^6 & \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} & \phi \\ \phi^7 & \phi^6 & \phi^5 & \phi^4 & \phi^3 & \phi^2 & \phi & \sigma^{2}_{e} \end{array}\right]\], and the between-person residuals structure is “fully unstructured” \[\left[\begin{array} {rrr} \sigma^{2}_{u0} & \sigma_{u0u1} & \sigma_{u0u2} \\ \sigma_{u1u0} & \sigma^{2}_{u1} & \sigma_{u1u2}\\ \sigma_{u2u0} & \sigma_{u2u1} & \sigma^{2}_{u2} \end{array}\right]\].

We redo a bit of data prep on the pss variable, flipping the scale and creating the state scores.

#Needed to re-read the file
#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

#reverse coding the pss variable into a new stresslevel variable
daily$stresslevel <- 4 - daily$pss
#calculating intraindividual means
daily.imeans <- ddply(daily, "id", summarize, imean.stresslevel=mean(stresslevel, na.rm=TRUE),
                                              imean.negaff=mean(negaff, na.rm=TRUE))
#centering the person-level variables
daily.imeans$imean.stresslevel.c <- scale(daily.imeans$imean.stresslevel,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
## id                     1 190 318.29 130.44 321.50  318.99 151.23 101.00
## imean.stresslevel      2 190   1.40   0.48   1.41    1.39   0.51   0.19
## imean.negaff           3 190   2.48   0.73   2.41    2.43   0.72   1.11
## imean.stresslevel.c    4 190   0.00   0.48   0.01    0.00   0.51  -1.21
## imean.negaff.c         5 190   0.00   0.73  -0.07   -0.05   0.72  -1.37
##                        max  range  skew kurtosis   se
## id                  532.00 431.00 -0.04    -1.09 9.46
## imean.stresslevel     2.56   2.38 -0.04    -0.23 0.03
## imean.negaff          5.09   3.98  0.68     0.45 0.05
## imean.stresslevel.c   1.17   2.38 -0.04    -0.23 0.03
## imean.negaff.c        2.61   3.98  0.68     0.45 0.05
#merging intraindividual means into long data
daily <- merge(daily,daily.imeans,by="id")

#calculating state variables
daily$stresslevel.state <- daily$stresslevel - daily$imean.stresslevel
daily$negaff.state <- daily$negaff - daily$imean.negaff

#plotting trait and state predictor scores
#Trait Scores
ggplot(data = daily[which(daily$id <=105),], aes(x = day, y = imean.stresslevel, group = id, color=factor(id))) +
  geom_line(data=daily[which(daily$id <= 105 & daily$imean.stresslevel !="NA"),]) +
  xlab("Day") + 
  ylab("Trait Stress Level") + ylim(0,4) +
  scale_x_continuous(breaks=seq(0,7,by=1)) 

#State Scores (note change of the ylim)
ggplot(data = daily[which(daily$id <=105),], aes(x = day, y = stresslevel.state, group = id, color=factor(id))) +
  geom_line(data=daily[which(daily$id <= 105 & daily$stresslevel.state !="NA"),]) +
  xlab("Day") + 
  ylab("State Negative Affect") + ylim(-3,3) +
  scale_x_continuous(breaks=seq(0,7,by=1)) 

Ok - now we are ready to run the model.

#Full intraindividual covariation model with detrending
lg.fit.q6 <- lme(fixed= negaff ~ 1 + day + stresslevel.state + imean.stresslevel.c + imean.stresslevel.c*day + imean.stresslevel.c*stresslevel.state, 
                 random= ~ 1 + day + stresslevel.state|id, 
                 correlation = corAR1(),
                 data=daily,
                 na.action=na.exclude)
summary(lg.fit.q6)
## Linear mixed-effects model fit by REML
##  Data: daily 
##        AIC      BIC    logLik
##   3163.956 3237.691 -1567.978
## 
## Random effects:
##  Formula: ~1 + day + stresslevel.state | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##                   StdDev    Corr         
## (Intercept)       0.4631407 (Intr) day   
## day               0.0569535 -0.265       
## stresslevel.state 0.3555564  0.331  0.418
## Residual          0.6274212              
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | id 
##  Parameter estimate(s):
##       Phi 
## 0.1030767 
## Fixed effects: negaff ~ 1 + day + stresslevel.state + imean.stresslevel.c +      imean.stresslevel.c * day + imean.stresslevel.c * stresslevel.state 
##                                            Value  Std.Error   DF  t-value
## (Intercept)                            2.6958345 0.04706531 1244 57.27859
## day                                   -0.0660158 0.00893997 1244 -7.38434
## stresslevel.state                      0.7513456 0.04557269 1244 16.48675
## imean.stresslevel.c                    1.2080205 0.09850104  188 12.26404
## day:imean.stresslevel.c               -0.0498111 0.01869773 1244 -2.66402
## stresslevel.state:imean.stresslevel.c  0.1419311 0.09714895 1244  1.46096
##                                       p-value
## (Intercept)                            0.0000
## day                                    0.0000
## stresslevel.state                      0.0000
## imean.stresslevel.c                    0.0000
## day:imean.stresslevel.c                0.0078
## stresslevel.state:imean.stresslevel.c  0.1443
##  Correlation: 
##                                       (Intr) day    strss. imn.s. dy:m..
## day                                   -0.601                            
## stresslevel.state                      0.124  0.130                     
## imean.stresslevel.c                    0.009 -0.009  0.034              
## day:imean.stresslevel.c               -0.010  0.024 -0.048 -0.593       
## stresslevel.state:imean.stresslevel.c  0.033 -0.046 -0.114  0.113  0.149
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.46990102 -0.59357636 -0.07534063  0.48062925  4.59342234 
## 
## Number of Observations: 1438
## Number of Groups: 190
intervals(lg.fit.q6)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                                             lower        est.       upper
## (Intercept)                            2.60349831  2.69583447  2.78817062
## day                                   -0.08355488 -0.06601579 -0.04847670
## stresslevel.state                      0.66193781  0.75134564  0.84075346
## imean.stresslevel.c                    1.01371122  1.20802055  1.40232988
## day:imean.stresslevel.c               -0.08649371 -0.04981114 -0.01312857
## stresslevel.state:imean.stresslevel.c -0.04866276  0.14193112  0.33252500
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                                          lower       est.      upper
## sd((Intercept))                     0.37443385  0.4631407 0.57286315
## sd(day)                             0.03272832  0.0569535 0.09910993
## sd(stresslevel.state)               0.26917795  0.3555564 0.46965340
## cor((Intercept),day)               -0.61119789 -0.2651138 0.16607698
## cor((Intercept),stresslevel.state) -0.03180671  0.3312140 0.61703307
## cor(day,stresslevel.state)         -0.16048528  0.4183131 0.78303804
## 
##  Correlation structure:
##          lower      est.     upper
## Phi 0.01817013 0.1030767 0.1865073
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.5940599 0.6274212 0.6626561

Interpretation of output.

And lets make those interaction plots.

#xlevels = is the list of values that we want to evaluate at.
describe(daily.imeans$imean.stresslevel.c)
##    vars   n mean   sd median trimmed  mad   min  max range  skew kurtosis
## X1    1 190    0 0.48   0.01       0 0.51 -1.21 1.17  2.38 -0.04    -0.23
##      se
## X1 0.03
describe(daily$stresslevel.state)
##    vars    n mean   sd median trimmed  mad   min  max range skew kurtosis
## X1    1 1445    0 0.49  -0.03   -0.02 0.46 -1.75 2.12  3.88 0.36     0.79
##      se
## X1 0.01
#calculate effect
ef6 <- effect(term="stresslevel.state*imean.stresslevel.c", mod=lg.fit.q6, 
              xlevels=list(imean.stresslevel.c=c(-0.48, +0.48), stresslevel.state=c(-0.49,+0.49)))
summary(ef6)
## 
##  stresslevel.state*imean.stresslevel.c effect
##                  imean.stresslevel.c
## stresslevel.state    -0.48     0.48
##             -0.49 1.634391 2.560593
##             0.49  2.303946 3.363676
## 
##  Lower 95 Percent Confidence Limits
##                  imean.stresslevel.c
## stresslevel.state    -0.48     0.48
##             -0.49 1.525135 2.454759
##             0.49  2.166950 3.228815
## 
##  Upper 95 Percent Confidence Limits
##                  imean.stresslevel.c
## stresslevel.state    -0.48     0.48
##             -0.49 1.743647 2.666427
##             0.49  2.440942 3.498537

Everything we need is in there!

We can convert to a data frame and plot it!

#convert to dataframe
preddata6 <- as.data.frame(ef6)
#plotting the effect evaluation (with standard error ribbon)
ggplot(data=preddata6, aes(x=stresslevel.state, y=fit, group=imean.stresslevel.c), legend=FALSE) + 
  geom_point() +
  geom_line() +
  #geom_ribbon(aes(ymin=lower, ymax=upper), alpha=.3) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.15) +
  xlab("Stress Level State") + xlim(-2,2) +
  ylab("Predicted Negative Affect") + ylim(1,7) +
  ggtitle("Prediction from Continuous Predictor Model")

There we go!

We can now do everything up through Chapter 5 in Bolger & Laurenceau, plus more!