day18-Nonbipartite.Rmd

---
title: |
 | Matching with more than two treatments: 
 | Non-bipartite matching
date: '`r format(Sys.Date(), "%B %d, %Y")`'
author: |
  | ICPSR 2023 Session 1
  | Jake Bowers, Ben Hansen, Tom Leavitt
bibliography:
 - 'BIB/MasterBibliography.bib'
fontsize: 10pt
geometry: margin=1in
graphics: yes
biblio-style: authoryear-comp
biblatexoptions:
  - natbib=true
output:
  beamer_presentation:
    slide_level: 2
    keep_tex: true
    latex_engine: xelatex
    citation_package: biblatex
    template: icpsr.beamer
    incremental: true
    includes:
        in_header:
           - defs-all.sty
    md_extensions: +raw_attribute-tex_math_single_backslash+autolink_bare_uris+ascii_identifiers+tex_math_dollars
    pandoc_args: [ "--csl", "chicago-author-date.csl" ]
---


<!-- To show notes  -->
<!-- https://stackoverflow.com/questions/44906264/add-speaker-notes-to-beamer-presentations-using-rmarkdown -->

```{r setup1_env, echo=FALSE, include=FALSE}
library(here)
source(here::here("rmd_setup.R"))
```

```{r setup2_loadlibs, echo=FALSE, include=FALSE}
## Load all of the libraries that we will use when we compile this file
## We are using the renv system. So these will all be loaded from a local library directory
library(tidyverse)
library(dplyr)
library(ggplot2)
library(coin)
library(RItools)
library(optmatch)
library(estimatr)
library(sensitivitymw)
library(sensitivitymult)
library(sensitivityfull)
library(senstrat)
library(rbounds)
library(nbpMatching)
library(lme4)
library(rstanarm)
```

## Today

 1. Agenda: Non-bipartite matching: How to created stratified comparisons if we
    have more than two groups to compare? Multiple treatments, continuous
    treatments.
 2. We are now in the "Non-bipartite Matching" day of the course (see the
    reading there).
 3. Questions arising from the reading or assignments or life?

# But first, review

## Due Diligence and Stratified Observational Designs

 - **Before looking at outcomes** we explain our designs to ourselves by
   comparing the design to our background subtantive understanding of the
   context for causality. (What are the drivers of the "treatment"? How
   **much** adjustment in substantive terms is required? What are the most
   compelling alternative explanations for the treatment$\rightarrow$outcome
   relationship? (Alternative to the explanation that we are
   exploring/assessing))
 - **Before looking at outcomes** we explain our designs to ourselves by
   comparing the design to an equivalently design randomized experiment.
 - **After estimating effects/testing hypotheses** we again engage with
   alternative explanations by modeling how *unobserved covariates* might
   confound the relationship (Sensivity Analysis).


## Sensitivity analysis as a formalized thought experiment {.fragile}


*All non-randomized studies have some bias from unobserved confounders.* The
question is not *whether bias* but about *how much bias* or even *what kinds of bias* would change substantive conclusions.

\begin{center}
\begin{tikzcd}[column sep=large,every arrow/.append style=-latex]
 u \arrow[from=1-1,to=2-2, bend right, "\Gamma" near start] \arrow[from=1-1,to=2-3, bend left, "\rho" near start] \\
& Z  \arrow[from=2-2,to=2-3, "\tau"] &  y \\
x_1 \arrow[from=3-1,to=2-2, "\beta_1" ] \arrow[from=3-1,to=2-3,grey] &
x_2 \arrow[from=3-2,to=2-2, "\beta_2" ] \arrow[from=3-2,to=2-3,grey]
\end{tikzcd}
\end{center}

What is the worst case $u=\{u_1, \ldots, u_N \}$ for tests of the sharp null of no effects?

