# kh.R R script for models discussed in Kenny & Hoyt (2009) nimh <- read.csv(file.choose()) # Select downloaded file "nimh.csv" str(nimh) # see variable names and types install.packages("nlme") # (skip this line if you have already installed nlme) library(nlme) # Load nlme package # Optional (and optimal): Create a "grouped data" object, specifying nesting factor(s) nimh2 <- groupedData(y ~ t | therapist, data=nimh, outer= ~t) # Model 1: Ignoring ypre (pre-scores) mod1 <- lme(y ~ 1, data=nimh2, random = ~ 1|therapist, na.action="na.omit") # (empty model) summary(mod1) # s_T = 2.7851 "(Intercept)" # s_E = 8.9712 "Residual" 2.7851^2/(2.78512^2 + 8.9712^2) # Simple computation of ICC = .0879 # A little R programming to compute the ICC VarCorr(mod1) # extracts both variances and sds from an object of class lme varests <- as.numeric(VarCorr(mod1)[1:2]) # vector of variance estimates ICC <- varests[1]/sum(varests) ICC # displays value # Simple function for computing ICC from lme() output ICClme <- function(out) { varests <- as.numeric(VarCorr(out)[1:2]) return(paste("ICC =", varests[1]/sum(varests))) } ICClme(mod1) # Model 2: Controlling for (centered) pre-scores nimh2$ypre.c <- with(nimh, ypre - mean(ypre)) # ypre.c = centered pre-score mod2 <- lme(y ~ ypre.c, data=nimh2, random = ~ 1 | therapist, na.action="na.omit") summary(mod2) # s_T = 2.3700 "(Intercept)" # s_E = 8.8238 "Residual" 2.3700^2/(2.3700^2 + 8.8238^2) # ICC = .067 VarCorr(mod2) ICClme(mod2) # Model 3: Pt's as repeated measures for each therapist # (Allows for possibility that ICC < 0) mod3 = gls(y ~ 1, data=nimh2, na.action="na.omit",verbose=TRUE, correlation=corCompSymm(form=~1|therapist)) summary(mod3) # ICC (called "Rho") = .0879 (like Model 1) # Model 4: Including treatment effect mod4 <- lme(y ~ t, data=nimh, random = ~ 1 | therapist, na.action = "na.omit") summary(mod4) intervals(mod4) # Model 5: Inclusion of patient gender mod5 <- lme(y ~ 1, data=nimh2, random = ~ p_gender | therapist, na.action = "na.omit") summary(mod5)