R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

> # 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
'data.frame':   120 obs. of  9 variables:
 $ X        : int  1 2 3 4 5 6 7 8 9 10 ...
 $ patid    : int  1329 1349 1428 1502 1520 1342 1368 1401 1424 1460 ...
 $ p_gender : int  1 2 2 2 2 1 1 2 2 1 ...
 $ ypre     : int  26 31 24 25 24 19 41 14 20 12 ...
 $ y        : int  1 4 NA 4 4 2 1 NA NA 13 ...
 $ therapist: int  101 101 101 101 101 102 102 102 102 102 ...
 $ t_gender : int  1 1 1 1 1 1 1 1 1 1 ...
 $ bdi_gain : num  -7.69 -6.22 NA -4.39 -4.08 ...
 $ t        : Factor w/ 2 levels "CBT","IPT": 2 2 2 2 2 2 2 2 2 2 ...
> 
> 
> install.packages("nlme")  # (skip this line if you have already installed nlme)
--- Please select a CRAN mirror for use in this session ---
trying URL 'http://cran.irsn.fr/bin/windows/contrib/2.15/nlme_3.1-108.zip'
Content type 'application/zip' length 2119493 bytes (2.0 Mb)
opened URL
downloaded 2.0 Mb

package ‘nlme’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
        C:\Users\Kenny\AppData\Local\Temp\RtmpIVukkA\downloaded_packages
> library(nlme)    # Load nlme package
Warning message:
package ‘nlme’ was built under R version 2.15.3 
> # 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)
Linear mixed-effects model fit by REML
 Data: nimh2 
       AIC      BIC    logLik
  616.3183 623.5748 -305.1591

Random effects:
 Formula: ~1 | therapist
        (Intercept) Residual
StdDev:    2.785076 8.971209

Fixed effects: y ~ 1 
               Value Std.Error DF  t-value p-value
(Intercept) 8.671625  1.200477 66 7.223484       0

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.4239908 -0.7030746 -0.2414293  0.4067050  4.3330192 

Number of Observations: 84
Number of Groups: 18 
> # s_T = 2.7851     "(Intercept)"
> # s_E = 8.9712     "Residual"
> 2.7851^2/(2.78512^2 + 8.9712^2)  # Simple computation of ICC = .0879
[1] 0.08790618
> 
> # A little R programming to compute the ICC
> VarCorr(mod1)  # extracts both variances and sds from an object of class lme
therapist = pdLogChol(1) 
            Variance  StdDev  
(Intercept)  7.756648 2.785076
Residual    80.482583 8.971209
> varests <- as.numeric(VarCorr(mod1)[1:2])  # vector of variance estimates
> ICC <- varests[1]/sum(varests)
> ICC    # displays value
[1] 0.08790476
> 
> # 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)
[1] "ICC = 0.0879047551989659"
> 
> # 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)
Linear mixed-effects model fit by REML
 Data: nimh2 
       AIC      BIC    logLik
  615.5249 625.1517 -303.7624

Random effects:
 Formula: ~1 | therapist
        (Intercept) Residual
StdDev:    2.370091 8.823759

Fixed effects: y ~ ypre.c 
               Value Std.Error DF  t-value p-value
(Intercept) 8.841337 1.1330053 65 7.803438  0.0000
ypre.c      0.291651 0.1273675 65 2.289840  0.0253
 Correlation: 
       (Intr)
ypre.c 0.063 

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.4292448 -0.6020348 -0.2625727  0.3880986  4.2534321 

Number of Observations: 84
Number of Groups: 18 
> # s_T = 2.3700     "(Intercept)"
> # s_E = 8.8238     "Residual"
> 2.3700^2/(2.3700^2 + 8.8238^2)  # ICC = .067
[1] 0.06728732
> VarCorr(mod2)
therapist = pdLogChol(1) 
            Variance  StdDev  
(Intercept)  5.617333 2.370091
Residual    77.858717 8.823759
> ICClme(mod2)
[1] "ICC = 0.067292750435604"
> 
> # 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)
Generalized least squares fit by REML
  Model: y ~ 1 
  Data: nimh2 
       AIC      BIC    logLik
  616.3183 623.5748 -305.1591

Correlation Structure: Compound symmetry
 Formula: ~1 | therapist 
 Parameter estimate(s):
       Rho 
0.08790543 

Coefficients:
               Value Std.Error  t-value p-value
(Intercept) 8.671625  1.200478 7.223476       0

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-0.9231442 -0.7102328 -0.2844098  0.4607804  4.2931871 

Residual standard error: 9.393575 
Degrees of freedom: 84 total; 83 residual
> 
> # Model 4:  Including treatment effect
> mod4 <- lme(y ~ t, data=nimh, random  = ~ 1 | therapist,
+             na.action = "na.omit") 
> summary(mod4)
Linear mixed-effects model fit by REML
 Data: nimh 
       AIC      BIC    logLik
  612.4669 622.0938 -302.2335

Random effects:
 Formula: ~1 | therapist
        (Intercept) Residual
StdDev:    2.366985 8.978599

Fixed effects: y ~ t 
                Value Std.Error DF   t-value p-value
(Intercept) 10.698848  1.728778 66  6.188676  0.0000
tIPT        -3.596908  2.307256 16 -1.558955  0.1386
 Correlation: 
     (Intr)
tIPT -0.749

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.4889673 -0.6507164 -0.2521612  0.5105497  4.2079855 

Number of Observations: 84
Number of Groups: 18 
> intervals(mod4)
Approximate 95% confidence intervals

 Fixed effects:
                lower      est.     upper
(Intercept)  7.247231 10.698848 14.150466
tIPT        -8.488072 -3.596908  1.294256
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: therapist 
                    lower     est.   upper
sd((Intercept)) 0.6218186 2.366985 9.01005

 Within-group standard error:
    lower      est.     upper 
 7.607386  8.978599 10.596969 
> 
> # Model 5:  Inclusion of patient gender
> mod5 <- lme(y ~ 1, data=nimh2, random  = ~ p_gender | therapist,
+             na.action = "na.omit")
> summary(mod5)
Linear mixed-effects model fit by REML
 Data: nimh2 
       AIC      BIC   logLik
  620.2919 632.3861 -305.146

Random effects:
 Formula: ~p_gender | therapist
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev       Corr  
(Intercept) 0.0008969124 (Intr)
p_gender    1.4922377434 0.119 
Residual    9.0056538954       

Fixed effects: y ~ 1 
               Value Std.Error DF  t-value p-value
(Intercept) 8.610599  1.164232 66 7.395949       0

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.4506811 -0.7062470 -0.2252426  0.4125560  4.4500517 

Number of Observations: 84
Number of Groups: 18