## Ben Hansen and Jake Bowers wrote a code snippet that I've found useful ## for calculating mahalanobis distances: library(optmatch) data(nuclearplants) nuclear.nopt <- nuclearplants myMH <- function(Tnms, Cnms, inv.cov, data) { stopifnot(!is.null(dimnames(inv.cov)[[1]]), dim(inv.cov)[1] > 1, all.equal(dimnames(inv.cov)[[1]], dimnames(inv.cov)[[2]]), all(dimnames(inv.cov)[[1]] %in% names(data))) covars <- dimnames(inv.cov)[[1]] xdiffs <- as.matrix(data[Tnms, covars]) xdiffs <- xdiffs - as.matrix(data[Cnms, covars]) rowSums((xdiffs %*% inv.cov) * xdiffs) } icv <- solve(cov(nuclear.nopt[, c("cap", "date")])) trtnms <- row.names(nuclear.nopt)[as.logical(nuclear.nopt$pr)] ctlnms <- row.names(nuclear.nopt)[!as.logical(nuclear.nopt$pr)] ## This line is where the magic happens mdist <- outer(trtnms, ctlnms, FUN = myMH, inv.cov = icv, data = nuclear.nopt) dimnames(mdist) <- list(trtnms, ctlnms) round(mdist, 2) sort(fullmatch(mdist))