day16.Rmd

---
title: Non-bipartite Matching
date: '`r format(Sys.Date(), "%B %d, %Y")`'
author: ICPSR 2023 Session 1
bibliography:
 - refs.bib
 - BIB/master.bib
 - BIB/misc.bib
fontsize: 10pt
geometry: margin=1in
graphics: yes
biblio-style: authoryear-comp
output:
  beamer_presentation:
    slide_level: 2
    keep_tex: true
    latex_engine: xelatex
    citation_package: biblatex
    template: icpsr.beamer
    includes:
        in_header:
           - defs-all.sty
---


<!-- Make this document using library(rmarkdown); render("day12.Rmd") -->


```{r include=FALSE, cache=FALSE}
# Some customization.  You can alter or delete as desired (if you know what you are doing).
# knitr settings to control how R chunks work.
rm(list=ls())

require(knitr)

## This plus size="\\scriptsize" from https://stackoverflow.com/questions/26372138/beamer-presentation-rstudio-change-font-size-for-chunk

knitr::knit_hooks$set(mysize = function(before, options, envir) {
  if (before)
    return(options$size)
})

knit_hooks$set(plotdefault = function(before, options, envir) {
    if (before) par(mar = c(3, 3, .1, .1),oma=rep(0,4),mgp=c(1.5,.5,0))
})

opts_chunk$set(
  tidy=FALSE,     # display code as typed
  echo=TRUE,
  results='markup',
  strip.white=TRUE,
  fig.path='figs/fig',
  cache=FALSE,
  highlight=TRUE,
  width.cutoff=132,
  size='\\scriptsize',
  out.width='.8\\textwidth',
  fig.retina=FALSE,
  message=FALSE,
  comment=NA,
  mysize=TRUE,
  plotdefault=TRUE)

if(!file.exists('figs')) dir.create('figs')

options(digits=4,
	scipen=8,
	width=132
	)
```

```{r eval=FALSE, include=FALSE, echo=FALSE}
## Run this only once and then not again until we want a new version from github
library('devtools')
library('withr')
with_libpaths('./lib', install_github("markmfredrickson/RItools"), 'pre')

## Having downloaded optmatch from Box OR http://jakebowers.org/ICPSR for mac (tgz) or windows (zip)
## For Mac
## with_libpaths('./lib',install.packages('optmatch_0.9-8.9003.tgz', repos=NULL),'pre')
## For Windows
## with_libpaths('./lib',install.packages('optmatch_0.9-8.9003.zip', repos=NULL),'pre')

## Or if you have all of the required libraries (like fortran and c++) for compilation use
with_libpaths('./lib', install_github("markmfredrickson/optmatch"), 'pre')
```

```{r echo=FALSE}
library(dplyr)
library(RItools,lib.loc="./lib")
library(optmatch,lib.loc="./lib")
library(nbpMatching)
library(lmtest)
library(sandwich)
```

## Today

\begin{enumerate}
\item Agenda: Non-bipartite matching.
\item Reading for this week: (1) Non-bipartite matching DOS Chap 11 and \autocite{lu2011optimal}  and DOS 12 (longitudinal applications of non-bipartite matching) and (2) Sensitivity analysis DOS Chap 3 \autocite{hhh2010}
\item Questions arising from the reading or assignments or life?
\item Recap: How to make a matched design that warrants causal interpretations of comparisons (decisions and strategies that are part of research design; matching on missingness and `fill.NAs`; `exactMatch`; `caliper`; `min.controls`; `effectiveSampleSize`); Modes of statistical inference for causal quantities (and the argument for the design-based / as-if-randomized approach after matching; the `optmatch` matching creates a design like a block-randomized study so analyze it as if it were a block-randomized study. ).
\item Tomorrow and Wed: Sensitivity Analysis
\item Thurs: Interference
\end{enumerate}


## How do perceptions of place influence attitudes?

\autocite{wong2012jop} set out to measure perceptions of environments using an
internet survey of Canadians during 2012 which each respondent drew a free hand
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 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 perceptions are 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 here matches people with similar neighborhood proportions of visible minorities.

```{r }
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, 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 }
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)
```

```{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

Now we cannot ask whether the treatment and control groups look appropriately similar, but we can still compare the **relationships** between the adjusted variable (`vmdaPct`) and other covariates.

```{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, to give
  special attention to `csd.pop` (so that we are not comparing people in small
  towns to people in large towns), and `community_area_km` (so that we are not
  comparing people who drew big maps to people who drew small maps). And we
  might also want to drop the 8 least well matched observations. (Choosing 8
  here arbitrarily.)

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

```{r 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 cache=TRUE}
vmdaPhPenMat <- distancematrix(vmdaPhPen)
nbp2match<-nonbimatch(vmdaPhPenMat)
nbp2<-get.sets(nbp2match$matches,remove.unpaired=TRUE)
wrkdat[names(nbp2),"nbp2"]<-nbp2
```

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

## Strengthening the treatment

The difference in "treatment" within sets varies:
```{r}
summary(nbp1vmdiffs)
summary(nbp1percdiffs)

percDist <- scalar.dist(wrkdat$vm.community.norm2*100)

da <- vmdaDist[1:5,1:5]
perc <- percDist[1:5,1:5]

da/perc

da + 1000*(perc < 2)

```

\normalsize
To prevent so 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.

\medskip

Class challenge: How would you do this? (use the code from above to do this or at least sketch it)

## Now assess hypotheses about effects

Now, test the hypothesis of no relationship between perceptions
`vm.community.norm2` and `social capital`.

```{r eval=FALSE,echo=FALSE,results="hide"}
## These are the same test in this case
library(coin)
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

One idea is to ask whether, within a pair, the person who perceives more
visible minorities in their community tends to be higher (or lower) in
`social.capital`. What do you think the following analysis shows?

```{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)
}

library(sandwich)
library(lmtest)
source(url("http://jakebowers.org/ICPSR/confintHC.R"))
```

```{r}
wrkdat$scMD<-with(wrkdat,align.by.block(social.capital01,nbp1))
wrkdat$vmcn2MD<-with(wrkdat,align.by.block(vm.community.norm2,nbp1))
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(scMD~vmcn2MD,data=wrkdat[!is.na(wrkdat$nbp1),])
theHC2ci<-confint.HC(lm1,level=.95,parm="vmcn2MD",thevcov=vcovHC(lm1,type="HC2"))
c(meandiff=coef(lm1)[[2]],ci=theHC2ci)
```

## Summarize mean differences within pairs

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

c(meandiff=coef(lm1)[[2]],ci=theHC2ci)
```

\normalsize
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.

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


## 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=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}).


## Anything Else?

## References