Three steps that "go beyond classical estimation:"
- Form a posterior distribution, represented by a set of simulated parameter values.
- Propogate uncertainty in this distribution/these parameters when making "simulation-based predictions for unobserved or future outcomes."
- Feed additional information into the mode using a prior distribution.
- "[S]tart by considering Bayesian simulation simply as a way of expressing uncertainty in inferences and predictions[...]."
- "The real advantage of summarizing inference by simulations is that we can
directly use these to propagate uncertainty." Especially for derived
quantities, they use
$a/b$ , intercept over slope, as an example. You just calculate that for all 1000 simulated$(a,b)$ pairs and report quantiles.
If you have the Bayesian regression fit (i.e., the simulation matrix of
- Collapse the simulations into a point estimate set of parameters, then use
that point estimate to generate a prediction from
$x$ ("point prediction"). "The point prediction ignores uncertainty." - Calculate predictions based on
$x$ for each of the simulated parameter sets, then report on the empirical distribution of those predictions ("linear predictor with uncertainty"). The specific uncertainty you're reporting is the "uncertainty about the expected or average value of$y$ ," emphasis mine. - Do what you did in (2), but now also rope in estimates of the error term(s).
They call this "the predictive distribution for a new observation." In our
case, this would be augmenting the uncertainty in the mean value of
$y$ with estimates of what we think$\sigma$ is.
These are actually very useful terms of art for me to keep in mind. I think the one to lead with every time is "predictive distribution" though.
The "by hand" implementation of the predictive distribution is cruder than I'd
thought it'd be: it looks like they're just turning each simulated
The "Propagating uncertainty" section is interesting. They feed a distinct
They provide some off-the-shelf formulas for predictive uncertainty
($\hat{\sigma}\text{linpred}$ and $\hat{sigma}\text{prediction}$), the
standard deviation of the point estimate for
Finally, when I look at their R code snippets: I have no idea where these functions are defined. You should care about stuff like that! Your code style should insist on making this easy to determine at the point of use!
They won't talk about Bayesian inference in general, but they do say we're going to have a compromise between the prior and the data we provide, when making our estimations.
Their formula for the parameter estimate, $\hat{\theta}\text{Bayes}$, might look nicer with a rewrite (multiply top and bottom by $\text{se}^2\text{prior}\text{se}^2_\text{data}$):
The standard error of the Bayes estimate goes to zero:
Typically,
It's important to highlight what they say about the standard error of the Bayes estimate: it must be less than the standard error you'd get for the data-only estimate. If the prior is proper, you are reducing uncertainty in your reported estimates. Maybe that's bad, if you pick a bad prior!
They note a key assumption about Bayesian information aggregation: the prior and the data are independent of each other. That's kind of implied by the name "prior," if you fix it prior to seeing the data, it's an independent source of data. Or it isn't definitely dependent. It may still be accidentally independent if whoever is providing the data is aware of, and responding to, the same prior info you're incorporating.
I take the closing subsection to mean there's no standard way of picking a prior distribution. Sometimes, like in the election and births examples, you use subject matter expertise. Sometimes, they allow that you just won't know enough about what you're estimating ("[multivariate regression] coefficients for which we have little knowledge or about which we do not want to make any strong assumptions").
They talk about a study where 3,000 parent-pairs were surveyed for their attractiveness and how many of their kids were girls.
This is obviously a sociological/bureaucratic process by which we declare
something is "well known," the way that "[i]t is well known that variation in
the human sex ratio occurs in a very narrow range." They have a set of
expectations for what a truly large change in girl-boy birth share would be, and
encode it as
They run a bunch of numbers on standard errors of proportions from Chapter 4,
on the whole 3,000. The data on "very attractive" parents
(
Estimating a comparison to an accuracy of 2 or 3 percentage points is pretty good -- for many purposes. But for studying differences in the sex ratio, we need much higher precision, more on the scale of 0.1 percentage point.
I don't think this quite sells their actual beef with the study. Yes, the Bayes estimate fails to reject the 50-50 null. That's cuz you set a prior to things "well known" to you! I guess the actual issues -- like, who set the rules for what counts as "very attractive" and did they set them so that they could reject the null -- are outside this particular book's remit. But as is, the section is light on arguments for why no, really, you need sub-0.1% precision for studying sex ratios.
They want you to know this could be a Type M error, but they can only compare the standard error to the prior that they insist upon: "the data standard error is more than 10 times the prior uncertainty." I think one more paragraph where they adjust the prior to once again admit a null-rejecting Bayes estimate would be interesting. They could point to how wide it is, and say, if sex ratios really can vary this much, how come we've never seen a (say) 66-33 split ever, anywhere, for any reason.
This section is about using stan_glm to set priors on regression coefficients.
By opening with the "uniform distribution" (actually an improper prior): it's
wild to me they don't even get into the API docs of stan_glm here. This
is the first time we see prior_intercept, prior, and prior_aux, and
there's not yet been a description of what those arguments expect. Integers?
Callbacks?
The default prior is norms for all coefficients. They word it weirdly -- they
make it sound like you need to know the
This part is worth going into:
If
$x$ and$y$ have both been standardized, then the coefficient is the expected difference in standard deviations in$y$ corresponding to a change of 1 standard deviation in$x$ . We typically expect such a difference to be less than 1 in absolute value, hence a normal prior with mean 0 and scale 2.5 will partially pool noisy coefficient estimates toward that range.
Why do we expect such a difference to be less than 1? Because we're doing statistics, and no one needs to do statistics to learn obvious effects. This loops back around to the sociology of where priors come from, where this is the biggest fact of all: all of this stuff is only useful for the narrow range of study between hopelessly noisy and blindingly obvious effects.
I'm still not loving their case against the beauty-sex-ratio study. "There's no steady, linear increase in girl-share as a function of attractiveness bucket" is different than "something's going on, there's lots of girls in the most attractive bucket." It's interesting that if you did think about the dataset in the first version ("girl-share is a linear function of parental hotness"), that the data rule that out. But that's a different thing! And the averages in each bucket are insufficient to check the second thing right now (but are sufficient to fit the regression, especially if you weight each bucket by corresponding number of parents).
Plots and computation powered by Chapter09.ipynb
A linear regression is fit on high school students modeling grade point average given household income. Write R code to compute the 90% predictive interval for the difference in grade point average comparing two students, one with household incomes of $40,000 and one with household income of $80,000.
Using Bambi as best I can right now:
rng = numpy.random.default_rng("Exercise 9.1")
data = pandas.DataFrame({"income": income, "gpa": gpa})
model = bambi.Model("gpa ~ income", data)
fitted = model.fit()
diffs = []
sigmas = []
for chain in range(4):
chain_df = (
fitted.posterior.sel(chain=chain).to_dataframe().groupby("draw").mean())
diffs.extend(inc * 40_000 for inc in chain_df["income"])
sigmas.extend(chain_df["sigma"].to_list())
noisy_diffs = [d + rng.normal(0, 2 * s) for d, s in zip(diffs, sigmas)]
print(numpy.percentiles(noisy_diffs, [5, 90]))Using data of interest to you, fit a linear regression. Use the output from this model to simulate a predictive distribution for observations with a particular combination of levels of all the predictors in the regression.
Bringing back my fantasy football dataset from Chapter 6, Exercise 6.7:
I fit a regression in Bambi to this, and evaluated the predictive posterior distribution for players at the 25th, 50th, and 75th percentile IDP level in 2023:
If I redo the histogram after transforming the sampled predictions to the linear scale, you can see that the upside of drafting a 75%-ile player is in that greater chance of a truly outstanding season:
Consider the economy and voting example from Section 7.1. Fit the linear regression model to the data through 2012; these are available in the folder
ElectionsEconomy. Make a forecast for the incumbent party’s share of the two-party vote in a future election where economic growth is 2%.(a) Compute the point forecast, the standard deviation of the predicted expectation from (9.1), and the standard deviation of the predicted value from (9.2).
(b) Now compute these using the relevant prediction functions discussed in Section 9.2. Check that you get the same values as in part (a) of this problem.
When I fit the (Bayesian) regression model to the data as vote ~ growth, I get
point estimates as:
That makes a point forecast of a 52.3% share of the vote if growth is 2%.
With
When I run this manually (and by sampling normal RVs according to the sampled
- Point forecast: 52.4%
- Standard deviation of point predictions: 1.0 percentage points
- Standard deviation of predictions (with noise/error added): 4.2 percentage points
Not bad!
Consider the example in Section 9.3 of combining prior information and survey data to forecast an election. Suppose your prior forecast for a given candidate’s share of the two-party vote is 42% with a forecast standard deviation of 5 percentage points, and you have survey data in which this candidate has 54% support. Assume that the survey is a simple random sample of voters with no nonresponse, and that voters will not change their mind between the time of the survey and the election.
(a) When you do the calculation, it turns out that your posterior mean estimate of the candidate’s vote share is 49%. What must the sample size of the survey be?
(b) Given this information, what is the posterior probability that the candidate wins the election?
Here's what we're given:
From Equation 9.3:
From Chapter 4:
We can derive:
Since we have a value for
So we can calculate the area under a stats.norm.sf(0.5, loc=0.49, scale=0.032). It's 38%.
A new job training program is being tested. Based on the successes and failures of previously proposed innovations, your prior distribution on the effect size on log(income) is normal with a mean of -0.02 and a standard deviation of 0.05. You then conduct an experiment which gives an unbiased estimate of the treatment effect of 0.16 with a standard deviation of 0.08. What is the posterior mean and standard deviation of the treatment effect?
Given:
We derive:
So our posterior can be modeled as a log-scale additive effect of
Perform the Bayesian analysis for the Jamaica experiment described on page 15. We shall work on the logarithmic scale:
(a) Do Exercise 4.8 to get the estimate and standard error of the log of the multiplicative treatment effect.
(b) Combine these with a normal prior distribution with mean 0 and standard deviation 0.10 to get the posterior distribution of the log of the multiplicative treatment effect.
(c) Exponentiate the mean of this distribution and check that it comes to 1.1, that is, an estimated treatment effect of +10%.
(d) Compute the 95% interval on the log scale and exponentiate it; check that this comes to [0.9, 1.3], that is, a range from -10% to +30% on the original scale.
From Exercise 4.8, the log-scale additive effect measured in
the data is distributed as
Combining this with the
If we use
Follow the steps of Section 9.5 for a different example, a regression of earnings on height using the data from the folder
Earnings. You will need to think what could be an informative prior distribution in this setting.
I transformed the Earnings data, creating a new variable by subtracting
sixty inches off each subject's height, so that the intercept term made more
sense.
For the strongly informative priors,
- Median earnings are 16k, and mean earnings are 21k. So using
$\mathcal{N}(20, 10)$ as the intercept prior will encompass the "no effect" outcome. - The quartiles (25, 50, 75%-ile) of earnings are spaced around $10k apart from
each other. And height quartiless are spaced around 2 inches apart from
each other. It would be very weird if height percentile entirely determined
earnings percentile, so
$5k/inch is an extreme overestimate. I'll go with a prior of $ \mathcal{N}(0, 0.3)$ to allow for at-most a quarter of that effect. - Like the book example in 9.5, I'll leave the prior on the noise level unspecified/flat.
Here are the inferences for the three parameters across the three choices of priors:
| Param | Flat | Weak | Strong |
|---|---|---|---|
| Intercept | |||
height |
|||
Huh! That strong prior on height didn't drastically affect the estimates.
I guess with
An experiment is performed to measure the efficacy of a television advertising program. The result is an estimate that each minute spent on a national advertising program will increase sales by $500 000, and this estimate has a standard error of $200,000. Assume the uncertainty in the treatment effect can be approximated by a normal distribution. Suppose ads cost $300,000 per minute. What is the expected net gain for purchasing 20 minutes of ads? What is the probability that the net gain is negative?
The gross benefit of one minute of advertising is distributed as
If you buy twenty minutes of ads, the net benefit is distributed as
(a) Expect a net gain of $4MM,
(b) but there's a 16% chance that the net gain winds up negative (that's the
share of a normal curve that's to the left of
Consider a regression predicting final exam score
$y$ from midterm exam score$x$ . Suppose that both exams are on the scale of 0–100, that typical scores range from 50–100, that the correlation of midterm and final is somewhere in the range 0.5–0.8, and that the average score on the final exam might be up to 20 points higher or lower than the average score on the midterm. Given the above information, set up reasonable priors for the slope and the intercept after centering.
I translate this as a prior of
Consider the model,
$y = a + bx + \text{error}$ , with predictors$x$ uniformly sampled from the range (-1, 1), independent prior distributions for the coefficients$a$ and$b$ , and the default prior distribution for the residual standard deviation$\sigma$ . For$a$ , assume a normal prior distribution with mean 0 and standard deviation 1; for$b$ , assume a normal prior distribution with mean 0 and standard deviation 0.2.(a) Simulate
$n = 100$ data points from this model, assuming the true parameter values are$a = 1, b = 0.1, \sigma = 0.5$ . Compute the least squares estimate of$a, b$ and compare to the Bayesian estimate obtained fromstan_glmusing the above priors.(b) Repeat the simulations with different values of
$n$ . Graph the Bayes estimates for$a$ and$b$ as a function of$n$ . What will these values be in the limit of$n = 0$ ?$n = \infty$ ?(c) For what value of
$n$ is the Bayes estimate for$a$ halfway between the prior mean and the least squares estimate? For what value of$n$ is the Bayes estimate for$b$ halfway between the prior mean and the least squares estimate?
With ordinary least squares, we recover the true parameters exactly, $\hat{a}\text{OLS} = 1, \hat{b}\text{OLS} = 0.1$.
When I sweep
| 2 | 0.96 | 0.12 |
| 4 | 1.05 | 0.16 |
| 10 | 1.22 | -0.13 |
| 20 | 0.80 | -0.05 |
| 35 | 0.94 | 0.05 |
| 50 | 1.10 | 0.18 |
| 67 | 1.05 | 0.01 |
| 83 | 0.94 | 0.16 |
| 100 | 1.14 | 0.28 |
| 150 | 1.02 | 0.12 |
| 200 | 0.99 | -0.04 |
| 300 | 1.04 | 0.13 |
| 400 | 1.00 | 0.07 |
If I greatly narrow the prior distribution on
| 2 | 0.51 | 0.02 |
| 4 | 0.89 | -0.00 |
| 10 | 1.00 | 0.00 |
| 20 | 0.73 | 0.00 |
| 35 | 1.08 | 0.00 |
| 50 | 0.88 | 0.01 |
| 67 | 0.93 | -0.00 |
| 83 | 1.00 | -0.00 |
| 100 | 0.95 | 0.00 |
| 150 | 0.97 | 0.01 |
| 200 | 1.02 | 0.01 |
| 300 | 0.97 | 0.02 |
| 400 | 1.01 | 0.01 |
And when I just use Bambi's default priors I get:
| 2 | -1.58 | -2.81 |
| 4 | 0.69 | -0.73 |
| 10 | 1.17 | 0.26 |
| 20 | 0.79 | 0.23 |
| 35 | 1.13 | 0.23 |
| 50 | 1.00 | 0.06 |
| 67 | 1.04 | 0.18 |
| 83 | 0.96 | 0.11 |
| 100 | 0.92 | 0.11 |
| 150 | 0.94 | 0.12 |
| 200 | 1.05 | 0.08 |
| 300 | 0.96 | 0.02 |
| 400 | 0.97 | 0.15 |
Maybe if I were extracting the sample median instead of the sample mean, this would be more stable, or produce the results this question is looking for.