Allen B. Downey's Blog: Probably Overthinking It, page 20

Internet use and religion, part two

In the previous article, I posted a preliminary exploration of the relationship between Internet use and religious affiliation in Europe.  In this article I clean up some data issues and present results broken by country.
Cleaning and resamplingHere are the steps of the data cleaning pipeline:

I replace sentinel values with NaNs.I recode some of the explanatory variables, and shift "year born" and "interview year" so their mean is 0.Within each data collection round and each country, I resample the respondents using their post stratification weights (pspwght). In Rounds 1 and 2, there are a few countries that were not asked about Internet use.  I remove these countries from those rounds.For the variables eduyrs and hinctnta, I replace each value with its rank (from 0-1) among respondents in the same round and country.I replace missing values with random samples from the same round and country.Finally, I merge rounds 1-5 into a single DataFrame and then group by country.This IPython notebook has the details, and summaries of the variables after processing.

One other difference compared to the previous notebook/article: I've added the variable invyr70, which is the year the respondent was interviewed (between 2002 and 2012), shifted by 2007 so the mean is near 0.

Aggregate resultsUsing variables with missing values filled (as indicated by the _f suffix), I get the following results from logistic regression on rlgblg_f (belonging to a religion), with sample size 233,856:

coefstd errzP>|z|[95.0% Conf. Int.]Intercept1.10960.01766.3810.0001.077 1.142inwyr07_f0.04770.00230.7820.0000.045 0.051yrbrn60_f-0.00900.000-32.6820.000-0.010 -0.008edurank_f0.01320.0170.7750.438-0.020 0.046hincrank_f0.07870.0164.9500.0000.048 0.110tvtot_f-0.01760.002-7.8880.000-0.022 -0.013rdtot_f-0.01300.002-7.7100.000-0.016 -0.010nwsptot_f-0.03560.004-9.9860.000-0.043 -0.029netuse_f-0.10820.002-61.5710.000-0.112 -0.105
Compared to the results from last time, there are a few changesInterview year has a substantial effect, but probably should not be taken too seriously in this model, since the set of countries included in each round varies.  The apparent effect of time might reflect the changing mix of countries.  I expect this variable to be more useful after we group by country.The effect of "year born" is similar to what we saw before.  Younger people are less likely to be affiliated.The effect of education, now expressed in relative terms within each country, is no longer statistically significant.  The apparent effect we saw before might have been due to variation across countries.The effect of income, now expressed in relative terms within each country, is now positive, which is more consistent with results in other studies.  Again, the apparent negative effect in the previous analysis might have been due to variation across countries (see Simpson's paradox).The effect of the media variables is similar to what we saw before:  Internet use has the strongest effect, 2-3 times bigger than newspapers, which are 2-3 times bigger than television or radio.  And all are negative.The inconsistent behavior of education and income as control variables is a minor concern, but I think the symptoms are most likely the result of combining countries, possibly made worse because I am not weighting countries by population, so smaller countries are overrepresented.
Here are the results from linear regression with rlgdgr_f (degree of religiosity) as the dependent variable:
coefstd errtP>|t|[95.0% Conf. Int.]Intercept6.01400.022270.4070.0005.970 6.058inwyr07_f0.02530.00212.0890.0000.021 0.029yrbrn60_f-0.01720.000-46.1210.000-0.018 -0.016edurank_f-0.24290.023-10.5450.000-0.288 -0.198hincrank_f-0.15410.022-7.1280.000-0.196 -0.112tvtot_f-0.07340.003-24.3990.000-0.079 -0.067rdtot_f-0.01990.002-8.7600.000-0.024 -0.015nwsptot_f-0.07620.005-15.6730.000-0.086 -0.067netuse_f-0.13740.002-57.5570.000-0.142 -0.133In this model, all parameters are statistically significant.  The effect of the media variables, including Internet use, is similar to what we saw before.
The effect of education and income is negative in this model, but I am not inclined to take it too seriously, again because we are combining countries in a way that doesn't mean much.
Breakdown by countryThe following table shows results for logistic regression, with rlgblg_f as the dependent variable, broken down by country; the columns are country code, number of observations, and the estimated parameter associated with Internet use:
Country Num      Coef of     code    obs.     netuse_f------- ----     --------AT 6918 -0.0795 ** BE 8939 -0.0299 ** BG 6064 0.0145 CH 9310 -0.0668 ** CY 3293 -0.229 ** CZ 8790 -0.0364 ** DE 11568 -0.0195 * DK 7684 -0.0406 ** EE 6960 -0.0205 ES 9729 -0.0741 ** FI 7969 -0.0228 FR 5787 -0.0185 GB 11117 -0.0262 ** GR 9759 -0.0245 HR 3133 -0.0375 HU 7806 -0.0175 IE 10472 -0.0276 * IL 7283 -0.0636 ** IS 579 0.0333 IT 1207 -0.107 ** LT 1677 -0.0576 * LU 3187 -0.0789 ** LV 1980 -0.00623 NL 9741 -0.0589 ** NO 8643 -0.0304 ** PL 8917 -0.108 ** PT 10302 -0.103 ** RO 2146 0.00855 RU 7544 0.00437 SE 9201 -0.0374 ** SI 7126 -0.0336 ** SK 6944 -0.0635 ** TR 4272 -0.0857 * UA 7809 -0.0422 **
** p < 0.01, * p < 0.05
In the majority of countries, there is a statistically significant relationship between Internet use and religious affiliation.  In all of those countries the relationship is negative, with the magnitude of most coefficients between 0.03 and 0.11 (with one exceptionally large value in Cyprus).
Degree of religiosityAnd here are the results of linear regression, with rlgdgr_f as the dependent variable:
Country Num      Coef of     code    obs.     netuse_f------- ----     --------AT 6918 -0.0151     ** BE 8939 -0.0072     ** BG 6064  0.0023     CH 9310 -0.0132     ** CY 3293 -0.00221     ** CZ 8790 -0.005     ** DE 11568 -0.0045     * DK 7684 -0.00909     ** EE 6960 -0.00363     ES 9729 -0.0165     ** FI 7969 -0.00501     FR 5787 -0.00429     GB 11117 -0.0061     ** GR 9759 -0.00362     ** HR 3133 -0.00559     HU 7806 -0.00478     * IE 10472 -0.00412     * IL 7283 -0.00419     ** IS 579  0.00752     IT 1207 -0.0212     ** LT 1677 -0.00746     LU 3187 -0.0152     ** LV 1980 -0.00147     NL 9741 -0.014     ** NO 8643 -0.00721     ** PL 8917 -0.00919     ** PT 10302 -0.0149     ** RO 2146  0.000303     RU 7544  0.00102     SE 9201 -0.00815     ** SI 7126 -0.00835     ** SK 6944 -0.0119     ** TR 4272 -0.00331     ** UA 7809 -0.00947     **
In most countries there is a negative and statistically significant relationship between Internet use and degree of religiosity.
In this model the effect of education is consistent: in most countries it is negative and statistically significant.  In the two countries where it is positive, it is not statistically significant.
The effect of income is less consistent: in most countries it is not statistically significant; when it is, it is positive as often as negative.
But education and income are in the model primarily as control variables; they are not the focus on this study.  If they are actually associated with religious affiliation, these variables should be effective controls; if not, they contribute some noise, but otherwise do no harm.
Next stepsFor now I am using StatsModels to estimate parameters and compute confidence intervals, but that's not quite right because I am using resampled data and filling missing values with random samples.  To account correctly for these sources of random error, I have to run the whole process repeatedly:Resample the data.Fill missing values.Estimate parameters.Collecting the estimated parameters from multiple runs, I can estimate the sampling distribution of the parameters and compute confidence intervals.
Once I have implemented that, I plan to translate the results into a form that is easier to interpret (rather than just estimated coefficients), and generate visualizations to make the results easier to explore.
I would also like to relate the effect of Internet use in each country with the average level of religiosity, to see whether, for example, the effect is bigger in more religious countries.
While I am working on that, I am open to suggestions for additional explorations people might be interested in.  You can explore the variables in the ESS using their "Cumulative Data Wizard"; let me know what you find!

Internet use and religious affiliation in Europe

A few years ago I wrote a paper about Internet use and religious affiliation using data from the General Social Survey (GSS).  After controlling for things like education, income, and religious upbringing, I found that people who use the Internet more are less likely to have a religious affiliation; that is, they are more likely to be "Nones".  I made an argument for why I think this effect is plausibly causative and, if so, I estimated that it accounts for about 20% of the decrease in religious affiliation in the U.S. between 1990 and 2010.  Last April it got some media attention (including this brief discussion on Fox News); here is the blog post I wrote about it at the time.

Since then, I have been planning to replicate the study using data from the European Social Survey, but I didn't get back to it until a few weeks ago.  I was reminded about it because they recently released data from Round 7, conducted in 2014.  I am always excited about new data, but in this case it doesn't help me: in Rounds 6 and 7 they did not ask about Internet use.

But they did ask in Rounds 1 through 5, collected between 2002 and 2010.  So that's something.  I downloaded the data, loaded it into Pandas dataframes, and started cleaning, validating, and doing some preliminary exploration.  This IPython notebook shows what the first pass looks like.

In the interest of open, replicable science, I am posting preliminary work here, but at this point we should not take the results too seriously.

Data inventoryThe dependent variables I plan to study are

rlgblg: Do you consider yourself as belonging to any particular religion or denomination?

rlgdgr: Regardless of whether you belong to a particular religion, how religious would you say you are?

The explanatory variables are

yrbrn: And in what year were you born?

hinctnta: Using this card, please tell me which letter describes your household's total income, after tax and compulsory deductions, from all sources? If you don't know the exact figure, please give an estimate. Use the part of the card that you know best: weekly, monthly or annual income.

eduyrs: About how many years of education have you completed, whether full-time or part-time? Please report these in full-time equivalents and include compulsory years of schooling.

tvtot: On an average weekday, how much time, in total, do you spend watching television?

rdtot: On an average weekday, how much time, in total, do you spend listening to the radio?

nwsptot: On an average weekday, how much time, in total, do you spend reading the newspapers?

netuse: Now, using this card, how often do you use the internet, the World Wide Web or e-mail - whether at home or at work - for your personal use?

RecodesIncome: I created a variable, hinctnta5, which subtracts 5 from hinctnta, so the mean is near 0.  This shift makes the parameters of the model easier to interpret.

Year born: Similarly, I created yrbrn60, which subtracts 1960 from yrbrn.

Years of education: The distribution of eduyrs includes some large values that might be errors, and the question was posed differently in the first few rounds.  I will investigate more carefully later, but for now I am replacing values greater than 25 years with 25, and subtracting off the mean, 12, to create eduyrs12.


ResultsJust to get a quick look at things, I ran a logistic regression with rlgblg as the dependent variable, using data Rounds 1-5 and including all countries.  The sample size is 229,307.  Here are the estimated parameters (computed by StatsModels):
coefstd errzP>|z|[95.0% Conf. Int.]Intercept0.98110.01470.0140.0000.954 1.009yrbrn60-0.00780.000-27.8260.000-0.008 -0.007eduyrs12-0.03760.001-29.6190.000-0.040 -0.035hinctnta5-0.02200.002-12.9340.000-0.025 -0.019tvtot-0.01610.002-7.2050.000-0.021 -0.012rdtot-0.01490.002-8.8260.000-0.018 -0.012nwsptot-0.03200.004-8.9240.000-0.039 -0.025netuse-0.07580.002-42.0620.000-0.079 -0.072
The parameters are all statistically significant with very small p-values.  And they are all negative, which indicates:

Younger people are less likely to report a religious affiliation.More educated people are less likely...People with higher income are less likely...People who consume more media (television, radio, newspaper) are less likely...People who use the Internet more are less likely...The effect of Internet use (per hour per week) appears to be about twice the effect of reading the newspaper, which is about twice the effect of television or radio.
The effect of the Internet is comparable to about a decade of age, two years of education, or 3 deciles of income.

Most of these results are consistent with what I saw in my previous study and what other studies have shown.  One exception is income: in other studies, the usual pattern is that people in the lowest income groups are less likely to be affiliated, and after that, income has no effect.  We'll see if this preliminary result holds up.

I ran a similar model using rlgdgr (degree of religiosity) as the dependent variable:
coefstd errtP>|t|[95.0% Conf. Int.]Intercept5.66680.019300.0710.0005.630 5.704yrbrn60-0.01800.000-47.3520.000-0.019 -0.017eduyrs12-0.06880.002-40.1590.000-0.072 -0.065hinctnta5-0.02660.002-11.3970.000-0.031 -0.022tvtot-0.08010.003-26.3340.000-0.086 -0.074rdtot-0.01790.002-7.7910.000-0.022 -0.013nwsptot-0.05310.005-10.8730.000-0.063 -0.044netuse-0.10200.002-40.9420.000-0.107 -0.097
The results are similar.  Again, this IPython notebook has the details.

LimitationsAgain, we should not take these results too seriously yet:

So far I am not taking into account the weights associated with respondents, either within or across countries.  So for now I am oversampling people in small countries, as well as some groups within countries.At this point I haven't done anything careful to fill missing values, so the results will change when I get that done.And I think it will be more meaningful to break the results down by country.Stay tuned.  More coming soon!

Recidivism and logistic regression

In my previous article, I presented the problem of estimating a criminal's risk of recidivism, focusing on the philosophical problems of attributing a probability to an individual.  In this article I turn to the more practical problem of doing the estimation.

My collaborator asked me to compare two methods and explain their assumptions, properties, pros and cons:

1) Logistic regression: This is the standard method in the field.  Researchers have designed a survey instrument that assigns each offender a score from -3 to 12.  Using data from previous offenders, they estimate the parameters of a logistic regression model and generate a predicted probability of recidivism for each score.

