---
title: "Prototyping with `representr`"
author: "Andee Kaplan and Rebecca C. Steorts"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    fig_caption: true
vignette: >
  %\VignetteIndexEntry{Prototyping with `representr`}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
references:
- id: blink
  title: "blink: Record Linkage for Empirically Motivated Priors"
  author:
  - family: Steorts
    given: Rebecca C.
  note: "R package version 0.1.0"
  URL: 'https://CRAN.R-project.org/package=blink'
  issued:
    year: ""
- id: steorts2015entity
  title: Entity resolution with empirically motivated priors
  author:
  - family: Steorts
    given: Rebecca C.
  container-title: Bayesian Analysis
  volume: 10
  issue: 4
  page: 849-875
  type: article-journal
  issued:
    year: 2015
  publisher: International Society for Bayesian Analysis
- id: marchant2020d
  title: "d-blink: Distributed end-to-end Bayesian entity resolution"
  author:
    - family: Marchant  
      given: Neil G
    - family: Kaplan
      given: Andee
    - family: Elazar
      given: Daniel N
    - family: Rubinstein
      given: Benjamin IP
    - family: Steorts
      given: Rebecca C
  container-title: Journal of Computational and Graphical Statistics
  volume: just-accepted
  page: 1--42
  type: article-journal
  issued:
    year: 2020
  publisher: Taylor & Francis
- id: kaplan2018posterior
  title: "Posterior Prototyping: Bridging the Gap between Record Linkage and Regression"
  author:
    - family: Kaplan
      given: Andee
    - family: Betancourt
      given: Brenda
    - family: Steorts
      given: Rebecca C.
  container-title: "arXiv preprint arXiv:1810.01538"
  type: article-journal
  URL: "https://arxiv.org/abs/1810.01538"
  issued:
    year: ""
- id: dblinkR
  title: "dblinkR: R interface for dblink"
  author:
    - family: Marchant
      given: Neil
  issued:
    year: "2021"
  container-title: "R package version 0.2.0"
  URL: "https://github.com/cleanzr/dblinkR"
- id: binette2020almost
  title: "(Almost) All of Entity Resolution"
  author:
    - family: Binette
      given: Olivier
    - family: Steorts
      given: Rebecca C
  container-title: "arXiv preprint arXiv:2008.04443"
  URL: "https://arxiv.org/abs/2008.04443"
  issued:
    year: ""
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE
)

library(ggplot2)
set.seed(1234)
```


# Background and Motivating Example

Record linkage (entity resolution or de-duplication) is used to join multiple databases to remove duplicate entities. While record linkage removes the duplicate entities from the data, many researchers are interested in performing inference, prediction, or post-linkage analysis on the linked data (e.g., regression or capture-recapture), which we call the *downstream task*. Depending on the downstream task, one may wish to find the most representative record before performing the post-linkage analysis. For example, when the values of features used in a downstream task differ for linked data, which values should be used? This is where `representr` comes in. Before introducing our new package `representr` from the paper @kaplan2018posterior, we first provide an introduction to record linkage. 

Throughout this vignette, we will use data that is available in the `representr` package, `rl_reg1` (rl = record linkage, reg = regression, 1 = amount of noisiness).

```{r ex-records}
# load library
library(representr)

# load data
data("rl_reg1") # data for record linkage and regression
data("identity.rl_reg1") # true identity of each record
```

```{r echo = FALSE}
knitr::kable(head(rl_reg1))
```

This is simulated data, which consists of 500 records with 30% duplication and the following attributes:

- `fname`: First name
- `lname`: Last name
- `bm`: Birth month (numeric)
- `bd`: Birth day
- `by`: Birth year
- `sex`: Sex ("M" or "F")
- `education`: Education level ("Less than a high school diploma", ""High school graduates, no college", "Some college or associate degree", "Bachelor's degree only", or "Advanced degree")
- `income`: Yearly income (in $1000s)
- `bp`: Systolic blood pressure

## Record linkage

Before we perform prototyping to get a representative data set using `representr`, we must first perform record linkage to remove duplication in the data set.  In the absence of unique identifier (such as a social security number), we can use probabilistic methods to perform record linkage. We recommend the use of clustering records to a latent entity, known in the literature as graphical entity resolution. See [@binette2020almost] for a review. 

<!-- For the examples in this vignette, we will use `blink` [@blink],a Bayesian record linkage model [@steorts2015entity]. For this example, we will run only a small number of iterations, which will certainly not have converged. For a faster (more scalable version of this), we recommend recent work of [@marchant2020d] using either dblink or dblinkR.  -->

<!-- For a real example, we would want to let this model run longer. -->

<!-- ```{r rl, results='hide'} -->
<!-- # load blink -->
<!-- library(blink) -->

