Plotting the accuracy of cutpoints in a logistic regression model

Question

I would like to perform a simple logistic regression and draw a plot showing how the accuracy (that is "(Σ True positive + Σ True negative) / Σ Total population") change through different cutpoints, that is, the X axis should be the cutpoint for the value of the explanatory variable (between its minimum and maximum values) and the Y axis should be the accuracy (between 0 and 1). I use R and the package "Epi" does something sort of similar (see its function called "ROC") but I would like to do it on my own to easily tweak the plot to fit my study needs (I am not interested in the AUROC of the model).

For example. imagine my logistic regression is:

status <- c(0,0,0,0,1,0,1,1,1,1,0,0,0,0,1,0,1,1,1,1,0,0,0,0,1,0,1,1,1,1,0,0,0,0,1,0,1,1,1,1,0,0,0,0,1,0,1,1,1,1)
score <- c(3,3.03314007581356,3.06698316386694,3.10154799359025,3.1368538939095,3.17292081544931,3.20976935366332,3.24742077293537,3.28589703169672,3.32522080860666,3.36541552984653,3.40650539757956,3.44851541963181,3.49147144045215,3.53540017341248,3.58032923451229,3.62628717755537,3.67330353086979,3.72140883564605,3.77063468597247,3.82101377065092,3.87257991688045,3.92536813590135,3.97941467069665,4.03475704585397,4.09143411969572,4.14948613879184,4.20895479497545,4.26988328498839,4.33231637289088,4.39630045537665,4.46188363014311,4.52911576747436,4.59804858520385,4.6687357272328,4.7412328457907,4.81559768763467,4.8918901843959,4.97017254729336,5.05050936644761,5.13296771504133,5.21761725858725,5.3045303695798,5.39378224782294,5.48545104674383,5.57961800602068,5.67636759087249,5.77578763837942,5.87796951122469,5.98300825927263)
model <- glm(status~score,family=binomial)

Theoretically, how can I draw the plot?

Answer 1

I'm assuming that by cutpoint you mean the threshold for the probability/confidence score given by the logistic function. \begin{equation} \text{prediction} = \begin{cases} 1 & \text{if } p \ge \theta \\ 0 & \text{if } p \lt \theta \end{cases} \end{equation} Here $\theta$ is the threshold.

To plot the graph, get the probability estimates for all the samples within the data set. (If you are plotting this graph for the training data set, get the probability estimates for all the samples within the training data set).

If the data set has $m$ samples, you have the probability estimates as - $ p_{1}, p_{2} ..., p_{m}$.

What thresholds do I need to consider when plotting the graph ?
thresholds = unique($ \{ p_{1}, p_{2} ..., p_{m} \}$)
where the function unique($\cdot$) returns unique values for the input array.

Next, you iterate over thresholds in ascending order -> classify the samples using this threshold -> find accuracy to get a point on graph as - (threshold, accuracy). Following is the pseudocode -

    X <- data set

    prob_estimates = model.get_prob_estimates(X)
    thresholds <- unique(prob_estimates)
    points_on_graph = []

    for threshold in sort(thresholds)
        predictions = classify({p1, p2, ... , pm}, threshold)
        accuracy = score(predictions, true_labels)

    plot(points_on_graph)