day11_Nonbipartite.Rmd

---
title: Non-bipartite Matching
date: '`r format(Sys.Date(), "%B %d, %Y")`'
author:  |
  | ICPSR 2023 Session 1
  | Jake Bowers \& Tom Leavitt
bibliography:
 - 'BIB/MasterBibliography.bib'
fontsize: 10pt
geometry: margin=1in
graphics: yes
biblio-style: authoryear-comp
colorlinks: true
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"))
opts_chunk$set(echo = TRUE, digits = 4)
```

```{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(nbpMatching)
library(lme4)
library(rstanarm)
library(slam)
```

## Today

 1. Agenda: Non-bipartite matching: How to created stratified comparisons if we
    have more than two groups to compare? Multiple treatments, continuous
    treatments. See Rabb et al (2022) as a published example.
  3. Questions arising from the reading or assignments or life?

# But first, review

## What is regression doing?

Following Berk (2004)'s *Regression Analysis: A Constructive Critique*.

 - Describing relationships: If $Z \in {0,1}$ then `lm(Y~Z)` gives you `mean(Y[Z==1]) - mean(Y[Z==0])`; If $Z$ has more values, then this describes a linear relationship. Very useful. Simple interpretation.
 - Statistical Inference
 - Causal Inference

## 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 theoretical explanation that we are
   exploring/assessing))
 - **Before looking at outcomes** we explain our designs to ourselves by
   comparing the design to an equivalently designed randomized experiment using
   the known distribution of the $d^2$ statistic under the null hypothesis of
   no covariate-to-treatment relationships across any covariates (see the
   Hansen and Bowers 2008 piece).
 - We estimate (average) effects and test hypotheses about effects **as if the
   research design was randomized**.
 - **After estimating effects/testing hypotheses** we again engage with
   alternative explanations by modeling how *unobserved covariates* might
   confound the relationship (Sensivity Analysis).

## But first review

 - Statistical inference for Causal Effects and Causal Hypotheses in Randomized
   Experiments
 - Adjustment by stratification;
   - Matching to generate optimal stratifications  (decisions and strategies
     that are part of research design; matching on missingness and `fill.NAs`;
     `exactMatch`; `caliper`; `min.controls`; `effectiveSampleSize`);
   - Assessing success of stratified research designs in adjustment;
   - The As-If-Randomized mode of statistical inference for stratified research
     designs (treat a matched design as a block-randomized experiment).

## Review 2: An Adjustment Strategy to Address Alternative Explanations Effectively

How to strengthen evidence about the claim that Metrocable caused a decrease in crime?

 1. **List the main alternative explanations** (could crime have caused
    Metrocable stations; socio-economic status differences; \ldots). Can we
    operationalize these explanations?
 2. **Stratify data to minimize heterogeneity within set.** If education does
    not vary within set, then we have "adjusted for" education by conditioning
    on the set. The `optmatch` package for R finds sets that minimize the
    weighted sum of distances across the sets. (See also `rcbalance`, `DiPs`,
    `bigmatch`, `designmatch`, `quickmatch`).
    1. Create distance matrices using `match_on` (and `caliper` and
       `exactMatch`) (Scalar distances on  especially important variables like
       baseline outcomes; Multivariate distances in terms of other covariates
       via Mahalanobis or Propensity scores distances.)
    2. Find stratifications using `fullmatch` etc (`bmatch` from `designmatch`,
       etc.).



```{r echo=FALSE, cache=TRUE}
load(url("http://jakebowers.org/Data/meddat.rda"))
meddat<- mutate(meddat,
		HomRate03=(HomCount2003/Pop2003)*1000,
		HomRate08=(HomCount2008/Pop2008)*1000)
row.names(meddat) <- meddat$nh
```


```{r echo=FALSE}
covs <- unique(c(names(meddat)[c(5:7,9:24)],"HomRate03"))
balfmla <- reformulate(covs,response="nhTrt")
mhdist <- match_on(balfmla,data=meddat, method="rank_mahalanobis")
psmod <- arm::bayesglm(balfmla,data=meddat,family=binomial(link="logit"))
stopifnot(any(abs(coef(psmod))<10))
psdist <- match_on(psmod,data=meddat)
tmp <- meddat$HomRate03
names(tmp) <- rownames(meddat)
absdist <- match_on(tmp, z = meddat$nhTrt,data=meddat)
```

## Example design  and workflow {.allowframebreaks}

```{r echo=TRUE}
## Inspect the distance matrices
quantile(as.vector(psdist),seq(.9,1,.01))
quantile(as.vector(mhdist),seq(.9,1,.01))
quantile(as.vector(absdist),seq(.9,1,.01))
matchdist <- psdist + caliper(psdist,9) + caliper(absdist,2) + caliper(mhdist,60)
summary(matchdist)
fm1 <- fullmatch(matchdist, min.controls= 1, max.controls=Inf, data=meddat,tol=.00001)
summary(fm1,min.controls=0,max.controls=Inf,propensity.model=psmod)
meddat$fm1 <- factor(fm1)
meddat$nhTrtF <- factor(meddat$nhTrt)
```


## An Adjustment Strategy to Address Alternative Explanations

 3. **Assess the stratification in substantive terms** If we look within the
    sets, are the differences we see substantively concerning or trivial?



```{r echo=FALSE}
meddat[names(fm1),"fm1"] <- fm1
setmeanDiffs <- meddat %>% filter(!is.na(fm1)) %>% group_by(fm1) %>%
  summarise(ateb=mean(HomRate08[nhTrt==1])-mean(HomRate08[nhTrt==0]),
            nb=n(),
            nTb = sum(nhTrt),
            nCb = sum(1-nhTrt),
            prob_trt = mean(nhTrt),
	    baselinediffs = mean(HomRate03[nhTrt==1])-mean(HomRate03[nhTrt==0]),
	    minbaselines = min(HomRate03),
	    maxbaseline = max(HomRate03)
            ) %>% arrange(abs(baselinediffs))
setmeanDiffs <- setmeanDiffs %>% mutate(nbwt = nb/sum(nb), hbwt0 = nbwt * prob_trt * ( 1- prob_trt))
setmeanDiffs$hbwt <- with(setmeanDiffs, hbwt0/sum(hbwt0))

setmeanDiffs %>% dplyr::select(fm1,ateb,nbwt,hbwt)
```

## An Adjustment Strategy to Address Alternative Explanations

 4. **Assess the stratification by comparison to a model of a block-randomized
    experiment** Does our research design look like a block-randomized
    experiment in terms of covariate balance? If so, move onto step 4.
    Otherwise, work to improve the research design by (a) changing scores; (b)
    combining scores (for example, using calipers); (c) excluding units (using
    calipers); (d) exact matching on subgroups; (e) reducing variation in
    set-size.

## An Adjustment Strategy to Address Alternative Explanations

 5. **Estimate effects and test hypothesis as-if-block-randomized** Estimators
    and tests refer to the finite "population" of the study pool and the fixed
    stratification in the same way common in the analysis of block-randomized
    experiments.

```{r  estandtest,  cache=TRUE, warning=FALSE}
meddat_new <-  meddat %>%  filter(!is.na(fm1)) %>% group_by(fm1) %>%
	mutate(trtprob=mean(nhTrt), nbwt=nhTrt/trtprob + (1-nhTrt)/(1-trtprob))
estate <-  lm_robust(HomRate08~nhTrt,data=meddat_new,weights=nbwt,subset=!is.na(fm1))
estate
with(setmeanDiffs,sum(ateb*nb/sum(nb)))
with(setmeanDiffs,sum(ateb*nbwt))

with(setmeanDiffs,sum(ateb*hbwt))

estate_fe <- lm_robust(HomRate08~nhTrt,fixed_effects = ~fm1, data=meddat_new,subset=!is.na(fm1))
estate_fe

estate_fe2 <- lm_robust(HomRate08~nhTrt+fm1, data=meddat_new,subset=!is.na(fm1))
coef(estate_fe2)["nhTrt"]

set.seed(12345)
ranktest1  <-  wilcox_test(HomRate08~nhTrtF|fm1,data=meddat_new,
		       distribution=approximate(nresample=1000))
## Notice that the "99 percent confidence interval" below refers to predicted differences in the p-value across difference simulations with different seeds.
coin::pvalue(ranktest1)
## Doing another 10000
blah<- wilcox_test(HomRate08~nhTrtF|fm1,data=meddat_new,
		       distribution=approximate(nresample=1000))
coin::pvalue(blah)
```

## An Adjustment Strategy to Address Alternative Explanations

 6. **Assess the sensitivity of the analysis to the assumptions of as-if-randomized** The design is not a randomized design. Is this likely to
    cause small or large changes in the substantive interpretation of our
    results? (R packages `sensitivityfull`, `sensitivitymv`, `sensitivitymw`,
    `sensitivitymult`, `rbounds`)

## What about data with multiple observations of each unit?

Panel, longitudinal, time-series cross-sectional, clustered, multilevel,
nested,\ldots

\medskip

\bh{Remember}:**A parametric model is a not a research design.**

 - \bh{If "treatment" occurred only once} ("birth of first child aka transition to
parenthood $\rightarrow$ political activity", "first seatbelt law $\rightarrow$
highway deaths", etc.): What unit(s) (observed at which point in time) is/are a good
counterfactual to the focal unit experiencing the change? How can we find those
units and focus comparisons on those units versus the focal units? (See
Rosenbaum 2010 and 2017 on "risk-set matching".)


 - \bh{If you find it difficult to make the as-if-randomized design} If you don't want to find units which were similar up until the moment of
treatment (say, you have few covariates), can you make other assumptions? (See
the difference-in-differences idea of parallel trends and/or lagged DV and/or
other assumptions in Tom's discussion.)


## What about data with multiple observations of each unit?


 - \bh{If "treatment" occurred more than once:} Are you estimating an effect averaging
over all occurrences? How do you want to weight each occurrence? Equally? What
are you assuming about SUTVA/interference across those treatments within a unit?
(Often these questions are hard to answer so you might prefer to break up the
problem into simpler pieces.)


- See papers on TSCS / panel / longitudinal data referred to here
  <https://imai.fas.harvard.edu/research/FEmatch.html> and
  <https://yiqingxu.org/research/> "A Practical Guide to Counterfactual
  Estimators for Causal Inference with Time-Series Cross-Sectional Data"
  <https://onlinelibrary.wiley.com/doi/10.1111/ajps.12723>.
 - And Ding and Li 2019, Blackwell and Glynn 2018 with a nice pre-post design
   application: Keele, L., Cubbison, W., and White, I. (2021). Suppressing
   black votes: A historical case study of voting restrictions in louisiana.
   American Political Science Review, pages 1–7.


# Non-bipartite Matching: The Medellin Data


## Hypothetical Setup {.allowframebreaks}

Imagine that there is a debate about whether housing insecurity is strongly
related to violence. We have neighborhoods in Medellin where
we have measured both violence scaled by the population of the place
(`HomRate08`), whether people own their own home (`nhOwn`), and potential confounders like the proportion of people who are employed (`nhEmp`).  However, we know that both housing insecurity as well as violence can be predicted from other background variables: maybe the relationships we would summarize between housing and violence  would be confounded by those other relationships.

## Designmatch setup {.allowframebreaks}

We will use an approach to adjustment called **non-bipartite** matching) which
doesn't require two groups. Rather it creates pairs of units (neighborhoods) in
this case, which are as similar as possible in regards to many covariates.

```{r echo=TRUE}
covs <- c("nhClass", "nhSisben","nhPopD",  "nhQP03",  "nhPV03",  "nhTP03",
    "nhBI03",  "nhCE03",  "nhNB03" , "nhMale",  "nhAgeYoung",
    "nhAgeMid","nhMarDom","nhSepDiv","nhAboveHS" , "nhHS", "HomRate03")

covmat <- dplyr::select(meddat,one_of(covs))

## Mahalanobis distances for each neighborhood
meddat$covmh <- mahalanobis(
  x = covmat ,
  center = slam::col_means(covmat),
  cov = cov(covmat)
)

## Absolute mahalanobis distances between neighborhoods
mhdist_mat <- outer(meddat$covmh, meddat$covmh, FUN = function(x, y){ abs(x - y) })
dimnames(mhdist_mat) <- list(meddat$nh,meddat$nh)
```

## Designmatch use {.allowframebreaks}

Now, we can match on those distances:


```{r echo=TRUE}
## Turns out that the designmatch software doesn't like too many decimals, and prefers
## mean-centered distances. This doesn't really matter in substantive terms but is important in
## regards to getting the software to work
matchdist_mat <- round(100*mhdist_mat / mean(mhdist_mat), 1)

## Restrict allowable matches. This is like a caliper but on two dimensions.
nearlist <- list(covs=as.matrix(meddat[,c("HomRate03","nhAboveHS")]),
		 pairs=c(HomRate03=5,nhAboveHS=.5))

## For larger problems you will want to install gurobi using an academic
## license. After installing the license, then I do something like the following
## where the details of the version numbers will differ
## install.packages("/Library/gurobi952/macos_universal2/R/gurobi_9.5-2_R_4.2.0.tgz",repos=NULL)
## also had to use a different version of designmatch for now:

## Only run this next one one time
### renv::install("bowers-illinois-edu/designmatch")
library(designmatch)
#library(slam)
library(highs)
#library(gurobi)
solverlist <- list(name = "highs", approximate = 0, t_max = 1000, trace = 1)
```

The function `nmatch` does the optimization. It is not full-matching, but is pair-matching. 

```{r}
mh_pairs <- nmatch(
  dist_mat = matchdist_mat,
  near = nearlist,
  subset_weight = 1,
  solver = solverlist
)
## mh_pairs
```

```{r def_fn, echo=FALSE}
#' Function to convert the output of nmatch into a factor variable for use in analysis
nmatch_to_df <- function(obj, origid) {
## We want a factor that we can merge onto our
## existing dataset. Here returning a data.frame so that
## we can merge --- seems less error prone than using
## rownames even if it is slower.
    matchesdat <- data.frame(
        bm = obj$group_id,
        match_id = c(obj$id_1, obj$id_2)
        )
      matchesdat$id <- origid[matchesdat$match_id]
      return(matchesdat)
  }
```

```{r convert_matches_to_factor, echo=TRUE}
mh_pairs_df <- nmatch_to_df(mh_pairs,origid=meddat$nh)
nrow(mh_pairs_df)

## So, in matched set 1 (bm==1) we see two neighborhoods:
mh_pairs_df %>% filter(bm==1)
mh_pairs_df$nh <- mh_pairs_df$id

# The nmatch_to_df function creates a column labeled "bm" which contains
meddat2 <- inner_join(meddat, mh_pairs_df, by = "nh")
meddat2 <- droplevels(meddat2)
stopifnot(nrow(meddat2) == nrow(mh_pairs_df))

## Number of matches:
# meddat2$bm is the matched set indicator.
stopifnot(length(unique(meddat2$bm)) == nrow(meddat2) / 2)
nrow(mh_pairs_df)
nrow(meddat2)
## Notice some observations were not matched:
nrow(meddat)
```

## Assessing the design {.allowframebreaks}

Now, what we are trying to do is break the relationship between covariates and
the main explanatory variables (just as we might in a pair randomized study):
the neighborhood higher on the explanatory variable shouldn't be systematically more or less likely to be the neighborhood higher on any given covariate in such a study. We assess this below:

```{r echo=TRUE}
## Make a new variable that is 1 for the neighborhood higher in home ownership
## and 0 for the neighborhood who is lower. (Similarly for Employment)
## We'd like to show that the covariates are not related to either home
## ownership or employment within pair.
meddat2 <- meddat2 %>%
  group_by(bm) %>%
  mutate(rank_own = rank(nhOwn) - 1,
      rank_emp = rank(nhEmp) - 1) %>%
  arrange(bm) %>%
  ungroup()

## Notice pair bm=1
meddat2 %>% dplyr::select(bm,nh, nhOwn,rank_own,nhEmp, rank_emp)

## Check for sets with a tie
table(meddat2$rank_own)

## Since balanceTest demands binary treatment, we remove ties for now.
meddat3 <- meddat2 %>% filter(rank_own!=.5)
table(meddat3$rank_own)

## We are trying to break the relationships between the covariates and the two
## explanatories. Let's look at one of them here.

## Since we have a smaller dataset, we need to use fewer covariates if we want to use the large sample approximation from balanceTest
newcovs <- c("nhClass","HomRate03","nhTP03","nhAgeYoung","nhAboveHS","nhEmp")

balfmla_new <- reformulate(newcovs, response = "rank_own")
## Using only the matched data and also conditional within sets
xb_own <- balanceTest(update(balfmla_new,.~.+strata(bm)), data = meddat3,p.adjust="none")
xb_own$overall
xb_own_vars <- data.frame(xb_own$results[, c("Control", "Treatment", "adj.diff", "std.diff", "p"), "bm"])
## xb_own_vars$padj <- p.adjust(xb_own_vars$p, method = "holm") ## already adjusted using holm adjustment by default in balanceTest
options(digits = 3)
arrange(xb_own_vars, p) %>% zapsmall(digits = 5)
stopifnot(xb_own$overall[, "p.value"] > .3)
```

An equivalent way to do what balanceTest is doing

```{r echo=TRUE}
library(formula.tools)
library(coin)
coin_fmla <- ~ rank_own | bmF
lhs(coin_fmla) <- rhs(balfmla_new)
meddat3$bmF <- factor(meddat3$bm)
coin_test <- independence_test(coin_fmla,data=meddat3,teststat="quadratic")
coin_test_perm <- independence_test(coin_fmla,data=meddat3,teststat="quadratic",distribution=approximate(nresample=1000))
```

## Outcome Analysis

Now, assuming we are happy with the design, we describe the relationships between home ownership and violence in 2008 at the neighborhood level.

```{r}
## Ways to assess the relationship between home ownership and the outcome
## conditional on sets. These are all the same.

## We will start with estimating the difference between the high and low home
## ownership neighborhoods and then move to estimating the smooth linear
## relationship between differences in proportion home ownership and the
## outcome.

## First, the most transparent way, but most typing is to convert the data
## into the strata level and create averages.
meddat2$bmF <- factor(meddat2$bm)
pair_diffs <- meddat2 %>% filter(rank_own!=.5) %>%
    group_by(bmF) %>%
    summarize(hr=mean(HomRate08),
    hr_diff=HomRate08[rank_own==1] - HomRate08[rank_own==0],
    own_diff=nhOwn[rank_own==1] - nhOwn[rank_own==0],
    own_diff_raw=diff(nhOwn),
    hr_diff_raw=diff(HomRate08),.groups="drop")

## Simply the mean of the differences within pair between the higher and lower
## home ownership neighborhoods. We will see that this is exactly the same as
## the other estimates.
est1 <- mean(pair_diffs$hr_diff)
est1

est2 <- difference_in_means(HomRate08~rank_own,blocks=bm,data=meddat2,subset=rank_own!=.5)
est3 <- lm_robust(HomRate08~rank_own,fixed_effects=~bm,data=meddat2,subset=rank_own!=.5)
est4 <- lm_robust(HomRate08~rank_own+bmF,data=meddat2,subset=rank_own!=.5)
## This next estimate is often called the group-mean centered or mean-deviated version
## it is what is happening the background of the fixed_effects approach
meddat2 <- meddat2 %>% group_by(bmF) %>%
    mutate(hr_md = ifelse(rank_own!=.5,HomRate08- mean(HomRate08),NA),
        rank_own_md = ifelse(rank_own!=.5,rank_own - mean(rank_own),NA))
est5 <- lm_robust(hr_md~rank_own_md,data=meddat2)

rbind(est1=est1,
    est2=coef(est2),
    est3=coef(est3),
    est4=coef(est4)[["rank_own"]],
    est5=coef(est5)[["rank_own_md"]])

all.equal(est1,coef(est4)[["rank_own"]])
all.equal(est1,coef(est2)[["rank_own"]])
all.equal(est1,coef(est3)[["rank_own"]])
all.equal(est1,coef(est5)[["rank_own_md"]])
```


```{r}
## More information about the mean-deviated approach to adjusting for pairs
meddat2 %>% dplyr::select(bmF,nhOwn,rank_own, rank_own_md, HomRate08, hr_md ) %>% head()
meddat2 %>% ungroup() %>% filter(rank_own!=.5) %>% summarize(mean(rank_own_md),mean(hr_md))
```

```{r}
## Notice exactly the same as the mean outcome within each pair
group_means <- lm_robust(HomRate08~bmF,data=meddat2,subset=rank_own!=.5)
coef(group_means)
rbind(pair_diffs$hr,
c(coef(group_means)[1],coef(group_means)[1]+coef(group_means)[2:length(coef(group_means))]))

## What about this?
coef(est4)

## Notice that all of the coefficients are the same.
coef(est4)[3:length(coef(est4))]
coef(group_means)[2:length(coef(group_means))]

## So what is happening with the intercept?
## Maybe this will help us understand:
## Create yhat for rank_own==1, 0, and .5 (even though rank_own==.5 is
## excluded), it turns out that the mean of rank_own is .5
mean(filter(meddat2,rank_own!=.5)$rank_own)
pred_est4 <- predict(est4,newdata=data.frame(rank_own=c(0,.5,1),bmF="1"))
pred_est4
all.equal(pred_est4[["2"]],coef(group_means)[[1]])
## So, again, the intercept is the **predicted** mean of the outcome in the first group (the
## excluded group) when the explanatory variable is 0. (Although, as we see
## here, this prediction is not exactly the same as the mean of the outcome in
## that group).
meddat2 %>% filter(bmF=="1") %>% dplyr::select( rank_own, nhOwn, HomRate08)
meddat2 %>% filter(bmF=="1") %>% dplyr::select( rank_own, nhOwn, HomRate08) %>%
    summarize(mean(HomRate08))
```

## Graphing the possibly non-linear/heterogeneous relationships

This next allows us to explore the within pair differences --- here we look at how differences in proportion home ownership within pair relate to differences in homocide rate within pair.

```{r}
## More exploring about the pair-differences
g1 <- ggplot(data=pair_diffs,aes(x=own_diff,y=hr_diff))+
    geom_point()+
    geom_smooth(method="loess",se = FALSE,method.args=list(family="gaussian",deg=2,span=.6)) +
    geom_smooth(method="loess",se =
        FALSE,method.args=list(family="symmetric",span=.8,deg=1),col="orange")

g1
```

## Outcome analysis 2: Size of the difference within pairs

So far our analysis asked, "Did the neighborhood in the pair with higher home ownership have less or more violence, on average, than the neighborhood in the pair with less home ownership." This ignores the *size* of the difference in proportion owning a home and in exchange allows us to simplify the question. That said, we can also look how the mean neighborhood violence differs given different magnitude of differences within pair. What about when we are looking at the difference in violence associated linearly
with continuous differences in home ownership? (i.e. looking at how differences
in violence are associated with differences in home ownership in proportions).
Notice below that we have the same methods as above (only that the
`difference_in_means` doesn't work because we don't have a binary explanatory
variable.)

## Outcome analysis 2: Size of the difference within pairs

In each case the interpretation is about average differences in outcome for a
one unit difference in the explanatory variable (which is really large, it is
the maximum difference between any two neighborhoods on the explanatory.)

```{r echo=TRUE}
## Still restricting attention to pairs that are not identical so that we can be
## using the same observations for both analyses.
est1cont <- lm_robust(hr_diff~own_diff-1,data=pair_diffs)

est3cont <- lm_robust(HomRate08~nhOwn,fixed_effects=~bmF,data=meddat2,subset=rank_own!=.5)
est4cont <- lm_robust(HomRate08~nhOwn+bmF,data=meddat2,subset=rank_own!=.5)

meddat2 <- meddat2 %>% group_by(bmF) %>% mutate(own_md=nhOwn - mean(nhOwn)) %>% ungroup()
est5cont <- lm_robust(hr_md~own_md,data=meddat2,subset=rank_own!=.5)

meddat2 %>% filter(bmF=="1") %>% dplyr::select(nhOwn,rank_own,own_md,HomRate08,hr_md) %>% head()
pair_diffs %>% filter(bmF=="1")

## Again, showing how all of these aproaches which appear different on their face are the same:
rbind(est1cont=coef(est1cont)[["own_diff"]],
    est3cont=coef(est3cont)[["nhOwn"]],
    est4cont=coef(est4cont)[["nhOwn"]],
    est5cont=coef(est5cont)[["own_md"]])

```

## Summary of non-bipartite matching

 - We can make pairs of units within which we can claim to have broken the
   relationship between many background covariates and another causal driver,
   intervention, or treatment even if that $Z$ variable has many values. This
   is called **non-bipartite matching**.
 - We can compare these relationships to (1) our substantive and contextual
   knowledge and (2) the kind of $X \rightarrow Z$ relationships we would see
   had $Z$ been randomly assigned within pair (imagine $Z$ having multiple
   values and the higher value being assigned at random within pair).
 - We can compare how $Z \rightarrow Y$ conditional on pair in a variety of
   ways: estimation and testing comparing the higher-vs-lower treatment value member of a pair or by averaging over the size of the higher-vs-lower treatment value differences (say, using OLS to focus on the linear relationship). We can also visualize the relationships to assess linearity and/or learn more.

# 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)
```

## 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
wrkdat0<-canadamapdat[!is.na(canadamapdat$vm.community.norm2) &
		      !is.na(canadamapdat$social.capital01),]
## Take a random sample so that the lecture compiles
set.seed(12345)
wrkdat <- droplevels(sample_n(wrkdat0,500))
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<-round(scalar.dist(wrkdat$vmdaPct),1)
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}