<!-- # params for running record linkage -->
<!-- a <- 1; b <- 99 # distortion hyperparams -->
<!-- c <- 1 # string density hyperparams -->
<!-- num.gs <- 10 # number of iterations -->
<!-- M <- nrow(rl_reg1) # upper bound on number of entities -->
<!-- str_idx <- c(1, 2) # string columns -->
<!-- cat_idx <- c(3, 4, 5) # categorical columns -->

<!-- # data prep  -->
<!-- # X.c contains the categorical variables -->
<!-- # X.s contains the string variables  -->
<!-- X.c <- apply(as.matrix(rl_reg1[, cat_idx]), 2, as.character) -->
<!-- X.s <- as.matrix(rl_reg1[, str_idx]) -->

<!-- # X.c and X.s include all files "stacked" on top of each other. -->
<!-- # The vector below keeps track of which rows of X.c and X.s are in which files. -->
<!-- file.num <- rep(1, nrow(rl_reg1)) -->

<!-- # perform record linkage   -->
<!-- linkage.rl <- rl.gibbs(file.num, X.s, X.c, num.gs=num.gs, a=a, b=b, c=c, d=utils::adist, M=M) -->
<!-- ``` -->

<!-- ```{r delete-blink-res, echo=FALSE, results="hide"} -->
<!-- file.remove(dir(path = ".", pattern="lambda*") ) -->
<!-- ``` -->

For the examples in this vignette, we have fit the model in recent work of [@marchant2020d] using `dblinkR` [@dblinkR]. Please see associated vignette for details on using `dblinkR`. We load the results of running this record linkage model for $100,000$ iterations, which have been stored in the package as a data object called `linkage.rl`.

```{r}
data("linkage.rl")
```


## The Downstream Task

After record linkage is complete, one may want to perform analyses of the linked data. This is what we call the "downstream task". As motivation, consider modeling blood pressure (bp) using the following two features (covariates): income and sex in our example data `rl_reg1`. We want to fit this model after performing record linkage using the following features: first and last name and full data of birth. Here is an example of four records that represent the same individual (based on the results from record linkage) using data that is in the `representr` package. 

```{r, echo=FALSE}
# take the last iteration of linkage for lambda value
lambda <- linkage.rl[nrow(linkage.rl),]

