## 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.

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

#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))