Cohen's d calculation in R I am using following code to find Cohen's d (means are m1 and m2, standard deviations are s1 and s2 while sample sizes are n1 and n2):
lx <- n1- 1
ly <- n2- 1
md  <- abs(m1-m2)        ## mean difference (numerator)
csd <- lx * s1^2 + ly * s2^2
csd <- csd/(lx + ly)
csd <- sqrt(csd)                     ## common sd computation

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

For the example given on http://en.wikipedia.org/wiki/Effect_size#Cohen.27s_d (men's and women's heights by Aaron et al) I am getting a value of 1.756 while example states the value to be 1.72. Although the difference is small, I want to know if there is some mistake I am making in the code above?
 A: It looks like the calculation on the website uses an unweighted estimate of the pooled variance. Try your calculation again with equal ns of any size over 1 and you'll see that the answer is 1.72. You might want to edit the wikipedia page.
As to whether you're doing anything wrong, you're not. However, Cohen's d isn't well specified in terms of equation. Always state the equation you used when reporting Cohen's d.
A: In a recent article, we made this R function available which computes Hedges's g (an unbiased version of Cohen's) along with confidence interval:
gethedgesg <-function( x1, x2, design = "between", coverage = 0.95) {
  # mandatory arguments are x1 and x2, both a vector of data

  require(psych) # for the functions SD and harmonic.mean.

  # store the columns in a dataframe: more convenient to handle one variable than two
  X <- data.frame(x1,x2)

  # get basic descriptive statistics
  ns  <- lengths(X)
  mns <- colMeans(X)
  sds <- SD(X)

  # 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(X)[1,2]
    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")

}

It can be used for example with:
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 )

I hope it helps.
A: To add to Denis's helpful answer and to hopefully make this more easily discoverable (I just spent half an hour finding a solution): gethedgesg() currently fails with large samples (N > 340ish) - replacing the calculation for J as follows makes the function work more robustly:
J <- exp ( lgamma(eta/2) - (log(sqrt(eta/2)) + lgamma((eta-1)/2) ) )