# split data by linkage results
clusters <- split(rl_reg1, lambda)
```

```{r, echo=FALSE}
knitr::kable(clusters[[names(which(table(lambda) == 4)[1])]])
```

Examination of this table raises important questions that need to be addressed before performing a particular downstream task, such as which values of bp, income, and sex should be used as the representative features (or covariates) in a regression model? In this vignette, we will provide multiple solutions to this question using a prototyping approach.

# Prototyping methods

We have four methods to choose or create the representative record from linked data included in `representr`. This process is a function of the data and the linkage structure, and we present both probabilistic and deterministic functions. The result in all cases is a representative data set to be passed on to the downstream task. The prototyping is completed using the `represent()` function.


## Random Prototyping

Our first proposal to choose a representative record (*prototype*) for a cluster is the simplest and serves as a baseline or benchmark. One simply chooses the representative record uniformly at random or using a more informed distribution. 

For demonstration purposes, we can create a representative dataset using the last iteration of the results from running the record linkage model using `blink`. This is accomplished using the `represent()` function, and passing through the type of prototyping to be `proto_random`.

```{r random}
# ids for representative records (random)
random_id <- represent(rl_reg1, lambda, "proto_random", parallel = FALSE)
rep_random <- rl_reg1[random_id,] # representative records (random)
```

We can have a look at a few records chosen as representative in this way.

```{r, echo = FALSE}
knitr::kable(rep_random[as.numeric(names(which(table(lambda) == 4)))[1:5],])
```

## Minimax Prototyping

Our second proposal to choose a representative record is to select the record that "most closely captures" that of the latent entity. Of course, this is quite subjective. We propose selecting the record whose farthest neighbors within the cluster is closest, where closeness is measured by a record distance function, $d_r(\cdot)$. We can write this as the record $r = (i, j)$ within each cluster $\Lambda_{j'}$ such that
$$
r = \arg\min\limits_{(i, j) \in \Lambda_{j'}} \max\limits_{(i^*, j^*) \in \Lambda_{j'}} d_r((i, j), (i^*, j^*)).
$$
The result is a set of representative records, one for each latent individual, that is closest to the other records in each cluster. When there is a tie within the cluster, we select a record uniformly at random.

There are many distance functions that can be used for $d_r(\cdot, \cdot)$. We define the distance function to be a weighted average of individual variable-level distances that depend on the column type. Given two records, $(i, j)$ and $(i*, j*)$, we use a weighted average of column-wise distances (based on the column type) to produce the following single distance metric:
$$
d_r((i, j), (i*, j*)) = \sum\limits_{\ell = 1}^p w_\ell d_{r\ell}((i, j), (i^*, j^*)),
$$
where $\sum\limits_{\ell = 1}^p w_\ell = 1$. The column-wise distance functions $d_{r\ell}(\cdot, \cdot)$ we use are presented below.

| Column      | $d_{r\ell}(\cdot, \cdot)$                                      |
|:------------|----------------------------------------------------------------|
| String      | Any string distance function, i.e. Jaro-Winkler string distance |
| Numeric     | Absolute distance, $d_{r\ell}((i, j), (i^*, j^*)) = \mid x_{ij\ell} - x_{i^*j^*\ell} \mid$ |
| Categorical | Binary distance, $d_{r\ell}((i, j), (i^*, j^*)) = \mathbb{I}(x_{ij\ell} != x_{i^*j^*\ell})$ |
| Ordinal     | Absolute distance between levels. Let $\gamma(x_{ij\ell})$ be the order of the value $x_{ij\ell}$, then $d_{r\ell}((i, j), (i^*, j^*)) = \mid \gamma(x_{ij\ell}) - \gamma(x_{i^*j^*\ell}) \mid$  |

The weighting of variable distances is used to place importance on individual features according to prior knowledge of the data set and to scale the feature distances to a common range. In this vignette, we scale all column-wise distances to be values between $0$ and $1$.

Again, we can create a representative dataset using the last iteration of the results from running the record linkage model using `blink`. But this time we need to specify some more parameters, like what types the columns are. This is accomplished using the `represent()` function, and passing through the type of prototyping to be `proto_minimax`.

```{r minimax}
# additional parameters for minimax prototyping
# need column types, the order levels for any ordinal variables, and column weights
col_type <- c("string", "string", "numeric", "numeric", "numeric", "categorical", "ordinal", "numeric", "numeric")
orders <- list(education = c("Less than a high school diploma", "High school graduates, no college", "Some college or associate degree", "Bachelor's degree only", "Advanced degree"))
weights <- c(.25, .25, .05, .05, .1, .15, .05, .05, .05)

# ids for representative records (minimax)
minimax_id <- represent(rl_reg1, linkage.rl[nrow(linkage.rl),], "proto_minimax",
                        distance = dist_col_type, col_type = col_type, 
                        weights = weights, orders = orders, scale = TRUE, parallel = FALSE)
rep_minimax <- rl_reg1[minimax_id,] # representative records (minimax)
```

We can have a look at some of the representative records chosen via minimax prototyping.

```{r, echo = FALSE}
knitr::kable(rep_minimax[as.numeric(names(which(table(lambda) == 4)))[1:5],])
```

## Composite Records

Our third proposal to choose a representative record is by aggregating the records (in each cluster) to form a composite record that includes information from each linked record. The form of aggregation can depend on the column type, and the aggregation itself can be weighted by some prior knowledge of the data sources or use the posterior information from the record linkage model. For quantitative variables, we use a weighted arithmetic mean to combine linked values, whereas for categorical variables, a weighted majority vote is used. For string variables, we use a weighted majority vote for each character, which allows for noisy strings to differ on a continuum. This is accomplished using the `represent()` function, and passing through the type of prototyping to be `composite`.

```{r composite}
# representative records (composite)
rep_composite <- represent(rl_reg1, linkage.rl[nrow(linkage.rl),], "composite", col_type = col_type, parallel = FALSE)
```

We can have a look at some of the representative records.

```{r, echo = FALSE}
knitr::kable(rep_composite[as.numeric(names(which(table(lambda) == 4)))[1:5],])
```

## Posterior Prototyping

Our fourth proposal to choose a representative record utilizes the minimax prototyping method in a fully Bayesian setting. This is desirable as the posterior distribution of the linkage is used to weight the downstream tasks, which allows the error from the record linkage task to be naturally propagated into the downstream task.

We propose two methods for utilizing the posterior prototyping (PP) weights --- a weighted downstream task and a thresholded representative data set based on the weights. As already mentioned, PP weights naturally propagate the linkage error into the downstream task, which we now explain. For each MCMC iteration from the Bayesian record linkage model, we obtain the most representative records using minimax prototyping and then compute the probability of each record being selected over all MCMC iterations. The posterior prototyping (PP) probabilities can then either be used as weights for each record in the regression or as a thresholded variant where we only include records whose PP weights are above $0.5$. Note that a record with PP weight above 0.5 has a posterior probability greater than 0.5 of being chosen as a prototype and should be included in the final data set.

```{r pp}
# Posterior prototyping weights
pp_weights <- pp_weights(rl_reg1, linkage.rl[seq(80000, 100000, by = 100), ], 
                         "proto_minimax", distance = dist_col_type, 
                         col_type = col_type, weights = weights, orders = orders, 
                         scale = TRUE, parallel = FALSE)
