Using a Fitted Model for Decision Analysis — Classification Competition

Summary

A short, self-contained case study that runs a Bayesian analysis all the way through to a decision: should you pay $10toentera$100,000 contest? The task is to classify 1000 length-135 time series as trendless, $+1°$C/century, or $-1°$C/century, with at least 900 correct required to win. A three-component mixture on the estimated slopes gives an expected 854 correct with sd 10.3 — the target is about 5 standard errors away, a 1-in-190,000 chance against the 1-in-10,000 needed for the bet to break even on expected value. The methodological point: the answer cannot be read off any parameter’s posterior.

Overview

The contest (Ch. 20.1, p. 317)

The organizer’s framing: “It has often been claimed that alarm about global warming is supported by observational evidence. I have argued that there is no observational evidence … A prize of $100,000 will be awarded to the first person who submits an entry that correctly identifies at least 900 series: which series were generated by a trendless process and which were generated by a trending process.”

The data generation, as described. 1000 series of length 135 (matching 1880-2014 global temperature records), generated “via a trendless statistical model fit for global temperatures,” after which “some randomly-selected series had a trend added to them. Some trends were positive; the others were negative. Each individual trend was $1°$C/century (in magnitude) — which is greater than the trend claimed for global temperatures.”

And the catch: “Each entry must be accompanied by a payment of $10.”

Main Content

Step 1 — data exploration and a deliberate simplification

Reducing each series to a slope (Figure 20.1-20.2, Ch. 20.2, pp. 317-318)

“Aha! The lines are fanning out from a common starting point. We’ll fit a regression to each line and then summarize each line by its average slope.”

for (n in 1:N) {
  fit <- lm(series[n,] ~ time)
  coefs <- summary(fit)$coefficients
  slope[n] <- 100 * coefs[2, "Estimate"]
  se[n]    <- 100 * coefs[2, "Std. Error"]
}

(Multiplied by 100 to put slopes on a per-century scale.)

The simplification acknowledged up front: “This is not necessarily the most efficient way to estimate a slope from correlated data, but we will not be concerned about that, as our purpose here is to demonstrate decision analysis in the context of a fitted model, rather than to fit the model in an optimal way.”

What the plots show: “There does not seem to be much information in the standard errors” (Figure 20.2a); the histogram of slopes (20.2b) shows what looks like three components — “Based on the problem description, we’d expect to see distributions centered at 0, $-1$, and 1. It looks like this might be the case.”

Step 2 — a three-component mixture with known means

The model, and what is fixed by the problem description (Ch. 20.3, p. 319)

data {
  int K; int N;
  array[N] real y;
  array[K] real mu;
}
parameters {
  simplex[K] theta;
  real<lower=0> sigma;
}
model {
  array[K] real ps;
  sigma ~ exponential(0.1);
  for (n in 1:N) {
    for (k in 1:K) {
      ps[k] = log(theta[k]) + normal_lpdf(y[n] | mu[k], sigma);
    }
    target += log_sum_exp(ps);
  }
}

• $\mu =(0,-1,1)$ supplied as data — “based on the problem description of there being a trendless process for which some series were added trends of $-1$ and $+1°$C/century.”
• A single shared $\sigma$ — “again based on the description of how the data were constructed.”
• simplex[K] theta “imposes the constraint that the mixture probabilities are all nonnegative and sum to 1.”

Fixing $\mu$ and sharing $\sigma$ is what makes this mixture trivially identified — the label-switching and unbounded-likelihood problems of Failure 4 — Label switching in mixture models simply do not arise. Exercise 20.3 asks the reader to relax both and watch it break.

The fit

 variable  mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
 theta[1]  0.54   0.54 0.02 0.02 0.50 0.57 1.00     2641     2668
 theta[2]  0.24   0.24 0.02 0.02 0.21 0.27 1.00     3466     2832
 theta[3]  0.22   0.22 0.02 0.02 0.20 0.25 1.00     2776     2454
 sigma     0.40   0.40 0.02 0.02 0.38 0.44 1.00     3028     2640

“The estimated weights of the three mixture components are approximately 0.5, 0.25, 0.25. Given that the problem was made up, we might guess that the weights of the underlying data-generation process were exactly 1/2, 1/4, and 1/4.”

A model expansion considered and declined: “We could also try fitting a model where the standard deviations of the three components differ, but we won’t, partly because the description given with the simulated data described the change as adding a trend, and partly because the histogram in Figure 20.2b doesn’t seem to show any varying of the widths.”

Step 3 — per-series membership probabilities

Computing classification probabilities in generated quantities

generated quantities {
  matrix[N,K] p;
  for (n in 1:N) {
    vector[K] p_raw;
    for (k in 1:K) {
      p_raw[k] = theta[k] * exp(normal_lpdf(y[n] | mu[k], sigma));
    }
    for (k in 1:K) {
      p[n,k] = p_raw[k] / sum(p_raw);
    }
  }
}

Then averaged over the posterior draws:

prob_sims <- posterior::as_draws_rvars(fit_mix$draws())
prob <- mean(prob_sims$p)

The first ten rows:

      [,1] [,2] [,3]
 [1,] 0.09 0.00 0.91
 [2,] 0.40 0.60 0.00
 [3,] 0.93 0.01 0.06
 ...

“the first series is probably drawn from the sloping-upward model; the second might be from the null model or it might be from the sloping-downward model; the third … through eighth are probably from the null model.”

Note that these are posterior means of the classification probabilities, averaged over uncertainty in $\theta$ and $\sigma$ — the “simulate first, summarize last” discipline of Simulation to Express Uncertainty applied to a derived quantity.

Step 4 — the decision, and its evaluation

Picking, counting, and computing the odds (Ch. 20.4, pp. 320-321)

