Example when using accuracy as an outcome measure will lead to a wrong conclusion

Question

I am looking into various performance measures for predictive models. A lot was written about problems of using accuracy, instead of something more continuous to evaluate model performance. Frank Harrell provides an example when adding an informative variable to a model will lead to a drop in accuracy, clearly counterintuitive and wrong conclusion.

However, in this case, this seems to be caused by having imbalanced classes, and thus it can be solved just by using balanced accuracy instead ((sens+spec)/2). Is there some example where using accuracy on a balanced dataset will lead to some clearly wrong or counterintuitive conclusions?

Edit

I am looking for something where accuracy will drop even when the model is clearly better, or that using accuracy will lead to a false positive selection of some features. It's easy to make false-negative examples, where accuracy is the same for two models where one is clearly better using other criteria.

Answer 1

I'll cheat.

Specifically, I have argued often (e.g., here) that the statistical part of modeling and prediction extends only to making probabilistic predictions for class memberships (or giving predictive densities, in the case of numerical forecasting). Treating a specific instance as if it belonged to a specific class (or point predictions in the numerical case), is not properly statistics any more. It is part of the decision theoretic aspect.

And decisions should not only be predicated on the probabilistic prediction, but also on costs of misclassifications, and on a host of other possible actions. For instance, even if you have only two possible classes, "sick" vs. "healthy", you could have a large range of possible actions depending on how likely it is that a patient suffers from the disease, from sending him home because he is almost certainly healthy, to giving him two aspirin, to running additional tests, to immediately calling an ambulance and putting him on life support.

Assessing accuracy presupposes such a decision. Accuracy as a evaluation metric for classification is a category error.

So, to answer your question, I will walk down the path of just such a category error. We will consider a simple scenario with balanced classes where classifying without regard for the costs of misclassification will indeed mislead us badly.

---

Suppose an epidemic of Malignant Gutrot runs rampant in the population. Happily, we can screen everybody easily for some trait $t$ ($0\leq t \leq 1$), and we know that the probability of developing MG depends linearly on $t$, $p=\gamma t$ for some parameter $\gamma$ ($0\leq \gamma \leq 1$). The trait $t$ is uniformly distributed in the population.

Fortunately, there is a vaccine. Unfortunately, it is expensive, and the side effects are very uncomfortable. (I'll let your imagination supply the details.) However, they are better than to suffer from MG.

In the interest of abstraction, I posit that there are indeed only two possible courses of action for any given patient, given their trait value $t$: either vaccinate, or do not vaccinate.

Thus, the question is: how should we decide who to vaccinate and who not to, given $t$? We will be utilitarian about this and aim at having the lowest total expected costs. It is obvious that this comes down to choosing a threshold $\theta$ and to vaccinate everyone with $t\geq\theta$.

---

Model-and-decision 1 are accuracy-driven. Fit a model. Fortunately, we already know the model. Pick the threshold $\theta$ that maximizes accuracy when classifying patients, and vaccinate everyone with $t\geq \theta$. We easily see that $\theta=\frac{1}{2\gamma}$ is the magic number - everyone with $t\geq \theta$ has a higher chance of contracting MG than not, and vice versa, so this classification probability threshold will maximize accuracy. Assuming balanced classes, $\gamma=1$, we will vaccinate half the population. Funnily enough, if $\gamma<\frac{1}{2}$, we will vaccinate nobody. (We are mostly interested in balanced classes, so let's disregard that we just let part of the population die a Horrible Painful Death.)

Needless to say, this does not take the differential costs of misclassification into account.

---

Model-and-decision 2 leverage both our probabilistic prediction ("given your trait $t$, your probability of contracting MG is $\gamma t$") and the cost structure.

First, here is a little graph. The horizontal axis gives the trait, the vertical axis the MG probability. The shaded triangle gives the proportion of the population who will contract MG. The vertical line gives some particular $\theta$. The horizontal dashed line at $\gamma\theta$ will make the calculations below a bit simpler to follow. We assume $\gamma>\frac{1}{2}$, just to make life easier.

Let's give our costs names and calculate their contributions to total expected costs, given $\theta$ and $\gamma$ (and the fact that the trait is uniformly distributed in the population).

• Let $c^+_+$ denote the cost for a patient who is vaccinated and would have contracted MG. Given $\theta$, the proportion of the population who incurs this cost is the shaded trapezoid at the bottom right with area $$ (1-\theta)\gamma\theta + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta). $$
• Let $c^-_+$ denote the cost for a patient who is vaccinated and would not have contracted MG. Given $\theta$, the proportion of the population who incurs this cost is the unshaded trapezoid at the top right with area $$ (1-\theta)(1-\gamma) + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta). $$
• Let $c^-_-$ denote the cost for a patient who is not vaccinated and would not have contracted MG. Given $\theta$, the proportion of the population who incurs this cost is the unshaded trapezoid at the top left with area $$ \theta(1-\gamma\theta) + \frac{1}{2}\theta\gamma\theta. $$
• Let $c^+_-$ denote the cost for a patient who is not vaccinated and would have contracted MG. Given $\theta$, the proportion of the population who incurs this cost is the shaded triangle at the bottom left with area $$ \frac{1}{2}\theta\gamma\theta. $$