**Other approaches:** calibrating the model of unobserved confounding using
*observed* covariates (@hosman2010) ("What if the unobserved confounder related
to $Z$ and $Y$ the way that one of our observed $x$ variable does? (What if it
is like education? Or income?)"); placebo tests ("Does our treatment predict
irrelevant differences within strata?")


## Confounding depends on both $\Gamma$ and $\rho$


  @rosenbaum2009amplification propose a way decompose $\Gamma$ into two
  pieces: one $\Delta$ gauges the relationship between an unobserved
  confounder at the outcome (it records the maximum effect of the unobserved
  confounder on the odds of a positive response (imagining a binary outcome))
  and the other $\Lambda$ gauges the maximum relationship between the unobserved
  confounder and treatment assignment.

```{r amplify}
library(sensitivitymv)
lambdas1 <- seq(2 + .01, 2 * 2, length = 100)
ampres1 <- amplify(2, lambda = lambdas1)
lambdas2 <- seq(3 + .01, 2 * 3, length = 100)
ampres2 <- amplify(3, lambda = lambdas2)
```
\centering
```{r echo=FALSE, out.width=".6\\textwidth"}
par(mar = c(3, 3, 1, 1), mgp = c(1.5, .5, 0))
plot(as.numeric(names(ampres1)), ampres1,
  xlab = "Lambda (maximum selection effect of confounder)",
  ylab = "Delta/rho (maximum outcome effect of confounder)",
  main = "Decomposition of Gamma=2 and 3",
  ylim=range(c(ampres1,ampres2)),
  xlim=c(1.9,4),
  type="b"
)
## Add a line for another Gamma
lines(as.numeric(names(ampres2)), ampres2, type = "b",pch=23)
```



# Non-bipartite Matching: An Application with the Study of Race and Place

## How do perceptions of place influence attitudes?

@wong2012jop set out to measure perceptions of environments using an
internet survey of Canadians during 2012 where each respondent drew a
map of their "local community" and then reported their understanding of the
demographic breakdown of this place.

```{r echo=FALSE, results='hide'}
## White English Speaking Canadians only
load(url("http://jakebowers.org/ICPSR/canadamapdat.rda"))
## summary(canadamapdat)
```
\centering
\igrphx{TwoMapsToronto.png}


## Capturing perceptions

Here are 50 maps drawn by people based in Toronto.

\centering
\igrphx{TorontoAllCommunities1.png}

## Capturing perceptions

And here is the question people were asked (groups in random order).

\centering
\igrphx{MLCCPerceptionsQuestion.pdf}

## Capturing perceptions

White, Engish-speaking, Canadian respondents' reports about "visible minorities" in their hand drawn "local communities".

\centering
```{r echo=FALSE}
par(mfrow = c(1, 2))
with(canadamapdat, scatter.smooth(vm.da, vm.community.norm2,
  col = "gray",
  ylab = "Perceptions", xlab = "Census Neighborhood (DA)",
  xlim = c(0, 1), ylim = c(0, 1), lpars = list(lwd = 2)
))
with(canadamapdat, scatter.smooth(vm.csd, vm.community.norm2,
  col = "gray",
  ylab = "Perceptions", xlab = "Census Municipality (CSD)",
  xlim = c(0, 1), ylim = c(0, 1), lpars = list(lwd = 2)
))
## summary(canadamapdat$vm.community.norm2)
```

\note{
Codebook: Mainly for Rmd file

The variables are: age in years, income as a scale, sex in categories, a
social.capital scale coded to run 0 to 1, country of ancestry in categories,
csd.pop is population of the Census Subdivision (like a municipality), vm.csd
is 2006 proportion visible minority in the CSD, vm.da is proportion visible
minority in the Census Dissemination Area (a small area containing 400--700
persons), and vm.community.norm2 is the proportion of visible minorities
reported by respondents in their map of their local community,
community_area_km is the area within their drawing in square km.
}

## How to make the case for perceptions?

If we could randomly assign different perceptions to people, we could claim
that differences of perceptions matter (above and beyond and independent of
objective characteristics of the context).

\medskip

What is an observational design that would do this? Match people on objective
context (and maybe covariates) who differ in perceptions.

\medskip

But objective context is continuous not binary: rather than matching $m$ "treated"
to $n-m$ "controls", we want to compare all $n$ with all $n$ respondents.

```{r echo=FALSE}
## Exclude people who did not offer a perception or an outcome
wrkdat <- canadamapdat[!is.na(canadamapdat$vm.community.norm2) &
  !is.na(canadamapdat$social.capital01), ]
wrkdat$vmdaPct <- wrkdat$vm.da * 100 ## express in pct
```

## Create $n \times n$ distance matrices

Our main design compares white, English-speaking, Canadians with similar
neighborhood proportions of visible minorities (as measured by the Canadian Census in 2006).

```{r echo=TRUE}
scalar.dist <- function(v) {
  ## Utility function to make n x n abs dist matrices
  outer(v, v, FUN = function(x, y) {
    abs(x - y)
  })
}

vmdaDist <- scalar.dist(wrkdat$vmdaPct)
dimnames(vmdaDist) <- list(row.names(wrkdat), row.names(wrkdat))
## The nbpmatching way (Mahalanobis \equiv standardized in one dimension) takes a while:
## obj.com.dist.mat2<-distancematrix(gendistance(wrkdat[,"vmdaPct",drop=FALSE]))
## compare to tmp<-scalar.dist(wrkdat$vmdaPct/sd(wrkdat$vmdaPct))
wrkdat$vmdaPct[1:4]
diff(wrkdat$vmdaPct[1:4])
vmdaDist[1:4, 1:4]
```

## Non-bipartite match

```{r nbp1, echo=TRUE, cache=TRUE}
vmdaDistMat <- distancematrix(vmdaDist)
nbp1match <- nonbimatch(vmdaDistMat)
nbp1 <- get.sets(nbp1match$matches, remove.unpaired = TRUE)
wrkdat[names(nbp1), "nbp1"] <- nbp1
nbp1[1:5]
table(is.na(wrkdat$nbp1)) ## recall the "ghost message"
```

## Inspect the solution

```{r nbpsol, echo=TRUE }
wrkdat[order(wrkdat$nbp1), c("nbp1", "vmdaPct", "vm.community.norm2")][1:6, ]
## table(wrkdat$nbp1)
nbp1vmdiffs <- tapply(wrkdat$vmdaPct, wrkdat$nbp1, function(x) {
  abs(diff(x))
})
nbp1percdiffs <- tapply(wrkdat$vm.community.norm2, wrkdat$nbp1, function(x) {
  abs(diff(x))
})
summary(nbp1vmdiffs)
summary(nbp1percdiffs)
```

```{r echo=FALSE, warning=FALSE, message=FALSE}
source(url("http://jakebowers.org/Matching/nonbimatchingfunctions.R"))
```

## Inspect the solution

\centering
```{r out.width=".8\\textwidth"}
nbmplot(wrkdat,
  yvar = "vmdaPct", xvar = "vm.community.norm2", strata = "nbp1", points = FALSE,
  ylim = range(wrkdat$vmdaPct)
)
```

## Assess balance

No treatment and control groups to compare. But we can still compare the **relationships** between the adjusted variable (`vmdaPct`) and other covariates conditional on pair.

```{r balnbp1, cache=TRUE }
thecovs <- c(
  "age", "income.coded", "education", "x.years", "sex",
  "csd.pop", "vm.csd", "community_area_km"
)
balfmla <- reformulate(thecovs, response = "vmdaPct")
xb1 <- xBalance(balfmla, strata = list(unstrat = NULL, nbp1 = ~nbp1), report = "all", data = wrkdat)
xb1$overall
xb1$results[, c("z", "p"), "nbp1"]

```

## Improve balance using penalties and dropping observations

  For example, we might want to:
  
   -  require matches within Province,
   -  avoid comparing people in small towns to people in large towns,
   -  avoid comparing people who drew big maps to people who drew small maps,
   -  drop the 8 least well matched observations. (Choosing 8 arbitrarily to
      demonstrate.)

```{r echo=FALSE}
rescale01 <- function(x) {
  (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
}
```

```{r nbpdists, echo=TRUE, results="hide"}
csdpopDist <- scalar.dist(wrkdat$csd.pop)
dimnames(csdpopDist) <- list(row.names(wrkdat), row.names(wrkdat))

## Since we have some missing values on community area, and we would like to
## match people who are both missing, we will give it a very large value.
wrkdat$commarea <- ifelse(is.na(wrkdat$community_area_km),
  max(wrkdat$community_area_km, na.rm = TRUE) * 10,
  wrkdat$community_area_km
)

areaDist <- scalar.dist(log(wrkdat$commarea))
dimnames(areaDist) <- list(row.names(wrkdat), row.names(wrkdat))

csdpopDist01 <- rescale01(csdpopDist)
areaDist01 <- rescale01(areaDist)

summary(as.vector(csdpopDist01))
summary(as.vector(areaDist01))
summary(as.vector(vmdaDist))

maxvmdaDist <- max(as.vector(vmdaDist))
```

## Improve balance using penalties and dropping observations

```{r cache=TRUE}
vmdaPen1 <- vmdaDist + (maxvmdaDist * csdpopDist01) + (maxvmdaDist * areaDist01)
vmdaDist[1:5, 1:5]
csdpopDist01[1:5, 1:5]
areaDist01[1:5, 1:5]
vmdaPen1[1:5, 1:5]
```

## Dropping some observations

```{r cache=TRUE}
## now decide how many to drop
vmdaPhPen <- make.phantoms(vmdaPen1, 8, maxval = max(as.vector(vmdaPen1)) * 10)
```

## Improve balance using penalties and dropping observations

```{r echo=TRUE,cache=TRUE}
vmdaPhPenMat <- distancematrix(vmdaPhPen)
nbp2match <- nonbimatch(vmdaPhPenMat)
nbp2 <- get.sets(nbp2match$matches, remove.unpaired = TRUE)
wrkdat[names(nbp2), "nbp2"] <- nbp2
nbp2vmdiffs <- tapply(wrkdat$vmdaPct, wrkdat$nbp2, function(x) {
  abs(diff(x))
})
```

## Assess this new match

Is this match better or worse (in terms of balance? in terms of within-set distances?)

```{r echo=FALSE,results="hide"}
xb2 <- xBalance(balfmla,
  strata = list(unstrat = NULL, nbp1 = ~nbp1, nbp2 = ~nbp2),
  report = "all", data = wrkdat
)
xb2$overall[2:3, ]
xb2$results[, "p", c("nbp1", "nbp2")]
```


## Assess this new match

```{r echo=FALSE,out.width=".8\\textwidth"}
nbmplot(wrkdat, yvar = "vmdaPct", xvar = "vm.community.norm2", strata = "nbp2", points = FALSE, ylim = range(wrkdat$vmdaPct))
```

## Strength of the treatment

The difference in "treatment" within sets varies --- and so we expect the size
of the effect to vary.

```{r treatmentstr, echo=TRUE}
summary(nbp1vmdiffs)
summary(nbp1percdiffs)
summary(nbp2vmdiffs)
```

## Strength of the treatment

To prevent many sets with no variation on "treatment" we could add a penalty
for pairs that are too close on treatment: to maximize size of treatment (i.e.
differences in perceptions) while minimizing differences on covariates. That
is, we should be able to get more precision / power about the perceptions
related differences if they are larger.

For example: if difference in  perceptions is  small,  make the distance very large.


```{r penaltytostrengthentreatment, echo=TRUE}
## Just showing the first 5 respondents as an example
percDist <- scalar.dist(wrkdat$vm.community.norm2 * 100)
da <- vmdaDist[1:5, 1:5]
perc <- percDist[1:5, 1:5]
da + 1000 * (perc < 2)
```


## Try designmatch to balance covs and strengthen treatment

The `designmatch` package also has the ability to produce non-bipartite stratifications. It has some nice features that allow for strenghening the effects of the intervention, too.

```{r designmatch, cache=TRUE,eval=TRUE}
library(designmatch)
library(gurobi)

wrkdat_new <- na.omit(wrkdat[,c("csd.pop","community_area_km","social.capital01","vm.community.norm2","vmdaPct")])

thecovs <- c("csd.pop","community_area_km")
mom_covs <- wrkdat_new[, thecovs]
## calculated absolute std mean diffs between trt and ctrl on covs
## differences should be at most .5 sds apart
mom_tols <- round(apply(mom_covs, 2, sd,na.rm=TRUE)*.5, 2)
momlist <- list(covs = mom_covs, tols = mom_tols)

vmdaDist_new <- scalar.dist(wrkdat_new$vmdaPct)
dimnames(vmdaDist_new) <- list(row.names(wrkdat_new), row.names(wrkdat_new))
vmDist_integers <- round(vmdaDist_new*100,1)

farlist <- list(covs=as.matrix(round(wrkdat_new[,"vm.community.norm2",drop=FALSE]*100,1)),far.pairs=c("vm.community.norm2"=1))

#solverlist <- list(name = "glpk", approximate = 1)
solverlist <- list(name = "gurobi", approximate = 1, trace=1)
resnb <- nmatch(dist_mat= vmDist_integers,
    solver = solverlist,
    mom=momlist,
    far=farlist
)

str(resnb)

```



## Assess hypotheses about effects

Test the hypothesis of no relationship between perceptions
`vm.community.norm2` and `social capital`.

```{r eval=TRUE,echo=TRUE}
library(coin)
## These are the same test in this case
test1 <- independence_test(social.capital01 ~ vm.community.norm2 | nbp1, data = wrkdat[!is.na(wrkdat$nbp1), ])
test2 <- xBalance(vm.community.norm2 ~ social.capital01,
  strata = list(nbp1 = ~nbp1), data = wrkdat, report = c("adj.mean.diffs", "chisquare.test", "z.scores", "p.values")
)
test1
test2$overall
```


## Describe the differences within pairs

Does the person who perceives more visible minorities in their community tend
to be higher (or lower) in `social.capital` than the other person in the pair?

```{r echo=FALSE}
rank.pairs <- function(x, block) { ## Identify the low and high subj in each pair
  unsplit(lapply(split(x, block), function(x) {
    rank(x)
  }), block)
}
```

```{r}
wrkdat$scRank <- with(wrkdat, rank.pairs(social.capital01, nbp1))
wrkdat$vmCRank <- with(wrkdat, rank.pairs(vm.community.norm2, nbp1))
wrkdat[order(wrkdat$nbp1), c("nbp1", "social.capital01", "scRank", "vm.community.norm2", "vmCRank")][1:6, ]
with(wrkdat, tapply(scRank, vmCRank, mean))
```

## Summarize mean differences within pairs

If perceptions matters for social capital then we would expect pairs differing
greatly in subjective context to display greater differences in social capital
than pairs that differ a little.


```{r echo=FALSE,results="hide"}
align.by.block <- function(x, block, fn = mean, thenames = NULL) { ## By default, this rescales each observation to be the distance from the group mean.
  newx <- unsplit(lapply(split(x, block), function(x) {
    x - fn(x)
  }), block)
  if (!is.null(names)) {
    names(newx) <- thenames
  }
  return(newx)
}
```

```{r}
wrkdat$scMD <- with(wrkdat, align.by.block(social.capital01, nbp2))
wrkdat$vmcn2MD <- with(wrkdat, align.by.block(vm.community.norm2, nbp2))
wrkdat[order(wrkdat$nbp1), c("social.capital01", "scMD", "vm.community.norm2", "vmcn2MD", "nbp1")][1:4, ]
## notice that aligning or pair-mean-centering the data preserves the within
## set relationships
## summary(tapply(wrkdat$scMD,wrkdat$nbp1,function(x){ abs(diff(x)) }))
## summary(tapply(wrkdat$social.capital01,wrkdat$nbp1,function(x){ abs(diff(x)) }))
lm1 <- lm_robust(scMD ~ vmcn2MD, data = wrkdat[!is.na(wrkdat$nbp2), ])
lm1
```
## Summarize mean differences within pairs

```{r warning=FALSE,cache=TRUE}
lm2 <- lm_robust(scMD ~ vmcn2MD, data = wrkdat[!is.na(wrkdat$nbp2), ])
lm2
lm3 <- lm_robust(social.capital01 ~ vm.community.norm2, fixed_effects = ~nbp2, data = wrkdat, subset=!is.na(wrkdat$nbp2))
lm3
table(wrkdat$vmCRank, exclude = c())
lm4 <- lm_robust(social.capital01 ~ I(vmCRank - 1), fixed_effects = ~nbp2, data = wrkdat,subset=!is.na(wrkdat$nbp2))
lm4
```

## Summarize mean differences within pairs

If perceptions matter for social capital above and beyond objective context
then we would expect pairs differing greatly in subjective context to display
greater differences in social capital than pairs that differ a little.

```{r}
lm1
lm2
lm3
pairdiffs <- wrkdat %>%
  filter(!is.na(vmCRank) & !is.na(social.capital01) & !is.na(nbp2)) %>%
  group_by(vmCRank) %>%
  summarize(mnsc = mean(social.capital01))
wrkdat[order(wrkdat$nbp2), c("social.capital01", "scRank", "scMD", "vm.community.norm2", "vmcn2MD", "vmCRank", "nbp2")][1:4, ]
lm4
```

## Summarize mean differences within pairs

```{r}
summary(wrkdat$vmcn2MD)
summary(wrkdat$scMD)
```

Within matched pair, the person who perceives more visible minorities within set tends to report
lower social capital than the person who perceives fewer visible minorities
within set.

\medskip

The largest difference is about `r round(max(wrkdat$vmcn2MD,na.rm=TRUE),2)`. The model
predicts that social capital would differ by about `r coef(lm1)[[2]]*.48` for such a difference. This is about
`r coef(lm1)[[2]]*.48/sd(wrkdat$scMD,na.rm=TRUE)` of a standard deviation
of the social capital scale. Or about
`r coef(lm1)[[2]]*.48/abs(diff(range(wrkdat$scMD,na.rm=TRUE)))` of the range.


## Summarize mean differences within pairs

Here is a look at the within-pair differences in perceptions of visible minorities as well as social capital.

```{r smoothplot, out.width=".7\\textwidth", echo=FALSE}
with(wrkdat, scatter.smooth(vmcn2MD, scMD, span = .3, cex = .7, col = "gray", pch = 19, lpars = list(lwd = 2)))
abline(h = 0, lwd = .5)
```


## Summary of matching without groups

 - Workflow in general is the same as matching with groups (covariates,
   distance matrices, optimization to select a stratification, assessment of
   the stratification by comparison to an experiment)
 - Estimation is more flexible --- could look simply at "higher versus lower"
   within  pair, or could average over scores.



## Another estimation approach

\autocite{smith:1997} presents a multi-level modelling approach to taking
matched sets into account. The weights implied here are a bit different from
the weights that we've discussed before (although with pairs they might be more
or less the same). What is the data model? What additional assumptions are involved
here?

```{r lmer, cache=TRUE, message=FALSE, warning=FALSE}
library(lme4)
wrkdat$vmCbi <- ave(wrkdat$vm.community.norm2, wrkdat$nbp1)
lmer1 <- lmer(social.capital01 ~ vm.community.norm2 + vmCbi + (1 | nbp1), data = wrkdat)
confint(lmer1)["vm.community.norm2", ]
## Notice that we may have less information than we thought (some clustered of observation by DA).
table(table(wrkdat$dauid))
## So maybe:
lmer2 <- lmer(social.capital01 ~ vm.community.norm2 + vmCbi + (1 | nbp1) + (1 | dauid), data = wrkdat)
confint(lmer2)["vm.community.norm2", ]
```

## Other applications of non-bipartite matching?

See: DOS Chapter 11.

Also: this has a lot of applications in experimental design (see `blockTools` and \autocite{moore2012blocktools,moore2012multivariate}).



## Summary and Questions:

What do you think? What questions arise?



## References