Overview

The MultiLevel Model (MLM) was developed to accommodate the interdependence in nested data. It is particularly useful for examining individual differences in bivariate associations - thus providing for more robust statistical inference about within-person associations (intraindividual covariation) and the between-person differences therein.

Outline

A. Overview of model

B. Set up the basic multilevel model

C. Interpretation of results

D. Display (data visualization) of the result

A. Introduction to the Generalized Multilevel Model

The basic linear multilevel model is written as

\[y_{it} = \beta_{0i} + \beta_{1i}x_{it}\] where \(\beta_{1i}\) is individual i’s level of … -ity, and where

\[\beta_{0i} = \gamma_{00} + \gamma_{01}z_{i} + u_{0i}\] \[\beta_{1i} = \gamma_{10} + \gamma_{11}z_{i} + u_{1i}\]

where \[e_{it} \tilde{} N(0,\mathbf{R})\], and \[\mathbf{U}_{i} \tilde{} N(0,\mathbf{G})\]

\[\mathbf{R} = \mathbf{I}\sigma^{2}_{e}\], where where \(I\) is the identity matrix (diagonal matrix of 1s) and \(\sigma^{2}_{e}\) is the residual (within-person) variance.

\[\mathbf{G} = \left[\begin{array} {rr} \sigma^{2}_{u0} & \sigma_{u0u1} \\ \sigma_{u1u0} & \sigma^{2}_{u1} \end{array}\right]\]

B. Set up the Basic Multilevel Model

Preliminaries

Loading Libraries

Loading libraries used in this script.

library(psych)      # for describing the data
library(plyr)       # for data manipulation
library(ggplot2)    # for data visualization
library(lme4)       # for multilevel models
library(lmerTest)   # for p-values

Loading Data

Our example makes use of the person-level and day-level (EMA-type) AMIB data files. Specifically, we make use of three varaibles:

Outcome: daily negative affect (continuous), obtained through negaff in the day-level file;

Predictor: daily stress, obtained through reverse coding of pss in the day-level file;
Person-level predictor: trait neuroticism, bfi_n in the person-level file

Loading person-level file and subsetting to variables of interest

#set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/AMIBshare_persons_2019_0501.csv"
#read in the .csv file using the url() function
AMIB_persons <- read.csv(file=url(filepath),header=TRUE)

#subsetting to variables of interest
AMIB_persons <- AMIB_persons[ ,c("id","bfi_n")]

Loading day-level file (T = 8) and subsetting to variables of interest.

#set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/AMIBshare_daily_2019_0501.csv"
#read in the .csv file using the url() function
AMIB_daily <- read.csv(file=url(filepath),header=TRUE)

#subsetting to variables of interest
AMIB_daily <- AMIB_daily[ ,c("id","day","negaff","pss")]

Data Preparation

Reverse code pss so higher values indicate higher perceived stress.

#reverse coding the pss variable into a new stress variable
AMIB_daily$stress <- 4 - AMIB_daily$pss
#describing new variable
describe(AMIB_daily$stress)
##    vars    n mean   sd median trimmed  mad min max range skew kurtosis
## X1    1 1445 1.39 0.68   1.25    1.36 0.74   0   4     4 0.35     0.13
##      se
## X1 0.02
#histogram
ggplot(data=AMIB_daily, aes(x=stress)) +
  geom_histogram(fill="white", color="black",bins=20) +
  labs(x = "Stress (high = more stressed)")

We now split the predictor into “trait” (between-person differences) and “state” (within-person deviations) components. Specifically, the daily variable stress is split into two varaibles: stress_trait is the sample-mean centered between-person component, and stress_state is the person-centered within-person component.

Making trait variables.

#calculating intraindividual means
#Alternative approach using dplyr
    #AMIB_imeans <- AMIB_daily %>% group_by(id) %>% summarize(stress_trait=mean(stress, na.rm=TRUE))
AMIB_imeans <- ddply(AMIB_daily, "id", summarize,
                       stress_trait = mean(stress, na.rm=TRUE),
                       negaff_trait = mean(negaff, na.rm=TRUE))
