within= argument for scoreCaliper()

[benthestatistician]

ee_s_c needs scoreCaliper() to gain an optional within= argument, accepting a InfinitySparseMatrix and incorporating into the InfinitySparseMatrix that it returns exclusions of whatever pairings were excluded by the within ISM.

If the within ISM is a BISM, then scoreCaliper() should use the grouping structure encoded within it to improve performance in large problems. In this case its return value should be a BISM as well.

[benthestatistician]

In adsy_mistem_school PR27 @kkbrum reports

I had to modify the match_on function in optmatch because it wasn't really using "within" when doing a "scoreCaliper", which resulted in a long runtime of "makeModelPSCal" for ~150,000 students.

The modification is that her block of match_on.numeric() corresponding to exclude=NULL ceases to be

 } else {
        allowed <- scoreCaliper(x, z, caliper)
      }

and instead becomes

 } else {
      # This did not use to take "within" into account; instead "within" gets added later
      # However, if the problem is very very large, this is very slow, and "within" needs to be considered
      if (!is.null(within)) {
        allowed_list <- list()
        for (blk in levels(within@groups)) {
          ix <- within@groups == blk
          if (sum(ix & z) > 0 & sum(ix & !z) > 0) {
            allowed_list[[blk]] <- optmatch:::scoreCaliper(x[ix], z[ix], caliper)
          }
        }
        allowed <- dbind(allowed_list) 
        allowed <- meerktjob:::addRows(allowed, within@rownames[!within@rownames %in% allowed@rownames], Inf)
        allowed <- t(meerktjob:::addRows(t(allowed), within@colnames[!within@colnames %in% allowed@colnames], Inf))
      } else {
        allowed <- scoreCaliper(x, z, caliper)
      }
    }

(see her file 99-match_on.R.) We should update optmatch:::scoreCaliper() so that the work of this block (at least its performance aspects) is instead done inside of there.

Other notes:

1. We may need to require that a within= passed to scoreCaliper() be a BISM, not an ISM or dense matrix. (Or perhaps we would allow ISM arguments, but only offer performance advantage in the case of ISMs.) Interested to hear @markmfredrickson's thoughts on this.
2. @nullsatz I'm likely to request your assistance with this.
3. This issue is similar in spirit to match_on's matrix & ISM methods ignore within arguments #181, with tests at test.match_on.R#L401-L412.

[nullsatz]

I'm working on this in the branch within-score-caliper

[benthestatistician]

Perhaps the addRows() function being imported from ee_s_c could be instead be called addEligibleTreatments(). Then there could be a partly analogous function for adding "columns" labeled addIneligibleControls(). (It could have val=Inf, as opposed to val=0 in the current addRows().)

[nullsatz]

Perhaps the addRows() function being imported from ee_s_c could be instead be called addEligibleTreatments(). Then there could be a partly analogous function for adding "columns" labeled addIneligibleControls(). (It could have val=Inf, as opposed to val=0 in the current addRows().)

Ok, I switched to addEligibleTreatments() and added an addIneligibleControls(). I think the next thing to do is to write a couple of tests for these functions.

[nullsatz]

Perhaps the addRows() function being imported from ee_s_c could be instead be called addEligibleTreatments(). Then there could be a partly analogous function for adding "columns" labeled addIneligibleControls(). (It could have val=Inf, as opposed to val=0 in the current addRows().)

Ok, I switched to addEligibleTreatments() and added an addEligibleTreatments(). I think the next thing to do is to write a couple of tests for these functions.

addEligibleTreatments

test

addIneligibleControls

test

[benthestatistician]

These helpers look great, Josh. How about giving these functions an @keywords internal line in the roxygen block? When they become a part of optmatch, ee_s_c should stop using addRows() and addCols() in favor of these.

Next step is to use them to enable scoreCaliper(...,within=...)

[benthestatistician]

How about using addIneligibleControls() instead of t() and addEligibleTreatments() at match_on.R#L773-L775:

allowed <- t(addEligibleTreatments(t(allowed),
                    within@colnames[!within@colnames %in% allowed@colnames],
                    Inf))

