Probably Overthinking It

Belly Button Biodiversity: The End Game

2013-05-30T11:43:00.004-07:00

In the previous installment of this saga, I admitted that my predictions had completely failed, and I outlined the debugging process I began. Then the semester happened, so I didn't get to work on it again until last week.

It turns out that there were several problems, but the algorithm is now calibrating and validating! Before I proceed, I should explain how I am using these words:

Calibrate: Generate fake data from the same model the analysis is based on. Run the analysis on fake data and generate predictive distributions. Check whether the predictive distributions are correct.
Validate: Using real data, generate a rarefied sample. Run the analysis on the sample and generate predictive distributions. Check whether the predictive distributions are correct.

If the analysis calibrates, but fails to validate, that suggests that there is some difference between the model and reality that is causing a problem. And that turned out to be the case.

Here are the problems I discovered, and what I had to do to fix them:

The prior distribution of prevalences

For the prior I used a Dirichlet distribution with all parameters set to 1. I neglected to consider the "concentration parameter," which represents the prior belief about how uniform or concentrated the prevalences are. As the concentration parameter approaches 0, prevalences tend to be close to 1 or 0; that is, there tends to be one dominant species and many species with small prevalences. As the concentration parameter gets larger, all species tend to have the same prevalence. It turns out that a concentration parameter near 0.1 yields a distribution of prevalences that resembles real data.

The prior distribution of n

With a smaller concentration parameter, there are more species with small prevalences, so I had to increase the range of n (the number of species). The prior distribution for n is uniform up to an upper bound, where I choose the upper bound to be big enough to avoid cutting off non-negligible probability. I had to increase this upper bound to 1000, which slows the analysis down a little, but it still takes only a few seconds per subject (on my not-very-fast computer).

Up to this point I hadn't discovered any real errors; it was just a matter of choosing appropriate prior distributions, which is ordinary work for Bayesian analysis.

But it turns out there were two legitimate errors.

Bias due to the definition of "unseen"

I was consistently underestimating the prevalence of unseen species because of a bias that underlies the definition of "unseen." To see the problem, consider a simple scenario where there are two species, A and B, with equal prevalence. If I only collect one sample, I get A or B with equal probability.

Suppose I am trying to estimate the prevalence of A. If my sample is A, the posterior marginal distribution for the prevalence of A is Beta(2, 1), which has mean 2/3. If the sample is B, the posterior is Beta(1, 2), which has mean 1/3. So the expected posterior mean is the average of 2/3 and 1/3, which is 1/2. That is the actual prevalence of A, so this analysis is unbiased.

But now suppose I am trying to estimate the prevalence of the unseen species. If I draw A, the unseen species is B and the posterior mean is 1/3. If I draw B, the unseen species is A and the posterior mean is 1/3. So either way I believe that the prevalence of the unseen species is 1/3, but it is actually 1/2. Since I did not specify in advance which species is unseen, the result is biased.

This seems obvious in retrospect. So that's embarrassing (the price I pay for this experiment in Open Science), but it is easy to fix:

a) The posterior distribution I generate has the right relative prevalences for the seen species (based on the data) and the right relative prevalences for the unseen species (all the same), but the total prevalence for the unseen species (which I call q) is too low.

b) Fortunately, there is only one part of the analysis where this bias is a problem: when I draw a sample from the posterior distribution. To fix it, I can draw a value of q from the correct posterior distribution (just by running a forward simulation) and then unbias the posterior distribution with the selected value of q.

Here's the code that generates q:

def RandomQ(self, n):

# generate random prevalences
dirichlet = thinkbayes.Dirichlet(n, conc=self.conc)
prevalences = dirichlet.Random()

# generate a simulated sample
pmf = thinkbayes.MakePmfFromItems(enumerate(prevalences))
cdf = pmf.MakeCdf()
sample = cdf.Sample(self.num_reads)
seen = set(sample)

# add up the prevalence of unseen species
q = 0
for species, prev in enumerate(prevalences):
if species not in seen:
q += prev

return q

n is the hypothetical number of species. conc is the concentration parameter. RandomQ creates a Dirichlet distribution, draws a set of prevalences from it, then draws a simulated sample with the appropriate number of reads, and adds up the total prevalence of the species that don't appear in the sample.

And here's the code that unbiases the posterior:

def Unbias(self, n, m, q_desired):

params = self.params[:n].copy()

x = sum(params[:m])
y = sum(params[m:])
a = x + y

g = q_desired * a / y
f = (a - g * y) / x
params[:m] *= f
params[m:] *= g

n is the hypothetical number of species; m is the number seen in the actual data.

x is the total prevalence of the seen species; y is the total prevalence of the unseen species. f and g are the factors we have to multiply by so that the corrected prevalence of unseen species is q_desired.

After fixing this error, I find that the analysis calibrates nicely.

From each predictive distribution I generate credible intervals with ideal percentages 10, 20, ... 90, and then count how often the actual value falls in each interval.

For example, the blue line is the calibration curve for n, the number of species. After 100 runs, the 10% credible interval contains the actual value 9.5% of of the time.The 50% credible interval contains the actual value 51.5% of the time. And the 90% credible interval contains the actual value 88% of the time. These results show that the posterior distribution for n is, in fact, the posterior distribution for n.

The results are similar for q, the prevalence of unseen species, and l, the predicted number of new species seen after additional sampling.

To check whether the analysis validates, I used the dataset collected by the Belly Button Biodiversity project. For each subject with more than 400 reads, I chose a random subset of 100 reads, ran the analysis, and checked the predictive distributions for q and n. I can't check the predictive distribution of n, because I don't know the actual value.

Sadly, the analysis does not validate with the collected data. The reason is:

The data do not fit the model

The data deviate substantially from the model that underlies the analysis. To see this, I tried this experiment:

a) Use the data to estimate the parameters of the model.
b) Generate fake samples from the model.
c) Compare the fake samples to the real data.

Here's a typical result:

The blue line is the CDF of prevalences, in order by rank. The top-ranked species makes up about 27% of the sample. The top 10 species make up about 75%, and the top 100 species make up about 90%.

The green lines show CDFs from 10 fake samples. The model is a good match for the data for the first 10-20 species, but then it deviates substantially. The prevalence of rare species is higher in the data than in the model.

The problem is that the real data seem to come from a mixture of two distributions, one for dominant species and one for rare species. Among the dominant species the concentration parameter is near 0.1. For the rare species, it is much higher; that is, the rare species all have about the same prevalence.

There are two possible explanations: this effect might be real or it might be an artifact of errors in identifying reads. If it's real, I would have to extend my model to account for it. If it is due to errors, it might be possible to clean the data.

I have heard from biologists that when a large sample yields only a single read of a particular species, that read is likely to be in error; that is, the identified species might not actually be present.

So I explored a simple error model with the following features:

1) If a species appears only once after r reads, the probability that the read is bogus is p = (1 - alpha/r), where alpha is a parameter.

2) If a species appears k times after n reads, the probability that all k reads are bogus is p^k.

To clean the data, I compute the probability that each observed species is bogus, and then delete it with the computed probability.

With cleaned data (alpha=50), the model fits very nicely. And since the model fits, and the analysis calibrates, we expect the analysis to validate. And it does.

For n there is no validation curve because we don't know the actual values.

For q the validation curve is a little off because we only have a lower bound for the prevalence of unseen species, so the actual values used for validation are too high.

But for l the validation curve is quite good, and that's what we are actually trying to predict, after all.

At this point the analysis depends on two free parameters, the concentration parameter and the cleaning parameter, alpha, which controls how much of the data gets discarded as erroneous.

So the next step is to check whether these parameters cross-validate. That is, if we tune the parameters based on a training set, how well do those values do on a test set?

Another next step is to improve the error model. I chose something very simple, and it does a nice job of getting the data to conform to the analysis model, but it is not well motivated. If I can get more information about where the errors are coming from, I could take a Bayesian approach (what else?) and compute the probability that each datum is legit or bogus.

Or if the data are legit and the prevalences are drawn from a mixture of Dirichlet distributions with different concentrations, I will have to extend the analysis accordingly.

Summary

There were four good reasons my predictions failed:

1) The prior distribution of prevalences had the wrong concentration parameter.

2) The prior distribution of n was too narrow.

3) I neglected an implicit bias due to the definition of "unseen species."

4) The data deviate from the model the analysis is based on. If we "clean" the data, it fits the model and the analysis validates, but the cleaning process is a bit of a hack.

I was able to solve these problems, but I had to introduce two free parameters, so my algorithm is not as versatile as I hoped. However, it seems like it should be possible to choose these parameters automatically, which would be an improvement.

And now I have to stop, incorporate these corrections into Think Bayes, and then finish the manuscript!

Python Epistemology at PyCon Taiwan

2013-05-28T13:09:00.000-07:00

This weekend I gave a talk entitled "Python Epistemology" for PyCon Taiwan 2013. I would have loved to be in Taipei for the talk, but sadly I was in an empty room in front of a teleconference screen.

Python Epistemology: PyCon Taiwan 2013

As I explained, the title is more grandiose than accurate. In general, epistemology is the theory of knowledge: how we know what we think we know, etc. This talk is mostly about what Python has taught me about programming, and how programming in Python has changed the way I learn and the way I think.

About programming, I wrote:

Programming is not about translating a well-known solution into code, it is about discovering solutions by writing code, and then creating the language to express them.

I gave an example using the Counter data structure to check for anagrams:

from collections import Counter

def is_anagram(word1, word2):

return Counter(word1) == Counter(word2)

This is a nice solution because it is concise and demonstrably correct, but I suggested that one limitation is that it does not extend easily to handle "The Scrabble Problem": given a set of tiles, check to see whether you can spell a given word.

We can define a new class, called Multiset, that extends Counter and provides

is_subset

class Multiset(Counter):
"""A set with repeated elements."""

def is_subset(self, other):

       for char, count in self.items():
           if other[char] < count:
               return False
       return True

Now we can write can_spell concisely:

def can_spell(word, tiles):
return Multiset(word).is_subset(Multiset(tiles))

I summarized by quoting Paul Graham:

"... you don't just write your program down toward the language, you also build the language up toward your program.

"In the end your program will look as if the language had been designed for it. And ... you end up with code which is clear, small, and efficient."

Paul Graham, "Programming Bottom Up," 1993.

In the second half of the talk, I suggested that Python and other modern programming languages provide a new approach to solving problems. Traditionally, we tend to think and explore using natural language, do analysis and solve problems using mathematical notation, and then translate solutions from math notation into programming languages.

In some sense, we are always doing two translations, from natural language to math and from math to a computer program. With the previous generation of programming languages, this process was probably necessary (for reasons I explained), but I claim that it is less necessary now, and that it might be possible and advantageous to skip the intermediate mathematics and do analysis and problem-solving directly in programming languages.

After the talk, I got two interesting questions by email. Yung-Cheng Lin suggested that although programming languages are more precise than natural language, they might not be sufficiently precise to replace mathematical notation, and he asked if I think that using programming to teach mathematical concepts might cause misunderstandings for students.

I replied:

I understand what you mean when you say that programming languages are less rigorous that mathematical notation. I think many people have the same impression, but I wonder if it is a real difference or a bias we have.

I would argue that programming languages and math notation are similar in the sense that they are both formal languages designed by people to express particular ideas concisely and precisely.

There are some kinds of work that are easier to do in math notation, like algebraic manipulation, but other kinds of work that are easier in programming languages, like specifying computations, especially computations with state.

You asked if there is a danger that students might misunderstand mathematical ideas if they come to them through programming, rather than mathematical instruction. I'm sure it's possible, but I don't think the danger is specific to the programming approach.

And on the other side, I think a computational approach to mathematical topics creates opportunities for deeper understanding by running experiments, and (as I said in the talk) by getting your ideas out of your head and into a program so that, by debugging the program, you are also debugging your own understanding.

In response to some of my comments about pseudocode, A. T. Cheng wrote:

When we do algorithms or pseudocodes in the traditional way, we used to figure out the time complexity at the same time. But the Python examples you showed us, it seems not so easy to learn the time complexity in the first place. So, does it mean that when we think Python, we don't really care about the time complexity that much?

I replied:

You are right that it can be more difficult to analyze a Python program; you have to know a lot about how the Python data structures are implemented. And there are some gotchas; for example, it takes constant time to add elements to the end of a list, but linear time to add elements in the beginning or the middle.

It would be better if Python made these performance characteristics part of the interface, but they are not. In fact, some implementations have changed over time; for example, the += operator on lists used to create a new list. Now it is equivalent to append.

Thanks to both of my correspondents for these questions (and for permission to quote them). And thanks to the organizers of PyCon Taiwan, especially Albert Chun-Chieh Huang, for inviting me to speak. I really enjoyed it.

The Red Line problem

2013-05-09T10:35:00.002-07:00

I've just added a new chapter to Think Bayes; it is a case study based on a class project two of my students worked on this semester. It presents "The Red Line Problem," which is the problem of predicting the time until the next train arrives, based on the number of passengers on the platform.

Here's the introduction:

In Boston, the Red Line is a subway that runs north-south from Cambridge to Boston. When I was working in Cambridge I took the Red Line from Kendall Square to South Station and caught the commuter rail to Needham. During rush hour Red Line trains run every 7--8 minutes, on average.

When I arrived at the station, I could estimate the time until the next train based on the number of passengers on the platform. If there were only a few people, I inferred that I just missed a train and expected to wait about 7 minutes. If there were more passengers, I expected the train to arrive sooner. But if there were a large number of passengers, I suspected that trains were not running on schedule, so I would go back to the street level and get a taxi.

While I was waiting for trains, I thought about how Bayesian estimation could help predict my wait time and decide when I should give up and take a taxi. This chapter presents the analysis I came up with.

Sadly, this problem has been overtaken by history: the Red Line now provides real-time estimates for the arrival of the next train. But I think the analysis is interesting, and still applies for subway systems that don't provide estimates.

One interesting tidbit:

As it turns out, the average time between trains, as seen by a random passenger, is substantially higher than the true average.

Why? Because a passenger is more like to arrive during a large interval than a small one. Consider a simple example: suppose that the time between trains is either 5 minutes or 10 minutes with equal probability. In that case the average time between trains is 7.5 minutes.
But a passenger is more likely to arrive during a 10 minute gap than a 5 minute gap; in fact, twice as likely. If we surveyed arriving passengers, we would find that 2/3 of them arrived during a 10 minute gap, and only 1/3 during a 5 minute gap. So the average time between trains, as seen by an arriving passenger, is 8.33 minutes.

This kind of observer bias appears in many contexts. Students think that classes are bigger than they are, because more of them are in the big classes. Airline passengers think that planes are fuller than they are, because more of them are on full flights.

In each case, values from the actual distribution are oversampled in proportion to their value. In the Red Line example, a gap that is twice as big is twice as likely to be observed.

The data for the Red Line are close to this example. The actual time between trains is 7.6 minutes (based on 45 trains that arrived at Kendall square between 4pm and 6pm so far this week). The average gap as seen by random passengers is 8.3 minutes.

Interestingly, the Red Line schedule reports that trains run every 9 minutes during peak times. This is close to the average seen by passengers, but higher than the true average. I wonder if they are deliberately reporting the mean as seen by passengers in order to forestall complaints.

You can read the rest of the chapter here. One of the figures there didn't render very well. Here is a prettier version:

This figure shows the predictive distribution of wait times if you arrive and find 15 passengers on the platform. Since we don't know the passenger arrival rate, we have to estimate it. Each possible arrival rate yields one of the light blue lines; the dark blue line is the weighted mixture of the light blue lines.

So in this scenario, you expect the next train in 5 minutes or less, with 80% confidence.

UPDATE 10 May 2013: I got the following note from developer@mbta.com, confirming that their reported gap between trains is deliberately conservative:

Thank you for writing to let us know about the Red Line case study in your book, and thank you for your question. You’re right that the scheduled time between trains listed on the subway schedule card for rush hour is greater than what you observed at Kendall Square. It’s meant as a slightly conservative simplification of the actual frequency of trains, which varies by time throughout rush hour – to provide maximum capacity during the very peak of rush hour when ridership is normally highest – as well as by location along the Red Line during those different times, since when trains begin to leave more frequently from Alewife it takes time for that frequency to “travel” down the line. So yes it is meant to be slightly conservative for that reason. We hope this information answers your question.

Sincerely,
developer@mbta.com

Software engineering practices for graduate students

2013-05-06T12:28:00.001-07:00