2) Post test probabilities: An alternative of interest to my collaborator is the idea of treating survey scores the way some tests are used in medicine, with a threshold that distinguishes positive and negative results.  Then for people who "test positive" we can compute the post-test probability of recidivism.

In addition to these two approaches, I considered three other models:

3) Independent estimates for each group: In the previous article, I cited the notion that "the probability of recidivism for an individual offender will be the same as the observed recidivism rate for the group to which he most closely belongs." (Harris and Hanson 2004).  If we take this idea to an extreme, the estimate for each individual should be based only on observations of the same risk group.  The logistic regression model is holistic in the sense that observations from each group influence the estimates for all groups.  An individual who scores a 6, for example, might reasonably object if the inclusion of offenders from other risk classes has the effect of increasing the estimated scores for his own class.  To evaluate this concern, I also considered a simple model where the risk in each group is estimated independently.

4) Logistic regression with a quadratic term: Simple logistic regression is based on the assumption that the scores are linear in the sense that each increment corresponds to the same increase in risk; for example, it assumes that the odds ratio between groups 1 and 2 is the same as the odds ratio between groups 9 and 10.  To test that assumption, I ran a model that includes the square of the scores as a predictive variable.  If risks are actually non-linear, this quadratic term might provide a better fit to the data.  But it doesn't, so the linear assumption holds up, at least to this simple test.