canada_nearlist <- list(covs = as.matrix(wrkdat[, c("csd.pop","community_area_km")]),
        pairs=c(csd.pop=100000,community_area_km=5))

## Try not to match two people with the same perceptions --- that doesn't add anything to our analysis
canada_farlist <- list(covs = as.matrix(wrkdat[, "vm.community.norm2"]),
    pairs = c(vm.community.norm2 = .1))

canada_pairs <- nmatch(
  dist_mat = vmdaDist,
  near = canada_nearlist,
  far = canada_farlist,
  subset_weight = 1,
  solver = solverlist
)

## Version using nonbimatch
## vmdaDistMat <- distancematrix(vmdaDist)
## nbp1match<-nonbimatch(vmdaDistMat)
## nbp1<-get.sets(nbp1match$matches,remove.unpaired=TRUE)
wrkdat$id <- row.names(wrkdat)
canada_pairs_df <- nmatch_to_df(canada_pairs,origid=wrkdat$id)
nrow(canada_pairs_df)

## So, in matched set 1 (bm==1) we see two neighborhoods:
canada_pairs_df %>% filter(bm==1)

# The nmatch_to_df function creates a column labeled "bm" which contains
wrkdat2 <- inner_join(wrkdat, canada_pairs_df, by = "id")
wrkdat2 <- droplevels(wrkdat2)
stopifnot(nrow(wrkdat2) == nrow(canada_pairs_df))