(In each trapezoid, I first calculate the area of the rectangle, then add the area of the triangle.)

Total expected costs are $$ c^+_+\bigg((1-\theta)\gamma\theta + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta)\bigg) + c^-_+\bigg((1-\theta)(1-\gamma) + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta)\bigg) + c^-_-\bigg(\theta(1-\gamma\theta) + \frac{1}{2}\theta\gamma\theta\bigg) + c^+_-\frac{1}{2}\theta\gamma\theta. $$

Differentiating and setting the derivative to zero, we obtain that expected costs are minimized by $$ \theta^\ast = \frac{c^-_+-c^-_-}{\gamma(c^+_-+c^-_+-c^+_+-c^-_-)}.$$

This is only equal to the accuracy maximizing value of $\theta$ for a very specific cost structure, namely if and only if $$ \frac{1}{2\gamma} = \frac{c^-_+-c^-_-}{\gamma(c^+_-+c^-_+-c^+_+-c^-_-)},$$ or $$ \frac{1}{2} = \frac{c^-_+-c^-_-}{c^+_-+c^-_+-c^+_+-c^-_-}.$$

As an example, suppose that $\gamma=1$ for balanced classes and that costs are $$ c^+_+ = 1, \quad c^-_+=2, \quad c^+_-=10, \quad c^-_-=0.$$ Then the accuracy maximizing $\theta=\frac{1}{2}$ will yield expected costs of $1.875$, whereas the cost minimizing $\theta=\frac{2}{11}$ will yield expected costs of $1.318$.

In this example, basing our decisions on non-probabilistic classifications that maximized accuracy led to more vaccinations and higher costs than using a decision rule that explicitly used the differential cost structures in the context of a probabilistic prediction.

---

Bottom line: accuracy is only a valid decision criterion if

• there is a one-to-one relationship between classes and possible actions
• and the costs of actions applied to classes follow a very specific structure.

In the general case, evaluating accuracy asks a wrong question, and maximizing accuracy is a so-called type III error: providing the correct answer to the wrong question.

---

R code:

rm(list=ls())
gamma <- 0.7

cost_treated_positive <- 1          # cost of treatment, side effects unimportant
cost_treated_negative <- 2          # cost of treatment, side effects unnecessary
cost_untreated_positive <- 10       # horrible, painful death
cost_untreated_negative <- 0        # nothing

expected_cost <- function ( theta ) {
    cost_treated_positive * ( (1-theta)*theta*gamma + (1-theta)*(gamma-gamma*theta)/2 ) +
    cost_treated_negative * ( (1-theta)*(1-gamma) + (1-theta)*(gamma-gamma*theta)/2 ) +
    cost_untreated_negative *( theta*(1-gamma*theta) + theta*gamma*theta/2 ) +
    cost_untreated_positive * theta*gamma*theta/2
}