describe(AMIB_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
## stress_trait    2 190   1.40   0.48   1.41    1.39   0.51   0.19   2.56
## negaff_trait    3 190   2.48   0.73   2.41    2.43   0.72   1.11   5.09
##               range  skew kurtosis   se
## id           431.00 -0.04    -1.09 9.46
## stress_trait   2.38 -0.04    -0.23 0.03
## negaff_trait   3.98  0.68     0.45 0.05
#merging into person-level file
AMIB_persons <- merge(AMIB_persons, AMIB_imeans, by="id")                                              
#make centered versions of the person-level scores
AMIB_persons$bfi_n_c <- scale(AMIB_persons$bfi_n,center=TRUE,scale=FALSE)
AMIB_persons$stress_trait_c <- scale(AMIB_persons$stress_trait,center=TRUE,scale=FALSE)
#describe person-level data
describe(AMIB_persons)
##                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
## bfi_n             2 190   2.98   0.96   3.00    3.00   1.48   1.00   5.00
## stress_trait      3 190   1.40   0.48   1.41    1.39   0.51   0.19   2.56
## negaff_trait      4 190   2.48   0.73   2.41    2.43   0.72   1.11   5.09
## bfi_n_c           5 190   0.00   0.96   0.02    0.02   1.48  -1.98   2.02
## stress_trait_c    6 190   0.00   0.48   0.01    0.00   0.51  -1.21   1.17
##                 range  skew kurtosis   se
## id             431.00 -0.04    -1.09 9.46
## bfi_n            4.00 -0.09    -0.82 0.07
## stress_trait     2.38 -0.04    -0.23 0.03
## negaff_trait     3.98  0.68     0.45 0.05
## bfi_n_c          4.00 -0.09    -0.82 0.07
## stress_trait_c   2.38 -0.04    -0.23 0.03

Making state variables in long data

#merging person-level data into daily data
daily_long <- merge(AMIB_daily,AMIB_persons,by="id")

#calculating state variables
daily_long$stress_state <- daily_long$stress - daily_long$stress_trait
#describing data
describe(daily_long)
##                vars    n   mean     sd median trimmed    mad    min    max
## id                1 1458 322.53 129.08 324.00  324.20 151.23 101.00 532.00
## day               2 1458   3.48   2.30   3.00    3.47   2.97   0.00   7.00
## negaff            3 1441   2.45   1.04   2.20    2.34   1.04   1.00   6.90
## pss               4 1445   2.61   0.68   2.75    2.64   0.74   0.00   4.00
## stress            5 1445   1.39   0.68   1.25    1.36   0.74   0.00   4.00
## bfi_n             6 1458   2.97   0.96   3.00    2.98   1.48   1.00   5.00
## stress_trait      7 1458   1.39   0.47   1.41    1.39   0.51   0.19   2.56
## negaff_trait      8 1458   2.45   0.71   2.38    2.41   0.67   1.11   5.09
## bfi_n_c           9 1458  -0.01   0.96   0.02    0.00   1.48  -1.98   2.02
## stress_trait_c   10 1458  -0.01   0.47   0.01   -0.01   0.51  -1.21   1.17
## stress_state     11 1445   0.00   0.49  -0.03   -0.02   0.46  -1.75   2.12
##                 range  skew kurtosis   se
## id             431.00 -0.07    -1.06 3.38
## day              7.00  0.01    -1.24 0.06
## negaff           5.90  0.96     0.77 0.03
## pss              4.00 -0.35     0.13 0.02
## stress           4.00  0.35     0.13 0.02
## bfi_n            4.00 -0.08    -0.79 0.03
## stress_trait     2.38 -0.08    -0.21 0.01
## negaff_trait     3.98  0.66     0.52 0.02
## bfi_n_c          4.00 -0.08    -0.79 0.03
## stress_trait_c   2.38 -0.08    -0.21 0.01
## stress_state     3.88  0.36     0.79 0.01

Great! Now the data are all prepared as needed.

Set-up for basic multilevel model with continuous outcome

Outlining the substantive inquiry.

We define stress reactivty (a person-level dynamic characteristic) as the extent to which a person’s level of daily negative affect is related to daily stress level. Stress reactivity is the intraindividual covariation of daily negative affect and daily stress.

Formally,
\[y_{it} = \beta_{0i} + \beta_{1i}x_{it}\] First check out the distribution of the outcome variable (negaff):

#pdf("NegativeAffectHistogram.pdf",height=3, width=3.5)
ggplot(data=AMIB_daily, aes(x=negaff)) +
  geom_histogram(aes(y=..density..),      
                   binwidth=.25,
                   colour="black", fill="white") +
  labs(x = "Negative Affect", y="Density") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=16))
## Warning: Removed 17 rows containing non-finite values (stat_bin).

#dev.off()

Now look at a subset of persons to see what stress reactivity looks like.

#faceted plot
ggplot(data=daily_long[which(daily_long$id <= 125),], aes(x=stress_state,y=negaff)) +
  geom_point() +
  stat_smooth(method="lm", fullrange=TRUE) +
  xlab("Stress") + ylab("Negative Affect (Continuous)") + 
  facet_wrap( ~ id) +
  theme(axis.title=element_text(size=16),
        axis.text=element_text(size=14),
        strip.text=element_text(size=14))