Overview

This session we work through the probing of interactions (moderation) in a growth modeling framework. In particular we provide exposure to the Johnson-Neyman method for identifying zones of significance. We used existing packages and plotting function.

Outline

  1. Data Preparation and Description
  2. Linear Growth Model
  3. Conditional Growth Model and Probing the “Interaction”

Prelim - Loading libraries used in this script.

library(psych)
library(ggplot2)
library(nlme) #for mixed effects models
library(reghelper)  #for probing interactions

1. Data Preparation and Description

For our example, we use 4-occasion WISC data.

Load the repeated measures data

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

Subsetting to the variables of interest. Specifically, we include the id variable; the repeated measures outcome variables verb1, verb2, verb4, verb6; and the predictors grad and momed variables.

#Creating a list of variables of interest
varnames <- c("id","verb1","verb2","verb4","verb6","grad","momed")
#sub-setting the columns of interest
wiscsub <- wisc3raw[ ,varnames]

#Looking at the sample-level descriptives for the chosen variables
describe(wiscsub)
##       vars   n   mean    sd median trimmed   mad   min    max  range  skew
## id       1 204 102.50 59.03 102.50  102.50 75.61  1.00 204.00 203.00  0.00
## verb1    2 204  19.59  5.81  19.34   19.50  5.41  3.33  35.15  31.82  0.13
## verb2    3 204  25.42  6.11  25.98   25.40  6.57  5.95  39.85  33.90 -0.06
## verb4    4 204  32.61  7.32  32.82   32.42  7.18 12.60  52.84  40.24  0.23
## verb6    5 204  43.75 10.67  42.55   43.46 11.30 17.35  72.59  55.24  0.24
## grad     6 204   0.23  0.42   0.00    0.16  0.00  0.00   1.00   1.00  1.30
## momed    7 204  10.81  2.70  11.50   11.00  2.97  5.50  18.00  12.50 -0.36
##       kurtosis   se
## id       -1.22 4.13
## verb1    -0.05 0.41
## verb2    -0.34 0.43
## verb4    -0.08 0.51
## verb6    -0.36 0.75
## grad     -0.30 0.03
## momed     0.01 0.19

We note the characteristics of our person-level predictors, grad and momed. We will want to retain information about the standard deviation, range, and plausible values of these variables (most easily done in a wide data set).

#descriptives of person-level predictors
table(wiscsub$grad)
## 
##   0   1 
## 158  46
describe(wiscsub$momed)
##    vars   n  mean  sd median trimmed  mad min max range  skew kurtosis
## X1    1 204 10.81 2.7   11.5      11 2.97 5.5  18  12.5 -0.36     0.01
##      se
## X1 0.19

Growth models (in the multilevel framework) require a long data set. So, we also reshape from wide to long in order to have a long data set.

#reshaping wide to long
verblong <- reshape(data=wiscsub, 
                     varying=c("verb1","verb2","verb4","verb6"), 
                     timevar="grade", idvar="id", 
                     direction="long", sep="")
#reordering for convenience
verblong <- verblong[order(verblong$id,verblong$grade),c("id","grade","verb","grad","momed")]
#looking at dataset
head(verblong,12)
##     id grade  verb grad momed
## 1.1  1     1 24.42    0   9.5
## 1.2  1     2 26.98    0   9.5
## 1.4  1     4 39.61    0   9.5
## 1.6  1     6 55.64    0   9.5
## 2.1  2     1 12.44    0   5.5
## 2.2  2     2 14.38    0   5.5
## 2.4  2     4 21.92    0   5.5
## 2.6  2     6 37.81    0   5.5
## 3.1  3     1 32.43    1  14.0
## 3.2  3     2 33.51    1  14.0
## 3.4  3     4 34.30    1  14.0
## 3.6  3     6 50.18    1  14.0
#plotting the individual trajectories

#plotting intraindividual change
ggplot(data = verblong, aes(x = grade, y = verb, group = id)) +
  geom_point() + 
  geom_line() +
  xlab("Grade") + 
  ylab("WISC Verbal Score") + ylim(0,100) +
  scale_x_continuous(breaks=seq(1,6,by=1))

2. Linear Growth Model (Multilevel Model of Change)

We look at the within-person linear relation beween the time variable (grade) and the outcome variable (verb).

Linear growth model (MLM): random intercepts and random slopes.

model0 <- lme(fixed = verb ~ 1 + grade, 
                 random = ~ 1 + grade|id, 
                 data=verblong,
                 na.action = na.exclude)

summary(model0)
## Linear mixed-effects model fit by REML
##  Data: verblong 
##        AIC      BIC    logLik
##   5053.632 5081.844 -2520.816
## 
## Random effects:
##  Formula: ~1 + grade | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 3.915536 (Intr)
## grade       1.241299 0.321 
## Residual    3.581590       
## 
## Fixed effects: verb ~ 1 + grade 
##                Value Std.Error  DF  t-value p-value
## (Intercept) 15.15099 0.3686513 611 41.09844       0
## grade        4.67339 0.1087023 611 42.99255       0
##  Correlation: 
##       (Intr)
## grade -0.155
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.61615621 -0.54270691 -0.02898115  0.52402090  3.16468334 
## 
## Number of Observations: 816
## Number of Groups: 204
VarCorr(model0)
## id = pdLogChol(1 + grade) 
##             Variance  StdDev   Corr  
## (Intercept) 15.331421 3.915536 (Intr)
## grade        1.540823 1.241299 0.321 
## Residual    12.827783 3.581590
intervals(model0)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                 lower     est.     upper
## (Intercept) 14.427016 15.15099 15.874970
## grade        4.459914  4.67339  4.886865
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                              lower      est.     upper
## sd((Intercept))        3.255506646 3.9155359 4.7093810
## sd(grade)              1.057964280 1.2412989 1.4564035
## cor((Intercept),grade) 0.006320221 0.3205334 0.5771394
## 
##  Within-group standard error:
##    lower     est.    upper 
## 3.343438 3.581590 3.836705

We can obtain predicted scores and plot them

#obtaining predicted scores for individuals
verblong$pred.model0 <- predict(model0)

#obtaining predicted scores for prototype
verblong$proto.model0 <- predict(model0, level=0)

#plotting predicted trajectories
ggplot(data = verblong, aes(x = grade, y = pred.model0, group = id)) +
  #geom_point() + 
  geom_line(color="gray", alpha=.8) +
  geom_line(aes(x = grade, y = proto.model0), color="red",size=2) +
  xlab("Grade") + 
  ylab("WISC Verbal Score") + ylim(0,100) +
  scale_x_continuous(breaks=seq(1,6,by=1))