```

We can look at the minimax PP weights distribution for the true and duplicated records in the data set as an example. Note that the true records consistently have higher PP weights and the proportion of duplicated records with high weights is relatively low.

```{r thresh-plot, fig.width = 5, echo = FALSE, fig.align="center", fig.cap = "The distribution of PP weights as generated from posterior draws of $\\boldsymbol \\Lambda$ colored by if they are true or duplicated records in the original data set. The dotted vertical line shows the threshold value of 0.5. The true records have consistently higher PP weights and the proportion of duplicated records with high weights is relatively low."}
data.frame(pp_weights = pp_weights,
           true = seq_len(nrow(rl_reg1)) %in% unique(identity.rl_reg1),
           included = pp_weights >= .5) -> threshold_df

ggplot(threshold_df) +
  geom_histogram(aes(pp_weights, fill = true), position = "dodge") +
  geom_vline(aes(xintercept = .5), lty = 2, alpha = .5) +
  xlab("Posterior Prototyping weights") +
  ylab("") +
  scale_fill_discrete("True record")  +
  theme(legend.position = "bottom")
```

We can make a representative dataset with these weights by using the cutoff of $0.5$, and look at some of the records.

```{r pp_thresh}
# representative records (PP threshold)
rep_pp_thresh <- rl_reg1[pp_weights > .5, ]
```
```{r, echo = FALSE}
knitr::kable(head(rep_pp_thresh))
```


These four proposed methods each have potential benefits. The goal of prototyping is to select the correct representations of latent entities as often as possible; however, uniform random selection has no means to achieve this goal. Turning to minimax selection, if a distance function can accurately reflect the distance between pairs of records in the data set, then this method may perform well. Alternatively, composite records necessarily alter the data for all entities with multiple copies in the data, affecting some downstream tasks (like linear regression) heavily. The ability of posterior prototyping to propagate record linkage error to the downstream task is an attractive feature and a great strength of the Bayesian paradigm. In addition, the ability to use the entire posterior distribution of the linkage structure also poses the potential for superior downstream performance.

# Evaluation

We can evaluate the performance of our methods by assessing the distributional closeness of the representative dataset to the true records. The distributional closeness of the representative datasets to the true records is useful because one of the benefits of using a two-stage approach to record linkage and downstream analyses is the ability to perform multiple analyses with the same data set. As such, downstream performance of representative records may be dependent on the type of downstream task that is being performed. In order to assess the distributional closeness of the representative data sets to the truth, we use an empirical Kullback-Leibler (KL) divergence metric. Let $\hat{F}_{rep}(\boldsymbol x)$ and $\hat{F}_{true}(\boldsymbol x)$ be the empirical distribution functions for the representative data set and true data set, respectively (with continuous variables transformed to categorical using a histogram approach with statistically equivalent data-dependent bins). The empirical KL divergence metric we use is then defined as
$$
\hat{D}_{KL}(\hat{F}_{rep} || \hat{F}_{true}) = \sum_{\boldsymbol x} \hat{F}_{rep}(\boldsymbol x) \log\left(\frac{\hat{F}_{rep}(\boldsymbol x)}{\hat{F}_{true}(\boldsymbol x)}\right).
$$

This metric is accessed in `representr` using the `emp_kl_div()` command.

```{r empkl_div, message = FALSE, warning = FALSE}
true_dat <- rl_reg1[unique(identity.rl_reg1),] # true records
emp_kl_div(true_dat, rep_random, c("sex"), c("income", "bp"))
emp_kl_div(true_dat, rep_minimax, c("sex"), c("income", "bp"))
emp_kl_div(true_dat, rep_composite, c("sex"), c("income", "bp"))
emp_kl_div(true_dat, rep_pp_thresh, c("sex"), c("income", "bp"))
```

The representative dataset based on the posterior prototyping weights is the closest to the truth using the three variables we might be interested in using for regression. This might indicate that we should use this representation in a downstream model, like linear regression.

# References