Recently I was talking with an Olin student who will start graduate school in the fall, and I suggested a few things I wish I had done in grad school. And then I thought I should write them down. So here is my list of Software Engineering Practices All Graduate Students Should Adopt:

Version Control

Every keystroke you type should be under version control from the time you initiate a project until you retire it. Here are the reasons:

1) Everything you do will be backed up. But instead of organizing your backups by date (which is what most backup systems do) they are organized by revision. So, for example, if you break something, you can roll back to an earlier working revision.

2) When you are collaborating with other people, you can share repositories. Version control systems are well designed for managing this kind of collaboration. If you are emailing documents back and forth, you are doing it wrong.

3) At various stages of the project, you can save a tagged copy of the repo. For example, when you submit a paper for publication, make a tagged copy. You can keep working on the trunk, and when you get reviewer comments (or a question 5 years later) you have something to refer back to.

I use Subversion (SVN) primarily, so I keep many of my projects on Google Code (if they are open source) or on my own SVN server. But these days it seems like all the cool kids are using git and keeping their repositories on github.

Either way, find a version control system you like, learn how to use it, and find someplace to host your repository.

Build Automation

This goes hand in hand with version control. If someone checks out your repository, they should be able to rebuild your project by running a single command. That means that everything someone needs to replicate your results should be in the repo, and you should have scripts that process the data, generate figures and tables, and integrate them into your papers, slides, and other documents.

One simple tool for automating the build is Make. Every directory in your project should contain a Makefile. The top-level directory should contain the Makefile that runs all the others.

If you use GUI-based tools to process data, it might not be easy to automate your build. But it will be worth it. The night before your paper is due, you will find a bug somewhere in your data flow. If you've done things right, you should be able to rebuild the paper with just five keystrokes (m-a-k-e, and Enter).

Also, put a README in the top-level directory that documents the directory structure and the build process. If your build depends on other software, include it in the repo if practical; otherwise provide a list of required packages.

Or, if your software environment is not easy to replicate, put your whole development environment in a virtual machine and ship the VM.

Agile Development

For many people, the most challenging part of grad school is time management. If you are an undergraduate taking 4-5 classes, you can do deadline-driven scheduling; that is, you can work on whatever task is due next and you will probably get everything done on time.

In grad school, you have more responsibility for how you spend your time and fewer deadlines to guide you. It is easy to lose track of what you are doing, waste time doing things that are not important (see Yak Shaving), and neglect the things that move you toward the goal of graduation.

One of the purposes of agile development tools is to help people decide what to do next. They provide several features that apply to grad school as well as software development:

1) They encourage planners to divide large tasks into smaller tasks that have a clearly-defined end condition.

2) They maintain a priority-ranking of tasks so that when you complete one you can start work on the next, or one of the next few.

3) They provide mechanisms for collaborating with a team and for getting feedback from an adviser.

4) They involve planning on at least two time scales. On a daily basis you decide what to work on by selecting tasks from the backlog. On a weekly (or longer) basis, you create and reorder tasks, and decide which ones you should work on during the next cycle.

If you use github or Google code for version control, you get an issue tracker as part of the deal. You can use issue trackers for agile planning, but there are other tools, like Pivotal Tracker, that have more of the agile methodology built in. I suggest you start with Pivotal Tracker because it has excellent documentation, but you might have to try out a few tools to find one you like.

Do these things -- Version Control, Build Automation, and Agile Development -- and you will get through grad school in less than the average time, with less than the average drama.

The Price is Right Problem: Part Two

2013-04-24T06:49:00.000-07:00

This article is an excerpt from Think Bayes, a book I am working on. The entire current draft is available from http://thinkbayes.com. I welcome comments and suggestions.

In the previous article, I described presented The Price is Right problem and a Bayesian approach to estimating the value of a showcase of prizes. This article picks up from there...

Optimal bidding

Now that we have a posterior distribution, we can use it to compute the optimal bid, which I define as the bid that maximizes expected gain.

To compute optimal bids, I wrote a class called GainCalculator:

class GainCalculator(object):

    def __init__(self, player, opponent):
        self.player = player
        self.opponent = opponent

player and opponent are Player objects.

GainCalculator provides ExpectedGains, which computes a sequence of bids and the expected gain for each bid:

    def ExpectedGains(self, low=0, high=75000, n=101):
        bids = numpy.linspace(low, high, n)

        gains = [self.ExpectedGain(bid) for bid in bids]

        return bids, gains

low and high specify the range of possible bids; n is the number of bids to try. Here is the function that computes expected gain for a given bid:

    def ExpectedGain(self, bid):
        suite = self.player.posterior
        total = 0
        for price, prob in suite.Items():
            gain = self.Gain(bid, price)
            total += prob * gain
        return total

ExpectedGain loops through the values in the posterior and computes the gain for each bid, given the actual prices of the showcase. It weights each gain with the corresponding probability and returns the total.
Gain takes a bid and an actual price and returns the expected gain:

    def Gain(self, bid, price):
        if bid > price:
            return 0

        diff = price - bid
        prob = self.ProbWin(diff)

        if diff <= 250:
            return 2 * price * prob
        else:
            return price * prob

If you overbid, you get nothing. Otherwise we compute the difference between your bid and the price, which determines your probability of winning.

If diff is less than $250, you win both showcases. For simplicity, I assume that both showcases have the same price. Since this outcome is rare, it doesn’t make much difference.

Finally, we have to compute the probability of winning based on diff:

    def ProbWin(self, diff):
        prob = (self.opponent.ProbOverbid() + 
                self.opponent.ProbWorseThan(diff))
        return prob

If your opponent overbids, you win. Otherwise, you have to hope that your opponent is off by more than diff. Player provides methods to compute both probabilities:

# class Player:

    def ProbOverbid(self):
        return self.cdf_diff.Prob(-1)

    def ProbWorseThan(self, diff):
        return 1 - self.cdf_diff.Prob(diff)

This code might be confusing because the computation is now from the point of view of the opponent, who is computing, “What is the probability that I overbid?” and “What is the probability that my bid is off by more than diff?”

Both answers are based on the CDF of diff [CDFs are described here]. If your opponent’s diff is less than or equal to -1, you win. If your opponent’s diff is worse than yours, you win. Otherwise you lose.

Finally, here’s the code that computes optimal bids:

# class Player:

    def OptimalBid(self, guess, opponent):
        self.MakeBeliefs(guess)
        calc = GainCalculator(self, opponent)
        bids, gains = calc.ExpectedGains()
        gain, bid = max(zip(gains, bids))
        return bid, gain

Given a guess and an opponent, OptimalBid computes the posterior distribution, instantiates a GainCalculator, computes expected gains for a range of bids and returns the optimal bid and expected gain. Whew!

Figure 6.4

Figure 6.4 shows the results for both players, based on a scenario where Player 1’s best guess is $20,000 and Player 2’s best guess is $40,000.

For Player 1 the optimal bid is $21,000, yielding an expected return of almost $16,700. This is a case (which turns out to be unusual) where the optimal bid is actually higher than the contestant’s best guess.

For Player 2 the optimal bid is $31,500, yielding an expected return of almost $19,400. This is the more typical case where the optimal bid is less than the best guess.

Discussion

One of the most useful features of Bayesian estimation is that the result comes in the form of a posterior distribution. Classical estimation usually generates a single point estimate or a confidence interval, which is sufficient if estimation is the last step in the process, but if you want to use an estimate as an input to a subsequent analysis, point estimates and intervals are often not much help.

In this example, the Bayesian analysis yields a posterior distribution we can use to compute an optimal bid. The gain function is asymmetric and discontinuous (if you overbid, you lose), so it would be hard to solve this problem analytically. But it is relatively simple to do computationally.

Newcomers to Bayesian thinking are often tempted to summarize the posterior distribution by computing the mean or the maximum likelihood estimate. These summaries can be useful, but if that’s all you need, then you probably don’t need Bayesian methods in the first place.

Bayesian methods are most useful when you can carry the posterior distribution into the next step of the process to perform some kind of optimization, as we did in this chapter, or some kind of prediction, as we will see in the next chapter [which you can read here].

The Price is Right Problem

2013-04-22T11:36:00.001-07:00

This article is an excerpt from Think Bayes, a book I am working on. The entire current draft is available from http://thinkbayes.com. I welcome comments and suggestions.

The Price is Right Problem

On November 1, 2007, contestants named Letia and Nathaniel appeared on The Price is Right, an American game show. They competed in a game called the Showcase, where the objective is to guess the price of a showcase of prizes. The contestant who comes closest to the actual price of the showcase, without going over, wins the prizes.

Nathaniel went first. His showcase included a dishwasher, a wine cabinet, a laptop computer, and a car. He bid $26,000.

Letia’s showcase included a pinball machine, a video arcade game, a pool table, and a cruise of the Bahamas. She bid $21,500.

The actual price of Nathaniel’s showcase was $25,347. His bid was too high, so he lost.

The actual price of Letia’s showcase was $21,578. She was only off by $78, so she won her showcase and, because her bid was off by less than $250, she also won Nathaniel’s showcase.

For a Bayesian thinker, this scenario suggests several questions:

Before seeing the prizes, what prior beliefs should the contestant have about the price of the showcase?
After seeing the prizes, how should the contestant update those prior beliefs?
Based on the posterior distribution, what should the contestant bid?

The third question demonstrates a common use of Bayesian analysis: optimization. Given a posterior distribution, we can choose the bid that maximizes the contestant’s expected return.

This problem is inspired by an example in Cameron Davidson-Pilon’s book, Bayesian Methods for Hackers.

The prior

To choose a prior distribution of prices, we can take advantage of data from previous episodes. Fortunately, fans of the show keep detailed records. When I corresponded with Mr. Davidson-Pilon about his book, he sent me data collected by Steve Gee at http://tpirsummaries.8m.com. It includes the price of each showcase from the 2011 and 2012 seasons, and the bids offered by the contestants.

Figure 6.1: Distribution of prices for showcases on The Price is Right, 2011-12.

Figure 6.1 shows the distribution of prices for these showcases. The most common value for both showcases is around $28,000, but the first showcase has a second mode near $50,000, and the second showcase is occasionally worth more than $70,000.

These distributions are based on actual data, but they have been smoothed by Gaussian kernel density estimation (KDE). So before we go on, I want to take a detour to talk about probability density functions and KDE.

Probability density functions

So far [in Think Bayes, that is] we have been working with probability mass functions, or PMFs. A PMF is a mapping from each possible value to its probability. In my implementation, a Pmf object provides a method named Prob that takes a value and returns a probability, also known as a “probability mass.”

A probability density function, or PDF, is the continuous version of a PMF, where the possible values make up a continuous range rather than a discrete set.

In mathematical notation, PDFs are usually written as functions; for example, here is the PDF of a Gaussian distribution with mean 0 and standard deviation 1:

f(x) = exp(−x²)

To represent this PDF in Python, I could define a class like this:

class StandardGaussianPdf(object):

    def Density(self, x):
        return math.exp(-x**2)

Density takes a value, x, and returns the probability density evaluated at x.

A probability density is similar to a probability mass in one way: higher density indicates that a value is more likely.

But a density is not a probability. If you integrate a density over a continuous range, the result is a probability. But for the applications in this book we seldom have to do that.

In this book we primarily use probability densities as part of a Likelihood function. We will see an example soon.

Representing Pdfs

Before we get back to The Price is Right, I want to present a more general way to represent PDFs.
thinkbayes.py provides a class named Pdf that defines two functions, Density and MakePmf:

class Pdf(object):

    def Density(self, x):
        raise UnimplementedMethodException()

    def MakePmf(self, xs):
        pmf = Pmf()
        for x in xs:
            pmf.Set(x, self.Density(x))
        pmf.Normalize()
        return pmf

Pdf is an abstract type, which means that it defines the interface a Pdf is supposed to have, but does not provide a complete implementation. Specifically, Pdf provides MakePmf but not Density. [PMFs are described in Chapter 2 of Think Bayes.]

A concrete type is a class that extends an abstract parent class and provides an implementation of the missing methods.

For example, GaussianPdf extends Pdf and provides Density:

class GaussianPdf(Pdf):

    def __init__(self, mu, sigma):
        self.mu = mu
        self.sigma = sigma
        
    def Density(self, x):
        density = scipy.stats.norm.pdf(x, 
                                       loc=self.mu, 
                                       scale=self.sigma)
        return density

__init__ takes mu and sigma, which are the mean and standard deviation of the distribution.
Density uses a function from scipy to evaluate the Gaussian PDF.

The Gaussian PDF is defined by a simple mathematical function, so it is easy to evaluate. And it is useful because many quantities in the real world have distributions that are approximately Gaussian.
But with real data, there is no guarantee that the PDF is Gaussian, or any other simple mathematical function. In that case we can use a data sample to estimate the PDF of the whole population.

For example, in The Price Is Right data, we have 313 prices for the first showcase. We can think of these values as a sample from the population of all possible showcase prices.

Near the middle of the distribution, we see the following values:

28800, 28868, 28941, 28957, 28958

In the sample, no values appear between 28801 and 28867, but there is no reason to think that these values are impossible. Based on our background information, we would expect all values in this range to be equally likely. In other words, we expect the PDF to be reasonably smooth.

Kernel density estimation (KDE) is an algorithm that takes a sample of values and finds an appropriately-smooth PDF that fits the data. You can read about the details at http://en.wikipedia.org/wiki/Kernel_density_estimation.

scipy provides an implementation of KDE. thinkbayes provides a class called EstimatedPdf that extends Pdf and uses KDE:

class EstimatedPdf(Pdf):

    def __init__(self, sample):
        xs = numpy.array(sample, dtype=numpy.double)
        self.kde = scipy.stats.gaussian_kde(xs)

    def Density(self, x):
        return self.kde.evaluate(x)

__init__ takes a sample, converts it to a NumPy array, and computes a kernel density estimate. The result is a gaussian_kde object that provides an evaluate method.

Density takes a value, calls gaussian_kde.evaluate, and returns the resulting density.
Finally, here’s an outline of the code I used to generate Figure 6.1:

    prices = ReadData()
    kde = thinkbayes.EstimatedPdf(prices)

    low, high = 0, 75000
    n = 101
    xs = numpy.linspace(low, high, n) 
    pmf = kde.MakePmf(xs)

    myplot.Pmf(pmf)

And now back to The Price is Right.

Modeling the contestants

The PDFs in Figure 6.1 estimate the distribution of possible prices for each showcase. If you were a contestant on the show, you could use this distribution to quantify your prior belief about the price of the showcases (before you see the prizes).

To update these priors, you have to answer these questions:

What data should we consider and how should we quantify it?
Can we compute a Likelihood function; that is, for each hypothetical value of price, can we compute the conditional likelihood of the data?

To answer these questions, I am going to model the contestant as a price-guessing instrument with known error characteristics. In other words, when the contestant sees the prizes, he or she guesses the price of each prize—ideally without taking into consideration the fact that the prize is part of a showcase—and adds up the prices. Let’s call this total guess.

Under this model, the question we have to answer is, “If the actual price is price, what is the likelihood that the contestant’s total estimate would be guess?”

Or if we define

    error = price - guess

then we could ask, “What is the likelihood that the contestant’s estimate is off by error?”
To answer this question, we can use the historical data again. Figure 6.2 shows the cumulative distribution of diff, the difference between the contestant’s bid and the actual price of the showcase.
The definition of diff is

    diff = price - bid

When diff is negative, the bid is too high. As an aside, we can use this CDF to compute the probability that the contestants overbid: the first contestant overbids 25% of the time; the second contestant overbids 29% of the time.

Figure 6.2: Cumulative distribution (CDF) of the difference between the contestant’s bid and the actual price.

We can also use this distribution to estimate the reliability of the contestants’ guesses. This step is a little tricky because we don’t actually know the contestant’s guesses; we only know what they bid.
In Figure 6.2 we can see that the bids are biased; that is, they are more likely to be too low than too high. And that makes sense, given the rules of the game.

So we’ll have to make some assumptions. Specifically, I assume that the distribution of error is Gaussian with mean 0 and the same variance as diff.

The Player class implements this model:

class Player(object):

    def __init__(self, price, bid, diff):
        self.price = price
        self.bid = bid
        self.diff = diff

        self.pdf_price = thinkbayes.EstimatedPdf(price)
        self.cdf_diff = thinkbayes.MakeCdfFromList(diff)

        mu = 0
        sigma = numpy.std(self.diff)
        self.pdf_error = thinkbayes.GaussianPdf(mu, sigma)

price is a sequence of showcase prices, bid is a sequence of bids, and diff is a sequence of diffs, where again diff = price - bid.