## Number of matches:
# wrkdat2$bm is the matched set indicator.
stopifnot(length(unique(wrkdat2$bm)) == nrow(wrkdat2) / 2)
nrow(canada_pairs_df)
nrow(wrkdat2)
## Notice some observations were not matched:
nrow(wrkdat)
wrkdat2$nbp1 <- wrkdat2$bm

##wrkdat[names(nbp1),"nbp1"]<-nbp1
##nbp1[1:5]
##table(is.na(wrkdat$nbp1)) ## recall the "ghost message"
```

## Inspect the solution

```{r nbpsol, echo=TRUE }
wrkdat2[order(wrkdat2$nbp1),c("nbp1","vmdaPct","vm.community.norm2")][1:6,]
## table(wrkdat2$nbp1)
nbp1vmdiffs <- tapply(wrkdat2$vmdaPct, wrkdat2$nbp1, function(x) {
  abs(diff(x))
})
nbp1percdiffs <- tapply(wrkdat2$vm.community.norm2, wrkdat2$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(wrkdat2,
  yvar = "vmdaPct", xvar = "vm.community.norm2", strata = "nbp1", points = FALSE,
  ylim = range(wrkdat2$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. Here using `xBalance` because it can handle continuous treatments.

```{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=wrkdat2)
xb1$overall
xb1$results[,c("z","p"),"nbp1"]
```

## Assess balance: Approach with higher-vs-lower

No treatment and control groups to compare. But we can still compare the
**relationships** between which person is higher versus lower on the adjusted
variable (`vmdaPct`) and other covariates conditional on 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 balnbp1_ranked, cache=TRUE }
wrkdat2$id <- row.names(wrkdat2)
wrkdat2 <- wrkdat2 %>% group_by(nbp1) %>%
    mutate(vmdaPct_ranked=rank(vmdaPct,ties="random")-1)
wrkdat2 <- as.data.frame(wrkdat2)
row.names(wrkdat2) <- wrkdat2$id

wrkdat2 %>% arrange(nbp1) %>% dplyr::select(nbp1,vmdaPct,vmdaPct_ranked) %>% head()

thecovs <- c(
  "age", "income.coded", "education", "x.years", "sex",
  "csd.pop", "vm.csd", "community_area_km"
)
balfmla_ranked<-reformulate(thecovs,response="vmdaPct_ranked")
xb1_ranked<-balanceTest(update(balfmla_ranked,.~.+strata(nbp1)),data=wrkdat2,p.adjust="none")
xb1_ranked$overall
xb1_ranked$results[,,"nbp1"]
```

## Strength of the treatment

The difference in "treatment" within sets varies --- and so we expect the size
of the effect to vary. For example, consider the ratio of  objective context
differences to  perceived context differences:

```{r treatmentstr, echo=TRUE}
summary(nbp1vmdiffs)
summary(nbp1percdiffs)
percDist <- scalar.dist(wrkdat2$vm.community.norm2*100)
da <- vmdaDist[1:5,1:5]
perc <- percDist[1:5,1:5]
da/perc
```


```{r summarizenbp}

## Size of causal driver differences: bigger is better for statistical power later
perc_diffs_by_nbp1 <- wrkdat2 %>% filter(!is.na(nbp1)) %>% group_by(nbp1) %>% summarize(perc_diff=diff(vm.community.norm2))

## Notice no pairs with 0 difference in the designmatch result
summary(abs(perc_diffs_by_nbp1$perc_diff))

```



## Assess hypotheses about effects

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

```{r eval=TRUE,echo=TRUE}
library(coin)
wrkdat2$nbp1F <- factor(wrkdat2$nbp1)
test1<-independence_test(social.capital01~vm.community.norm2|nbp1F,data=wrkdat2[!is.na(wrkdat2$nbp1F),])
test1
```


## 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}
wrkdat2$scRank<-with(wrkdat2,rank.pairs(social.capital01,nbp1))
wrkdat2$vmCRank<-with(wrkdat2,rank.pairs(vm.community.norm2,nbp1))
wrkdat2[order(wrkdat2$nbp1),c("nbp1","social.capital01","scRank","vm.community.norm2","vmCRank")][1:6,]
with(wrkdat2,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"}
## By default, this rescales each observation to be the distance from the group mean.
align.by.block<-function (x, block, fn = mean, thenames=NULL) {
	newx<-unsplit(lapply(split(x, block), function(x) {
				     x - fn(x) 
				 }), block)
	if (!is.null(names)) {
		names(newx) <- thenames
	}
	return(newx)
}
```

```{r}
wrkdat2$scMD <- with(wrkdat2, align.by.block(social.capital01, nbp1))
wrkdat2$vmcn2MD <- with(wrkdat2, align.by.block(vm.community.norm2, nbp1))
wrkdat2[order(wrkdat2$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(wrkdat2$scMD,wrkdat2$nbp1,function(x){ abs(diff(x)) }))
## summary(tapply(wrkdat2$social.capital01,wrkdat2$nbp1,function(x){ abs(diff(x)) }))
lm1 <- lm_robust(scMD ~ vmcn2MD, data = wrkdat2[!is.na(wrkdat2$nbp1), ])
lm1
lm1_fe <- lm_robust(social.capital01~vm.community.norm2,fixed_effects=~nbp1,data=wrkdat2[!is.na(wrkdat2$nbp1), ])
lm1_fe
##library(fixest) ## for more speed with strata-by-strata estimation
```
## Summarize mean differences within pairs

```{r warning=FALSE,cache=TRUE}
lm2 <- lm_robust(scMD~vmcn2MD,data=wrkdat2[!is.na(wrkdat2$nbp1),])
lm2
lm3 <- lm_robust(social.capital01 ~ vm.community.norm2, fixed_effects = ~nbp1, data = wrkdat2, subset = !is.na(wrkdat2$nbp1))
lm3
table(wrkdat2$vmCRank,exclude=c())
lm4 <- lm_robust(social.capital01 ~ I(vmCRank - 1), fixed_effects = ~nbp1, data = wrkdat2, subset = !is.na(wrkdat2$nbp1))
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}
lm2
lm3
pairdiffs <- wrkdat2 %>%
  filter(!is.na(vmCRank) & !is.na(social.capital01) & !is.na(nbp1)) %>%
  group_by(vmCRank) %>%
  summarize(mnsc = mean(social.capital01))
wrkdat2[order(wrkdat2$nbp1),c("social.capital01","scRank","scMD","vm.community.norm2","vmcn2MD","vmCRank","nbp1")][1:4,]
lm4
```

## Summarize mean differences within pairs

```{r}
summary(wrkdat2$vmcn2MD)
summary(wrkdat2$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(wrkdat2$vmcn2MD,na.rm=TRUE),2)`. 

The model predicts that social capital would differ by about `r
coef(lm1)[[2]]*.4` for such a difference. This is about `r
coef(lm1)[[2]]*.4/sd(wrkdat2$scMD,na.rm=TRUE)` of a standard deviation of the
social capital scale. Or about `r
coef(lm1)[[2]]*.4/abs(diff(range(wrkdat2$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(wrkdat2,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)
wrkdat2$vmCbi<-ave(wrkdat2$vm.community.norm2,wrkdat2$nbp1)
lmer1<-lmer(social.capital01~vm.community.norm2+vmCbi+(1|nbp1),data=wrkdat2)
confint(lmer1)["vm.community.norm2",]
## Notice that we may have less information than we thought (some clustered of observation by DA).
table(table(wrkdat2$dauid))
## So maybe:
lmer2<-lmer(social.capital01~vm.community.norm2+vmCbi+(1|nbp1)+(1|dauid),data=wrkdat2)
confint(lmer2)["vm.community.norm2",]
```

## Other applications of non-bipartite matching?

See: DOS Chapter 11.

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

## Summary and Questions:

 - We can make pairs of people within which we can claim to have broken the
   relationship between many background covariates and another causal driver,
   intervention, or treatment even if that $Z$ variable has many values. This
   is called non-bipartite matching.
 - We can compare these relationships to (1) our substantive and contextual
   knowledge and (2) the kind of $X \rightarrow Z$ relationships we would see
   had $Z$ been randomly assigned within pair (imagine $Z$ having multiple
   values and the higher value being assigned at random within pair).
 - We can compare how $Z \rightarrow Y$ conditional on pair in a variety of
   ways: estimation and testing.


## References