0.1 Overview

This session we work through basics of growth modeling. We describe the data, run some individual level models, and work through no-growth and linear growth modeling examples. We then expand on that model in two ways. We add a predictor - to get the conditional growth model; and we change the time-metric - to illustrate how alternative time metrics facilitate different interpretations.

0.2 Outline

  1. Data Preparation and Description
  2. Individual Growth Models
  3. Growth Model (Multilevel Model of Change)
  1. No-Growth Model
  2. Linear Growth Model
  1. Comparing MLM and individual growth results.
  2. Conditional Growth Model
  3. Growth Models with Alternative Time Metrics

0.2.0.1 Prelim - Loading libraries used in this script.

library(psych)
library(ggplot2)
library(data.table) #for fast data management
library(nlme) #for mixed effects models
library(plyr) #for data management

0.3 1. Data Preparation and Description

For our examples, 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 kurtosis se
id 1 204 102.5000000 59.0338886 102.500 102.5000000 75.612600 1.00 204.00 203.00 0.0000000 -1.2176609 4.1331989
verb1 2 204 19.5850490 5.8076950 19.335 19.4976829 5.411490 3.33 35.15 31.82 0.1301071 -0.0528376 0.4066200
verb2 3 204 25.4153431 6.1063772 25.980 25.4026220 6.567918 5.95 39.85 33.90 -0.0639307 -0.3445494 0.4275319
verb4 4 204 32.6077451 7.3198841 32.820 32.4240244 7.183197 12.60 52.84 40.24 0.2317998 -0.0776524 0.5124944
verb6 5 204 43.7499020 10.6650511 42.545 43.4560976 11.297412 17.35 72.59 55.24 0.2356459 -0.3605210 0.7467029
grad 6 204 0.2254902 0.4189328 0.000 0.1585366 0.000000 0.00 1.00 1.00 1.3040954 -0.3007374 0.0293312
momed 7 204 10.8112745 2.6982790 11.500 10.9969512 2.965200 5.50 18.00 12.50 -0.3586143 0.0056095 0.1889173

Multilevel modeling analyses typically 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

0.4 2. Individual Growth Models

Lets make a dataset with just 1 person of interest, id = 23

#Subsetting persons
verb_id23 <- verblong[which(verblong$id == 23), ]

Lets make a plot of this person’s data

#plotting intraindividual change
ggplot(data = verb_id23, 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))

Lets do an individual regression with time as a predictor. Conceptually, this is a model of intraindividual change.

#regress verb on grade 
linear_id23 <- lm(formula = verb ~ 1 + grade, 
                  data = verb_id23, 
                  na.action=na.exclude)
#show results
summary(linear_id23) 
## 
## Call:
## lm(formula = verb ~ 1 + grade, data = verb_id23, na.action = na.exclude)
## 
## Residuals:
##    23.1    23.2    23.4    23.6 
## -2.6556  3.5536 -0.4681 -0.4298 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  21.3247     3.1147   6.846   0.0207 *
## grade         1.6308     0.8251   1.977   0.1867  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.169 on 2 degrees of freedom
## Multiple R-squared:  0.6614, Adjusted R-squared:  0.4921 
## F-statistic: 3.907 on 1 and 2 DF,  p-value: 0.1867

Lets save the 3 parameters into objects and look at them.

id23_reg_linear <- as.list(coef(linear_id23))
id23_rse_linear <- summary(linear_id23)$sigma #to obtain residual standard error

id23_reg_linear
## $`(Intercept)`
## [1] 21.32475
## 
## $grade
## [1] 1.630847
id23_rse_linear
## [1] 3.168901

Lets add a regression line to the plot

#plotting intraindividual change
ggplot(data = verb_id23, aes(x = grade, y = verb, group = id)) +
  geom_point() + 
  geom_line() +
  geom_smooth(method=lm, se=FALSE,colour="red", size=1) +
  xlab("Grade") + 
  ylab("WISC Verbal Score") + ylim(0,100) +
  scale_x_continuous(breaks=seq(1,6,by=1))

#note that geom_smooth() uses its own calculation, not the one from the model above

Now lets do the same thing for all the persons!

We do this in a speedy way using the data.table package.

#converting to a data.table object
verblong_dt <- data.table(verblong)
#collecting regression output by id 
indiv_reg <- verblong_dt[,c(reg_1 = as.list(coef(lm(verb ~ grade))), 
                            reg_1_sigma = as.list(summary(lm(verb ~ grade))$sigma)), by=id]

Lets look at the distributions/correlation of the parameters

#converting back to data.frame 
names(indiv_reg)
## [1] "id"                "reg_1.(Intercept)" "reg_1.grade"      
## [4] "reg_1_sigma"
indiv_reg_data <- as.data.frame(indiv_reg)

#descriptives
describe(indiv_reg_data[-1])
vars n mean sd median trimmed mad min max range skew kurtosis se
reg_1.(Intercept) 1 204 15.150993 5.265393 15.138644 15.231877 5.119996 0.5020339 29.413729 28.911695 -0.0786938 0.0565420 0.3686513
reg_1.grade 2 204 4.673390 1.552579 4.599831 4.613264 1.525068 1.4115254 9.380678 7.969152 0.3764694 0.0229344 0.1087023
reg_1_sigma 3 204 3.146516 1.715122 2.888048 3.004732 1.500142 0.0630711 8.483826 8.420754 0.7937263 0.4596089 0.1200826
#correlations among parameters
cor(indiv_reg_data[-1], use="complete.obs",method="pearson")
##                   reg_1.(Intercept) reg_1.grade reg_1_sigma
## reg_1.(Intercept)         1.0000000  -0.1551763   0.1060735
## reg_1.grade              -0.1551763   1.0000000   0.1927477
## reg_1_sigma               0.1060735   0.1927477   1.0000000
#pairs in the psych library
pairs.panels(indiv_reg_data[-1])