Modified McCarthy Data for the Effect Size Puzzler

12 subjects (Measured four times, twice per condition)

Condition one is measured at early time points 30 and 36 months
Condition two is measured at later time points 42 and 48 months

Method #1 – Classic Cohen’s d

library(readr)
dat <- read_csv("~/R/DataHouse/Delaney/SimulateMLM/McCarthyMod.csv", col_names = FALSE)
dat$mo <- c(rep(30, 12), rep(36, 12), rep(42, 12), rep(48, 12))
dat$ID <- rep(1:12, 4)
colnames(dat) <- c("AgeGRP", "Score", "Month", "Id")

condition.means <- with(dat, tapply(Score, AgeGRP, mean))
condition.sd <- with(dat, sqrt(mean(tapply(Score, AgeGRP, var))))


Cohen.d <- diff(condition.means)/condition.sd
names(Cohen.d) <- "Cohen.d"
print(Cohen.d)

  Cohen.d 
0.4350369

Method 2 – Average over replicates

sub_means <- with(dat, tapply(Score, list(Id, AgeGRP), mean))
sd.a <- sqrt(mean(apply(sub_means, 2, var)))
d_a <- diff(condition.means)/sd.a
names(d_a) <- "d_a"
print(d_a)

      d_a 
0.4488602

Method 3 – Standerdized difference scores

sd.z <- sd(apply(sub_means, 1, diff))
d_z <- diff(condition.means)/sd.z
names(d_z) <- "d_z"
print(d_z)

      d_z 
0.5655712

Method 4 – Naive conversion from the test statistic

t <- t.test(sub_means[, 2] - sub_means[, 1])$statistic
d_t <- t * sqrt(2/nrow(sub_means))
names(d_t) <- "d_t"
print(d_t)

      d_t 
0.7998384

Method 5 – Using the variance within individuals (residual) – between removed

library(nlme)
mem1 <- lme(fixed = Score ~ AgeGRP, random = ~1 | Id, method = "ML", data = dat)
summary(mem1)

Linear mixed-effects model fit by maximum likelihood
 Data: dat 
       AIC      BIC   logLik
  366.8059 374.2907 -179.403

Random effects:
 Formula: ~1 | Id
        (Intercept) Residual
StdDev:    11.10118 7.684761

Fixed effects: Score ~ AgeGRP 
            Value Std.Error DF   t-value p-value
(Intercept)    99  4.853282 35 20.398567  0.0000
AgeGRP          6  2.266112 35  2.647706  0.0121
 Correlation: 
       (Intr)
AgeGRP -0.7  

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-1.789111263 -0.653198682 -0.009377154  0.726143382  1.538521222 

Number of Observations: 48
Number of Groups: 12

sd.r <- VarCorr(mem1)
d_r <- diff(condition.means)/as.numeric(sd.r[2, 2])
names(d_r) <- "d_r"
print(d_r)

     d_r 
0.780766

Method 6 – Using the sum of the variance components

sd.vc <- sqrt(sum(as.numeric(sd.r[, 1])))
d_vc <- diff(condition.means)/sd.vc
names(d_vc) <- "d_vc"
print(d_vc)

     d_vc 
0.4443936

Method 7 – Solution advocated for in Maxwell and Delaney

grp_var <- sqrt(sum(with(dat, tapply(Score, Month, var)))/4)
d_md <- diff(condition.means)/grp_var
names(d_md) <- "d_md"
d_md

     d_md 
0.4284224