`compute_rank_mahalanobis` ignores index argument

[benthestatistician]

...in the event that there's an exact matching constraint in force. The function presently computes and returns a distance for each treatment/control combination, even when provided a list of specific treatment-control matching possibilities for which a distances is needed.

to do's:

• If possible, fix it in compute_rank_mahananobis
• If not possible, workaround on R side
• Add a test case to test.rank.mahal.R

[benthestatistician]

In the example in #122 , we have

> length(within@rows)
[1] 4889857

whereas compute_rank_mahalanobis has give a result of length

length(dists)
[1] 22887076

matching the product of the overall number of rows and columns.

> within@dimension

[1] 2252 10163

[benthestatistician]

@jwbowers since you've caught the rank_mahalanobis bug, would you be willing to help out w/ writing of test cases?

I figured I'd start by writing out some tests where you do a rank_mahalanobis calculation by pre-computing the ranks in R, then applying ordinary Mahalanobis. That led me to try the following, which does not work, but which is more likely to represent a misunderstanding on my part of what rank_mahalanobis is trying to do than a problem with the code. I.e., maybe you can help me improve it w/o much effort.

test_that("compute_rank.mahal results match ordinary Mahalanobis's", {
   nr <- 10L
   z <- integer(nr)
    z[sample(1:nr, nr / 2L)] <- 1L

   X <- as.matrix(1L:nr)
   df <- data.frame(z = z, X)
   expect_equivalent(match_on(z~., data=df, method="rank_mahalanobis"),
                     match_on(z~., data=df, method="mahalanobis")
                     )
}
)

(Once I get something like this going, should be easy to write test cases for the combo of exact matching and rank Mahalanobis.)

[benthestatistician]

In tests/testthat/test.rank.mahal.R, compute_rank_mahalanobis is tested against the rank Mahalanobis function that Rosenbaum has posted to his site. I confirmed that their results match up in this example also, even as they differ from those of compute_mahalanobis.

The likely cause of the difference is underlying differences in how the covariance is being computed. optmatch has always calculated the covariances needed for Mahalanobis distance computations by pooling covariances in treatment and control groups. In defining rank-based Mahalanobis, Rosenbaum uses an unpooled covariance. Anyway there's no indication of anything amiss, other than the problem with compute_rank_mahalanobis noted at the top of this ticket. We'll need another strategy for test cases.

[kkbrum]

@benthestatistician FYI, this "ignoring" the "within" argument during the actual computation process is causing issues for the new version of student matching for ADSY (since the data is really too large when "within" is ignored). It won't help me (I'll work around it), but I would advocate for actually fixing this issue. At the very least, a note in the "match_on" documentation page about the futility of "within" for "method = rank_mahalanobis" would be nice.

[benthestatistician]

Thanks for taking the trouble to track this down and post a comment on it, @kkbrum. I think I resorted to a shortcut because I wasn't confident I'd be able to update @nullsatz 's code in src/r_smahal.* and src/smahal.*.

Now that we've got Josh back on the team, though: Josh, could you look into enabling the index argument of the rank mahanalobis computation, so that when a BISM is requested it only has to compute and store the within-block distances? My workaround was to use the existing rank mahalanobis routines to compute a distance for each treatment control pair. In sufficiently large matching problems this is no good.

[nullsatz]

I was working on making the smallest change possible.

Simply replacing the C code version,
sqrt(r_smahal(NULL, data, z)),
with the Rcpp version that does support an index argument,
compute_mahalanobis(index, apply(data, 2, rank), z),
and I don't know why.

Failure (test.rank.mahal.R:74:5): compute_rank.mahal returns results similar to the Rosenbaum code
as.matrix(match_on(z ~ X, method = "rank")) not equivalent to `reference_rankmahal`.
25/25 mismatches (average diff: 0.658)
[1] 4.45 - 3.49 == 0.962
[2] 4.24 - 3.45 == 0.790
[3] 5.57 - 3.54 == 2.034
[4] 3.64 - 3.07 == 0.568
[5] 4.38 - 3.19 == 1.191
[6] 4.34 - 3.78 == 0.561
[7] 4.28 - 3.89 == 0.388
[8] 5.42 - 3.93 == 1.481
[9] 3.85 - 3.68 == 0.178
...