pdf_price is the smoothed PDF of prices, estimated by KDE. cdf_diff is the cumulative distribution of diff, which we saw in Figure 6.2. And pdf_error is the PDF that characterizes the distribution of errors; where error = price - guess.

Again, we use the variance of diff to estimate the variance of error. But contestant’s bids are sometimes strategic; for example, if Player 2 thinks that Player 1 has overbid, Player 2 might make a very low bid. In that case diff does not reflect error. If this strategy is common, the observed variance in diff might overestimate the variance in error. Nevertheless, I think it is a reasonable modeling decision.

As an alternative, someone preparing to appear on the show could estimate their own distribution of error by watching previous shows and recording their guesses and the actual prices.

Likelihood

Now we are ready to write the likelihood function. As usual, I define a new class that extends thinkbayes.Suite:

class Price(thinkbayes.Suite):

    def __init__(self, pmf, player):
        thinkbayes.Suite.__init__(self)

        for price, prob in pmf.Items():
            self.Set(price, prob)

        self.player = player

pmf represents the prior distribution. The for loop copies the values and probabilities from pmf into the new Suite.

player is a Player object as described in the previous section. And here’s Likelihood:

    def Likelihood(self, hypo, data):
        price = hypo
        guess = data

        error = price - guess
        like = self.player.ErrorDensity(error)

        return like

hypo is the hypothetical price of the showcase. data is the contestant’s best guess at the price. error is the difference, and like is the likelihood of the data, given the hypothesis.

ErrorDensity is defined in Player:

# class Player:

    def ErrorDensity(self, error):
        return self.pdf_error.Density(error)

ErrorDensity works by evaluating pdf_error at the given value of error.

The result is a probability density, which means we can’t treat it as a probability. But remember that Likelihood does not really need to compute a probability; it only has to compute something proportional to a probability. As long as the constant of proportionality is the same for all likelihoods, it gets cancelled out when we normalize the posterior distribution.

And therefore, a probability density is a perfectly good likelihood.

Update

Player provides a method that takes the contestant’s guess and computes the posterior distribution:

    def MakeBeliefs(self, guess):
        pmf = self.PmfPrice()
        self.prior = Price(pmf, self, name='prior')
        self.posterior = self.prior.Copy(name='posterior')
        self.posterior.Update(guess)

PmfPrice evaluates pdf_price at an equally-spaced series of values:

    def PmfPrice(self):
        return self.pdf_price.MakePmf(self.price_xs)

The result is a new Pmf object, which we use to construct the prior. To construct the posterior, we make a copy of the prior and then invoke Update, which invokes Likelihood for each hypothesis, multiplies the priors by the likelihoods, and then renormalizes.

So let’s get back to the original scenario. Suppose you are Player 1 and when you see your showcase, your best guess is that the total price of the prizes is $20,000.

The following code constructs and plots your prior and posterior beliefs about the actual price:

    player1.MakeBeliefs(20000)
    myplot.Pmf(player1.prior)
    myplot.Pmf(player2.prior)









Figure 6.3: Prior and posterior distributions for Player 1, based on a best guess of $20,000.

Figure 6.3 shows the results. The value of your guess is on the low end of the prior range, so the posterior is shifted to the left. The mean of the posterior is $25,096; the most likely value is $24,000.
On one level, this result makes sense. The most likely value in the prior is $27,750. Your best guess is $20,000. And the most likely value in the posterior is about half way in between.

On another level, you might find this result bizarre, because it suggests that if you think the price is $20,000, then you should believe the price is $24,000.

To resolve this apparent paradox, remember that you are combining two sources of information, historical data about past showcases and guesses about the prizes you see.

We are treating the historical data as the prior and updating it based on your guesses, but we could equivalently use your guess as a prior and update it based on historical data.

If you think of it that way, maybe it is less surprising that the most likely value in the posterior is not your original guess.

In the next installment, we'll use the posterior distribution to compute the optimal bid for each player.

The price is right

2013-04-11T08:55:00.001-07:00

On the reddit statistics forum, I recently posted a link to my tutorial on Bayesian statistics. One of my fellow redditors drew my attention to this article, which uses pymc to do a Bayesian analysis of The Price is Right. He or she asked how I would solve this problem using the framework in Think Bayes. So, here goes.

First, we have to define a Suite of hypotheses. In this example, each hypothesis represents a belief about the total price of the showcase. We are told, based on data from previous shows, that a reasonable prior distribution is normal with mean 35000 dollars and standard deviation 7500.

So here's the code that creates the Suite:

class Price(thinkbayes.Suite):

def __init__(self, error_sigma):

"""Constructs the suite.

error_sigma: standard deviation of the distribution of error

"""

thinkbayes.Suite.__init__(self)

pmf = thinkbayes.MakeGaussianPmf(35000, 7500, num_sigmas=4)

# copy items from pmf to self

for val, prob in pmf.Items():

self.Set(val, prob)

# store error_sigma for use in Likelihood

self.error_sigma = error_sigma

MakeGaussianPmf makes a Pmf (probability mass function) that approximates a normal distribution. It truncates the range of the distribution 4 sigmas in each direction from the mean.

I'll explain error_sigma in a minute.

Now that we have a suite of hypotheses, we think about how to represent the data. In this case, the "data" is my best guess about the total value of the showcase, which we can think of as a measurement produced by a somewhat unreliable instrument, my brain.

According to the problem statement, my best guess is 3000 for the price of the snowmobile, 12000 for the price of the trip, so 15000 for the total price.

Now we need a Likelihood function that takes a hypothetical price for the showcase, and my best guess, and returns the Likelihood of the data given the hypothesis. That is, if the actual price of the showcase is X, what is the probability that I would guess 12000?

To answer that question, we need some information about how good I am at guessing prices, and that's where error_sigma comes in.

We are told that my uncertainty about the price of the snowmobile can be captured by a normal distribution with sigma=500. And my uncertainty about the price of the trip is normal with sigma=3000. So let's compute the standard deviation of my total error:

error_snowmobile = 500

error_trip = 3000

error_total = math.sqrt(error_snowmobile**2 + error_trip**2)

Now we can create the Suite

suite = Price(error_total)

and update it with my guess

my_guess = 15000

suite.Update(my_guess)

When we invoke Update, it invokes Likelihood once for each hypothetical showcase price. Here's the Likelihood function:

def Likelihood(self, hypo, data):

"""Computes the likelihood of the data under the hypothesis.

hypo: actual price

data: my guess

"""

actual_price = hypo

my_guess = data

error = my_guess - actual_price

like = thinkbayes.EvalGaussianPdf(

mu=0,

sigma=self.error_sigma,

x=error)

return like

x is the error; that is, how much my guess is off by. like is the likelihood that I would be off by that much, computing by evaluating the density function of the Gaussian with sigma=error_sigma.

And we're done. Let's see what the results look like.

The mean of the posterior distribution is 17800, substantially less than the prior mean (35000). It is also substantially less than the result reported in the original article, near 28000. So we are left with a suite of three hypotheses:

1) There is a mistake in my implementation.

2) There is a mistake in the other author's implementation.

3) We are actually making subtly different modeling assumptions.

Option #3 is possible because I am not sure I understand the model in the original article. The author describes my beliefs about the snowmobile and the trip as "priors", which suggests that they are going to get updated. In contrast, I am treating my guess about the prices as data (that is, a summary of what I learned by seeing the contents of the showcase), but I am also modeling myself as a measurement instrument with a characteristic distribution of errors.

Under my interpretation, the posterior shown above makes sense. For example, if my guess is 15000, and the standard deviation of my guesses is 3050, then it is very unlikely that I am off by 4 standard deviations, so the upper bound of the posterior should be around 15000 + 4 * 3050 = 27200. That makes the mean reported in the original article (28000) seem too high.

But maybe I am not interpreting the statement of the problem (or the model) as intended. I will check in with my correspondent and update this article when we have an explanation!

UPDATE April 18, 2013: I exchanged a few emails with the author of the original article, Cameron Davidson-Pilon. He found a bug in his code that explains at least part of the difference between his results and mine. So I think he is planning to update his article and the book he is working on.

He also sent me some data on the value of recent showcases on The Price is Right and the bids offered by the contestants. The data were collected by Steve Gee and posted at The Price is Right Stats.

I have written code that uses this data to form the prior distribution of prices, and also to estimate the distribution of errors for the contestants. And I am writing it up as a chapter in Think Bayes. I'll post the new chapter here when it is done!

UPDATE April 22, 2013: I have added a chapter to Think Bayes, and I am publishing it as a two-part series, starting here.

Freshman hordes regress to the mean

2013-04-09T10:53:00.000-07:00

More nones, no nuns

For several years I have been following one of the most under-reported stories of the decade: the fraction of college freshmen who report no religious preference has tripled since 1985, from 8% to 24%, and the trend is accelerating.

Two years ago I wrote Freshman hordes more godless than ever; last year I updated it with Freshman hordes even more godless. Each year, the number of students with no religious preference increased, and the number attending religious services decreased.

In last year's installment, I made the bold prediction that the trend would continue, and that the students starting college in 2012 would again, be the most godless ever. It turns out I was wrong: attendance went up slightly, and the fraction of "Nones" dropped slightly, in both cases reverting toward long-term trends.

My analysis is based on survey results from the Cooperative Institutional Research Program (CIRP) of the Higher Education Research Insitute (HERI). In 2012, more than 190,000 students at 283 colleges and universities completed the CIRP Freshman Survey, which includes questions about students’ backgrounds, activities, and attitudes.

In one question, students select their “current religious preference,” from a choice of seventeen common religions, “Other religion,” or “None.”

Another question asks students how often they “attended a religious service” in the last year. The choices are “Frequently,” “Occasionally,” and “Not at all.” Students are instructed to select “Occasionally” if they attended one or more times.

The following figure shows the fraction of Nones over more than 40 years of the survey:

The blue line shows actual data through 2011; the red line shows a quadratic fit to the data. The dark gray region shows a 90% confidence interval, taking into account sampling error, so it reflects uncertainty about the parameters of the fit.

The light gray region shows a 90% confidence interval taking into account both sampling error and residual error. So it reflects total uncertainty about the predicted value, including uncertainty due to random variation from year to year.

We expect the new data point from 2012, shown as a blue square, to fall within the light gray interval, and it does. In fact, at 23.8% it falls almost exactly on the fitted curve.

Here is the corresponding plot for attendance at religious services:

Again, the new data point for 2012, 26.8%, falls comfortably in the predicted range. Don't listen to Nate Silver; prediction is easy :)

Predictions for 2013

Using the new 2012 data, we can generate predictions for 2013. Here is the revised plot for "Nones":

The prediction for next year is that the fraction of Nones will hit a new all-time high at 25% (up from 23.8%).

And here is the prediction for "No attendance":

The prediction for 2013 is a small decrease to 26.6% (from 26.8%). I'll be back next year to check on these predictions.

Other updates

1) This year the survey repeated two questions from 2010, asking students if they consider themselves "Born again Christian" or "Evangelical". The fraction reporting "Born again" dropped from 22.8% to 20.2%. The fraction who consider themselves Evangelical dropped from 8.9% to 8.5%. But it's too early to declare a trend.

2) As always, more males than females report no religious preference. The gender gap increased this year, but still falls in the predicted range, as shown in the following plot:

Evidence that the gender gap is increasing is strong. The p-value of the slope of the fitted curve is less than 10e-5.

Data Source

The American Freshman: National Norms Fall 2012
Pryor, J.H., Eagan, K., Palucki Blake, L., Hurtado, S., Berdan, J., Case, M.H.
ISBN: 978-1-878477-22-4 90 pages.
Jan 2013

This and all previous reports are available from the HERI publications page.

Belly Button Biodiversity: Part Four

2013-03-22T11:40:00.002-07:00

March 22, 2013

Well, I've started testing the predictions I made in my previous post, and exactly as I expected and deserved, I am getting killed. The actual results pretty consistently show more species than I predicted, sometimes way more.

I have started the process of debugging the problem. Of course, now that I have looked at the right answers, I can no longer use this data set for validation, especially since I plan to bang on my algorithm until it produces the right answers.

But in the interest of transparent science, I will at least document the debugging process. My first step was to review the most recent (and least-tested) code for obvious bugs, and I found one. I made an error parsing one of the data files, which had the effect of double-counting the total reads for each subject. I fixed that, but it didn't help the results much.

To start the debugging process, I am looking at the various places where my predictions could go wrong:

0) The data could be wrong. In particular, I assume that the rarefacted data I got from one data file is consistent with the complete dataset I got from another.

1) The posterior distribution could be right, but the predictive distribution could be wrong.

2) The posterior distribution might be wrong because of modeling errors.

3) The posterior distribution might be wrong because of implementation errors.

To check (0), I used the complete dataset to generate a few re-rarefacted datasets to see if my rarefaction process looks like theirs. It does, so I accept the data, at least for now.

To check (1), I used the posterior distribution to generate a 90% credible interval for the total number of species. Since the number of observed species in the complete dataset is necessarily less than the total number of species, the actual values should fall in or below the CIs, but in fact they often exceed the CIs, meaning that there are just more species than my algorithm expects.

While investigating (1) I discovered one problem. In my prior distribution on the number of species, I was setting the upper bound too low, cutting off some values of n with non-negligible probability. So I cranked it up high enough that any additional increase has no further effect on the results. That helps, but my predictions are still too low.

The next step is to test (2). I will generate simulated datasets, generate predictions and then validation them. Since the simulated data come straight from the model, there can be no modeling errors. If I can't validate on simulated data, the problem has to be the algorithm or the implementation, not the model.

Of course, I should have done all this first, before blowing my testing data.

I won't have a chance to get back to this for a little while, but I'll update this post when I do.

Belly Button Biodiversity: Part Three

2013-02-18T11:35:00.000-08:00

This is part three of a series of articles about a Bayesian solution to the Unseen Species problem, applied to data from the Belly Button Biodiversity project.

In Part One I presented the simplest version of the algorithm, which I think is easy to understand, but slow. In Think Bayes I present some ways to optimize it. In Part Two I apply the algorithm to real data and generate predictive distributions. Now in Part Three, as promised, I publish the predictions the algorithm generates. In Part Four I will compare the predictions to actual data.

Background: Belly Button Biodiversity 2.0 (BBB2) is a nation-wide citizen science project with the goal of identifying bacterial species that can be found in human navels (http://bbdata.yourwildlife.org).

Transparent science

In an effort to explore the limits of transparent science, I have started publishing my research in this blog as I go along. This past summer I wrote a series of articles exploring the relationship between Internet use and religious disaffiliation. This "publish as you go" model should help keep researchers honest. Among other things, it might mitigate publication bias due to the "file drawer effect." And if the data and code are published along with the results, that should help make experiments more reproducible.

Toward that end, I will now subject myself to public humiliation by generating a set of predictions using my almost-entirely-unvalidated solution to the Unseen Species problem. In the next installment I will publish the correct answers and score my predictions. Here are the details:

1) I am working with data from the Belly Button Biodiversity project; this data was used in a paper published in PLOS ONE and made available on the web pages of the journal and the researchers. The data consists of rDNA "reads" from 60 subjects. In order to facilitate comparisons between subjects, the researchers chose subjects with at least 400 reads, and for each subject they chose a random subset of 400 reads. The data for the other reads was not published.

2) For each subject, I know the results of the 400 selected reads, and the total number of reads. I will use my algorithm to generate a "prediction" for each subject, which is the number of additional species in the complete dataset.

3) Specifically, for each subject I will generate 9 posterior credible intervals (CIs) for the number of additional species: a 10% CI, a 20% CI, and so on up to a 90% CI.

4) To validate my predictions, I will count the number of CIs that contain the actual, correct value. Ideally, 10% of the correct values should fall in the 10% CIs, 20% should fall in the 20% CIs, and so on. Since the predictions and actual values are integers, a value that hits one end of a predicted CI counts as a half-hit.

Predictions

And now, without further ado, here are my predictions. The columns labelled 10, 20, etc. are 10% credible intervals, 20% CIs, and so on.

