= R syntax for bootstrap confidence intervals for effect sizes =

The bootstrap can be used to obtain confidence intervals for effect sizes by repeatedly sampling from the data a large number of times and obtaining percentiles corresponding to a given confidence interval.

The SPSS spreadsheet file, [attachment:skewd.sav skewed data.sav] is passed into R and a 95% confidence interval for Cohen's d evaluated and plotted representing the standardised difference between males and females in the sample.
This also will work on EXCEL spreadsheet files as they can be entered  into R via SPSS. 

{{{
library(foreign)
x <- read.spss("skewed data.sav")
x1 <- data.frame(x)
y <- x1$time_o2
y <- na.omit(y)

# Adapted syntax in Gilbert MacKenzie's Crash course on R handout

#
nb<-1000
n <- length(y)
boot<-matrix(NA,nb,6)
attributes(boot)
for (i in 1:nb)
{
bs<-sample(y,n,replace =T)
boot[i,1]<-mean(bs)
boot[i,2]<-var(bs)
}
boot[,3]<-sqrt(boot[,2])
#
#look at shape of bootstrap distributions
#
hist(boot[,1])
hist(boot[,2])
hist(boot[,3])

#
# Do separately for Males and Females as in RW macro
# This is R equivalent of the SPSS macro of RW (2004)

nb <- 1000
ymale <- matrix(NA,length(y),1)

for (i in 1:length(y)) {
if (x1$sex[i] == 'Male') {
 ymale[i] <- y[i] }
}
ymale <- na.omit(ymale)

nm <- length(ymale)

nm <- length(ymale)

for (i in 1:nb)
{
bs<-sample(ymale,nm,replace =T)
boot[i,1]<-mean(bs)
boot[i,2]<-var(bs)
}
boot[,3]<-sqrt(boot[,2])
#
#look at shape of bootstrap distributions
#
hist(boot[,1])
hist(boot[,2])
hist(boot[,3])
}

yfemale <- matrix(NA,length(y),1)

for (i in 1:length(y)) {
if (x1$sex[i] == 'Female') {
 yfemale[i] <- y[i] }
}
yfemale <- na.omit(yfemale)

nf <- length(yfemale)

for (i in 1:nb)
{
bs<-sample(yfemale,nf,replace =T)
boot[i,4]<-mean(bs)
boot[i,5]<-var(bs)
}
boot[,6]<-sqrt(boot[,5])
#
#look at shape of bootstrap distributions
#
hist(boot[,4])
hist(boot[,5])
hist(boot[,6])
}

poolsd <- matrix(NA,1000,1)
cohensd <- matrix(NA,1000,1)

for (i in 1:nb) {
poolsd[i] <- ( (nm-1)*boot[i,3]+(nf-1)*boot[i,6] )/(nm+nf-2)

cohensd[i] <- (boot[i,1] - boot[i,4])/poolsd[i]
}
hist(cohensd)
quantile(cohensd,0.05)
quantile(cohensd,0.95)
}}}