(theta <- optim(par=0.5,fn=expected_cost,lower=0,upper=1,method="L-BFGS-B")$par)
(cost_treated_negative-cost_untreated_negative)/
    (gamma*(cost_treated_negative+cost_untreated_positive-cost_treated_positive-cost_untreated_negative))

plot(c(0,1),c(0,1),type="n",bty="n",xaxt="n",xlab="Trait t",yaxt="n",ylab="MG probability")
rect(0,0,1,1)
axis(1,c(0,theta,1),c(0,"theta",1),lty=0,line=-1)
axis(2,c(0,1),lty=0,line=-1,las=1)
axis(4,c(0,gamma,1),c(0,"gamma",1),lty=0,line=-1.8,las=1)
polygon(c(0,1,1),c(0,0,gamma),col="lightgray")
abline(v=theta,col="red",lwd=2)
abline(h=gamma*theta,lty=2,col="red",lwd=2)

expected_cost(1/(2*gamma))
expected_cost(theta)

Answer 2

It might worth adding another, perhaps more straightforward example to Stephen's excellent answer.

Let's consider a medical test, the result of which is normally distributed, both in sick and in healthy people, with different parameters of course (but for simplicity, let's assume homoscedasticity, i.e., that the variance is the same): $$\begin{gather*}T \mid D \ominus \sim \mathcal{N}\left(\mu_{-},\sigma^2\right)\\T \mid D \oplus \sim \mathcal{N}\left(\mu_{+},\sigma^2\right)\end{gather*}.$$ Let's denote the prevalence of the disease with $p$ (i.e. $D\oplus\sim Bern\left(p\right)$), so this, together with the above, which are essentially conditional distributions, fully specifies the joint distribution.

Thus the confusion matrix with threshold $b$ (i.e., those with test results above $b$ are classified as sick) is $$\begin{pmatrix} & D\oplus & D\ominus\\ T\oplus & p\left(1-\Phi_{+}\left(b\right)\right) & \left(1-p\right)\left(1-\Phi_{-}\left(b\right)\right)\\ T\ominus & p\Phi_{+}\left(b\right) & \left(1-p\right)\Phi_{-}\left(b\right)\\ \end{pmatrix}.$$

---

Accuracy-based approach

The accuracy is $$p\left(1-\Phi_{+}\left(b\right)\right)+\left(1-p\right)\Phi_{-}\left(b\right),$$

we take its derivative w.r.t. $b$, set it equal to 0, multiply with $\sqrt{1\pi\sigma^2}$ and rearrange a bit: $$\begin{gather*} -p\varphi_{+}\left(b\right)+\varphi_{-}\left(b\right)-p\varphi_{-}\left(b\right)=0\\ e^{-\frac{\left(b-\mu_{-}\right)^2}{2\sigma^2}}\left[\left(1-p\right)-pe^{-\frac{2b\left(\mu_{-}-\mu_{+}\right)+\left(\mu_{+}^2-\mu_{-}^2\right)}{2\sigma^2}}\right]=0\end{gather*}$$ The first term can't be zero, so the only way the product can be zero is if the second term is zero: $$\begin{gather*}\left(1-p\right)-pe^{-\frac{2b\left(\mu_{-}-\mu_{+}\right)+\left(\mu_{+}^2-\mu_{-}^2\right)}{2\sigma^2}}=0\\-\frac{2b\left(\mu_{-}-\mu_{+}\right)+\left(\mu_{+}^2-\mu_{-}^2\right)}{2\sigma^2}=\log\frac{1-p}{p}\\ 2b\left(\mu_{+}-\mu_{-}\right)+\left(\mu_{-}^2-\mu_{+}^2\right)=2\sigma^2\log\frac{1-p}{p}\\ \end{gather*}$$ So the solution is $$b^{\ast}=\frac{\left(\mu_{+}^2-\mu_{-}^2\right)+2\sigma^2\log\frac{1-p}{p}}{2\left(\mu_{+}-\mu_{-}\right)}=\frac{\mu_{+}+\mu_{-}}{2}+\frac{\sigma^2}{\mu_{+}-\mu_{-}}\log\frac{1-p}{p}.$$