5) Bayesian logistic regression: This is a natural extension of logistic regression that takes as inputs prior distributions for the parameters, and generates posterior distributions for the parameters and predictive distributions for the risks in each group.  For this application, I didn't expect this approach to provide any great advantage over conventional logistic regression, but I think the Bayesian version is a powerful and surprisingly simple idea, so I used this project as an excuse to explore it.

Details of these methods are in this IPython notebook; I summarize the results below.

Logistic regressionMy collaborator provided a dataset with 10-year recidivism outcomes for 703 offenders.  The definition of recidivism in this case is that the individual was charged with another offense within 10 years of release (but not necessarily convicted).  The range of scores for offenders in the dataset is from -2 to 11.

I estimated parameters using the StatsModels module.  The following figure shows the estimated probability of recidivism for each score and the predictive 95% confidence interval.



A striking feature of this graph is the range of risks, from less than 10% to more than 60%.  This range suggests that it is important to get this right, both for the individuals involved and for society.

[One note: the logistic regression model is linear when expressed in log-odds, so it is not a straight line when expressed in probabilities.]

Post test probabilitiesIn the context of medical testing, a "post-test probability/positive" (PTPP) is the answer to the question, "Given that a patient has tested positive for a condition of interest (COI), what is the probability that the patient actually has the condition, as opposed to a false positive test".  The advantage of PTPP in medicine is that it addresses the question that is arguably most relevant to the patient.

In the context of recidivism, if we want to treat risk scores as a binary test, we have to choose a threshold that splits the range into low and high risk groups.  As an example, if I choose "6 or more" as the threshold, I can compute the test's sensitivity, specificity, and PTPP:
With threshold 6, the sensitivity of the test is 42%, meaning that the test successfully identifies (or predicts) 42% of recidivists.The specificity is 76%, which is the fraction of non-recidivists who correctly test negative.The PTPP is 42%, meaning that 42% of the people who test positive will reoffend.  [Note: it is only a coincidence in this case that sensitivity and PTPP are approximately equal.]One of the questions my collaborator asked is whether it's possible to extend this analysis to generate a posterior distribution for PTPP, rather than a point estimate.  In the notebook I show how to do that using beta distributions for the priors and posteriors.  For this example, the 95% posterior credible interval is 35% to 49%.

Now suppose we customize the test for each individual.  So if someone scores 3, we would set the threshold to 3 and compute a PTPP and a credible interval.  We could treat the result as the probability of recidivism for people who test 3 or higher.  The following figure shows the results (in blue) superimposed on the results from logistic regression (in grey).


A few features are immediately apparent:
At the low end of the range, PTPP is substantially higher than the risk estimated by logistic regression.At the high end, PTPP is lower.With this dataset, it's not possible to compute PTPP for the lowest and highest scores (as contrasted with logistic regression, which can extrapolate).For the lowest risk scores, the credible intervals on PTPP are a little smaller than the confidence intervals of logistic regression.For the high risk scores, the credible intervals are very large, due to the small amount of data.So there are practical challenges in estimating PTPPs with a small dataset.  There are also philosophical challenges.  If the goal is to estimate the risk of recidivism for an individual, it's not clear that PTPP is the answer to the right question.
This way of using PTPP has the effect of putting each individual in a reference class with everyone who scored the same or higher.  Someone who is subject to this test could reasonably object to being lumped with higher-risk offenders.  And they might reasonably ask why it's not the other way around: why is the reference class "my risk or higher" and not "my risk or lower"?

I'm not sure there is a principled reason to that question.  In general, we should prefer a smaller reference class, and one that does not systematically bias the estimates.  Using PTPP (as proposed) doesn't do well by these criteria.

But logistic regression is vulnerable to the same objection: outcomes from high-risk offenders have some effect on the estimates for lower-risk groups.  Is that fair?

Independent estimatesTo answer that question, I also ran a simple model that estimates the risk for each score independently.  I used a uniform distribution as a prior for each group, computed a Bayesian posterior distribution, and plotted the estimated risks and credible intervals.  Here are the results, again superimposed on the logistic regression model.

People who score 0, 1, or 8 or more would like this model.  People who score -2, 2, or 7 would hate it.  But given the size of the credible intervals, it's clear that we shouldn't take these variations too seriously.

With a bigger dataset, the credible intervals would be smaller, and we expect the results to be ordered so that groups with higher scores have higher rates of recidivism.  But with this small dataset, we have some practical problems.

We could smooth the data by combining adjacent scores into risk ranges.  But then we're faced again with the objection that this kind of lumping is unfair to people at the low end of each range, and dangerously generous to people at the high end.

I think logistic regression balances these conflicting criteria reasonably: by taking advantage of background information -- the assumption that higher scores correspond to higher risks -- it uses data efficiently, generating sensible estimates for all scores and confidence intervals that represent the precision of the estimates.

However, it is based on an assumption that might be too strong, that the relationship between log-odds and risk score is linear.  If that's not true, the model might be unfair to some groups and too generous to others.

Linear regression with a quadratic termTo see whether the linear assumption is valid, we can run logistic regression again with two explanatory variables: score and the square of score.  If the relationship is actually non-linear, this quadratic model might fit the data better.

The estimated coefficient for the quadratic term is negative, so the line of best fit is a parabola with downward curvature.  That's consistent with the shape of the independent estimates in the previous section.  But the estimate for the quadratic parameter is not statistically significant (p=0.11), which suggests that it might be due to randomness.

Here is what the results look like graphically, again compared to the logistic regression model:


The estimates are lower for the lowest score groups, about the same in the middle, and lower again for the highest scores.  At the high end, it seems unlikely that the risk of recidivism actually declines for the highest scores, and more likely that we are overfitting the data.

This model suggests that the linear logistic regression might overestimate risk for the lowest-scoring groups.  It might be worth exploring this possibility with more data.

Bayesian logistic regressionOne of the best features of logistic regression is that the results are in the form of a probability, which is useful for all kinds of analysis and especially for the kind of risk-benefit analysis that underlies parole decisions.

But conventional logistic regression does not admit prior information about the parameters, and the result is not a posterior distribution, but just a point estimate and confidence interval.  In general, confidence intervals are less useful as part of a decision making process.  And if you have ever tried to explain a confidence interval to a general audience, you know they can be problematic.