(Seems like it'd be faster, but of course I haven't checked.)

[nullsatz]

(Seems like it'd be faster, but of course I haven't checked.)

I used microbenchmark to try to see if our changes are affecting the running time at all.

my test

library(optmatch)
library(microbenchmark)

n <- 1600
Z <- rep(c(0,1), n/2)
X1 <- rnorm(n, mean = 5)
X2 <- rnorm(n, mean = -2, sd = 2)
B <- as.factor(rep(c(0,1), each = n/2))
stratify <- exactMatch(Z ~ B)
res.strat <- match_on(Z ~ X1 + X2, within = stratify)
microbenchmark(match_on(Z ~ X1 + X2, within = stratify))

output

	min	lq	mean	median	uq	max	neval
master	72.9891	74.70754	80.75358	75.8132	79.37901	105.3166	100
t()	72.63367	73.99143	80.07324	75.41831	78.98814	103.3892	100
aIC	74.13087	75.63893	81.31143	76.54536	78.58816	106.8389	100

It seems like addIneligibleControls is a bit slower. It could probably be tinkered with to speed it up.

[benthestatistician]

Aha. Nothing like actual data!

If you see a promising shortcut to improve performance for addIneligibleControls(), do go ahead with it. But if you don't, or if the only promising shortcuts you see would be time consuming to code, then I'm happy to have us use the version with t(addEligibleTreatments()).

[nullsatz]

(Seems like it'd be faster, but of course I haven't checked.)

I used microbenchmark to try to see if our changes are affecting the running time at all.

I had a hunch about my script and decided to test it.

The script I was using wasn't actually calling addIneligibleControls. I put a print statement in there to verify. I found a test testthat test that does exercise that part of the code and increased the data size significantly. The function addIneligibleControls does perform significantly better than the previous t(addEligibleControls... method.

I stuck some calls to nanotime inside scoreCaliper and got a runtime of
00:00:00.000_657_081 for addIneligibleControls and a runtime of 00:00:00.013_576_985 for the previous t(addEligibleControls... method.

I used the script below.

library(optmatch)
library(testthat)

n <- 1000
n2 <- n * 2
n4 <- n * 4

scores <- rep(c(1,2,2,3), each = n)
z <- rep(c(0,1), each = n2)

names(scores) <- paste("A", 1:n4, sep = "")

## get the caliper only results
res.cal <- match_on(scores, z = z, caliper = 1)
expect_equal(max(res.cal), 1)

## now make up a within arg
names(z) <- names(scores)
www  <- exactMatch(x = as.factor(rep(c(0,1), n2)),
                   treatment = z)

res.w  <- match_on(scores, z = z, within = www)
expect_true(max(res.w) > 1)

## they should safely interact
res.cal.w <- match_on(scores, z = z, within = www, caliper = 1)
expect_equal(max(res.cal.w), 1)

[nullsatz]

I copied a small test from the other testthat tests to see whether I get a BISM returned from scoreCaliper when I pass scoreCaliper a BISM via the within argument. See below.

  scores <- rep(1:3, each = 4)
  z <- rep(c(0,1), 6)
  names(z) <- names(scores) <- letters[1:12]
  b <- rep(1:3, 4)

  ez <- exactMatch(z ~ b)
  a <- scoreCaliper(scores, z, caliper=1, within=ez)
  expect_true(is(a, "BlockedInfinitySparseMatrix")) # passes
  expect_equal(ez@groups, a@groups) # fails

I do get a BISM back. The BISM does have a groups member with the same names as the within argument, but the values are different. See below.

  ## currently getting:
  ## > ez@groups
  ## a b c d e f g h i j k l
  ## 1 2 3 1 2 3 1 2 3 1 2 3
  ## Levels: 1 2 3
  ## > a@groups
  ## a g d j e k b h c i f l
  ## 0 0 0 0 1 1 1 1 2 2 2 2
  ## Levels: 0 1 2

It looks like levels have been mapped 1 => 0, 2 => 1, 2 => 3 while the names have remained the same. Not sure this is a problem.