Generating the Toy Data

library(ggplot2)
toy_data <- data.frame(x1=runif(1000), x2=runif(1000))
toy_data$label <- as.factor((toy_data$x2 + toy_data$x1 > 1) * 1)
levels(toy_data$label) <- c("R", "NR")
toy_data$label
##    [1] R  NR R  NR NR R  R  R  NR NR R  R  R  R  R  NR R  R  NR R  R  NR R 
##   [24] R  NR NR R  NR R  R  NR NR NR NR NR R  NR NR R  NR NR R  NR NR NR NR
##   [47] NR NR R  R  NR NR R  NR R  R  NR R  NR NR R  NR R  NR R  R  NR NR NR
##   [70] NR R  R  NR NR R  NR R  R  NR R  R  NR NR R  R  R  NR NR NR R  R  NR
##   [93] R  NR NR NR NR NR NR R  R  R  R  R  R  NR NR NR R  NR R  NR NR NR R 
##  [116] R  R  R  R  NR NR R  NR R  NR R  NR NR R  R  R  R  NR R  R  R  R  R 
##  [139] R  R  R  R  R  NR NR R  NR NR NR R  NR NR NR NR R  NR NR NR R  NR R 
##  [162] R  R  NR NR NR R  NR R  R  NR NR R  NR R  NR R  R  R  NR NR NR R  NR
##  [185] NR R  R  NR NR NR NR R  NR R  NR R  R  R  R  NR NR R  R  NR NR NR R 
##  [208] NR NR R  R  NR NR R  R  NR R  NR R  NR NR NR R  R  NR NR NR R  R  NR
##  [231] R  R  NR NR NR R  R  NR NR R  NR NR NR NR R  R  NR NR NR NR R  R  R 
##  [254] R  R  R  NR R  NR R  R  R  R  R  R  NR NR R  NR NR R  R  R  NR NR R 
##  [277] R  NR R  NR R  R  R  R  R  R  R  NR NR R  R  NR R  NR R  NR NR R  NR
##  [300] R  NR R  R  R  NR R  R  NR R  NR NR NR NR NR NR R  NR R  NR R  NR NR
##  [323] R  NR NR R  R  R  NR NR R  R  NR R  NR NR R  R  NR R  NR R  R  R  R 
##  [346] NR NR R  R  NR NR R  R  R  NR R  NR NR R  NR NR R  R  R  NR R  NR R 
##  [369] R  R  NR NR NR R  R  R  R  NR R  NR R  R  R  R  R  R  R  R  NR NR R 
##  [392] NR NR R  R  NR R  R  R  NR R  R  R  R  R  R  NR NR R  NR R  R  R  NR
##  [415] NR NR R  R  R  NR NR NR NR R  R  R  R  NR NR R  R  NR R  NR R  NR R 
##  [438] R  R  NR NR R  R  NR R  NR R  R  R  NR R  R  NR R  R  NR R  NR NR R 
##  [461] NR R  R  R  R  NR NR R  NR NR NR R  NR R  R  NR NR NR NR NR NR R  R 
##  [484] R  NR NR R  R  R  NR R  NR R  NR R  R  R  R  NR NR R  R  R  NR R  NR
##  [507] R  R  R  R  R  NR R  NR NR R  R  NR R  R  NR R  R  NR NR NR R  R  NR
##  [530] NR NR NR R  NR NR R  NR R  NR NR R  NR NR NR NR R  R  NR NR NR R  NR
##  [553] NR NR R  R  NR NR NR NR R  R  NR R  R  R  R  NR NR NR R  NR NR NR NR
##  [576] NR R  R  NR R  R  R  NR NR R  NR R  NR R  NR R  R  NR NR R  R  R  NR
##  [599] NR R  NR R  R  R  R  NR R  R  R  R  R  R  R  R  R  NR R  NR NR NR NR
##  [622] NR NR R  R  NR NR NR R  R  NR R  NR R  NR NR NR NR NR R  NR NR R  NR
##  [645] NR NR R  R  NR R  R  R  NR R  R  NR R  R  NR R  NR NR R  R  R  R  R 
##  [668] NR NR NR NR NR NR NR R  R  R  R  NR R  NR NR R  R  NR R  R  NR NR NR
##  [691] R  NR NR R  R  R  NR R  NR R  NR NR R  R  R  R  NR R  NR NR NR R  NR
##  [714] R  NR NR R  R  NR R  NR NR R  R  NR R  NR R  R  R  R  R  R  NR NR R 
##  [737] NR NR R  NR NR R  R  R  NR R  R  NR R  R  NR R  R  NR R  NR R  R  R 
##  [760] NR NR NR NR NR R  R  NR NR R  NR NR NR NR NR NR R  R  R  R  NR R  R 
##  [783] NR R  R  NR R  R  NR R  NR R  R  R  NR R  R  NR NR NR NR R  R  R  NR
##  [806] NR R  R  NR NR NR R  R  R  R  R  R  NR R  R  NR R  NR R  R  NR R  NR
##  [829] NR NR NR NR R  R  R  R  NR NR R  R  R  NR R  R  NR R  R  NR NR NR R 
##  [852] R  NR NR NR NR R  NR NR R  R  NR NR R  R  NR NR NR NR NR R  R  NR R 
##  [875] NR NR R  R  NR NR NR R  NR R  R  NR R  R  NR NR NR NR R  NR R  NR NR
##  [898] R  R  R  R  R  R  R  NR R  R  NR R  NR R  R  NR R  NR NR NR NR R  NR
##  [921] R  R  R  NR NR NR R  NR NR NR R  R  NR R  R  R  NR NR R  NR NR NR R 
##  [944] NR NR R  R  NR NR R  NR R  NR NR R  NR NR NR NR R  R  R  R  NR NR R 
##  [967] R  R  NR NR R  NR NR R  NR R  NR R  R  R  NR R  R  R  NR NR NR R  R 
##  [990] NR NR R  R  R  NR NR R  R  NR NR
## Levels: R NR