So, mostly for my own interest, I applied Bayesian logistic regression to the same data.  The Bayesian version of logistic regression is conceptually simple and easy to implement, in part because at the core of logistic regression there is a tight connection to Bayes's theorem.  I wrote about this connection in a previous article.

My implementation, with explanatory comments, is in the IPython notebook; here are the results:


With a uniform prior for the parameters, the results are almost identical to the results from (conventional) logistic regression.  At first glance it seems like the Bayesian version has no advantage, but there are a few reasons it might:

For each score, we can generate a posterior distribution that represents not just a point estimate, but all possible estimates and their probabilities.  If the posterior mean is 55%, for example, we might want to know the probability that the correct value is more than 65%, or less than 45%.  The posterior distribution can answer this question; conventional logistic regression cannot.If we use estimates from the model as part of another analysis, we could carry the posterior distribution through the subsequent computation, producing a posterior predictive distribution for whatever results the analysis generates.If we have background information that guides the choice of the prior distributions, we could use that information to produce better estimates.But in this example, I don't have a compelling reason to choose a different prior, and there is no obvious need for a posterior distribution: a point estimate is sufficient.  So (although it hurts me to admit it) the Bayesian version of logistic regression has no practical advantage for this application.
SummaryUsing logistic regression to estimate the risk of recidivism for each score group makes good use of the data and background information.  Although the estimate for each group is influenced by outcomes in other groups, this influence is not unfair in any obvious way.
The alternatives I considered are problematic:  PTPP is vulnerable to both practical and philosophical objections.  Estimating risks for each group independently is impractical with a small dataset, but might work well with more data.
The quadratic logistic model provides no apparent benefit, and probably overfits the data.  The Bayesian version of logistic regression is consistent with the conventional version, and has no obvious advantages for this application.

Recidivism and single-case probabilities

I am collaborating with a criminologist who studies recidivism.  In the context of crime statistics, a recidivist is a convicted criminal who commits a new crime post conviction.  Statistical studies of recidivism are used in parole hearings to assess the risks of releasing a prisoner.  This application of statistics raises questions that go to the foundations of probability theory.

Specifically, assessing risk for an individual is an example of a "single-case probability", a well-known problem in the philosophy of probability.  For example, we would like to be able to make a statement like, "If released, there is a 55% probability that Mr. Smith will be charged with another crime within 10 years."  But how could we estimate a probability like that, and what would it mean?

I suggest we answer these questions in three steps.  The first is to choose a "reference class"; the second is to estimate the probability of recidivism in the reference class; the third is to interpret the result as it applies to Mr. Smith.

For example, if the reference class includes all people convicted of the same crime as Mr. Smith, we could find a sample of that population and compute the rate of recidivism in the sample.  If 55% of the sample were recidivists, we might claim that Mr. Smith has a 55% probability of reoffending.

Let's look at those steps in more detail:

1) The reference class problem  For any individual offender, there are an unbounded number of reference classes we might choose.  For example, if we consider characteristics like age, marital status, and severity of crime, we might form a reference class using any one of those characteristics, any two, or all three.  As the number of characteristics increases, the number of possible classes grows exponentially (and I mean that literally, not as a sloppy metaphor for "very fast").

Some reference classes are preferable to others; in general, we would like the people in each class to be as similar as possible to each other, and the classes to be as different as possible from each other.  But there is no objective procedure for choosing the "right" reference class.

2) Sampling and estimation  Assuming we have chosen a reference class, we would like to estimate the proportion of recidivists in the class.  First, we have to define recidivism in a way that's measurable.  Ideally we would like to know whether each member of the class commits another crime, but that's not possible.  Instead, recidivism is usually defined to mean that the person is either charged or convicted of another crime.

If we could enumerate the members of class and count recidivists and non-, we would know the true proportion, but normally we can't do that.  Instead we choose a sampling process intended to select a representative sample, meaning that every member of the class has the same probability of appearing in the sample, and then use the observed proportion in the sample to estimate the true proportion in the class.

3) Interpretation  Suppose we agree on a reference class, C, a sampling process, and an estimation process, and estimate that the true proportion of recidivists in C is 55%.  What can we say about Mr. Smith?

As my collaborator has demonstrated, this question is a topic of ongoing debate.  Among practitioners in the field, some take the position that "the probability of recidivism for an individual offender will be the same as the observed recidivism rate for the group to which he most closely belongs." (Harris and Hanson 2004).  On this view, we would conclude that Mr. Smith has a 55% chance of reoffending.

Others take the position that this claim is nonsense because it could never be confirmed or denied; whether Mr. Smith goes on to commit another crime or not, neither outcome supports or refutes the claim that the probability was 55%.  On this view, probability cannot be meaningfully applied to a single event.

To understand this view, consider an analogy suggested by my colleague Rehana Patel:  Suppose you estimate that the median height of people in class C is six feet.  You could not meaningfully say that the median height of Mr. Smith is six feet.  Only the class has a median height, individuals do not.  Similarly, only the class has a proportion of recidivists; individuals do not.

And the answer is...At this point I have framed the problem and tried to state the opposing views clearly.  Now I will explain my own view and try to justify it.

I think it is both meaningful and useful to say, in the example, that Mr. Smith has a 55% chance of offending again.  Contrary to the view that no possible outcome supports or refutes this claim, Bayes's theorem tells us otherwise.  Suppose we consider two hypotheses:

H55: Mr. Smith has a 55% chance of reoffending.
H45: Mr. Smith has a 45% chance of reoffending.

If Mr. Smith does in fact commit another crime, this datum supports H55 with a Bayes factor of (55)/(45) = 1.22.  And if he does not, that datum supports H45 by the same factor.  In either case it is weak evidence, but nevertheless it is evidence, which means that H55 is a meaningful claim that can be supported or refuted by data.  The same argument applies if there are more than two discrete hypotheses or a continuous set of hypotheses.

Furthermore, there is a natural operational interpretation of the claim.  If we consider some number, n, of individuals from class C, and estimate that each of them has probability, p, of reoffending, we can compute a predictive distribution for k, the number who reoffend.  It is just the binomial distribution of k with parameters p and n:



For example, if n=100 and p=55, the most likely value of k is 55, and the probability that k=55 is about 8%.  As far as I know, no one has a problem with that conclusion.

But if there is no philosophical problem with the claim, "Of these 100 members of class C, the probability that 55 of them will reoffend is 8%", there should be no special problem when n=1.  In that case we would say, "Of these 1 members of class C, the probability that 1 will reoffend is 55%."  Granted, that's a funny way to say it, but that's a grammatical problem, not a philosophical one.

Now let me address the challenge of the height analogy.  I agree that Mr. Smith does not have a median height; only groups have medians.  However, if the median height in C is six feet, I think it is correct to say that Mr. Smith has a 50% chance of being taller than six feet.

That might sound strange; you might reply, "Mr. Smith's height is a deterministic physical quantity.  If we measure it repeatedly, the result is either taller than six feet every time, or shorter every time.  It is not a random quantity, so you can't talk about its probability."

I think that reply is mistaken, because it leaves out a crucial step:

1) If we choose a random member of C and measure his height, the result is a random quantity, and the probability that it exceeds six feet is 50%.

