---
title: "Using prim to estimate highest density difference regions"
output: rmarkdown::html_vignette
description: Blah
date: "`r format(Sys.time(), '%d %B %Y')`"
author: Tarn Duong https://mvstat.net/tduong/
vignette: >
  %\VignetteIndexEntry{Using prim to estimate highest density difference regions}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r, echo=FALSE, message=FALSE}
knitr::opts_chunk$set(collapse=TRUE, comment="#>", fig.width=5, fig.height=5, fig.align="center", global.par=TRUE, dpi=96) 
options(tibble.print_min=4L, tibble.print_max=4L)
```


The Patient Rule Induction Method (PRIM) was introduced by Friedman and Fisher (1999). It is a technique from data mining for finding `interesting' regions in high-dimensional data. We start with regression-type data *(**X**~1~, Y~1~), ..., (**X**~n~, Y~n~)* where ***X**~i~* is *d*-dimensional and *Y~i~* is a scalar response variable. We are interested in the conditional expectation function

*m(**x**) = **E** (Y | **X** = **x**)*.

In the case where we have 2 samples, we can label the response as a binary variable with 

*Y~i~ = 1* if ***X**~i~* is in sample 1

or 

*Y~i~ = -1* if ***X**~i~* is in sample 2.

Then PRIM finds the regions where the samples are most different.  Here we have a positive HDR (where sample 1 points dominate) and a negative HDR (where sample 2 points dominate).

We look at a 3-dimensional data set `quasiflow` included in the `prim` library. It is a randomly generated data set from 
two normal mixture distributions whose structure mimics some light scattering data, taken from a machine known as a flow cytometer. 

```{r}
library(prim)
data(quasiflow)
yflow <- quasiflow[,4]
xflow <- quasiflow[,1:3]
xflowp <- quasiflow[yflow==1,1:3]
xflown <- quasiflow[yflow==-1,1:3]
```

We can think of `xflowp` as flow cytometric measurements from an HIV+ patient, and `xflown` from an HIV-- patient.
```{r}
pairs(xflowp, cex=0.5, pch=16, col=grey(0,0.1), xlim=c(0,1), ylim=c(0,1))
pairs(xflown, cex=0.5, pch=16, col=grey(0,0.1), xlim=c(0,1), ylim=c(0,1))
```

There are two ways of using `prim.box` to estimate where the two samples are most different (or equivalently 
to estimate the HDRs of the difference of the density functions). In the first way, we assume that we have suitable values for the thresholds. Then we can use
```{r}=
qflow.thr <- c(0.38, -0.23)
qflow.prim <- prim.box(x=xflow, y=yflow, threshold=qflow.thr, threshold.type=0)
```

An alternative is compute PRIM box sequences which cover the range of the response variable `y`, and then use `prim.hdr` to experiment with different threshold values. This two-step process is more efficient and faster than calling `prim.box` for each different threshold. We are happy with the positive HDR threshold so we can compute the positive HDR directly: 
``` {r}
qflow.hdr.pos <- prim.box(x=xflow, y=yflow, threshold=0.38, threshold.type=1)
summary(qflow.hdr.pos)
```

On the other hand, we are not sure about the negative HDR thresholds. So we try several different values for `threshold`. 
```{r}
qflow.neg <- prim.box(x=xflow, y=yflow, threshold.type=-1)
qflow.hdr.neg1 <- prim.hdr(qflow.neg, threshold=-0.23, threshold.type=-1)
qflow.hdr.neg2 <- prim.hdr(qflow.neg, threshold=-0.43, threshold.type=-1)
qflow.hdr.neg3 <- prim.hdr(qflow.neg, threshold=-0.63, threshold.type=-1)
```

After examining the summaries and plots,
```{r}
summary(qflow.hdr.neg1)
```

we choose `qflow.hdr.neg1` to combine with `qflow.hdr.pos`.
```{r}
qflow.prim2 <- prim.combine(qflow.hdr.pos, qflow.hdr.neg1)
summary(qflow.prim2)
```

In the plot below, the positive HDR is coloured orange (where the HIV+ sample is more prevalent), and the  negative HDR is coloured blue (where the HIV- sample is more prevalent)  
```{r, fig.asp=1, fig.height=10}
plot(qflow.prim2, x.pt=xflow, pch=16, cex=0.5, alpha=0.1)
```

or equivalently as a 3D scatter plot. 
```{r, fig.asp=1, fig.height=10}
plot(qflow.prim2, x.pt=xflow, pch=16, cex=0.5, alpha=0.1, splom=FALSE, colkey=FALSE, ticktype="detailed")
```


## References

Friedman, J. H. and Fisher, N. I. (1999). Bump-hunting for high dimensional data. *Statistics and Computing*, **9**, 123-143.