Note that this - of course - doesn't depend on the costs.

If the classes are balanced, the optimum is the average of the mean test values in sick and healthy people, otherwise it is displaced based on the imbalance.

---

Cost-based approach

Using Stephen's notation, the expected overall cost is $$c_{+}^{+}p\left(1-\Phi_{+}\left(b\right)\right) + c_{+}^{-}\left(1-p\right)\left(1-\Phi_{-}\left(b\right)\right) + c_{-}^{+} p\Phi_{+}\left(b\right) + c_{-}^{-} \left(1-p\right)\Phi_{-}\left(b\right).$$ Take its derivate w.r.t $b$ and set it equal to zero: $$\begin{gather*} -c_{+}^{+} p \varphi_{+}\left(b\right)-c_{+}^{-}\left(1-p\right)\varphi_{-}\left(b\right)+c_{-}^{+}p\varphi_{+}\left(b\right)+c_{-}^{-}\left(1-p\right)\varphi_{-}\left(b\right)=\\ =\varphi_{+}\left(b\right)p\left(c_{-}^{+}-c_{+}^{+}\right)+\varphi_{-}\left(b\right)\left(1-p\right)\left(c_{-}^{-}-c_{+}^{-}\right)=\\ = \varphi_{+}\left(b\right)pc_d^{+}-\varphi_{-}\left(b\right)\left(1-p\right)c_d^{-}= 0,\end{gather*}$$ using the notation I introduced in my comments below Stephen's answer, i.e., $c_d^{+}=c_{-}^{+}-c_{+}^{+}$ and $c_d^{-}=c_{+}^{-}-c_{-}^{-}$.

The optimal threshold is therefore given by the solution of the equation $$\boxed{\frac{\varphi_{+}\left(b\right)}{\varphi_{-}\left(b\right)}=\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}}.$$ Two things should be noted here:

1. This results is totally generic and works for any distribution of the test results, not only normal. ($\varphi$ in that case of course means the probability density function of the distribution, not the normal density.)
2. Whatever the solution for $b$ is, it is surely a function of $\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}$. (I.e., we immediately see how costs matter - in addition to class imbalance!)

I'd be really interested to see if this equation has a generic solution for $b$ (parametrized by the $\varphi$s), but I would be surprised.

Nevertheless, we can work it out for normal! $\sqrt{2\pi\sigma^2}$s cancel on the left hand side, so we have $$\begin{gather*} e^{-\frac{1}{2}\left(\frac{\left(b-\mu_{+}\right)^2}{\sigma^2}-\frac{\left(b-\mu_{-}\right)^2}{\sigma^2}\right)}=\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}} \\ \left(b-\mu_{-}\right)^2-\left(b-\mu_{+}\right)^2 =2\sigma^2 \log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}} \\ 2b\left(\mu_{+}-\mu_{-}\right)+\left(\mu_{-}^2-\mu_{+}^2\right) =2\sigma^2 \log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}\end{gather*}$$ therefore the solution is $$b^{\ast}=\frac{\left(\mu_{+}^2-\mu_{-}^2\right)+2\sigma^2 \log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}}{2\left(\mu_{+}-\mu_{-}\right)}=\frac{\mu_{+}+\mu_{-}}{2}+\frac{\sigma^2}{\mu_{+}-\mu_{-}}\log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}.$$

(Compare it the the previous result! We see that they are equal if and only if $c_d^{-}=c_d^{+}$, i.e. the differences in misclassification cost compared to the cost of correct classification is the same in sick and healthy people.)

---

A short demonstration

Let's say $c_{-}^{-}=0$ (it is quite natural medically), and that $c_{+}^{+}=1$ (we can always obtain it by dividing the costs with $c_{+}^{+}$, i.e., by measuring every cost in $c_{+}^{+}$ units). Let's say that the prevalence is $p=0.2$. Also, let's say that $\mu_{-}=9.5$, $\mu_{+}=10.5$ and $\sigma=1$.