2) We can consider Mr. Smith to be a randomly-chosen member of C, because part of the notion of a reference class is that we consider the members to be interchangeable.

3) Therefore Mr. Smith's height is a random quantity and we can make probabilistic claims about it.

My use of the word "consider" in the second step is meant to highlight that this step is a modeling decision: if we choose to regard Mr. Smith as a random selection from C, we can treat his characteristics as random variables.  The decision not to distinguish among the members of the class is part of what it means to choose a reference class.

Finally, if the proportion of recidivists in C is 55%, the probability that a random member of C will reoffend is 55%.  If we consider Mr. Smith to be a random member of C, his probability of reoffending is 55%.

Is this useful?I have argued that it is meaningful to claim that Mr. Smith has a 55% probability of recidivism, and addressed some of the challenges to that position.

I also think that claims like this are useful because they guide decision making under uncertainty.  For example, if Mr. Smith is being considered for parole, we have to balance the costs and harms of keeping him in prison with the possible costs and harms of releasing him.  It is probably not possible to quantify all of the factors that should be taken into consideration in this decision, but it seems clear that the probability of recidivism is an important factor.

Furthermore, this probability is most useful if expressed in absolute, as opposed to relative, terms.  If we know that one prisoner has a higher risk than another, that provides almost no guidance about whether either should be released.  But if one has a probability of 10% and another 55% (and those are realistic numbers for some crimes) those values could be used as part of a risk-benefit analysis, which would usefully inform individual decisions, and systematically yield better outcomes.


What about the Bayesians and the Frequentists?When I started this article, I thought the central issue was going to be the difference between the Frequentist and Bayesian interpretations of probability.  But I came to the conclusion that this distinction is mostly irrelevant.

Considering the three steps of the process again:

1) Reference class: Choosing a reference class is equally problematic under either interpretation of probability; there is no objective process that identifies the right, or best, reference class.  The choice is subjective, but that is not to say arbitrary.  There are reasons to prefer one model over another, but because there are multiple relevant criteria, there is no unique best choice.

2) Estimation: The estimation step can be done using frequentist or Bayesian methods, and there are reasons to prefer one or the other.  Some people argue that Bayesian methods are more subjective because they depend on a prior distribution, but I think both approaches are equally subjective; the only difference is whether the priors are explicit.  Regardless, the method you use to estimate the proportion of recidivists in the reference class has no bearing on the third step.

3) Interpretation: In my defense of the claim that the proportion of recidivists in the reference class is the probability of recidivism for the individual, I used Bayes's theorem, which is a universally accepted law of probability, but I did not invoke any particular interpretation of probability.

I argued that we could treat an unknown quantity as a random variable.  That idea is entirely unproblematic under Bayesianism, but less obviously compatible with frequentism.  Some sources claim that frequentism is specifically incompatible with single-case probabilities.

For example the Wikipedia article on probability interpretations claims that "Frequentists are unable to take this approach [a propensity interpretation], since relative frequencies do not exist for single tosses of a coin, but only for large ensembles or collectives."

I don't agree that assigning a probability to a single case is a special problem for frequentism.  Single case probabilities seem hard because they make the choice of the reference class more difficult.  But choosing a reference class is hard under any interpretation of probability; it is not a special problem for frequentism.

And once you have chosen a reference class, you can estimate parameters of the class and generate predictions for individuals, or groups, without commitment to a particular interpretation of probability.

[For more about alternative interpretations of probability, and how they handle single-case probabilities, see this article in the Stanford Encyclopedia of Philosophy, especially the section on frequency interpretations.  As I read it, the author agrees with me that (1) the problem of the single case relates to choosing a reference class, not attributing probabilities to individuals, and (2) in choosing a reference class, there is no special problem with the single case that is not also a problem in other cases.  If there is any difference, it is one of degree, not kind.]


Summary

In summary, the assignment of a probability to an individual depends on three subjective choices:

1) The reference class,
2) The prior used for estimation,
3) The modeling decision to regard an individual as a random sample with n=1.

You can think of (3) as an additional choice or as part of the definition of the reference class.

These choices are subjective but not arbitrary; that is, there are justified reasons to prefer one over others, and to accept some as good enough for practical purposes.

Finally, subject to those choices, it is meaningful and useful to make claims like "Mr. Smith has a 55% probability of recidivism".


Q&A

1) Isn't it dehumanizing to treat an individual as if he were an interchangeable, identical member of a reference class?  Every human is unique; you can't treat a person like a number!

I conjecture that everyone who applies quantitative methods to human endeavors has heard a complaint like this.  If the intended point is, "You should never apply quantitative models to humans," I disagree.  Everything we know about the world, including the people in it, is mediated by models, and all models are based on simplifications.  You have to decide what to include and what to leave out.  If your model includes the factors that matter and omits the ones that don't, it will be useful for practical purposes.  If you make bad modeling decisions, it won't.

But if the intent of this question is to say, "Think carefully about your modeling decisions, validate your models as much as practically possible, and consider the consequences of getting it wrong," then I completely agree.

Model selection has consequences.  If we fail to include a relevant factor (that is, one that has predictive power), we will treat some individuals unfairly and systematically make decisions that are less good for society.  If we include factors that are not predictive, our predictions will be more random and possibly less fair.

And those are not the only criteria.  For example, if it turns out that a factor, like race, has predictive power, we might choose to exclude it from the model anyway, possibly decreasing accuracy in order to avoid a particularly unacceptable kind of injustice.

So yes, we should think carefully about model selection, but no, we should not exclude humans from the domain of statistical modeling.

2) Are you saying that everyone in a reference class has the same probability of recidivism?  That can't be true; there must be variation within the group.

I'm saying that an individual only has a probability AS a member of a reference class (or, for the philosophers, qua a member of a reference class).  You can't choose a reference class, estimate the proportion for the class, and then speculate about different probabilities for different members of the class.  If you do, you are just re-opening the reference class question.

To make that concrete, suppose there is a standard model of recidivism that includes factors A, B, and C, and excludes factors D, E, and F.  And suppose that the model estimates that Mr. Smith's probability of recidivism is 55%.

You might be tempted to think that 55% is the average probability in the group, and the probability for Mr. Smith might be higher or lower.  And you might be tempted to adjust the estimate for Mr. Smith based on factor D, E, or F.

But if you do that, you are effectively replacing the standard model with a new model that includes additional factors.  In effect, you are saying that you think the standard model leaves out an important factor, and would be better if it included more factors.

That might be true, but it is a question of model selection, and should be resolved by considering model selection criteria.

It is not meaningful or useful to talk about differences in probability among members of a reference class.  Part of the definition of reference class is the commitment to treat the members as equivalent.

That commitment is a modeling decision, not a claim about the world.  In other words, when we choose a model, we are not saying that we think the model captures everything about the world.  On the contrary, we are explicitly acknowledging that it does not.  But the claim (or sometimes the hope) is that it captures enough about the world to be useful.

3) What happened to the problem of the single case?  Is it a special problem for frequentism?  Is it a special problem at all?

There seems to be some confusion about whether the problem of the single case relates to choosing a reference class (my step 1) or attributing a probability to an individual (step 3).

