Cohen's d effect size

gistfile1.r

# usage
# R --vanilla < effect-size.r 
#
# requires ggplot2
# aptitude install r-cran-ggplot2
#
# remixed from http://rpsychologist.com/short-r-script-to-plot-effect-sizes-cohens-d-and-shade-overlapping-area

require("ggplot2")

mean0 <- 1

stddev <- 1

# Standardized Mean Difference = Cohen's d)
cohen_theory <- .05
mean_diff <- cohen_theory *stddev

# get mean2 depending on value of ES from d = (u1 - u2)/sd
mean1 <- mean_diff + mean0
# create x sequence
scale = 2
x <- seq(mean0 - scale*stddev, mean1 + scale*stddev, .01)
# generate normal dist #1
y1 <- dnorm(x, mean0, stddev)
# generate normal dist #2
y2 <- dnorm(x, mean1, stddev)

# put in data frame
df1 <- data.frame("x" = x, "y" = y1)

# put in data frame
df2 <- data.frame("x" = x, "y" = y2)


# http://stackoverflow.com/questions/15436702/estimate-cohens-d-for-effect-size
cohens_d <- function(x, y) {
    lx <- length(x)- 1
    ly <- length(y)- 1
    md  <- abs(mean(x) - mean(y))        ## mean difference (numerator)
    csd <- lx * var(x) + ly * var(y)
    csd <- csd/(lx + ly)
    csd <- sqrt(csd)                     ## common sd computation

    cd  <- md/csd                        ## cohen's d
}

a<-rnorm(3000,mean=mean0,sd=stddev)
b<-rnorm(3000,mean=mean1,sd=stddev)
print(mean0)
print(mean(a))
print(mean1)
print(mean(b))

d <-cohens_d(a,b)
print(d)
chart_title <-paste("Cohen's d = ",cohen_theory,"", sep="")



# plot with ggplot2
foo <- ggplot(df1, aes(x,y, color="treatment")) +
  # add line for treatment group
  geom_line(size=1) + 
  # add line for control group
  geom_line(data=df2, aes(color="control"),size=1) +
  # add vlines for group means
  geom_vline(xintercept = 1, linetype="dotted") + 
  geom_vline(xintercept = mean1, linetype="dotted") + 
  # add plot title
  labs(title=(chart_title)) +
  # change colors and legend annotation
  # scale_color_manual("Group", values= c("treatment" = "black","control" = "red")) +
  # remove axis labels
  ylab(NULL) + xlab(NULL)

ggsave(paste("cohen-d-",cohen_theory,"-effect-size.png", sep=""), height=4, width=6, dpi = 80)