Code	# reads	# species	10	20	30	40	50	60	70	80	90
B1234	1392	48	(4, 4)	(4, 5)	(3, 5)	(3, 5)	(3, 6)	(2, 6)	(2, 7)	(1, 7)	(1, 9)
B1235	2452	69	(11, 12)	(10, 12)	(10, 13)	(9, 13)	(8, 14)	(8, 15)	(7, 16)	(6, 17)	(5, 19)
B1236	2964	45	(4, 5)	(4, 5)	(4, 6)	(4, 6)	(3, 7)	(3, 7)	(3, 8)	(2, 9)	(1, 10)
B1237	3090	62	(9, 10)	(9, 11)	(8, 11)	(8, 11)	(7, 12)	(7, 12)	(6, 13)	(5, 14)	(4, 16)
B1242	3056	61	(9, 9)	(8, 10)	(8, 10)	(7, 11)	(7, 11)	(6, 12)	(6, 14)	(5, 15)	(5, 16)
B1243	1518	71	(10, 11)	(10, 12)	(9, 12)	(8, 13)	(8, 13)	(8, 14)	(7, 15)	(6, 16)	(5, 17)
B1246	4230	91	(23, 24)	(22, 25)	(21, 26)	(21, 27)	(19, 28)	(18, 29)	(17, 30)	(16, 33)	(14, 35)
B1253	1928	86	(16, 17)	(15, 18)	(14, 18)	(14, 20)	(13, 20)	(13, 21)	(12, 23)	(11, 24)	(10, 26)
B1254	918	58	(5, 5)	(4, 6)	(4, 6)	(3, 6)	(3, 7)	(3, 7)	(2, 8)	(2, 9)	(1, 10)
B1258	1350	87	(15, 16)	(14, 17)	(14, 17)	(13, 18)	(12, 19)	(11, 19)	(11, 20)	(10, 21)	(8, 24)
B1259	1002	80	(10, 11)	(10, 12)	(10, 12)	(9, 13)	(9, 14)	(8, 14)	(7, 15)	(6, 16)	(6, 18)
B1260	1944	96	(22, 23)	(21, 24)	(20, 25)	(19, 25)	(19, 26)	(18, 27)	(17, 29)	(15, 30)	(14, 32)
B1264	1122	83	(12, 13)	(12, 14)	(11, 14)	(10, 15)	(10, 15)	(9, 16)	(8, 17)	(7, 18)	(6, 20)
B1265	2928	85	(18, 19)	(17, 20)	(16, 21)	(16, 22)	(15, 23)	(14, 24)	(13, 25)	(12, 26)	(11, 28)
B1273	2980	61	(9, 9)	(8, 10)	(8, 10)	(7, 11)	(7, 12)	(6, 12)	(6, 13)	(5, 14)	(4, 16)
B1275	1672	85	(16, 17)	(15, 18)	(15, 19)	(14, 19)	(13, 20)	(13, 21)	(12, 22)	(11, 24)	(9, 25)
B1278	1242	47	(4, 4)	(3, 4)	(3, 5)	(3, 5)	(2, 6)	(2, 6)	(2, 6)	(2, 7)	(1, 8)
B1280	1772	46	(4, 4)	(4, 5)	(3, 5)	(3, 5)	(3, 6)	(2, 6)	(2, 7)	(2, 8)	(1, 9)
B1282	1132	67	(8, 9)	(7, 9)	(7, 10)	(6, 10)	(6, 11)	(6, 11)	(5, 12)	(5, 13)	(3, 15)
B1283	1414	67	(8, 9)	(8, 10)	(7, 10)	(7, 11)	(7, 11)	(6, 12)	(5, 13)	(4, 14)	(3, 16)
B1284	1158	91	(15, 16)	(14, 17)	(14, 17)	(13, 18)	(13, 19)	(12, 19)	(12, 20)	(10, 21)	(9, 23)
B1285	2340	55	(7, 7)	(6, 8)	(6, 8)	(5, 8)	(5, 9)	(4, 9)	(4, 10)	(3, 12)	(2, 13)
B1286	1728	66	(9, 10)	(9, 11)	(8, 11)	(8, 12)	(8, 12)	(7, 13)	(6, 14)	(6, 14)	(4, 16)
B1288	1280	107	(23, 24)	(22, 25)	(21, 25)	(21, 26)	(20, 27)	(19, 27)	(18, 29)	(17, 31)	(15, 32)
B1289	2054	103	(26, 27)	(25, 28)	(24, 29)	(23, 30)	(23, 30)	(22, 32)	(21, 33)	(20, 34)	(17, 36)
B1291	1248	94	(17, 18)	(16, 19)	(16, 20)	(15, 20)	(15, 21)	(13, 22)	(13, 23)	(12, 25)	(10, 27)
B1292	1864	82	(15, 16)	(14, 16)	(13, 17)	(13, 18)	(13, 19)	(12, 20)	(11, 21)	(10, 22)	(9, 24)
B1293	1904	76	(13, 14)	(12, 14)	(12, 15)	(11, 16)	(11, 16)	(10, 17)	(9, 18)	(8, 19)	(7, 22)
B1294	1784	78	(14, 15)	(13, 16)	(12, 16)	(12, 17)	(11, 18)	(11, 19)	(10, 19)	(9, 20)	(8, 23)
B1295	1408	70	(10, 10)	(9, 11)	(9, 12)	(8, 12)	(8, 12)	(7, 13)	(7, 14)	(6, 15)	(4, 17)
B1296	2034	55	(7, 7)	(6, 8)	(6, 8)	(6, 8)	(5, 9)	(4, 9)	(4, 10)	(4, 11)	(3, 12)
B1298	1478	72	(10, 11)	(9, 12)	(9, 12)	(9, 13)	(8, 13)	(8, 14)	(7, 15)	(6, 16)	(5, 18)
B1308	1160	58	(6, 6)	(5, 7)	(5, 7)	(5, 7)	(4, 8)	(4, 8)	(3, 9)	(3, 10)	(2, 11)
B1310	1066	80	(11, 12)	(11, 13)	(10, 13)	(9, 14)	(9, 15)	(8, 15)	(7, 16)	(7, 17)	(5, 19)
B1374	2364	48	(5, 5)	(4, 6)	(4, 6)	(4, 6)	(3, 7)	(3, 7)	(3, 8)	(2, 9)	(2, 10)
B940	2874	93	(22, 24)	(21, 25)	(21, 25)	(20, 26)	(19, 27)	(19, 28)	(18, 30)	(16, 32)	(14, 33)
B941	2154	48	(5, 6)	(5, 6)	(4, 6)	(4, 7)	(4, 7)	(3, 7)	(3, 8)	(2, 9)	(2, 11)
B944	954	52	(4, 4)	(4, 5)	(3, 5)	(3, 5)	(3, 6)	(2, 6)	(2, 6)	(2, 7)	(1, 9)
B945	2390	67	(10, 11)	(10, 12)	(9, 12)	(9, 13)	(8, 13)	(8, 14)	(7, 15)	(7, 16)	(5, 17)
B946	5012	85	(20, 21)	(19, 22)	(19, 23)	(18, 24)	(18, 24)	(17, 26)	(16, 27)	(15, 28)	(12, 31)
B947	3356	62	(10, 11)	(9, 11)	(9, 12)	(8, 12)	(7, 13)	(7, 14)	(6, 14)	(5, 15)	(5, 17)
B948	2384	80	(16, 17)	(15, 18)	(14, 18)	(14, 19)	(13, 20)	(12, 21)	(11, 22)	(10, 23)	(9, 25)
B950	1560	63	(8, 9)	(8, 10)	(8, 10)	(7, 10)	(7, 11)	(6, 11)	(5, 12)	(5, 13)	(4, 15)
B952	1648	57	(7, 7)	(6, 8)	(6, 8)	(6, 8)	(5, 9)	(5, 9)	(4, 10)	(3, 11)	(3, 12)
B953	1452	32	(2, 2)	(1, 2)	(1, 2)	(1, 3)	(1, 3)	(1, 3)	(0, 3)	(0, 4)	(0, 5)
B954	1996	29	(2, 2)	(1, 2)	(1, 2)	(1, 2)	(1, 3)	(1, 3)	(0, 3)	(0, 4)	(0, 4)
B955	1474	65	(8, 9)	(8, 9)	(7, 9)	(7, 10)	(7, 10)	(6, 11)	(5, 12)	(5, 13)	(4, 14)
B956	1482	71	(10, 11)	(10, 12)	(10, 12)	(9, 13)	(8, 14)	(8, 14)	(7, 15)	(6, 16)	(5, 18)
B957	2604	36	(3, 3)	(3, 3)	(2, 4)	(2, 4)	(2, 5)	(1, 5)	(1, 6)	(1, 6)	(1, 7)
B958	2840	29	(2, 2)	(2, 2)	(1, 2)	(1, 3)	(1, 3)	(1, 3)	(1, 4)	(0, 4)	(0, 5)
B961	1214	36	(2, 3)	(2, 3)	(2, 3)	(2, 4)	(1, 4)	(1, 4)	(1, 5)	(1, 5)	(0, 6)
B962	1138	41	(3, 3)	(2, 3)	(2, 3)	(2, 4)	(2, 4)	(1, 4)	(1, 5)	(1, 6)	(0, 7)
B963	1600	71	(10, 12)	(10, 12)	(9, 12)	(9, 13)	(8, 14)	(8, 14)	(7, 15)	(5, 16)	(4, 19)
B966	1950	80	(15, 16)	(14, 16)	(14, 17)	(13, 17)	(12, 18)	(11, 18)	(11, 20)	(10, 22)	(9, 23)
B967	1108	47	(3, 4)	(3, 4)	(3, 4)	(2, 5)	(2, 5)	(2, 5)	(2, 6)	(1, 7)	(1, 7)
B968	2432	52	(6, 7)	(6, 7)	(6, 7)	(5, 8)	(5, 8)	(4, 9)	(3, 9)	(3, 10)	(2, 11)
B971	1462	49	(4, 5)	(4, 5)	(4, 5)	(3, 6)	(3, 6)	(3, 7)	(2, 7)	(2, 8)	(1, 9)
B972	1438	75	(11, 12)	(11, 13)	(10, 13)	(10, 14)	(9, 14)	(8, 15)	(8, 16)	(7, 17)	(6, 18)
B974	5072	54	(7, 8)	(7, 9)	(6, 9)	(6, 9)	(5, 10)	(5, 10)	(4, 11)	(4, 12)	(2, 14)
B975	1542	63	(7, 8)	(7, 9)	(6, 9)	(6, 9)	(6, 10)	(5, 11)	(5, 11)	(4, 12)	(3, 14)

Reading the last row, subject B975 yielded 1542 reads; in a random subset of 400 reads, there were 63 different species (or more precisely, OTUs). My algorithm predicts that if we look at all 1542 reads, the number of additional species we'll find is between 3 and 14, with 90% confidence.

I have to say that this table fills me with dread. The intervals seem quite small, which is to say that the algorithm is more confident than I am. The 90% CIs seem especially narrow to me; it is hard for me to believe that 90% of them will contain the correct values. Well, I guess that's why Karl Popper called them "bold hypotheses". We'll find out soon whether they are bold, or just reckless.

I want to thank Rob Dunn at BBB2 for his help with this project. The code and data I used to generate these results are available from this SVN repository.

EDIT 2-22-13: I ran the predictions again with more simulations. The results are not substantially different. I still haven't looked at the answers.

Belly Button Biodiversity: Part Two

2013-02-08T09:44:00.001-08:00

This is part two of a series of articles about a Bayesian solution to the Unseen Species problem, applied to data from the Belly Button Biodiversity project.

In Part One I presented the simplest version of the algorithm, which I think is easy to understand, but slow. In Think Bayes I present some ways to optimize it. Now in Part Two I apply the algorithm to real data and generate predictive distributions. In Part Three I will publish the predictions the algorithm generates, and in Part Four I will compare the predictions to actual data.