Visualizing the Toy Data

# Visualize the data labels and the separating hyperplane
ggplot(toy_data) + geom_point(aes(x=x1, y=x2, colour=label)) + geom_line(aes(x=x1, y=1-x1))

Density estimates of the Toy Data

# Look at the distribution of the first feature, broken down by label
ggplot(toy_data) + geom_density(aes(x=x1, fill=label), alpha=0.3)

# Look at the distribution of the second feature, broken down by label
ggplot(toy_data) + geom_density(aes(x=x2, fill=label), alpha=0.3) + coord_flip()

Finding the Cutpoint Using a Reasonable Approximation of the Distribution of Responders Given \(X_1\)

We want to find the value of \(X_1\) that gives us \(f(R | X_1) = 0.5\). One reasonable way to do that might be to just do a sliding window over \(X_1\) for all of our data, calculating the proportion of observations that are \(R\) vs. \(NR\) in each window. We keep track of the center of the bin and use that as our point estimate for the value of \(X_1\). More formally, we approximate the probability \(f(R | x_1)\) that a participant with a feature value \(X_1 = x_1\) will be a responder as follows:

\[ \begin{eqnarray} f(R | x_1) = \lim_{\Delta x \rightarrow \infty} \int_{x_1 - \Delta x}^{x_1 + \Delta x} f(R | x) dx &\approx \sum_{i=1}^{m} \frac{\mathbf{1}\left(\left(x_1 - \Delta x < x_1^{(i)} < x_1 + \Delta x\right) \cap \left(y^{(i)} = R\right)\right)}{\mathbf{1} \left(x_1 - \Delta x < x_1^{(i)} < x_1 + \Delta x\right)} \\ &= \frac{\text{# Responders in our window around } x_1}{\text{Total number of observations in window around } x_1} \end{eqnarray} \]

\[ \]

where \(\mathbf{1} \left( f(x^{(i)} \right) = 1\) if \(f(x^{(i)})\) is true and 0 otherwise.

We can see what the calculation/approximation looks like in practice as follows:

delta <- 0.01
bin_edges <- seq(0, 1, delta)  # Define a range for the left hand side of your sliding window
calib.dat <- data.frame(LHS=bin_edges, RHS=bin_edges + 20*delta)  # Define the window width for a number of windows
calib.dat$prop_R <- rep(0, length(bin_edges))  # Initialize your estimates of those proportions 
for(i in 1:(nrow(calib.dat))) {
  bin_labels <- toy_data$label[toy_data$x1 > calib.dat$LHS[i] & toy_data$x1 < calib.dat$RHS[i]]
  calib.dat[i, "x"] <- (calib.dat$LHS[i] + calib.dat$RHS[i])/2
  calib.dat[i, "prop_R"] <- mean(bin_labels == "R")
}

plot(calib.dat$x, calib.dat$prop_R)

Then we take a rolling average of our estimate for the marginal density of \(R\) given \(X_1\) (i.e. \(f(R | X_1)\)) to smooth our original findings and plot the result.

library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
my_zoo <- zoo(calib.dat[complete.cases(calib.dat),])
my_zoo <- rollapply(my_zoo, 5, mean)
plot(my_zoo$x, my_zoo$prop_R)

We choose our best guess to be \(x^*\), the point with an associated probability \(p(y^{(i)} = R | x - \Delta x < X_1 < x + \Delta x)\) that is closest to 0.5, i.e. \(x^*\) is the value of \(x\) that minimizes

\[ \left| \sum_{i=1}^{m} \frac{\mathbf{1}\left(\left(x - \Delta x < x_1^{(i)} < x + \Delta x\right) \cap \left(y^{(i)} = R\right)\right)}{\mathbf{1} \left(x - \Delta x < x_1^{(i)} < x + \Delta x\right)} - 0.5\right|\]

smoothed_props <- data.frame(my_zoo)
head(smoothed_props)
##    LHS  RHS    prop_R    x
## 3 0.02 0.22 0.8802190 0.12
## 4 0.03 0.23 0.8723186 0.13
## 5 0.04 0.24 0.8636487 0.14
## 6 0.05 0.25 0.8526662 0.15
## 7 0.06 0.26 0.8394023 0.16
## 8 0.07 0.27 0.8250218 0.17
smoothed_props$new_prop_R <- abs(smoothed_props$prop_R - 0.5)
head(smoothed_props)
##    LHS  RHS    prop_R    x new_prop_R
## 3 0.02 0.22 0.8802190 0.12  0.3802190
## 4 0.03 0.23 0.8723186 0.13  0.3723186
## 5 0.04 0.24 0.8636487 0.14  0.3636487
## 6 0.05 0.25 0.8526662 0.15  0.3526662
## 7 0.06 0.26 0.8394023 0.16  0.3394023
## 8 0.07 0.27 0.8250218 0.17  0.3250218
best_guess <- smoothed_props[which.min(smoothed_props$new_prop_R), "x"]

Just as a sanity check, let’s calculate what proportion of our observations fell between \(x^* - \Delta x\) and \(x^* + \Delta x\)…

mean(toy_data$label[toy_data$x1 > best_guess - 0.1 & toy_data$x1 < best_guess + 0.1] == "R")
## [1] 0.4951923

That looks about right! And if we look at the two density plots, we see that it about matches the point of intersection between the two estimated densities.

ggplot(toy_data) + geom_density(aes(x=x1, fill=label), alpha=0.3) + geom_vline(xintercept = best_guess)

Finding the Cutpoint by Directly Comparing Gaussian Kernel Density Estimates

Okay, but can we just find the actual point of intersection? It turns out we can. Check out https://stackoverflow.com/questions/25453706/how-to-find-the-intersection-of-two-densities-with-ggplot2-in-r for some additional insight. Basically, the geom_density function from ggplot2 uses the density function from R’s stats package (included in base R, you don’t have to install or load anything). If we use that function directly, we can find a more precise point of intersection.

lower.limit <- min(toy_data$x1)
upper.limit <- max(toy_data$x2)
responder.density <- density(toy_data[toy_data$label == "R", "x1"], from=lower.limit, to=upper.limit, n=1000)
nonresponder.density <- density(toy_data[toy_data$label == "NR", "x1"], from=lower.limit, to=upper.limit, n=1000)

# Take the difference between the densities for responders and nonresponders
density.difference <- responder.density$y - nonresponder.density$y

# Use as our best guess of intersection point the point associated with the
# minimum absolute valued difference between the two densities
intersection.point <- responder.density$x[which.min(abs(density.difference))] 

# And let's again verify that that gives us the point of intersection
ggplot(toy_data) + geom_density(aes(x=x1, fill=label), alpha=0.3) + geom_vline(xintercept = intersection.point)

Brilliant. Let’s use that as our cutpoint. We’ll say that an individual is a responder of $x_1^{(i)} < $ intersection.point.

Reference Scott (1992), Sheather et al (1991), and Silverman (1986) (or type ?density into the R console) to see how those kernel density estimates are created.

Finding the Cut Point Using Naive Bayes

Recall that we’re interested in finding a reasonable cutpoint of the data as a starting point for generating predictions based on each feature. If \(f(X_1)\) is the marginal distribution of our data along the first feature, then a reasonable cutpoint to try might be the point where \(f(R | X_1) = f(NR | X_1)\). As a heuristic, we would say that individual \(i\) is a Responder (i.e. \(y^{(i)} = R\)) if \(f(R | x_1^{(i)}) > f(NR | x_1^{(i)})\). From Bayes’ rule, we know that

\[f(R | X_1) = \frac{f(X_1 | R) f(R)}{f(X_1)}\]

and

\[f(NR | X_1) = \frac{f(X_1 | NR) f(NR)}{f(X_1)}\]

It follows that \(f(R | X_1) > f(NR | X_1)\) implies: \[\frac{f(X_1 | R) f(R)}{f(X_1)} > \frac{f(X_1 | NR) f(NR)}{f(X_1)}\] But notice that to solve this inequality we don’t actually have to explicitly calculate \(f(X_1)\), because it just cancels out. So our heuristic becomes

\[ f(X_1 | R) f(R) > f(X_1 | NR) f(NR)\]

This is the basis of Naive Bayes. \(f(R)\) and \(f(NR)\) are easy enough to calculate. We assume some sort of distribution for \(f(X_1 | R)\) and \(f(X_1 | NR)\) and then find where those two densities would overlap. In practice you could find that point of intersection by building a Naive Bayes classifier based on just a single feature, predicting R vs. NR for a series of points and find the point where the NB prediction switches from R to NR. I’ll implement that later, pending interest.