This tutorial extends examination of*intraindividual 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.

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

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.

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.

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