I have argued that it does not relate to step 3.  Once you have chosen a reference class and estimated its parameters, there is no special problem in applying the result to the case of n=1, not under frequentism or any other interpretation I am aware of.

During step 1, single case predictions might be more difficult, because the choice of reference class is less obvious or people might be less likely to agree.  But there are exceptions of both kinds: for some single cases, there might be an easy consensus on an obvious good reference class; for some non-single cases, there are many choices and no consensus.  In all cases, the choice of reference class is subjective, but guided by the criteria of model choice.

So I think the single case problem is just an instance of the reference class problem, and not a special problem at all.

Learning to Love Bayesian Statistics

At Strata NYC 2015, O'Reilly Media's data science conference, I gave a talk called "Learning to Love Bayesian Statistics".  The video is available now:


The sound quality is not great, and the clicker to advance the slides was a little wonky, but other than that, it went pretty well.

Here are the slides, if you want to follow along at home.


I borrowed illustrations from The Phantom Tollbooth, in a way I think it consistent with fair use.  I think they work remarkably well.

Thanks to the folks at Strata for inviting me to present, and for allowing me to make the video freely available.  It's actually the first video I have posted on YouTube.  I'm a late adopter.

One million is a lot

When I was in third grade, the principal of my elementary school announced a bottle cap drive with the goal of collecting one million bottle caps.  The point, as I recall, was to demonstrate that one million is a very large number.  After a few months, we ran out of storage space, the drive was cancelled, and we had to settle for the lesson that 100,000 is a lot of bottle caps.

So it is a special pleasure for me to announce that, early Sunday morning (November 1, 2015), this blog reached one million page views.  I am celebrating the occasion with a review of some of my favorite articles and, of course, some analysis of the page view statistics.

Here's a screenshot of my Blogger stats page to make it official:


And here are links to the 10 most read articles:
PostsEntryPageviewsAre first babies more likely to be late?Feb 7, 2011, 9 comments130446







All your Bayes are belong to us!Oct 27, 2011, 56 comments47773







My favorite Bayes's Theorem problemsOct 20, 2011, 13 comments33020







Bayesian survival analysis for "Game of Thrones"Mar 25, 2015, 5 comments32210







The Inspection Paradox is EverywhereAug 18, 2015, 23 comments30330







Bayesian statistics made simpleMar 14, 201221718







Are your data normal? Hint: no.Aug 7, 2013, 13 comments15035







Yet another reason SAT scores are non-predictiveFeb 2, 2011, 3 comments13806







Freshman hordes even more godless!Jan 29, 2012, 6 comments7904







Bayesian analysis of match rates on TinderFeb 10, 20157096
By far the most popular is my article about whether first babies are more likely to be late.  It turns out they are, but only by a couple of days.

Two of the top 10 are articles written by students in my Bayesian statistics class: "Bayesian survival analysis for Game of Thrones" by Erin Pierce and Ben Kahle, and "Bayesian analysis of match rates on Tinder", by Ankur Das and Mason del Rosario.  So congratulations, and thanks, to them!

Five of the top 10 are explicitly Bayesian, which is clearly the intersection of my interests and popular curiosity.  But the other common theme is the application of statistical methods (of any kind) to questions people are interested in.