max_prob <- apply(prob, 1, max)
choice   <- apply(prob, 1, which.max)
print(table(choice))

choice
  1   2   3
560 231 209

Why the guesses aren’t in proportion to the weights: “There seem to be too many guesses of zero slope and not enough of positive and negative slopes. But that makes sense given the decision problem: we want to maximize the number of correct guesses so we end up disproportionately guessing the most common category.”

This is the decision-theoretic point in miniature: the optimal action distribution is not the posterior distribution. Maximizing expected correct classifications is a different objective from reproducing the mixture weights.

expected_correct <- sum(max_prob)
sd_correct <- sqrt(sum(max_prob * (1 - max_prob)))

$$\text{expected correct}=854.0,\text{sd}=10.3$$

“An expectation of 854 does not reach the 900 that’s needed to win; indeed, with a standard error of 10.3, the target of 900 is about 5 standard errors away.”

The sd formula assumes “the reasonable approximation of independence of the 1000 series conditional on the model” — a Poisson-binomial variance, $\sum p_i(1-p_i)$.

Converting to odds, and comparing to the entry fee

1 / pnorm(expected_correct, 900, sd_correct)

$$\mathrm{Pr}(\ge 900)\approx 4.2\times 10^{-6}⟹\text{about 1 in }240,000$$

With a continuity correction (evaluating at 899.5): about 1 in 190,000.

“But a probability of 1 in 190,000 is still not enough. For the bet to be worth it, even in the crudest sense of expected monetary value, the probability of winning would have to be at least 1 in 10,000. Recall that the prize is $100,000butthecostofentryis$10.”

The bet loses by a factor of about 20.

Note the reporting choice: the raw probability $4.2\times 10^{-6}$ is “a number that is difficult to interpret, so we compute its reciprocal” — a small but real point about communicating extreme probabilities.

Two honest caveats about the analysis

1. The adversarial one. “That’s all conditional on the designer of the study doing everything exactly as he said, not playing with multiple seeds for the random number generator … He could well have first chosen a seed and generated the series, then performed something like the above analysis and checked that the most natural estimate gave only 850 correct or so, and in the very unlikely event that the natural estimate gave 900 or close to it, just re-running with a new seed.

We have no reason to think that the creator of this challenge did anything like that; our point here is only that, even if he did his simulation in a completely clean way, our odds of winning are about 1 in 200,000.”

2. The statistical one, which points at the remaining opportunity. “The data are highly autocorrelated, and least squares regression is not the most efficient way to estimate these slopes. If we can estimate the slopes more precisely, we can get more discrimination in our predictions. Maybe there is a way to win the game by extracting more information from each series, but it won’t be easy.”

Refusing to draw the designer's conclusion

“You could say that the above all demonstrates the designer’s point, that you can’t easily identify a trend in a time series of this length. But we don’t think it would make sense to draw that conclusion from this exercise.

After all, you can just tweak the parameters in the problem a bit — or simply set the goal to 800 correct instead of 900 — and the game becomes easy to win. Or, had the game been winnable as initially set up, you could just up the threshold to 950, and again it would become essentially impossible.

Conversely, if the designer had messed up his calculations and set the threshold to 800, and someone had sent in a winning entry, it wouldn’t disprove his claims about climate science; it would just mean that he hadn’t been careful enough in setting up his bet.”

A clean example of the reasoning in Statistical and Scientific Inference: the estimand must be connected to the scientific question, and here the contest’s threshold is an arbitrary parameter of the game, not a fact about climate.

Examples

General lessons (Ch. 20.5, pp. 321-322)

“We were able to perform a Bayesian analysis and pipe it all the way through to solve a decision problem. This is an example where the full posterior distribution contains important information that cannot be extracted from the posterior inferences for the parameters in the model.

It was necessary to keep the uncertainty throughout the analysis, which is done by performing all the steps separately on each posterior simulation draw and then only evaluating the decision at the end.”

Exercises 20.1-20.3 — three ways to improve the analysis (Ch. 20.6, p. 322)

• 20.1 — Integrate the two stages. “Write a Stan program to estimate the slopes and fit the mixture models all at once … For how many of the 1000 series does the best-estimated category change, and how much benefit is there in the expected number of correct classifications?”
• 20.2 — Model the autocorrelation, both in the two-stage and integrated forms, and compare “from parameter estimates through decision recommendations.”
• 20.3 — Break the mixture on purpose. Make $\mu$ and a vector $\sigma$ into parameters: “You should see some convergence problems.” Then fix it: “To avoid aliasing, the components of $\mu$ need to be constrained in some way, most simply by using the ordered type in Stan, and the relative values of $\sigma$ also need to be constrained, for example by giving them a lognormal prior, which will keep them away from zero.”

Exercise 20.3 is a compact, guided reproduction of both mixture pathologies from Failure Modes and Steps Forward: label switching (fixed by ordered) and the unbounded likelihood at $\sigma \to 0$ (fixed by a zero-avoiding lognormal).

Connections

• This is the book’s cleanest end-to-end instance of From Inference to Decision: a posterior, a loss (number correct), an action rule (argmax), and an expected-monetary-value comparison against a stated cost.
• The classification probabilities are computed in generated quantities and averaged over draws — the same post-processing pattern as LD50 in Bioassay - A First Probabilistic Program and SATE/PATE in Causal Inference as Generalization.
• The mixture here is well behaved only because the problem statement fixed $\mu$ and shared $\sigma$, which makes it a useful contrast case for Failure Modes and Steps Forward.

See Also

• From Inference to Decision — expected monetary value and when to use it
• Simulation to Express Uncertainty — propagating uncertainty to a derived decision
• Failure Modes and Steps Forward — mixture-model pathologies avoided here and reproduced in Exercise 20.3
• Decision Analysis — BDA3 background