In this case:

library( data.table )
library( lattice )

cminusminus <- 0
cplusplus <- 1
p <- 0.2
muminus <- 9.5
muplus <- 10.5
sigma <- 1

res <- data.table( expand.grid( b = seq( 6, 17, 0.1 ),
                                cplusminus = c( 1, 5, 10, 50, 100 ),
                                cminusplus = c( 2, 5, 10, 50, 100 ) ) )
res$cost <- cplusplus*p*( 1-pnorm( res$b, muplus, sigma ) ) +
  res$cplusminus*(1-p)*(1-pnorm( res$b, muminus, sigma ) ) +
  res$cminusplus*p*pnorm( res$b, muplus, sigma ) +
  cminusminus*(1-p)*pnorm( res$b, muminus, sigma )

xyplot( cost ~ b | factor( cminusplus ), groups = cplusminus, ylim = c( -1, 22 ),
        data = res, type = "l", xlab = "Threshold",
        ylab = "Expected overall cost", as.table = TRUE,
        abline = list( v = (muplus+muminus)/2+
                         sigma^2/(muplus-muminus)*log((1-p)/p) ),
        strip = strip.custom( var.name = expression( {"c"^{"+"}}["-"] ),
                              strip.names = c( TRUE, TRUE ) ),
        auto.key = list( space = "right", points = FALSE, lines = TRUE,
                         title = expression( {"c"^{"-"}}["+"] ) ),
        panel = panel.superpose, panel.groups = function( x, y, col.line, ... ) {
          panel.xyplot( x, y, col.line = col.line, ... )
          panel.points( x[ which.min( y ) ], min( y ), pch = 19, col = col.line )
        } )

The result is (points depict the minimum cost, and the vertical line shows the optimal threshold with the accuracy-based approach):

We can very nicely see how cost-based optimum can be different than the accuracy-based optimum. It is instructive to think over why: if it is more costly to classify a sick people erroneously healthy than the other way around ($c_{-}^{+}$ is high, $c_{+}^{-}$ is low) than the threshold goes down, as we prefer to classify more easily into the category sick, on the other hand, if it is more costly to classify a healthy people erroneously sick than the other way around ($c_{-}^{+}$ is low, $c_{+}^{-}$ is high) than the threshold goes up, as we prefer to classify more easily into the category healthy. (Check these on the figure!)

---

A real-life example

Let's have a look at an empirical example, instead of a theoretical derivation. This example will be different basically from two aspects:

• Instead of assuming normality, we will simply use the empirical data without any such assumption.
• Instead of using one single test, and its results in its own units, we will use several tests (and combine them with a logistic regression). Threshold will be given to the final predicted probability. This is actually the preferred approach, see Chapter 19 - Diagnosis - in Frank Harrell's BBR.

The dataset (acath from the package Hmisc) is from the Duke University Cardiovascular Disease Databank, and contains whether the patient had significant coronary disease, as assessed by cardiac catheterization, this will be our gold standard, i.e., the true disease status, and the "test" will be the combination of the subject's age, sex, cholesterol level and duration of symptoms:

library( rms )
library( lattice )
library( latticeExtra )
library( data.table )

getHdata( "acath" )
acath <- acath[ !is.na( acath$choleste ), ]
dd <- datadist( acath )
options( datadist = "dd" )

fit <- lrm( sigdz ~ rcs( age )*sex + rcs( choleste ) + cad.dur, data = acath )

It worth plotting the predicted risks on logit-scale, to see how normal they are (essentially, that was what we assumed previously, with one single test!):

densityplot( ~predict( fit ), groups = acath$sigdz, plot.points = FALSE, ref = TRUE,
             auto.key = list( columns = 2 ) )

Well, they're hardly normal...

Let's go on and calculate the expected overall cost:

ExpectedOverallCost <- function( b, p, y, cplusminus, cminusplus,
                                 cplusplus = 1, cminusminus = 0 ) {
  sum( table( factor( p>b, levels = c( FALSE, TRUE ) ), y )*matrix(
    c( cminusminus, cplusminus, cminusplus, cplusplus ), nc = 2 ) )
}

table( predict( fit, type = "fitted" )>0.5, acath$sigdz )

ExpectedOverallCost( 0.5, predict( fit, type = "fitted" ), acath$sigdz, 2, 4 )

And let's plot it for all possible costs (a computational note: we don't need to mindlessly iterate through numbers from 0 to 1, we can perfectly reconstruct the curve by calculating it for all unique values of predicted probabilities):

ps <- sort( unique( c( 0, 1, predict( fit, type = "fitted" ) ) ) )

xyplot( sapply( ps, ExpectedOverallCost,
                p = predict( fit, type = "fitted" ), y = acath$sigdz,
                cplusminus = 2, cminusplus = 4 ) ~ ps, type = "l", xlab = "Threshold",
        ylab = "Expected overall cost", panel = function( x, y, ... ) {
          panel.xyplot( x, y, ... )
          panel.points( x[ which.min( y ) ], min( y ), pch = 19, cex = 1.1 )
          panel.text( x[ which.min( y ) ], min( y ), round( x[ which.min( y ) ], 3 ),
                      pos = 3 )
        } )

We can very well see where we should put the threshold to optimize the expected overall cost (without using sensitivity, specificity or predictive values anywhere!). This is the correct approach.

It is especially instructive to contrast these metrics:

ExpectedOverallCost2 <- function( b, p, y, cplusminus, cminusplus,
                                  cplusplus = 1, cminusminus = 0 ) {
  tab <- table( factor( p>b, levels = c( FALSE, TRUE ) ), y )
  sens <- tab[ 2, 2 ] / sum( tab[ , 2 ] )
  spec <- tab[ 1, 1 ] / sum( tab[ , 1 ] )
  c( `Expected overall cost` = sum( tab*matrix( c( cminusminus, cplusminus, cminusplus,
                                                   cplusplus ), nc = 2 ) ),
     Sensitivity = sens,
     Specificity = spec,
     PPV = tab[ 2, 2 ] / sum( tab[ 2, ] ),
     NPV = tab[ 1, 1 ] / sum( tab[ 1, ] ),
     Accuracy = 1 - ( tab[ 1, 1 ] + tab[ 2, 2 ] )/sum( tab ),
     Youden = 1 - ( sens + spec - 1 ),
     Topleft = ( 1-sens )^2 + ( 1-spec )^2
  )
}

ExpectedOverallCost2( 0.5, predict( fit, type = "fitted" ), acath$sigdz, 2, 4 )

res <- melt( data.table( ps, t( sapply( ps, ExpectedOverallCost2,
                                        p = predict( fit, type = "fitted" ),
                                        y = acath$sigdz,
                                        cplusminus = 2, cminusplus = 4 ) ) ),
             id.vars = "ps" )

p1 <- xyplot( value ~ ps, data = res, subset = variable=="Expected overall cost",
              type = "l", xlab = "Threshold", ylab = "Expected overall cost",
              panel=function( x, y, ... ) {
                panel.xyplot( x, y,  ... )
                panel.abline( v = x[ which.min( y ) ],
                              col = trellis.par.get()$plot.line$col )
                panel.points( x[ which.min( y ) ], min( y ), pch = 19 )
              }  )
p2 <- xyplot( value ~ ps, groups = variable,
              data = droplevels( res[ variable%in%c( "Expected overall cost",
                                                     "Sensitivity",
                                                     "Specificity", "PPV", "NPV" ) ] ),
              subset = variable%in%c( "Sensitivity", "Specificity", "PPV", "NPV" ),
              type = "l", xlab = "Threshold", ylab = "Sensitivity/Specificity/PPV/NPV",
              auto.key = list( columns = 3, points = FALSE, lines = TRUE ) )
doubleYScale( p1, p2, use.style = FALSE, add.ylab2 = TRUE )