According to Blogger stats, my readers are mostly in the U.S., with most of the rest in Europe.  No surprises there, with the exception of Ukraine, which is higher in the rankings than expected.  Some of those views are probably bogus; anecdotally, Blogger does not do a great job of filtering robots and fake clicks (I don't have ads on my blog, so I am not sure how anyone benefits from fake clicks, but I have to conclude that some of my million are not genuine readers).


Most of my traffic comes from Google, Reddit, Twitter, and Green Tea Press, which is the home of my free books.  It looks like a lot of people find me through "organic" search, as opposed to my attempts at publicity.  And what are those people looking for?


People who find my blog are looking for Bayesian statistics, apparently, and the answer to the eternal question, "Are first babies more likely to be late?"

Those are all the reports you can get from Blogger (unless you are interested in which browsers my readers use).  But if I let it go at that, this blog wouldn't be called "Probably Overthinking It".

I used the Chrome extension SingleFile to grab the stats for each article in a form I can process, then used the Pandas read_html function to get it all into a table.  The results, and my analysis, are in this IPython notebook.

My first post, "Proofiness and Elections", was on January 4, 2011.  I've published 115 posts since then; the average time between posts is 15 days, but that includes a 180 day hiatus after "Secularization in America: Part Seven" in July 2012.  I spent the fall working on Think Bayes, and got back to blogging in January 2013.

Blogger provides stats on the most popular posts; I had to do my own analysis to extract the least popular posts:


Some of these deserve their obscurity, but not all.  "Will Millennials Ever Get Married?" is one of my favorite projects, and I think the video from the talk is pretty good.  And "When will I win the Great Bear Run?" is one of the best statistical models I've developed, albeit applied to a problem that is mostly silly.

Measures of popularity often follow Zipf's law, and my blog is no different.   As I suggest in Chapter 5 of Think Complexity, the most robust way to check for Zipf-like behavior is to plot the complementary CDF of frequency (for example, page views) on a log-log scale:

[image error]

For articles with more than 1000 page views, the CCDF is approximately straight, in compliance with Zipf's law.

The posts that elicited the most comments are:


Apparently, people like their veridical paradoxes!  The A few of my posts have attracted attention on the social network of Google employees, Google+:

I'm glad someone appreciates The Inspection Paradox.  I submitted it for publication in CHANCE magazine, but they didn't want it.  Thirty thousand readers, 909 Google employees, and I think they blew it.

One thing I have learned from this blog is that I can never predict whether an article will be popular.  One of the most technically challenging articles, "Bayes meets Fourier", apparently found an audience of people interested in Bayesian statistics and signal processing.  At the same time, some of my favorites, like The Rock Hyrax Problem and Belly Button Biodiversity have landed flat.  I've given up trying to predict what will hit.

I have posts coming up in the next few weeks that I am excited about, including an analysis of Internet use and religion using data from the European Social Survey.  Watch this space.

Thanks to everyone who contributed to the first million page views.  I hope you found it interesting and learned something, and I hope you'll be back for the next million!

When will I win the Great Bear Run?

Almost every year since 2008 I have participated in the Great Bear Run, a 5K road race in Needham MA. I usually finish in the top 40 or so, and in my age group I have come in 4th, 6th, 4th, 3rd, 2nd, 4th and 4th. In 2015 I didn't run because of a scheduling conflict, but based on the results I estimate that I would have come 4th again.

Having come close in 2012, I have to wonder what my chances are of winning my age group.  I've developed two models to predict the results.

Binomial modelAccording to a simple binomial model (the details are in this IPython notebook), the predictive distribution for the number of people who finish ahead of me is:

[image error]

And my chances of winning my age group are about 5%.

An improved modelThe binomial model ignores an important piece of information: some of the people who have displaced me from my rightful victory have come back year after year.  Here are the people who have finished ahead of me in my age group:

    2008: Gardiner, McNatt, Terry
    2009: McNatt, Ryan, Partridge, Turner, Demers
    2010: Gardiner, Barrett, Partridge
    2011: Barrett, Partridge
    2012: Sagar
    2013: Hammer, Wang, Hahn
    2014: Partridge, Hughes, Smith
    2015: Barrett, Sagar, Fernandez

Several of them are repeat interlopers.  I have developed a model that takes this information into account and predicts my chances of winning.  It is based on the assumption that there is a population of n runners who could displace me, each with some probability p.

In order to displace me, a runner has to

Show up,Outrun me, andBe in my age group.
For each runner, the probability of displacing me is a product of these factors:

pi=SOB

Some runners have a higher SOB factor than others; we can use previous results to estimate it.

First we have to think about an appropriate prior.   Again, the details are in this IPython notebook.

Based on my experience, I conjecture that the prior distribution of S is an increasing function, with many people who run nearly every year, and fewer who run only occasionally.The prior distribution of O is biased toward high values. Of the people who have the potential to beat me, many of them will beat me every time. I am only competitive with a few of them.  (For example, of the 15 people who have beat me, I have only ever beat 2).The probability that a runner is in my age group depends on the difference between his age and mine. Someone exactly my age will always be in my age group. Someone 4 years older will be in my age group only once every 5 years (the Great Bear run uses 5-year age groups).  So the distribution of B is uniform.I used Beta distributions for each of the three factors, so each piis the product of three Beta-distributed variates. In general, the result is not a Beta distribution, but it turns out there is a Beta distribution that is a good approximation of the actual distribution.
[image error]
Using this distribution as a prior, we can compute the posterior distribution of p for each runner who has displaced me, and another posterior for any hypothetical runner who has not displaced me in any of 8 races.  As an example, for Rich Partridge, who has displaced me in 4 of 8 years, the mean of the posterior distribution is 42%.  For a runner who has only displaced me once, it is 17%.
Then, for any hypothetical number of runners, n, we can draw samples from these distributions of p and compute the conditional distribution of k, the number of runners who finish ahead of me.  Here's what that looks like for a few values of n:

[image error]
If there are 18 runners, the mostly likely value of k is 3, so I would come in 4th.  As the number of runners increases, my prospects look a little worse.
These represent distributions of k conditioned on n, so we can use them to compute the likelihood of the observed values of k.  Then, using Bayes theorem, we can compute the posterior distribution of n, which looks like this:

[image error]
It's noisy because I used random sampling to estimate the conditional distributions of k. But that's ok because we don't really care about n; we care about the predictive distribution of k. And noise in the distribution of n has very little effect on k.
The predictive distribution of k is a weighted mixture of the conditional distributions we already computed, and it looks like this:
[image error]
Sadly, according to this model, my chance of winning my age group is less than 2% (compared to the binomial model, which predicts that my chances are more than 5%).
And one more time, the details are in this IPython notebook.

Bayes meets Fourier

Joseph Fourier never met Thomas Bayes—Fourier was born in 1768, seven years after Bayes died.  But recently I have been exploring connections between the Bayes filter and the Fourier transform.

By "Bayes filter", I don't mean spam filtering using a Bayesian classifier, but rather recursive Bayesian estimation, which is used in robotics and other domains to estimate the state of a system that evolves over time, for example, the position of a moving robot.  My interest in this topic was prompted by Roger Labbe's book, Kalman and Bayesian Filters in Python, which I am reading with my book club.

The Python code for this article is in this IPython notebook.  If you get sick of the math, you might want to jump into the code.

The Bayes filterA Bayes filter starts with a distribution that represents probabilistic beliefs about the initial position of the robot.  It proceeds by alternately executing predict and update steps:

The predict step uses a physical model of the system and its current state to predict the state in the future.  For example, given the position and velocity of the robot, you could predict its location at the next time step.The update step uses a measurement and a model of measurement error to update beliefs about the state of the system.  For example, if we use GPS to measure the location of the robot, the measurement would include some error, and we could characterize the distribution of errors.
The update step uses Bayes theorem, so computationally it involves multiplying the prior distribution by the likelihood distribution and then renormalizing.

The predict step involves adding random variables, like position and velocity, so computationally it involves the convolution of two distributions.

The Convolution TheoremThe Convolution Theorem suggests a symmetry between the predict and update steps, which leads to an efficient algorithm.  I present the Convolution Theorem in Chapter 8 of Think DSP.  For discrete signals, it says:

DFT(f ∗ g) = DFT(f) · DFT(g)
That is, if you want to compute the discrete Fourier transform (DFT) of a convolution, you can compute the DFT of the signals separately and then multiply the results.  In words, convolution in the time domain corresponds to multiplication in the frequency domain.  And conversely, multiplication in the time domain corresponds to convolution in the frequency domain.

In the context of Bayes filters, we are working with spatial distributions rather than temporal signals, but the convolution theorem applies: instead of the "time domain", we have the "space domain", and instead of the DFT, we have the characteristic function.
The characteristic functionThis section shows that the characteristic function and the Fourier transform are the same thing.  If you already knew that, you can skip to the next section.
If you took mathematical statistics, you might have seen this definition of the characteristic function:

Where X is a random variable, t is a transform variable that corresponds to frequency (temporal, spatial, or otherwise), and E is the expectation operator.
If you took signals and systems, you saw this definition of the Fourier transform:

where f is a function,  is its Fourier transform, and ξ is a transform variable that corresponds to a frequency.
The two definitions look different, but if you write out the expectation operator, you get 

where p(x) is the probability density function of X, and P(t) is its Fourier transform.  The only difference between the characteristic function and the Fourier transform is the sign of the exponent, which is just a convention choice.  (The bar above P(t) indicates the complex conjugate, which is there because of the sign of the exponent.)
And that brings us to the whole point of the characteristic function: to add two random variables, you have to convolve their distributions.  By the Convolution Theorem, convolution in the space domain corresponds to multiplication in the frequency domain, so you can add random variables like this:Compute the characteristic function of both distributions.Multiply the characteristic functions elementwise.Apply the inverse transform to compute the distribution of the sum.Back to BayesNow we can see the symmetry between the predict and update steps more clearly:The predict step involves convolution in the space domain, which corresponds to multiplication in the frequency domain.The update step involves multiplication in the space domain, which corresponds to convolution in the frequency domain.This observation is interesting, at least to me, and possibly useful, because it suggests an efficient algorithm based on the Fast Fourier Transform (FFT).  The simple implementation of convolution takes time proportional to N², where N is the number of values in the distributions.  Using the FFT, we can get that down to N log N.
Using this algorithm, a complete predict-update step looks like this:Compute the FFT of the distributions for position and velocity.Multiply them elementwise; the result is the characteristic function of the convolved distributions.Compute the inverse FFT of the result, which is the predictive distribution.Multiply the predictive distribution by the likelihood and renormalize.Steps 1 and 3 are N log N; steps 2 and 4 are linear.
I demonstrate this algorithm in this IPython notebook.
For more on the Bayesian part of the Bayes filter, you might want to read Chapter 6 of Think Bayes.  For more on the Convolution Theorem, see Chapter 8 of Think DSP.  And for more about adding random variables, see Chapter 6 of Think Stats.
AcknowledgementsThanks to Paul Ruvolo, the other member of my very exclusive reading club, and Wikipedia for the typeset equations I borrowed.

First babies are more likely to be late

If you are pregnant with your first child, you might have heard that first babies are more likely to be late.  This turns out to be true, although the difference is small.

Averaged across all live births, the mean gestation for first babies is 38.6 weeks, compared to 38.5 weeks for other babies.  This difference is about 16 hours.

Those means include pre-term babies, which affect the averages in a way that understates the difference.  For full-term babies, the differences are a little bigger.

For example, if you are at the beginning of week 36, the average time until delivery is 3.4 weeks for first babies and 3.1 weeks for others, a difference of 1.8 days.  The gap is about the same for weeks 37 through 40.  After that, there is no consistent difference between first babies and others.

The following figure shows average remaining duration in weeks, for first babies and others, computed for weeks 36 through 43.


The gap between first babies and others is consistent until Week 41.  As an aside, this figure also shows a surprising pattern: after Week 38, the expected remaining duration levels off at about one week.  For more than a month, the finish line is always a week away!

Looking at the probability of delivering in the next week, we see a similar pattern: from Week 38 on, the probability is almost the same, with some increase after Week 41.


The difference between first babies and others is highest in Weeks 39 and 40; for example, in Week 39, the chance of delivering in the next week is 52% for first babies, compared to 64% for others.  By Week 41, this gap has closed.

In summary, among full-term pregnancies, first babies arrive a little later than others, by about two days.  After Week 38, the expected remaining duration is about one week.

Methods

The code I used to generate these results is in this IPython Notebook.  I used data from the National Survey of Family Growth (NSFG).  During the last three survey cycles, they interviewed more than 25,000 women and collected data about more than 48,000 pregnancies.  Of those, I selected the 30,110 pregnancies whose outcome was a live birth.

Of those, there were 13,864 first babies and 16,246 others.  The mean gestation period for first babies is 38.61, with SE 0.024; for others it is 38.52 with SE 0.019.  The difference is statistically significant with p < 0.001.

However, those means could be misleading for two reasons: they include pre-term babies, which bring down the averages for both groups.  Also, they do not take into account the stratified survey design.

To address the second point, I use weighted resampling, running each analysis 101 times and selecting the 10th, 50th, and 90th percentile of the results.  The lines in the figure above show median values (50th percentile).  The gray areas show an 80% confidence interval (between the 10th and 90th percentiles).

Background

I use this question—whether first babies are more likely to be late—as a case study in my book, Think Stats .  There, I used data from only one cycle of the NSFG.  I report a small difference between first babies and others, but it is not statistically significant.

I also wrote about this question in a previous blog article, "Are first babies more likely to be late?", which has been viewed more than 100,000 times, more than any other article on this blog.

I am reviewing the question now for two reasons:

1) I worked on another project that required me to load data from other cycles of the NSFG.  Having done that work, I saw an opportunity to run my analysis again with more data.

