# see example of the use of this function at the end. gethedgesg <-function( x1, x2, design = "between", coverage = 0.95) { # mandatory arguments are x1 and x2, both a vector of data require(psych) # for the function harmonic.mean. # get basic descriptive statistics ns <- c(length(x1),length(x2)) mns <- c(mean(x1),mean(x2)) sds <- c(sd(x1),sd(x2)) # get pairwise statistics ntilde <- harmonic.mean(ns) dmn <- abs(mns[2]-mns[1]) sdp <- sqrt( (ns[1]-1) *sds[1]^2 + (ns[2]-1)*sds[2]^2) / sqrt(ns[1]+ns[2]-2) # compute biased Cohen's d (equation 1) cohend <- dmn / sdp # compute unbiased Hedges' g (equations 2a and 3) eta <- ns[1] + ns[2] - 2 J <- gamma(eta/2) / (sqrt(eta/2) * gamma((eta-1)/2) ) hedgesg <- cohend * J # compute noncentrality parameter (equation 5a or 5b depending on the design) lambda <- if(design == "between") { hedgesg * sqrt( ntilde/2) } else { r <- cor(x1,x2) hedgesg * sqrt( ntilde/(2 * (1-r)) ) } # confidence interval of the hedges g (equations 6 and 7) tlow <- qt(1/2 - coverage/2, df = eta, ncp = lambda ) thig <- qt(1/2 + coverage/2, df = eta, ncp = lambda ) dlow <- tlow / lambda * hedgesg dhig <- thig / lambda * hedgesg # all done! display the results cat("Hedges'g = ", hedgesg, "\n", coverage*100, "% CI = [", dlow, dhig, "]\n") } x1 <- c(53, 68, 66, 69, 83, 91) x2 <- c(49, 60, 67, 75, 78, 89) # using the defaults: between design and 95% coverage gethedgesg(x1, x2) # changing the defaults explicitely gethedgesg(x1, x2, design = "within", coverage = 0.90 )