We can now analyze those metrics that are sometimes specifically advertised as being able to come up with an optimal cutoff without costs, and contrast it with our cost-based approach! Let's use the three most often used metrics:

• Accuracy (maximize accuracy)
• Youden rule (maximize $Sens+Spec-1$)
• Topleft rule (minimize $\left(1-Sens\right)^2+\left(1-Spec\right)^2$)

(For simplicity, we will subtract the above values from 1 for the Youden and the Accuracy rule so that we have a minimization problem everywhere.)

Let's see the results:

p3 <- xyplot( value ~ ps, groups = variable,
              data = droplevels( res[ variable%in%c( "Expected overall cost", "Accuracy",
                                                     "Youden", "Topleft"  ) ] ),
              subset = variable%in%c( "Accuracy", "Youden", "Topleft"  ),
              type = "l", xlab = "Threshold", ylab = "Accuracy/Youden/Topleft",
              auto.key = list( columns = 3, points = FALSE, lines = TRUE ),
              panel = panel.superpose, panel.groups = function( x, y, col.line, ... ) {
                panel.xyplot( x, y, col.line = col.line, ... )
                panel.abline( v = x[ which.min( y ) ], col = col.line )
                panel.points( x[ which.min( y ) ], min( y ), pch = 19, col = col.line )
              } )
doubleYScale( p1, p3, use.style = FALSE, add.ylab2 = TRUE )

This of course pertains to one specific cost structure, $c_{-}^{-}=0$, $c_{+}^{+}=1$, $c_{+}^{-}=2$, $c_{-}^{+}=4$ (this obviously matters only for the optimal cost decision). To investigate the effect of cost structure, let's pick just the optimal threshold (instead of tracing the whole curve), but plot it as a function of costs. More specifically, as we have already seen, the optimal threshold depends on the four costs only through the $c_d^{-}/c_d^{+}$ ratio, so let's plot the optimal cutoff as a function of this, along with the typically used metrics that don't use costs:

res2 <- data.frame( rat = 10^( seq( log10( 0.02 ), log10( 50 ), length.out = 500 ) ) )
res2$OptThreshold <- sapply( res2$rat,
                             function( rat ) ps[ which.min(
                               sapply( ps, Vectorize( ExpectedOverallCost, "b" ),
                                       p = predict( fit, type = "fitted" ),
                                       y = acath$sigdz,
                                       cplusminus = rat,
                                       cminusplus = 1,
                                       cplusplus = 0 ) ) ] )

xyplot( OptThreshold ~ rat, data = res2, type = "l", ylim = c( -0.1, 1.1 ),
        xlab = expression( {"c"^{"-"}}["d"]/{"c"^{"+"}}["d"] ), ylab = "Optimal threshold",
        scales = list( x = list( log = 10, at = c( 0.02, 0.05, 0.1, 0.2, 0.5, 1,
                                                   2, 5, 10, 20, 50 ) ) ),
        panel = function( x, y, resin = res[ ,.( ps[ which.min( value ) ] ),
                                             .( variable ) ], ... ) {
          panel.xyplot( x, y, ... )
          panel.abline( h = resin[variable=="Youden"] )
          panel.text( log10( 0.02 ), resin[variable=="Youden"], "Y", pos = 3 )
          panel.abline( h = resin[variable=="Accuracy"] )
          panel.text( log10( 0.02 ), resin[variable=="Accuracy"], "A", pos = 3 )
          panel.abline( h = resin[variable=="Topleft"] )
          panel.text( log10( 0.02 ), resin[variable=="Topleft"], "TL", pos = 1 )
        } )

Horizontal lines indicate the approaches that don't use costs (and are therefore constant).

Again, we nicely see that as the additional cost of misclassification in the healthy group rises compared to that of the diseased group, the optimal threshold increases: if we really don't want healthy people to be classified as sick, we will use higher cutoff (and the other way around, of course!).

And, finally, we yet again see why those methods that don't use costs are not (and can't!) be always optimal.