2) Since my previous articles were intended partly for statistics education, I kept the analysis simple.  In particular, I ignored the stratified design of the survey, which made the results suspect.  Fortunately, it turns out that the effect is small; the new results are consistent with what I saw before.

Since I've been writing about this topic and using it as a teaching example for more than 5 years, I hope the question is settled now.

Bayesian analysis of gluten sensitivity

Last week a new study showed that many subjects diagnosed with non-celiac gluten sensitivity (NCGS) were not able to distinguish gluten flour from non-gluten flour in a blind challenge.

In this article, I review the the study and use a simple Bayesian model to show that the results support the hypothesis that none of the subjects are sensitive to gluten.  But there are complications in the design of the study that might invalidate the model.

Here is a description of the study:

"We studied 35 non-CD subjects (31 females) that were on a gluten-free diet (GFD), in a double-blind challenge study. Participants were randomised to receive either gluten-containing flour or gluten-free flour for 10 days, followed by a 2-week washout period and were then crossed over. The main outcome measure was their ability to identify which flour contained gluten.
"The gluten-containing flour was correctly identified by 12 participants (34%)..."

Since 12 out of 35 participants were able to identify the gluten flour, the authors conclude "Double-blind gluten challenge induces symptom recurrence in just one-third of patients fulfilling the clinical diagnostic criteria for non-coeliac gluten sensitivity."

This conclusion seems odd to me, because if none of the patients were sensitive to gluten, we would expect some of them to identify the gluten flour by chance.  So the results are consistent with the hypothesis that none of the subjects are actually gluten sensitive.

We can use a Bayesian approach to interpret the results more precisely.  But first, as always, we have to make some modeling decisions.

First, of the 35 subjects, 12 identified the gluten flour based on resumption of symptoms while they were eating it.  Another 17 subjects wrongly identified the gluten-free flour based on their symptoms, and 6 subjects were unable to distinguish.  So each subject gave one of three responses.  To keep things simple I follow the authors of the study and lump together the second two groups; that is, I consider two groups: those who identified the gluten flour and those who did not.

Second, I assume (1) people who are actually gluten sensitive have a 95% chance of correctly identifying gluten flour under the challenge conditions, and (2) subjects who are not gluten sensitive have only a 40% chance of identifying the gluten flour by chance (and a 60% chance of either choosing the other flour or failing to distinguish).

Under this model, we can estimate the actual number of subjects who are gluten sensitive, gs.  I chose a uniform prior for gs, from 0 to 35.  To perform the Bayesian analysis, we have to compute the likelihood of the data under each hypothetical value of gs.  Here is the likelihood function in Python

    def Likelihood(self, data, hypo):
        gs = hypo
        yes, no = data
        n = yes + no
        ngs = n - gs
        
        pmf1 = thinkbayes2.MakeBinomialPmf(gs, 0.95)
        pmf2 = thinkbayes2.MakeBinomialPmf(ngs, 0.4)
        pmf = pmf1 + pmf2
        return pmf[yes]

The function works by computing the PMF of the number of gluten identifications conditioned on gs, and then selecting the actual number of identifications, yes, from the PMF.  The details of the computation are in this IPython notebook.

And here is the posterior distribution:

The most likely value of gs is 0, so it is quite possible that none of the respondents are gluten sensitive.  The 95% credible interval for gs is (0, 8), so a reasonable upper bound on the number of gluten-sensitive respondents is 8 out of 35, or 23%.

We can also use this analysis to compare two hypotheses:

A) Some of the respondents are gluten sensitive (equally likely from 0 to 35).
B) None of the respondents are gluten sensitive.

The Bayes factor in support of B turns out to be about 8.4, which is moderately strong.  If you believed, before reading this study, that the probability of B was 50%, you should now believe it is about 90%.

However, there are complications in the design of the study that might invalidate this simple model.  In particular, the gluten free flour in the study contained corn starch, which some people may be sensitive to.  And several subjects reported symptoms under both challenge conditions; that is, when eating both gluten flour and gluten-free flour.  So it is possible that misidentification of the gluten flour, as well as failure to distinguish, indicates sensitivity to both gluten and corn starch.

But if we limit ourselves to the question of whether people diagnosed with non-celiac gluten sensitivity are specifically sensitive to gluten, this study suggests that they are not.

Thank yous: I heard about this study in this blog post.  And I want to thank the authors of the study and their publisher for making the entire paper available on the web, which made my analysis possible.