Background: Belly Button Biodiversity 2.0 (BBB2) is a nation-wide citizen science project with the goal of identifying bacterial species that can be found in human navels (http://bbdata.yourwildlife.org).

The belly button data

To get a sense of what the data look like, consider subject B1242, whose sample of 400 reads yielded 61 species with the following counts:

92, 53, 47, 38, 15, 14, 12, 10, 8, 7, 7, 5, 5, 
4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

There are a few dominant species that make up a substantial fraction of the whole, but many species that yielded only a single read. The number of these “singletons” suggests that there are likely to be at least a few unseen species.

In the example with lions and tigers, we assume that each animal in the preserve is equally likely to be observed. Similarly, for the belly button data, we assume that each bacterium is equally likely to yield a read.

In reality, it is possible that each step in the data-collection process might introduce consistent biases. Some species might be more likely to be picked up by a swab, or to yield identifiable amplicons. So when we talk about the prevalence of each species, we should remember this source of error.

I should also acknowledge that I am using the term “species” loosely. First, bacterial species are not well-defined. Second, some reads identify a particular species, others only identify a genus. To be more precise, I should say “operational taxonomic unit”, or OTU.

Now let’s process some of the belly button data. I defined a class called Subject to represent information about each subject in the study:

class Subject(object):

    def __init__(self, code):
        self.code = code
        self.species = []

Each subject has a string code, like “B1242”, and a list of (count, species name) pairs, sorted in increasing order by count. Subject provides several methods to make it easy to these counts and species names. You can see the details in http://thinkbayes.com/species.py.

Figure 12.3: Distribution of n for subject B1242.

In addition, Subject.Process creates a suite, specifically a suite of type Species5, which represents the distribution of n and the prevalences after processing the data.

It also provides PlotDistOfN, which plots the posterior distribution of n. Figure 12.3 shows this distribution for subject B1242. The probability that there are exactly 61 species, and no unseen species, is nearly zero. The most likely value is 72, with 90% credible interval 66 to 79. At the high end, it is unlikely that there are as many as 87 species.

Next we compute the posterior distribution of prevalence for each species. Species2 providesDistOfPrevalence:

# class Species2

    def DistOfPrevalence(self, index):
        pmfs = thinkbayes.Pmf()

        for n, prob in zip(self.ns, self.probs):
            beta = self.MarginalBeta(n, index)
            pmf = beta.MakePmf()
            pmfs.Set(pmf, prob)

        mix = thinkbayes.MakeMixture(pmfs)
        return pmfs, mix

index indicates which species we want. For each value of n, we have a different posterior distribution of prevalence.

Figure 12.4: Distribution of prevalences for subject B1242.

So the loop iterates through the possible values of n and their probabilities. For each value of n it gets a Beta object representing the marginal distribution for the indicated species. Remember that Beta objects contain the parameters alpha and beta; they don’t have values and probabilities like a Pmf, but they provide MakePmf which generates a discrete approximation to the continuous beta distribution.

pmfs is a MetaPmf that contains the distributions of prevalence, conditioned on n. MakeMixturecombines the MetaPmf into mix, which combines the conditional distributions into the answer, a single distribution of prevalence.

Figure 12.4 shows these distributions for the five species with the most reads. The most prevalent species accounts for 23% of the 400 reads, but since there are almost certainly unseen species, the most likely estimate for its prevalence is 20%, with 90% credible interval between 17% and 23%.

Predictive distributions

Figure 12.5: Simulated rarefaction curves for subject B1242.

I introduced the hidden species problem in the form of four related questions. We have answered the first two by computing the posterior distribution for n and the prevalence of each species.

The other two questions are:

If we are planning to collect additional samples, can we predict how many new species we are likely to discover?
How many additional reads are needed to increase the fraction of observed species to a given threshold?

To answer predictive questions like this we can use the posterior distributions to simulate possible future events and compute predictive distributions for the number of species, and fraction of the total, we are likely to see.

The kernel of these simulations looks like this:

Choose n from its posterior distribution.
Choose a prevalence for each species, including possible unseen species, using the Dirichlet distribution.
Generate a random sequence of future observations.
Compute the number of new species, num_new, as a function of the number of additional samples, k.
Repeat the previous steps and accumulate the joint distribution of num_new and k.

And here’s the code. RunSimulation runs a single simulation:

# class Subject

    def RunSimulation(self, num_samples):
        m, seen = self.GetSeenSpecies()
        n, observations = self.GenerateObservations(num_samples)

        curve = []
        for k, obs in enumerate(observations):
            seen.add(obs)

            num_new = len(seen) - m
            curve.append((k+1, num_new))

        return curve

num_samples is the number of additional samples to simulate. m is the number of seen species, and seen is a set of strings with a unique name for each species. n is a random value from the posterior distribution, and observations is a random sequence of species names.

The result of RunSimulation is a “rarefaction curve”, represented as a list of pairs with the number of samples and the number of new species seen.

Before we see the results, let’s look at GetSeenSpecies and GenerateObservations.

#class Subject

    def GetSeenSpecies(self):
        names = self.GetNames()
        m = len(names)
        seen = set(SpeciesGenerator(names, m))
        return m, seen

GetNames returns the list of species names that appear in the data files, but for many subjects these names are not unique. So I use SpeciesGenerator to extend each name with a serial number:

def SpeciesGenerator(names, num):
    i = 0
    for name in names:
        yield '%s-%d' % (name, i)
        i += 1

    while i < num:
        yield 'unseen-%d' % i
        i += 1

Given a name like Corynebacterium, SpeciesGenerator yields Corynebacterium-1. When the list of names is exhausted, it yields names like unseen-62.

Here is GenerateObservations:

# class Subject

    def GenerateObservations(self, num_samples):
        n, prevalences = self.suite.Sample()

        names = self.GetNames()
        name_iter = SpeciesGenerator(names, n)

        d = dict(zip(name_iter, prevalences))
        cdf = thinkbayes.MakeCdfFromDict(d)
        observations = cdf.Sample(num_samples)

        return n, observations

Again, num_samples is the number of additional samples to generate. n and prevalences are samples from the posterior distribution.

cdf is a Cdf object that maps species names, including the unseen, to cumulative probabilities. Using a Cdf makes it efficient to generate a random sequence of species names.

Finally, here is Species2.Sample:

    def Sample(self):
        pmf = self.DistOfN()
        n = pmf.Random()
        prevalences = self.SampleConditional(n)
        return n, prevalences

And SampleConditional, which generates a sample of prevalences conditioned on n:

# class Species2

    def SampleConditional(self, n):
        params = self.params[:n]
        gammas = numpy.random.gamma(params)
        gammas /= gammas.sum()
        return gammas

We saw this algorithm for generating prevalences previously in Species2.SampleLikelihood.

Figure 12.5 shows 100 simulated rarefaction curves for subject B1242. I shifted each curve by a random offset so they would not all overlap. By inspection we can estimate that after 400 more samples we are likely to find 2–6 new species.

Joint posterior

Figure 12.6: Distributions of the number of new species conditioned on the number of additional samples.

To be more precise, we can use the simulations to estimate the joint distribution of num_new andk, and from that we can get the distribution of num_new conditioned on any value of k.

# class Subject

    def MakeJointPredictive(self, curves):
        joint = thinkbayes.Joint()
        for curve in curves:
            for k, num_new in curve:
                joint.Incr((k, num_new))
        joint.Normalize()
        return joint

MakeJointPredictive makes a Joint object, which is a Pmf whose values are tuples.

curves is a list of rarefaction curves created by RunSimulation. Each curve contains a list of pairs of k and num_new.

The resulting joint distribution is a map from each pair to its probability of occurring. Given the joint distribution, we can get the distribution of num_new conditioned on k:

# class Joint

    def Conditional(self, i, j, val):
        pmf = Pmf()
        for vs, prob in self.Items():
            if vs[j] != val: continue
            pmf.Incr(vs[i], prob)

        pmf.Normalize()
        return pmf

i is the index of the variable whose distribution we want; j is the index of the conditional variables, and val is the value the jth variable has to have. You can think of this operation as taking vertical slices out of Figure 12.5.

Subject.MakeConditionals takes a list of ks and computes the conditional distribution of num_new for each k. The result is a list of Cdf objects.

# class Subject

    def MakeConditionals(self, curves, ks):
        joint = self.MakeJointPredictive(curves)

        cdfs = []
        for k in ks:
            pmf = joint.Conditional(1, 0, k)
            pmf.name = 'k=%d' % k
            cdf = thinkbayes.MakeCdfFromPmf(pmf)
            cdfs.append(cdf)

        return cdfs

Figure 12.6 shows the results. After 100 samples, the median predicted number of new species is 2; the 90% credible interval is 0 to 5. After 800 samples, we expect to see 3 to 12 new species.

Coverage

Figure 12.7: Complementary CDF of coverage for a range of additional samples.

The last question we want to answer is, “How many additional reads are needed to increase the fraction of observed species to a given threshold?”

To answer this question, we’ll need a version of RunSimulation that computes the fraction of observed species rather than the number of new species.

# class Subject

    def RunSimulation(self, num_samples):
        m, seen = self.GetSeenSpecies()
        n, observations = self.GenerateObservations(num_samples)

        curve = []
        for k, obs in enumerate(observations):
            seen.add(obs)

            frac_seen = len(seen) / float(n)
            curve.append((k+1, frac_seen))

        return curve

Next we loop through each curve and make a dictionary, d, that maps from the number of additional samples, k, to a list of fracs; that is, a list of values for the coverage achieved after ksamples.

    def MakeFracCdfs(self, curves):
        d = {}
        for curve in curves:
            for k, frac in curve:
                d.setdefault(k, []).append(frac)

        cdfs = {}
        for k, fracs in d.iteritems():
            cdf = thinkbayes.MakeCdfFromList(fracs)
            cdfs[k] = cdf

        return cdfs

Then for each value of k we make a Cdf of fracs; this Cdf represents the distribution of coverage after k samples.

Remember that the CDF tells you the probability of falling below a given threshold, so thecomplementary CDF tells you the probability of exceeding it. Figure 12.7 shows complementary CDFs for a range of values of k.

To read this figure, select the level of coverage you want to achieve along the x-axis. As an example, choose 90%.

Now you can read up the chart to find the probability of achieving 90% coverage after k samples. For example, with 300 samples, you have about a 60% of getting 90% coverage. With 700 samples, you have a 90% chance of getting 90% coverage.

With that, we have answered the four questions that make up the unseen species problem. Next time: validation!

Belly Button Biodiversity: Part One

2013-02-05T07:03:00.000-08:00

This post is a excerpt from Think Bayes: Bayesian Statistics Made Simple, the book I am working on now. You can read the entire current draft at http://thinkbayes.com.

Belly button bacteria

Belly Button Biodiversity 2.0 (BBB2) is a nation-wide citizen science project with the goal of identifying bacterial species that can be found in human navels (http://bbdata.yourwildlife.org).

The project might seem whimsical, but it is part of an increasing interest in the human microbiome, the set of microorganisms that live on human skin and other surfaces that contact the environment.

In their pilot study, BBB2 researchers collected swabs from the navels of 60 volunteers, used multiplex pyrosequencing to extract and sequence fragments of 16S rDNA, then identified the species or genus the fragments came from. Each identified fragment is called a “read.”

We can use these data to answer several related questions:

Based on the number of species observed, can we estimate the total number of species in the environment?
Can we estimate the prevalence of each species; that is, the fraction of the total population belonging to each species?
If we are planning to collect additional samples, can we predict how many new species we are likely to discover?
How many additional reads are needed to increase the fraction of observed species to a given threshold?

These questions make up what is called the “unseen species problem.”

Lions and tigers and bears

I’ll start with a simplified version of the problem where we know that there are exactly three species. Let’s call them lions, tigers and bears. Suppose we visit a wild animal preserve and see 3 lions, 2 tigers and one bear.

If we have an equal chance of observing any animal in the preserve then the number of each species we see is governed by the multinomial distribution. If the prevalence of lions and tigers and bears is p_lion and p_tiger and p_bear, the likelihood of seeing 3 lions, 2 tigers and one bear is

p_lion**3 * p_tiger**2 * p_bear**1

An approach that is tempting, but not correct, is to use beta distributions, as in Section 4.6, to describe the prevalence of each species separately. For example, we saw 3 lions and 3 non-lions; if we think of that as 3 “heads” and 3 “tails,” then the posterior distribution of p_lion is:

    beta = thinkbayes.Beta()
    beta.Update((3, 3))
    print beta.MaximumLikelihood()

The maximum likelihood estimate for p_lion is the observed rate, 50%. Similarly the MLEs for p_tiger and p_bear are 33% and 17%.
But there are two problems:

We have implicitly used a prior for each species that is uniform from 0 to 1, but since we know that there are three species, that prior is not correct. The right prior should have a mean of 1/3, and there should be zero likelihood that any species has a prevalence of 100%.
The distributions for each species are not independent, because the prevalences have to add up to 1. To capture this dependence, we need a joint distribution for the three prevalences.

We can use a Dirichlet distribution to solve both of these problems (see http://en.wikipedia.org/wiki/Dirichlet_distribution). In the same way we used the beta distribution to describe the distribution of bias for a coin, we can use a Dirichlet distribution to describe the joint distribution of p_lion, p_tiger and p_bear.

The Dirichlet distribution is the multi-dimensional generalization of the beta distribution. Instead of two possible outcomes, like heads and tails, the Dirichlet distribution handles any number of outcomes: in this example, three species.

If there are n outcomes, the Dirichlet distribution is described by n parameters, written α_i.
Here’s the definition, from thinkbayes.py, of a class that represents a Dirichlet distribution:

class Dirichlet(object):

    def __init__(self, n):
        self.n = n
        self.params = numpy.ones(n, dtype=numpy.int)

n is the number of dimensions; initially the parameters are all 1. I use a numpy array to store the parameters so I can take advantage of array operations.
Given a Dirichlet distribution, the marginal distribution for each prevalence is a beta distribution, which we can compute like this:

    def MarginalBeta(self, i):
        alpha0 = self.params.sum()
        alpha = self.params[i]
        return Beta(alpha, alpha0-alpha)

i is the index of the marginal distribution we want. alpha0 is the sum of the parameters; alpha is the parameter for the given species.
In the example, the prior marginal distribution for each species is Beta(1, 2). We can compute the prior means like this:

    dirichlet = thinkbayes.Dirichlet(3)
    for i in range(3):
        beta = dirichlet.MarginalBeta(i)
        print beta.Mean()

As expected, the prior mean prevalence for each species is 1/3.

To update the Dirichlet distribution, we add the number of observations to each parameter, like this:

    def Update(self, data):
        m = len(data)
        self.params[:m] += data

Here data is a sequence of counts in the same order as params, so in this example, it should be the number of lions, tigers and bears.

But data can be shorter than params; in that case there are some hypothetical species that have not been observed.

Here’s code that updates dirichlet with the observed data and computes the posterior marginal distributions.

    data = [3, 2, 1]
    dirichlet.Update(data)

    for i in range(3):
        beta = dirichlet.MarginalBeta(i)
        pmf = beta.MakePmf()
        print i, pmf.Mean()

This figure shows the results. The posterior mean prevalences are 44%, 33% and 22%.

A hierarchical model

We have solved a simplified version of the problem: if we know how many species there are, we can estimate the prevalence of each.

Now let’s get back to the original problem, estimating the total number of species. To solve this problem I’ll define a metasuite, which is a Suite that contains other Suites as hypotheses. In this case, the top-level Suite contains hypotheses about the number of species; the bottom level contains hypotheses about prevalences. A multi-level model like this is called “hierarchical.”
Here’s the class definition:

class Species(thinkbayes.Suite):

    def __init__(self, ns):
        hypos = [thinkbayes.Dirichlet(n) for n in ns]
        thinkbayes.Suite.__init__(self, hypos)

__init__ takes a list of possible values for n and makes a list of Dirichlet objects.

Here’s the code that creates the top-level suite:

    ns = range(3, 30)
    suite = Species(ns)

ns is the list of possible values for n. We have seen 3 species, so there have to be at least that many. I chose an upper bound that seemed reasonable, but we will have to check later that the probability of exceeding this bound is low. And at least initially we assume that any value in this range is equally likely.

To update a hierarchical model, you have to update all levels. Sometimes it is necessary or more efficient to update the bottom level first and work up. In this case it doesn’t matter, so I update the top level first:

#class Species

    def Update(self, data):
        thinkbayes.Suite.Update(self, data)
        for hypo in self.Values():
            hypo.Update(data)

Species.Update invokes Update in the parent class, then loops through the sub-hypotheses and updates them.

Now all we need is a likelihood function. As usual, Likelihood gets a hypothesis and a dataset as arguments:

# class Species

    def Likelihood(self, hypo, data):
        dirichlet = hypo
        like = 0
        for i in range(1000):
            like += dirichlet.Likelihood(data)

        return like

hypo is a Dirichlet object; data is a sequence of observed counts. Species.Likelihood calls Dirichlet.Likelihood 1000 times and returns the total.

Why do we have to call it 1000 times? Because Dirichlet.Likelihood doesn’t actually compute the likelihood of the data under the whole Dirichlet distribution. Instead, it draws one sample from the hypothetical distribution and computes the likelihood of the data under the sampled set of prevalences.
Here’s what it looks like:

# class Dirichlet

    def Likelihood(self, data):
        m = len(data)
        if self.n < m:
            return 0

        x = data
        p = self.Random()
        q = p[:m]**x
        return q.prod()

The length of data is the number of species observed. If we see more species than we thought existed, the likelihood is 0.

Otherwise we select a random set of prevalences, p, and compute the multinomial PDF, which is

c(x)

∏

p_i^x_i

p_i is the prevalence of the ith species, and x_i is the observed number. The first term, c(x), is the multinomial coefficient; I left it out of the computation because it is a multiplicative factor that depends only on the data, not the hypothesis, so it gets normalized away (see http://en.wikipedia.org/wiki/Multinomial_distribution).

Also, I truncated p at m, which is the number of observed species. For the unseen species, x_i is 0, so p_i^x_i is 1, so we can leave them out of the product.

Random sampling

There are two ways to generate a random sample from a Dirichlet distribution. One is to use the marginal beta distributions, but in that case you have to select one at a time and scale the rest so they add up to 1 (see http://en.wikipedia.org/wiki/Dirichlet_distribution#Random_number_generation).

A less obvious, but faster, way is to select values from n gamma distributions, then normalize by dividing through by the total. Here’s the code:

# class Dirichlet

    def Random(self):
        p = numpy.random.gamma(self.params)
        return p / p.sum()

Now we’re ready to look at some results. Here is the code that updates the top-level suite and extracts the posterior PMF of n:

    data = [3, 2, 1]
    suite.Update(data)
    pmf = suite.DistOfN()

To get the posterior distribution of n, DistOfN iterates through the top-level hypotheses:

    def DistOfN(self):
        pmf = thinkbayes.Pmf()
        for hypo, prob in self.Items():
            pmf.Set(hypo.n, prob)
        return pmf

This figure shows the result:

The most likely value is 5. Values from 3 to 8 are all likely; after that the probabilities drop off quickly. The probability that there are 29 species is low enough to be negligible; if we chose a higher bound, we would get the same result.

But remember that we started with a uniform prior for n. If we have background information about the number of species in the environment, we might choose a different prior.

Optimization

I have to admit that I am proud of this example. The unseen species problem is not easy, and I think this solution is simple and clear, and takes surprisingly few lines of code (about 50 so far).

The only problem is that it is slow. It’s good enough for the example with only 3 observed species, but not good enough for the belly button data, with more than 100 species in some samples.

In Think Bayes I present a series of optimizations we need to make this solution scale. Here’s a road map of the steps:

The first step is to recognize that if we update the Dirichlet distributions with the same data, the first m parameters are the same for all of them. The only difference is the number of hypothetical unseen species. So we don’t really need n Dirichlet objects; we can store the parameters in the top level of the hierarchy. Species2 implements this optimization.
Species2 also uses the same set of random values for all of the hypotheses. This saves time generating random values, but it has a second benefit that turns out to be more important: by giving all hypothesis the same selection from the sample space, we make the comparison between the hypotheses more fair, so it takes fewer iterations to converge.
But there is still a major performance problem. As the number of observed species increases, the array of random prevalences gets bigger, and the chance of choosing one that is approximately right becomes small. So the vast majority of iterations yield small likelihoods that don’t contribute much to the total, and don’t discriminate between hypotheses.The solution is to do the updates one species at a time. Species4 is a simple implementation of this strategy using Dirichlet objects to represent the sub-hypotheses.
Finally, Species5 combines the sub-hypothesis into the top level and uses numpy array operations to speed things up.

I won't present the details here. In the next part of this series, I will present results from the Belly Button Biodiversity project.

Bayesian Statistics Made Simple

2013-01-23T12:40:00.003-08:00

I am happy to announce that I will offer an updated and revised version of my tutorial, Bayesian Statistics Made Simple, at PyCon 2013.

Registration is open now. Here are the details:

PyCon 2013
Santa Clara, CA

Wednesday 13 March, 1:20 p.m.–4:40 p.m.

Bayesian statistics made simple
Allen Downey

Audience level: Intermediate

DESCRIPTION

An introduction to Bayesian statistics using Python. Bayesian statistics are usually presented mathematically, but many of the ideas are easier to understand computationally. People who know some Python have a head start.

We will use material from Think Stats: Probability and Statistics for Programmers (O’Reilly Media), and Think Bayes, a work in progress at http://thinkbayes.com.

ABSTRACT

Bayesian statistical methods are becoming more common and more important, but there are not many resources to help beginners get started. People who know Python can use their programming skills to get a head start.

I will present simple programs that demonstrate the concepts of Bayesian statistics, and apply them to a range of example problems. Participants will work hands-on with example code and practice on example problems.

Students should have at least basic Python and basic statistics. If you learned about Bayes’s Theorem and probability distributions at some time, that’s enough, even if you don’t remember it!

Students should bring a laptop with Python 2.x and matplotlib. You can work in any environment; you just need to be able to download a Python program and run it.

Are first babies more likely to be late, revisited.

2013-01-08T08:04:00.000-08:00

Two years ago I wrote an article called Are first babies more likely to be late?, based on a question that came up when my wife and I were expecting our first child. I compared the pool of first babies to the pool of all other babies, and found:

There is a small difference in the mean pregnancy length for the two groups, about 13 hours, but it is not practically or statistically significant.
If we group babies into Early, On Time, or Late (where On Time is 38, 39 or 40 weeks), first babies are a little more likely to be Early or Late, and less likely to be On Time.

Then yesterday I got the following question from an Unknown correspondent:

While interesting, I can't help but think you need to compare the first and others for the same woman. While may be unlikely it could still be that a tendency exists for a woman's second, third, etc, child comes earlier.

This is an excellent suggestion. It is possible that the variability between people is masking some of the variability between first and later babies. By pairing first and second babies with the same mother, we can control for variation between mothers.

So I ran that experiment, selecting all mothers with at least two children and computing the difference in pregnancy length between the second and first child (so a positive value means the second child was later). Here is the distribution of these value for 4387 women in the NSFG (National Survey of Family Growth):

Visually the distribution looks symmetric, and the summary statistics support that conclusion. The mean is -0.034, which means that (if anything) the second baby is born about 6 hours earlier, but this difference is not statistically significant.

Conclusion: good question, definitely worth running the experiment, but the primary result is the same as what we saw before: no significant difference in the means.

Call for Bayesian case studies

2013-01-07T08:50:00.000-08:00

It's been a while since the last post because I have been hard at work on Think Bayes. As always, I have been posting drafts as I go along, so you can read the current version at thinkbayes.com.

I am teaching Computational Bayesian Statistics in the spring, using the draft edition of the book. The students will work on case studies, some of which will be included in the book. And then I hope the book will be published as part of the Think X series (for all X). At least, that's the plan.

In the next couple of weeks, students will be looking for ideas for case studies. An ideal project has at least some of these characteristics:

An interesting real-world application (preferably not a toy problem).
Data that is either public or can be made available for use in the case study.
Permission to publish the case study!
A problem that lends itself to Bayesian analysis, in particular if there is a practical advantage to generating a posterior distribution rather than a point or interval estimate.

Examples in the book include:

The hockey problem: estimating the rate of goals scored by two hockey teams in order to predict the outcome of a seven-game series.
The paintball problem, a version of the lighthouse problem. This one verges on being a toy problem, but recasting it in the context of paintball got it over the bar for me.
The kidney problem. This one is as real as it gets -- it was prompted by a question posted by a cancer patient who needed a statistical estimate of when a tumor formed.
The unseen species problem: a nice Bayesian solution to a standard problem in ecology.

So far I have a couple of ideas prompted by questions on Reddit:

But I would love to get more ideas. If you have a problem you would like to contribute, let me know!

Secularization in America: part seven

2012-07-11T09:29:00.001-07:00

Abstract

Based on 2000-2010 data from the General Social Survey (GSS), I present results of a logistic regression that measures the relationship between Internet use and religious affiliation, controlling for religious upbringing, income and socioeconomic index, year born (age), and education.

I find that moderate Internet use reduces the chance of religious affiliation by 2 percentage points (odds ratio 0.8); heavier Internet use reduces affiliation by an additional 5 percentage points (odds ratio 0.7). Four years of college reduces affiliation by an additional 2 percentage points (odds ratio 0.8).

All reported effects are statistically significant with N=8960 respondents.

Results of logistic regression can be difficult to interpret; it might help to imagine the following progression:

Start with a hypothetical baseline person raised in any religion, with moderate or high household income ($25,000 per year or more), born in 1960, with high school education but no college, and low Internet use (less than 2 hours per week): in the GSS survey, 91% of people in this category have a religious affiliation. Now we change one variable at a time.
If this person were born 10 years later (in 1970) the fraction would drop to 89%.
If this person went to college, the fraction would drop to 87%
If this person used the Internet 2 or more hours per week, the fraction would drop to 85%.
If this person used the Internet 8 or more hours per week, the fraction would drop to 80%.

Taken together, college education and Internet use are associated with a decrease in religious affiliation of 9 percentage points.

Introduction

From 1990 to 2010 the fraction of Protestants in the U.S. population dropped from 62% to 51%; at the same time the fraction of people with no religious preference increased from 8% to 18%. The following graph shows these trends:

In a previous article I presented evidence that something happened in the 1990s, continuing through the 2000s, that is causing disaffiliation from religion across all generations, with the largest effect on the youngest generations in the survey, people born in the 1960s and 1970s.

There are many possible explanations, but for me, the Internet pops to the top of this list. First, the timing is at least approximately right. Here is data from the World Bank, showing number of Internet users per hundred people in the U.S.

Internet use increased rapidly from 1995 to 2010, which is the interval of steepest change in religious affiliation.

Regressions

To identify factors that contribute to disaffiliation, I ran logistic regressions with the following dependent variable:

has_relig: 1 if the respondent reported any religious affiliation when interviewed as an adult, or 0 if the respondent reported "None" (based on the GSS variable RELIG)

And these explanatory variables:

had_relig: 1 if the respondent reported being raised in a religion, 0 otherwise (based on RELIG16)

born_from_1960: year the respondent was born minus 1960 (based on AGE and survey year). Subtracting 1960 makes it easier to interpret the results of the regression.

educ_from_12: number of years of school completed, minus 12 (based on EDUC).

somewww: 1 if the respondent reported using the Internet 2 of more hours per week, 0 otherwise (based on WWWHR, with the threshold chosen near the median)

heavywww: 1 if the respondent uses the Internet more than 8 hours per week, 0 otherwise (threshold chosen near the 75th percentile)

SEI: Socioeconomic index (a measure of occupational prestige developed by the GSS).

high_income: 1 if the respondent reports annual household income of $25,000 or more, which includes 62% of respondents who answered the question.

I used data from GSS survey years 2000, 2002, 2004, 2006, and 2010 (the relevant questions were not asked in 2008). I excluded respondents who were not asked or did not answer one or more of the questions I used in my analysis.

It turns out that SEI does not make a contribution that is either statistically or practically significant, so I omit it from the model.

Here are the results of the model as reported by R:

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -0.164434 0.094978 -1.731 0.0834 .

had_relig 2.318141 0.087372 26.532 < 2e-16 ***

high_income 0.166673 0.072345 2.304 0.0212 *

born_from_1960 -0.020161 0.002128 -9.474 < 2e-16 ***

educ_from_12 -0.051850 0.012228 -4.240 2.23e-05 ***

somewww -0.178409 0.078490 -2.273 0.0230 *

heavywww -0.336658 0.080546 -4.180 2.92e-05 ***

---

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 7860.3 on 8959 degrees of freedom

Residual deviance: 6872.5 on 8953 degrees of freedom

AIC: 6886.5

Number of Fisher Scoring iterations: 5

All explanatory variables are statistically significant: high_income and somewww are borderline, both at p=0.02.

The odds ratios and cumulative probabilities are:

odds cumulative

ratio probability

(Intercept) 0.85 46

had_relig 10.16 90

high_income 1.18 91

born_from_1960 0.82 89

educ_from_12 0.81 87

somewww 0.84 85

heavywww 0.71 80

These results are summarized and interpreted in the Abstract, above.

Discussion

As always, statistical association does not prove causation, but in this case I think there are reasons to believe that Internet use causes disaffiliation from religion:

It is easy to imagine how Internet use could allow a person in a homogeneous community to find information about people of other religions (and none), and to interact with them personally. And there is anecdotal evidence that those interactions contribute to religious disaffiliation (for example, numerous personal reports on reddit.com/r/atheism).
Conversely it is harder to imagine plausible reasons why disaffiliation might cause increased Internet use (except possibly on Sunday mornings).
Although it is possible that a third factor causes both disaffiliation and Internet use, that factor would also have to be new, coincidentally rising in prevalence, like the Internet, during the 1990s and 2000s.
Whatever causes disaffiliation has the strongest effect on the youngest generations, which is consistent with the hypothesis that Internet use during adolescence and young adulthood has the strongest effect on religious affiliation.

So with appropriate caution, I think there is a strong case here for causation, and not just statistical association.

Furthermore, the magnitude of the effect is large enough to explain a substantial part of the observed changes in religious affiliation. In my next article I will incorporate this regression model into the generational model I presented in Part Six, in order to estimate the effect of Internet use on these trends.

Summary of previous reports

In Part One I described trends in market share of major religions in the U.S.: since 1988, the fraction of Protestants dropped from 60% to 51%, and the fraction of people with no religious affiliation increased from 8% to 18%.

In Part Two I used data from the 1988 General Social Survey (GSS) to model transmission of religion from parent to child, and found that the model failed to predict the decrease in Protestants and the increase in Nones that occurred between 1988 and 2010.

In Part Three I looked at changes, between 1988 and 2008, in the spouse tables (which describe the tendencies of people to marry within their religions), the environment table (which describes parents' decisions about their children's religious upbringing), and the transmission table (which describes the likely outcomes for children raised within each religion). I found that the transmission table has changed substantially since 1988, and accounts for a large part of the observed increase in Nones, but not the decrease in Protestants.

In Part Four I looked at changes in religiosity over the lifetime of respondents. I tentatively concluded that the differences between generations were larger than changes in affiliation, within generations, over time.

But in Part Five I looked more closely and saw that all generations were becoming more religious, or staying the same, prior to 1990, and that all generations began to disaffiliate during the 1990s, continuing into the 2000s.

In Part Six I presented a generational model that retroactively "predicts" the changes we have seen since 1988, and used it to predict how those changes are likely to continue in the next 30 years. I expect the fraction of Protestants to continue to decrease, and the fraction of Nones to increase and overtake Catholic as the second-largest affiliation by 2030.

Secularization in America: part six

2012-07-10T06:51:00.003-07:00

Summary so far

Generational Model

Now I am ready to get back to the generational model I have been working up to. The goal of the generational model is to separate these three effects:

Changes in religious preference from one generation to the next.
Changes in religious affiliation over the lifetime of respondents.
Changes in the composition of the GSS cohort over time.

The model works by simulation. Assuming that we are starting in 1988, here are the steps:

Read the survey data from 1988 and resample it. Compute and store the distribution of ages.
For each respondent, generate a hypothetical child. Use the BirthModel to determine year of birth, the UpbringingModel to determine what religion the child is raised in, and the TransmissionModel to determine what affiliation the child will have as an adult. Details of these models follow.
Form a combined cohort of parents and simulated children. Since the cohort of parents is a representative sample of the US population, the cohort of simulated children is a representative sample of the population one generation later (based, for now, on the simplifying assumptions that all groups have the same number of children on average, and there is no immigration).
In order to generate a cohort from a future survey year, draw a sample from the combined cohort, weighted so that the distribution of ages in the future year is the same as the original distribution of ages in 1988. As the simulation goes forward in time, this generated cohort contains fewer of the parents and more of the simulated children. After 20 years, about 25% of the "real" respondents have been replaced with "fake" respondents.

Now, where do all these auxiliary models come from?

BirthModel: This is just the distribution of parent's age when each child is born. It is based on data from the 1994 GSS, which includes questions about children. I had to do some work to correct for an obvious bias due to the ages of the respondents; I will skip the details here.

UpbringingModel: This is a combination of the SpouseTable and the EnvironmentTable, described in Part Three. It is a map from the parent's religion to the distribution of possible religions the child might be raised in.

TransmissionModel: This is the TransmissionTable described in Part Three. It is a map from the religious environment of the child to the distribution of religious affiliation reported by the child as an adult.

The Upbringing and Transmission models come in two flavors:

Time invariant: We use all respondents to estimate the parameters of the model, and apply the same model to generate all simulated children.

Time variant: We estimate different parameters for each generation (partitioned by decade born) and use different models to generate simulated children, depending on what year they are born.

For the time variant model, we have to extrapolate from observed data into the future. To keep this simple we copy the latest reliable data (based on sample size) and apply it to people born in later decades.

Ok, that's enough methodology for now. Let's take a look at some...

Results

The first step is to validate the model by showing that it can predict the observed changes using past data. Here I mean "predict" in a peculiar sense, which is that I will use the entire dataset (including data after 1988) to build the auxiliary models, then use the simulator to generate trends from 1988 to 2010.

Here is what the results look like:

The thick lines are the observed data; the thin lines are simulations. Here are my observations:

For Jews and Catholics, the observed data falls within the bounds of the simulations, so the model validates.
For Other, the observed data sometimes exceeds the bounds of the simulations, which may be due to immigration (not included in this model).
For None, the observed data is at the high end of the range, and for Prot it is at the low end of the range. This is most likely due to the disaffiliation we saw in Part Five, which is only partly captured in this model.

I conclude that the model is capturing a large part of the observed changes since 1988, but of course I am cheating by using data from after 1988. So these results validate my modeling decisions (what to include and what to leave out) but they don't test the predictive power of the model.

Predictive power

To make an honest test, we have to restrict ourselves to data from before 1988. That way we can tell what part of the observed changes would have been predictable in 1988.

Here's what the result looks like:

So if we had used this model in 1988, we would have predicted a small decrease in the fraction of Protestants and a small increase in None, but we would have underestimated both trends.

This supports my conclusion in Part Five that something happened in the 1990s that changed trends in religious affiliation, and suggests that these changes were unpredictable based on data observable before 1988.

Predictions

Finally, we can use all data to build the models, use 2010 as the starting place for the simulations, and make some predictions for the next 30 years:

So what should we expect?

The decline in fraction of Protestants will continue. The fraction of Catholics will also decrease, but more slowly.
The fraction of Nones will increase, overtaking Catholics as the second-largest religious affiliation around 2030.
The fraction of Others will increase slowly, about 1 percentage point in 30 years. If immigration from Asia continues at current rates, that would add another percentage point, bringing the total to 6%.
The fraction of Jews will decrease, possibly by half by 2040.

These predictions are likely to be conservative; that is, the rate of secularization will almost certainly be faster. Why?

Over the last several generations, the UpbringingModel and the TransmissionModel have changed substantially. Parents are less likely to raise their children with religion, and those children are less likely to adopt the religion they are raised with. The model captures these trends, but assumes that they will level off in 2010. It would probably be more accurate to assume that they will continue.
Rates of disaffiliation among adults are also increasing. Again, the model includes trends that have already occurred, but it assumes that they will level off rather than continue.

So there are reasons to expect the fraction of Nones to accelerate.

Conversely, it is hard to imagine that the trends will be any slower than these predictions. To a large extent, these results are not predictions about things that will happen in the future; rather, they are the future consequences of things that have already happened. For example, in 2020, the GSS survey will include a cohort of people in their 40s. What will they be like? They will be a lot like the people in the 2010 survey who are in their 30s. But they will be older. Changes in the general population are slow because is takes a long time to replace each generation with the next; but as a result, they are also predictable.

Next time: Was Rick Santorum right? Is college the #1 enemy of religious belief? (Hint: no.) I will look more closely at the TransmissionModel to see what factors make vertical transmission of religion more (or less) likely.

Secularization in America: part five

2012-07-09T09:45:00.000-07:00

Summary so far

Part Four revisited

In Part Four I looked at changes in religiosity over the lifetime of respondents. The GSS is not a longitudinal survey, so we can't follow individuals, but we can follow generations (which I partition by decade of birth) over time.

Last time I presented this figure, which shows religiosity (the fraction of respondents with any religious preference) as a function of respondent's age, partitioned by decade of birth, for people who were raised Protestant:

Each line represents a different generation. For example, the red line shows that people born in the 1920s were about 96% likely to report a religious preference when they were interviewed in their 40s, 50s, and 60s, and possibly less likely to be religious when they were in their 80s.

The conclusion I drew from this figure is that the differences between generations are larger than the changes, over time, within each generation. For purposes of modeling I concluded that religious disaffiliation accounts for only a small part of the observed changes in religious identity.

But I was bothered by one feature of these curves: many of them are concave down, and the maximum point in the curves is apparently shifting toward younger ages. I came to suspect that this picture of the data is "out of focus".

We can refocus the image by plotting the date of the survey (rather than the respondent's age) on the x-axis. Here's what that looks like:

In this figure, two trends are more apparent: before 1990, most generations were becoming more religious; after 1990, they all became less religious. So it seems clear that the explanation is something that affected all generations at a particular interval in time, not something that affects all people as they age.

We can see these changes more clearly by normalizing each curve with its 1990 value:

Again, most generation were becoming more religious before 1990; after 1990, all of them became less religious. Among people born in the 1960s, more than 10% lost their religion between 1990 and 2010 (when they were in their 30s and 40s).

Here's the same graph for people raised Catholic:

The general shape is the same: religious affiliation was flat or increasing prior to 1990, and decreasing for almost all generations after 1990.

Since the trends are similar for Catholics and Protestants, we can get a less noisy picture by combining them. Here is the same graph for respondents raised with any religion.

This figures makes it easier to compare across generations. It appears that more recent generations (born in the 1960s and 1970s) are disaffiliating at higher rates than earlier generations.

[As an aside, this result contradicts one of the primary (and widely-reported) claims of this article: Schwadel, Period and Cohort Effects on Religious Nonaffiliation and Religious Disaffiliation. Schwadel reports that people born in the 1960s and 1970s were disaffiliating at a slower rate than the previous generations. Some reasons my results might be different: Schwadel only had GSS data up to 2006, and he discards people under 30 years of age. So very little data about the youngest generations is included. Also, his results are based on statistical models that (if I understand correctly) don't include time as an explanatory variable, so they cannot account for an event that affects all generations during a particular interval.]

All right, it's audience participation time. What happened in the 1990s that caused widespread religious disaffiliation? Remember, idle speculations only. No evidence, please!

Secularization in America: part four

2012-06-28T09:02:00.001-07:00

Summary so far

Religiosity curves

Respondents in the GSS are surveyed at different ages, so we can get a sense of when people lose their religion (or acquire one). I collected all GSS respondents and partitioned them by the religion they were raised in and the decade they were born. For each of these subgroups, I plotted religiosity (the fraction with some religious preference) as a function of age when surveyed.

Here are the curves for people raised Protestant:

In the top right, we see that people born between 1900 and 1910 and raised Protestant were likely to be religious when they were interviewed in their 70s and 80s. In the lower left, we see that people born in the 1980s were less likely to be religious when they were interviewed in their 20s.

For the middle generations, we have a better sense of changes in religiosity over a respondent's lifetime. Several of the curves have an apparent peak in middle age; if this apparent effect is real, the location of the peak may be shifting left.

Overall, these curves are relatively flat, which suggests that respondents are not changing substantially after adulthood (everyone in the GSS is 18 or older).

The curves for Catholics are similar:

Again, there is a substantial differences between generations, but within each generation, little change over the respondents' lifetimes. People born in the 50s, 60s and 70s might be leaving the church as they age, but it is hard to tell in this plot whether these trends are statistically significant.

Finally, here are the curves for people raised with no religion:

There are only enough respondents in this category to plot curves for a few generations, and even then, the curves are noisy. Not surprisingly, people raised without religion are less likely to be religious, and recent generation are less religious than their elders. Again, the curves are generally flat, suggesting that people generally do not change religious affiliation as adults.

A possible exception is that people born in the 1970s and raised without religion might be finding religion in their 30s. But this data point is based on a small number of respondents, so it is probably too early to tell.

Why people switch

In 1988 the GSS asked respondents questions about changes in religious affiliation and the reasons for the change. Unfortunately, it looks like this data won't do me much good, because:

In many cases where a respondent switched from a religious preference to None, they were not asked why.
There are so many inconsistencies in the data, I wonder if it might have been mangled.
Because these questions were only asked once, we can't track trends.

So that's disappointing.

Modeling a mixed-age cohort

One of the challenges of working with GSS data is that the respondents each year are a mixture of people of all ages. From year to year, the oldest generation drops out of the cohort and the youngest generation joins the mix.

So when there is a trend from each generation to the next, as with religious behavior, there is a lag before the trend appears in a GSS time series, and the slope of the trend is much slower.

However, for purposes of prediction, this lag is actually useful. For example, 18 years after a baby boom, there is likely to be a spike in college enrollment; that's not really a prediction about the future; it's just a consequence of something that has already happened.

Similarly, we already know what most of the GSS cohort will look like next year. It will look like the cohort this year, one year older. The difference is that a few of the oldest respondents are replaced by the next group of 18 year olds.

In the 2010 cohort, the age range is roughly 20-80. To predict the 2020 cohort, we can:

Remove respondents older than 80.
Age the rest of the respondents by 10 years.
Add a new batch of respondents in their 20s.

Step 2 might be hard if people were changing religious affiliation as they age, but as we saw above, they generally do not. Step 3 is harder, but there are two reasonable options:

Conservatively, we can assume that the next generation will be like their immediate predecessors.
Alternatively, we can extrapolate from current trends. This option is probably better for prediction, but in some ways unsatisfying because it does not explain the cause of the trends, or why we should expect them to continue.

If we use this method to predict 20 years into the future, we replace about 25% of the cohort with simulated respondents. But since 75% of the prediction is based on simple population aging, it is likely to reliable.

The falling slinky problem

2012-06-26T07:31:00.000-07:00

Let's take a break from statistics and do some physics!

My friend Ted Bunn recently wrote about the falling slinky problem in his blog. He points to this video, which shows a falling slinky in slow motion. After the top of the slinky is released, the bottom seems to hover until the top reaches it. The effect is particularly strange because if you look carefully, the top of the slinky does not accelerate as we expect for an object in free fall. Rather, it falls at a constant rate.

Ted explains:

...the information that the top end has been dropped can’t propagate down the slinky any faster than the speed of sound in the slinky (i.e., the speed at which waves propagate down it), so there’s a delay before the bottom end “knows” it’s been dropped. But it’s surprising (at least to me) to see how long the delay is.

This explains why there is a delay, but to me it doesn't explain why the delay is the same as the time it takes for the top of the slinky to reach the bottom. There are lots of models out there that explain parts of this behavior, but the ones I found are either complicated or wrong.

Here's my take on it. First, let's assume that what we see in the video is correct: the slinky collapses from top to bottom, so that each coil doesn't move until the one above it comes down and (nearly) hits it.

Let's call the initial length L and the mass m. After some time, a fraction of the slinky, x, has collapsed. At that point, the collapsed part of the slinky has mass xm at height (1-x)L. The rest of the slinky is spread uniformly [EDIT: this assumption is not right...see Ted's comment below] between height 0 and (1-x)L. So the center of mass is

x(1-x)L + (1-x)(1-x)L/2

Since the slinky is in free fall, we know the center of mass as a function of time:

L/2 - g/2 t^2

If we set those equal and type them into WolframAlpha, we get

x = sqrt(g/L) t

Which means that the top of the slinky is moving at constant speed. Remember that x is the fraction of the slinky that collapsed; to get the distance traveled, we multiply by L:

d = xL = sqrt(gL) t

So the speed of the top of the slinky is sqrt(gL).

We can get to the same result a different way by using the formula for wave speed in a vibrating string: sqrt(T/μ), where T is tension and μ is mass per linear measure. In this case T=mg and μ=m/L. Plug that in and get wave speed sqrt(gL).

I think this analysis is useful, but to be rigorous, I haven't really explained why the slinky behaves the way it does. I have only shown that if the slinky collapses from top to bottom (as it appears to), then the top moves at a constant speed (as it appears to).

[UPDATE: Provoked by my amateurish attempts at Physics, Ted Bunn wrote up a version of this model that deals correctly with the change in the density of the spring from top to bottom. The result is that the speed of the top of the slinky is almost constant -- it slows down a bit at the end. ]

Secularization in America: part three

2012-06-22T12:47:00.001-07:00

The spouse tables are based on the parents of 1988 respondents. People from later generations might be increasingly likely to marry outside their religion.
The environment table is also based on the previous generation; again, later parents might be making different decisions about the religious environment of their children.
The transmission table is based on 1988 respondents; it's possible that after 1988, children were less likely to adopt the religion they were raised in. Anecdotally, the culprits most often blamed for this effect are college and the Internet.
Finally, I have not considered adult conversions from one religious identity to another. The GSS has data on these switches, so I could add them to the model.

I will investigate each possibility in turn, starting with the prevalence of mixed-religion marriages. In Secularization, Steve Bruce presents results from a study of intermarriage in the UK that found that the rate of vertical transmission:

"is halved if the parents are of different faiths (even when the differences are just Methodist-Anglican). Even if the parents agree on which faith they wish to pass on, the fact of disagreement makes the child aware that there are good people in other churches and introduces the relativism that weakens conviction. [page 71]"

So if the rate of mixed marriages is increasing, that could contribute to the increasing number of Nones.

To measure this effect, I used these GSS variables:

RELIG: What is your religous preference?
SPREL: What is your husband's/wife's religious preference?

In cases where one partner converts to the other's religion before marriage, that would count (for this model) as a same-religion marriage, since we are interested in the decision the couple makes about the religious environment they raise children in.

The Spouse Tables

The following graph shows the fraction of same-religion marriages over the history of the survey (data for SPREL were not collected every year):

Before 1988, the fraction of same-religion marriages was around 84%; after 1988 it fell to 78%. The abruptness of the change makes me worry that it may be an artifact; for example, a chance in the wording of the question. Also, these results only include respondents who are married, so they are biased toward older people and socio-economic groups that are more likely to be married.

But as it turns out, even if we take the data at face value, it has a small effect on the model's predictions.

I used the respondents from 2004-2010 to build spouse tables for men and women (see Part Two), then ran the 1988 model again with the anachronistic data. The results are almost identical to what we saw last time:

The only noticeable effect is that the prediction for Other got worse. I conclude:

It's possible that people are more likely now to marry outside their religion than in 1988, but the difference is small, and
Even if we cheat by using the 2004-2010 data in 1988, this change does not explain the subsequent changes in the fractions of Protestants and Nones.

The Environment Table

It seems unlikely that parents now are making different decisions about what religious environment to raise their children in, but just to rule it out, I compared the environment tables for 1988 and 2008.

prot cath jew other none change N excess

prot-prot 97 1 0 0 2 +1 634 8.7
prot-cath 43 46 0 1 10 +4 49 1.9
prot- jew 0 0 0 0 0 +0 1 0.0
prot-othe 44 21 0 36 0 +0 5 0.0
prot-none 85 4 0 0 11 +5 90 4.2
cath-prot 30 56 0 0 14 +14 37 5.1
cath-cath 1 98 0 0 1 +0 294 0.8
cath- jew 0 0 0 0 0 +0 0 0.0
cath-othe 0 0 0 0 0 +0 2 0.0
cath-none 2 82 0 0 16 +7 20 1.5
jew-prot 0 0 0 0 0 +0 1 0.0
jew-cath 0 100 0 0 0 +0 0 0.0
jew- jew 0 0 96 0 4 +0 27 0.0
jew-othe 0 0 0 0 0 +0 0 0.0
jew-none 0 0 67 0 33 +33 0 0.0
othe-prot 27 0 0 54 18 +18 4 0.7
othe-cath 0 43 0 57 0 +0 3 0.0
othe- jew 0 0 0 0 0 +0 0 0.0
othe-othe 7 0 0 84 9 +3 28 0.7
othe-none 25 25 0 25 25 +25 4 1.0
none-prot 83 0 0 0 17 +17 6 1.0
none-cath 14 68 0 0 18 -82 2 -1.6
none- jew 0 0 100 0 0 +0 0 0.0
none-othe 0 0 0 0 0 -100 1 -1.0
none-none 23 4 0 0 74 +18 34 6.1

The left column is the mother's-father's religion. The next five columns show the religious environments those parents chose, as reported by their children in 2008. For example, the second row shows that if the mother is Protestant and the father Catholic, 43% of the children were raised Protestant, 46% Catholic, and 10% None.

The next column shows the change in the None column, in percentage points, since the 1988 survey. N is the number of families in 1988 that fell into each category. Finally, Nones is the product of change and N, an estimate of the number of additional Nones in the 1988 survey that could be explained by changes in the environment table. The total of this column is 29, which is not nearly enough to explain the actual excess of 177.

Of course, most of the numbers in the change column are based on small samples, so we should not take them too seriously. By running simulations with resampled survey data, we can take account of these sample sizes.

Using the tables from 1988 to predict the fractions of Nones in 2008, we expect only 8.0% (compared to the actual 16.8%). If we used the environment table from 2008, the prediction goes to 8.5%. If we also use the spouse table, it goes up to 8.7%. So clearly the changes in these tables were not enough to explain the observed changes.

The transmission table

The transmission table is a cross-tabulation of the religion the respondent was brought up in and the religion reported when surveyed. It shows the outcome, after some years, of parents' decisions about their children's religious upbringing and the effect of the environment.

The following is the transmission table for 2008, with changes since 1988:

prot cath jew other none change N excess
prot 82 3 1 2 13 +7 951 67.7
cath 17 70 0 2 12 +6 414 24.6
jew 9 2 73 1 14 +9 31 2.7
other 12 2 0 75 10 +1 31 0.4
none 31 5 1 2 62 +5 53 2.6

Each row corresponds to a religious upbringing; each column shows a possible outcome. For example, the first row shows that of children raised Protestant, 82% report that their religious preference is Protestant, and 13% report None. The fraction of Nones has increased 7 percentage points since 1988. Since there are 951 people in this row, this increase accounts of 68 excess Nones in the 2008 survey.

Overall, the changes in the transmission table account for 98 excess Nones, which is a little more than half of the observed increase.

If we run the simulations again, applying the transmission table from 2008 in 1988, we get the following predictions:

The prediction for Nones is better, but it's clear that this model still misses the mark: it predicts that the fraction of Catholics should be going down, and fails to predict the decrease in the fraction of Protestants.

The problem is that I am treating everyone interviewed in 1988 as a cohort, but they represent people of all ages, who were raised in different environments. Also, I am using data from 2008 to predict what will happen in 2008, so I have got away from the original goal, to see whether the changes that occurred between 1988 and 2008 could have been predicted in 1988.

However, this model has given me some leads. It looks like a large part of the increase in Nones is due to changes in the transmission table, possibly a small part due to the environment table, and little or none due to the spouse tables.

Next time I will present a different model that reorganizes respondents into cohorts by age of birth, which will make it possible to compare people raised over the same time span. It will also allow me to look for trends that began prior to 1988.

Secularization in America: part two

2012-06-21T12:14:00.001-07:00

In Part One I described some trends in market share of the major religions in the U.S.; in particular, since 1988, the fraction of Protestants dropped from 60% to 51%, and the fraction of people with no religious affiliation increased from 8% to 18%.

I would like to know if something happened after 1988 to cause these changes, or if they could have been predicted based on patterns occurring before 1988. As a first step, I will use data from 1988 to model vertical transmission (from parent to child) and see if it predicts the observed changes

My model of vertical transmission works like this:

Each respondent chooses a spouse,
Each pair decides what religion to bring their children up in,
Each child chooses a religion.

I model each step of this process using data from the General Social Survey (GSS); specifically, I used these variables.

RELIG: What is your religous preference?
RELIG16: In what religion were you raised?
MARELIG: What was your mother's religious preference when you were growing up?
PARELIG: What was your fathers's religious preference when you were growing up?

The first two questions were asked every year, but questions about parents' religion were only asked in 1988 and 2008. I will use the data from 1988 to build and validate models, then use the data from 2008 to make predictions.

I used MARELIG and PARELIG to build two "Spouse tables", one for men and one for women. Here is the table for men:

Spouse Table (men)

	prot	cath	jew	other	none
prot	93	6	0	0	1
cath	14	85	0	0	1
jew	4	0	96	0	0
other	6	4	0	90	0
none	59	13	0	3	24

Each row indicates the religion of a male respondent; each column is the religion of a possible spouse; the numbers are percents. For example, the first row indicates that 93% of male Protestants married other Protestants, and another 6% married Catholics.

Here is the spouse table for women:

Spouse Table (women)

	prot	cath	jew	other	none
prot	82	7	0	0	12
cath	10	85	0	0	5
jew	4	0	96	0	0
other	8	3	0	74	15
none	13	7	0	0	80

In general, women are more likely to marry out of their religion than men, but still the great majority marry a co-religionist. One asymmetry is apparent: men with no religion seldom marry another None (24%), but women with no religion usually do (80%). This effect is partly due to the gender gap: 11% of male respondents are Nones, but only 5% of the women are (there is a similar, possibly smaller, gender gap in the CIRP data).

Once the respondents have paired up, they decide what religion to raise the children in. The following table shows results from the 1988 data. The rows enumerate all pairs of mother's and father's religion; the columns indicate the religious environment they chose. For example, the second row indicates that if a Protestant woman marries a Catholic man, they raise the children Protestant 58% of the time, Catholic 36% of the time, and None 6%.

Environment table

parents	prot	cath	jew	other	none
prot-prot	99	1	0	0	1
prot-cath	58	36	0	0	6
prot-jew	100	0	0	0	0
prot-other	100	0	0	0	0
prot-none	89	4	0	1	7
cath-prot	39	61	0	0	0
cath-cath	1	99	0	0	0
cath-jew	0	0	0	0	0
cath-other	100	0	0	0	0
cath-none	17	69	0	6	8
jew-prot	0	0	100	0	0
jew-cath	0	0	0	0	0
jew-jew	0	0	96	0	4
jew-other	0	0	0	0	0
jew-none	0	0	0	0	0
other-prot	60	0	0	40	0
other-cath	0	0	0	100	0
other-jew	0	0	0	0	0
other-other	7	2	0	89	2
other-none	33	0	0	67	0
none-prot	100	0	0	0	0
none-cath	0	0	0	0	100
none-jew	0	0	0	0	0
none-other	0	0	0	0	0
none-none	40	4	0	0	56

One surprise in this table is the last row: when two people with no religion marry, 40% of the time they apparently choose to raise their children Protestant. This seems unlikely, but there are several possible explanations: (1) the parents might have chosen to raise their children in the prevalent religion of their community, (2) a respondent might not have been raised by his parents, (3) a respondent might not be reporting his parents' religion accurately. For purposes of modeling I take these responses at face value.

Children raised with a religion usually adopt that religion, but not always. The following "transition table" shows possible outcomes for each religious environment. For example, 89% of respondents who say they were raised Protestant also report that their religious preference is Protestant, but 3% are Catholic and 6% have no religious preference. More people convert from Catholic to Protestant than the other way around.

Transition table

	prot	cath	jew	other	none
prot	89	3	0	1	6
cath	11	83	0	0	6
jew	0	0	95	0	5
other	5	3	0	83	9
none	32	11	0	0	57

As expected, the majority of people raised with no religion report no religious preference, but 32% of them identify as Protestant and 11% identify as Catholic. I found that surprising. I will look more closely later, but for now, again, I will take it at face value.

Finally, we can combine these results into a single "Generation table" that shows the transitions from one generation to the next. I ran simulations with following steps.

For each respondent, choose a spouse's religion from the Spouse Table.
For each parent pair, choose a religious environment from the Environment Table.
For each hypothetical child, choose a religious identity from the Transition Table.
For each parent-child pair, make an entry in the Generation Table, below.

Since this computation is based on random simulations, it varies from run to run, but here is a typical outcome:

Generation table

	prot	cath	jew	other	none
prot	86	6	0	1	7
cath	19	72	1	1	7
jew	0	0	95	0	5
other	29	10	0	55	6
none	67	9	0	1	23

Assuming that a generation time is about 22 years, we can use this model to predict the distribution of religions in 2010 (using only data from 1988). This figure shows the actual time series and the model predictions for each group:

On the right side of the plot, the vertical lines show the 90% confidence interval; the boxes show the mean of 20 simulation runs. [One technical note: each simulation is based on tables from resampled survey data, so the confidence intervals reflect both the sampling error of the survey and random variation of the simulations.]

The actual values for Catholics, Jews and Other fall within the prediction intervals, but the model fails to predict the decrease in Protestants or the increase in None.

So, what's missing from this model that could account for the observed changes?

The spouse tables are based on the parents of 1988 respondents. People from later generations are increasingly likely to marry outside their religion.
The environment table is also based on the previous generation; again, later parents might be making different decisions about the religious environment of their children.
The transition table is based on 1988 respondents; it's possible that after 1988, children were less likely to adopt the religion they were raised in. Anecdotally, the culprits most often blamed for this effect are college and the Internet.
Finally, I have not considered adult conversions from one religious identity to another. The GSS has data on these switches, so I could add them to the model.

Over the next few installments, I will investigate each of these factors to see which, if any, account for the observed changes.

Secularization in America, part one.

2012-06-19T07:43:00.000-07:00

In the last year or so I have written several articles about trends in religion among college students:

All of these are based on data from the Cooperative Institutional Research Program (CIRP) which runs the Freshman Survey, an annual survey of more than 200,000 incoming students at 270 colleges and universities in the U.S.

More recently, I read Secularization: In Defence of an Unfashionable Theory, by Steve Bruce. Bruce presents the "unfashionable theory" that as societies modernize, they secularize. In his formulation, modernization includes trends toward individualism, industrial capitalism, science and technology; and secularization means "decline in the social significance of religion."

The poster child for secularization is Western Europe, where the social influence of religion has been in decline for centuries, and where in every country the fraction of people with no religious affiliation has been increasing for decades.

But skeptics have suggested that countries where people are still religious, like the United States and many countries in the Middle East, are exceptions that disprove the theory. Bruce replies that religious countries in the Middle East are not exceptions because they are not modern, and the United States is not an exception because it is, in fact, secularizing.

The data from the Freshman Survey are consistent with secularization. The number of incoming college students with no religious affiliation has been climbing consistently since 1978, and the number of students reporting participation in religious service has fallen at about the same rate.

Of course, college students are not a random sample of the population; for that, we can use data from the General Social Survey (GSS), which is (according to the GSS) "widely regarded as the single best source of data on societal trends." It has run since 1972; each year (or every other year since 1994) it surveys a sample of about 2000 adults randomly sampled from the U.S. population. Respondents answer hundreds of questions about their background, life history, and beliefs. Many questions are repeated from year to year for trend analysis.

I will use this dataset to answer several questions:

Is there evidence of secularization in the U.S. (Hint: yes).
Can we explain the causes?
Can we predict how these trends will continue over the next few decades.

To get started, I tracked responses to the question, "What is your current religious preference?" The original set of options was Protestant, Catholic, Jewish, some other religion, or no religion. After 1994, the set of options was expanded, but for my purposes the original options are enough to describe large-scale trends. The following graph shows the fraction of the population in each group over time.

A few trends are apparent: the percentage of Protestants is declining; the percentages of Other and None are increasing. These trends are clearer in the following figures, broken into two intervals:

From 1972 to 1988, the fraction of Protestants and Catholics was unchanged, but the fraction of Nones may have increased.

From 1988 to 2010 (the most recent survey year), the fraction of Protestants and Jews declined, and the fraction of Nones increased by almost 250%. The number of Others increased during both intervals, with more variability.

This dataset shows signs of secularization in the U.S., at least since 1972. But religious affiliation is just one aspect of religious identity; there is a lot more data in the GSS to look at.

My particular interest is in explaining the trends we have seen so far, and predicting what's coming next. It is tempting to think that something happened in 1988 to cause the inflections in these curves, but I think it is more likely that the origin of these changes goes back farther.

To test that idea, let's pretend that it's 1988. We have see some changes in the market share of different religions since 1972, but nothing bigger than a few percentage points, and no indication of acceleration. Could we have predicted the much larger changes coming between 1988 and 2010?

In the next few articles, I develop several models intended to answer that question. Then I turn to prediction: using the data up to 2010 (and 2012 when it is available) what can we expect in the next 20 years?

Are religious colleges getting more religious?

2012-05-07T07:47:00.001-07:00

In response to my article about the increasing numbers of students entering college with no religious affiliation, a reader wrote:

A potentially interesting trend to look at is the religious participation of students attending Catholic or other religious institutions. I wonder if the trend is toward religious students (those who report both religious affiliation and participation in students) being more likely to go to religiously-affiliated colleges. This would mean a reduction of religious students at state/non-sectarian schools and an increase in the religiosity of students at affiliated schools (this might even be skewed because the survey doesn't include some of the most religious schools such as Liberty University). This would be a reflection of the increasing polarization of our society.

On another subject, I tend to distrust steadily increasing social trends. From a complexity theory perspective, I would expect more of a cycle in religious (dis-)belief, so I wouldn't be surprised to see a crash in the number of non-believers sometime in the near future, although predicting exactly when that will happen is virtually impossible with the relatively small amount of data available (40 years is not enough).

There are several interesting questions here. The first is whether religious colleges are getting more religious. This one is relatively easy to investigate: the HERI survey provides data broken down by several types of colleges, including Nonsectarian, Catholic, and Other Religious Affiliation.

The good news is that their reports are available in PDF now. The bad news is that most of the older ones are scanned and not OCRed, so they are not easy to search. I went through them by hand, but I only extracted the data at 5-year intervals. Here is what the trends look like for students responding "None" for religious preference at Private nonsectarian 4-year colleges, Catholic colleges, and Other religious colleges:

The percentage of Nones is increasing in all categories, more slowly at religious colleges than at other private colleges. Nevertheless, the fraction of students at religious colleges with no religious preference has nearly tripled in the last 35 years.

As an aside, I also plotted data for historically black colleges and universities (HBCU). Here's what that looks like:

Clearly the trend is slower; in fact, it is not obvious that it is statistically significant. And since about 1990, the percentage of Nones is higher at religious colleges than at HBCUs, by more than a factor of two (13%, compared with 7%).

Getting back to the reader's question, it doesn't look like the religious schools are getting more religious. In fact, they are getting less religious at almost the same rate as other schools. But maybe the fraction of students going to religious colleges is increasing?

Here is a table from the National Center For Educational Statistics, which publishes the Digest of Education Statistics. It shows "Fall enrollment and number of degree-granting institutions, by control and affiliation of institution: Selected years, 1980 through 2009."

Here's what that data looks like:

Enrollments have been increasing for all college types, with religious colleges growing faster than private nonsectarian colleges until 1995. Here's what these data look like expressed as a percent of the total:

Before 1995, religious colleges were gaining market share, at the expense of nonsectarian colleges; other than that, there is not much going on.

So, to address the reader's questions:

Are religious students more likely to attend religious schools? Maybe. It's hard to tell with the data I have.
Are more students going to religious schools? Not lately. Since 1995, the fraction of students at religious colleges has been flat.
Are students at religious schools increasingly religious? No. The percentage of Nones has increased in all college types, including religious colleges.
Do these trends introduce a bias in the results I presented? Not that I can see.

As for distrusting steadily increasing social trends, I agree that some caution is needed. The percentage of Nones can't keep accelerating forever. But it can keep growing forever.

Steven Pinker presents several one-way trends in his new book, The Better Angels of Our Nature. And Peter Singer, who reviewed Pinker's book in the New York Times, discussed related ideas in The Expanding Circle.

But who knows? I guess we'll see what next year's data point looks like.

Many thanks to the reader who posted the comments that prompted this update.

Here is my original article from March 2011. Here is the update from January 2012.

Fog warning system: part three

2012-04-25T11:08:00.001-07:00

Background: I am trying to evaluate the effect on traffic safety of a fog warning system deployed in California in November 1996. The system was installed by CalTrans on a section of I-5 and SR-120 near Stockton where the accident rate is generally high, particularly during the morning commute when ground fog is common. The warning system consists of (1) weather monitoring stations that detect fog and (2) changeable message signs that warn drivers to reduce speed.

I will post my findings as I go in order to solicit comments from professionals and demonstrate methods for students. If I can get permission, I will also post my data and code so you can follow along at home.

Previously: In the first installment I reviewed the first batch of data I am working with, and ran some tests to confirm that Poisson regression is appropriate for modeling the number of accidents in a given day. In part two I ran Poisson regressions to identify factors that influence the number of accidents per day.

Critical events

I have been waiting to get more details about several events that affected traffic safety during the observation period. I was able to get in touch with a Transportation Engineer in the Traffic Safety Branch of Caltrans District 10, which includes the study area. According to Caltrans records, the speed limit on the relevant section of I-5 was increased from 55 to 70 mph on March 25, 1996. The speed limit on SR-120 was increased from 55 to 65 mph about a month later, on April 22, 1996. Many thanks to my correspondent for this information!

The automated warning system was activated in November 1996. My collaborator has collected data on weather measurements made by the system and the warning it displayed. I hope to get this data processed soon.

Accidents per million vehicles

In the previous article, I ran models with raw accident counts as the dependent variable, and found that traffic volume is a significant explanatory variable. Not surprisingly, more cars yield more accidents.

Rather than use volume as an explanatory variable, an alternative is to express the dependent variable in terms of accidents per million vehicles. As a reminder, here's what the traffic volume (in thousands of cars per day) looks like during the observation period:

And here are the raw accident counts:

I divided counts by volume and converted to accidents per million cars. At the same time I smoothed the curves by aggregating quarterly. Here's what that looks like:

The vertical red lines show major events expected to affect traffic safety: increased speed limits in March and April 1996, and the activation of the warning system in November 1996.

This graph suggests several observations:

In the control directions, the accident rate was flat from 1992 through 1994, increased quickly in 1995 (before the speed limits were increased) and has been flat every since.
In the treatment directions, the accident rate was trending down until late 1996, including three quarters after the speed limit was increased. The accident rate increased sharply in 1997 and possibly again in 2000.
The accident rate in both directions was unusually low during the third quarter of 1996, when the warning system was activated. Other than that, there is no obvious relationship between accident rates and the events of 1996.

Since we don't expect the warning system to have much effect on the control directions (that's why they're called "control"), the speed limit changes are by far the most likely explanation for the accident rate changes. But it is puzzling that a large part of the change occurred before the new speed limits went into effect. One possibility is that as new speed limits were rolled out throughout California, drivers became accustomed to higher speeds and drove faster even on roads where the new limits were not in effect. But if that's true, it doesn't explain the continuing decline in the treatment directions.

My collaborator has some data on actual driving speeds before and after 1996. Once I process that data, I will be able to get back to this puzzle.

Injuries and fatal accidents

In response to a previous post, a reader suggested that if the warning system causes drivers to slow down, it might affect the severity of accidents more than the raw number. To investigate that possibility, I also plotted the rates for injury accidents (including fatalities) and fatal accidents.

Here is the graph for injury accidents:

The patterns we saw in the previous graph appear here, too. In addition, this graph suggests, more strongly, the possibility of a second changepoint in late 1999 or 2000.

And here is the graph for fatal accidents:

The number of fatal accidents is, fortunately, small. During more than 10 years of observation, there were only 26 in the study area. The trends in the other graphs are not apparent here, other than the general increase in the rate of fatal accidents in the second half of the observation period.

Summary

Accident rates in the control and treatment directions increased sharply around 1996, but neither effect is related in an obvious way to increased speed limits or deployment of the warning system.
Accident rates were unusually low in the quarter the warning system was activated; other than that, no effect of the warning system is apparent.
It looks like there was a second increase in accident rates in late 1999 or 2000. I will ask my correspondent at Caltrans if he has an explanation.

Next steps

There's not much more I want to do with this data. Now I need more numbers! In particular, I will be able to get data from the warning system itself, including:

Conditions measured at roadside weather stations, which should be better than the data I have from the airport 8 miles away, and
Messages displayed when the warning system was active.

If the warning system has an effect, it should be apparent on the days it is active. By comparing the treatment and control directions, it should be possible to quantify the effect.

Also, I have permission now to share the data; I will try to get it posted, along with my code, before the next update.

[UPDATE April 26, 2012]

A reader asked

I can think of two ways that overall traffic volume affects accident rates: (1) more cars = more accidents overall, which you control for by measuring accident rates, and now you're seeing rising accident rates per car. So this raises the next thought, (2) more cars = more traffic density, which raises accident rates per car for each car on the road.

What happens if you regress on traffic volume squared, or include traffic volume as an independent variable in the accident rate regression? The density effect is likely nonlinear but it's a thought.

This is a great question. If there is a non-linear relationship between traffic volume and the raw number of accidents, then even after we switch to accident rates, there might still be a positive relationship between traffic volume and accident rates.

I ran these regressions, and in fact there is a relationship, but with the limitations of the data I have, I don't think it means much. Specifically, I only have annual estimates for traffic volume, so there's no fluctuation over time; traffic volume increases at a nearly constant rate for the entire observation period (see the figure above).

So traffic volume will have a positive relationship with anything else that's increasing, and a negative relationship with anything decreasing. And that's what I see in the regressions:

All of the relationships are statistically significant, but notice that in the treatment directions, before 1996 when the accident rate was declining, the relationship with traffic volume is negative!

I don't think this variable has any explanatory content; any other ramp function would behave the same way. If I can get finer-grain data on traffic volume, I might be able to look for a more meaningful effect.