Carl Bialik and Andrew Flowers at fivethirtyeight.com (Nate Silver’s site) ran a story following up on our birthdays example—that time series decomposition of births by day, which is on the cover of the third edition of Bayesian Data Analysis using data from 1968-1988, and which then Aki redid using a new dataset from 2000-2014.
Friday the 13th: The Final Chapter
Bialik and Flowers start with Friday the 13th. For each “13th” (Monday the 13th, Tuesday the 13th, etc.), they took the average number of babies born on this day over the years in the dataset, minus the average number of babies born one week before and one week after (that is, comparing to the 6th and 20th of the same month). Here’s what they found:
The drop on the 13th is biggest on Fridays, of course. It’s smaller on the weekends, which makes sense. Weekend births are already suppressed by selection so there is not much room for there to be more selection because of being on the 13th.
All 366 days
Here’s the fivethirtyeight’s art department’s redrawing of a graph Aki made for the 2000-2014 data:
You can see dips on the 13th of each month, also a slight bump on Valentine’s, and drops on several other special days, including February 29th, Halloween, and September 11th, along with a bunch of major holidays.
Here’s the full decomposition:
If you look carefully on the bottom graph you’ll again see the little dips on the 13th of each month, but the blow-up above is better. The reason I’m showing all four graphs here is so you can see the decomposition with all the pieces of the model.
Comparing to the results from the earlier time period (just look at your copy of BDA3!), we see the new pattern on September 11th, also in general the day-of-week effects and special-day effects are much higher. This makes sense given that the rate of scheduled births has been going up. Indeed if the birth rate on Christmas is only 50% on other days, that suggests that something like 50% of all births in the U.S. are scheduled. Wow. Times have changed.
Checking and improving the model
There are some funny things about the fitted model. Look either at the graph above with the day-of-year patterns (the very second plot of this post) or, equivalently, the bottom of the four plots just above.
The first thing I noticed, after seeing the big dips on New Year’s, July 4th, and Christmas, were the blurs around Labor Day, Memorial Day, and Thanksgiving. These are floating days in the U.S. calendar, changing from year to year. (For example, Thanksgiving is the fourth Thursday in November.) At first I thought this was a problem, that the model was trying to fit these floating holidays with individual date effects, but it turns out that, no, Aki did it right: he looked up the floating holidays and added a term for each of them, in addition to the 366 individual date effects.
The second thing is that the day-of-week effects clearly go up over time, and we know from our separate analyses of the 1968-2008 and 2000-2014 data that the individual-date and special-day effects increase over time too. So really we’d want these to increase for each year in the model. Things are starting to get hard to fit if we let these parameters vary from year to year, but a simple linear trend would seem reasonable, I’d think.
The other problem with the model we fit can be seen, for example, around Christmas and New Year’s. There’s a lot more minus than plus. Similarly on April Fools and Halloween: the number of extra births on the days before and after these unlucky days is not enough to make up for the drop on the days themselves.
But this can’t be right. Once a baby’s in there, it’s gotta come out some time! The problem is that our individual-date and special-day effects are modeled as delta functions, so they won’t cancel out unless they just happen to do so. I think it would make more sense for us to explicitly fit the model with little ringing impulse-response functions that sum to 0. For example, something like this: (-1, -2, 6, -2, -1), on the theory that most of the missing babies on day X will appear between days X-2 and X+2. It’s an assumption, and it’s not perfect, but I think it makes more sense than our current (0, 0, 1, 0, 0).
You can also see this problem arising in the last half of September, where there seem to be “too many babies” for two entire weeks! Obviously what’s happening here is that the individual-date effects, unconstrained to sum locally to zero, are sucking up some of the variation that should be going into the time-of-year term. And, indeed, if you look at the estimated time-of-year series—the third in the set of four graphs above—you’ll see that the curves are too smooth. The fitted Gaussian process for the time-of-year series got estimated with a timescale hyperparameter that was too large, thus the estimated curve was too smooth, and the individual-date effects took up the slack. Put zero-summing impulse-response functions in there, and I think the time-of-year curve will be forced to do better.
This is a great example of posterior model checking, and a reminder that fitting the model is never the end of the story. As I said a few years ago, just as a chicken is an egg’s way of making another egg, Bayesian inference is just a theory’s way of uncovering problems with can lead to a better theory.
Results from U.S. CDC data 1994-2003 look similar to the results from U.S. SSA data 2000-2014, but for some reason in the overlapping dates in years 2000-2003 there are on average 2% more births per day in the SSA data:
It’s easy enough to adjust for this in the model by just including an indicator for the dataset.
We (that is, Aki) fit a hierarchical Gaussian process to the data, with several additive terms (on the log scale, I assume; at least, it should all be on the log scale). Details are in the appropriate chapter of BDA3.
Code or it didn’t happen
I’d love to give you the Stan code for this model—but it’s too big to fit in Stan. It’s not that the number of data points is too big, it’s that the Gaussian process has this big-ass covariance matrix which itself depends on hyperparameters. So lots of computing. We have plans to add some features to Stan to speed up some matrix operations so this sort of model can be fit more efficiently. Also there should be some ways to strip down the model to reduce its dimensionality.
But for now, the model was fit in GPstuff. And here’s the GPstuff code for various versions of the model (maybe not quite the one fit above).
P.S. The movie Friday the 13th and its sequels were supposed to be really bad—I never actually saw any of them myself, despite being in the target audience—but I just read the wikipedia article for the series, and it’s fascinating. Just for example:
When Harry Manfredini began working on the musical score for the 1980 film, the decision was made to play the music only alongside the killer so as not to trick the audience into believing that the killer was around during moments that they were not supposed to be. Manfredini explains that the lack of music for certain scenes was deliberate: “There’s a scene where one of the girls […] is setting up the archery area […] One of the guys shoots an arrow into the target and just misses her. It’s a huge scare, but if you notice, there’s no music. That was a choice.” Manfredini also noted that when something was about to happen, the music would cut off so that the audience would relax a bit, which allowed the scare to become more effective.
Since Mrs. Voorhees, the killer in the original Friday the 13th, does not show up until the final reel of the film, Manfredini had the job of creating a score that would represent the killer in her absence. Manfredini was inspired by the 1975 film Jaws, where the shark is not seen for the majority of the film, but the motif created by John Williams cued the audience as to when the shark was present during scenes and unseen. While listening to a piece of Krzysztof Penderecki music, which contained a chorus with “striking pronunciations”, Manfredini was inspired to recreate a similar sound for Friday the 13th. He came up with the sound “ki ki ki, ma ma ma”, based on the line “Kill her mommy!”, which Mrs. Voorhees recites repeatedly in the final reel. The “ki” comes from “kill”, and the “ma” from “mommy”. To achieve the unique sound he wanted for the film, Manfredini spoke the two words “harshly, distinctly, and rhythmically into a microphone” and ran them into an echo reverberation machine. Manfredini finished the original score after a few weeks and recorded it in a friend’s basement. Victor Miller and assistant editor Jay Keuper have commented on how memorable the music is, with Keuper describing it as “iconographic”. Manfredini makes note of the mispronunciation of the sounds: “Everybody thinks it’s cha, cha, cha. I’m like, ‘Cha, cha, cha’? What are you talking about?”