Monday, February 29, 2016

Epidemiology

So, what if you read yesterday's post and thought to yourself, "that's great, but I don't know about a field called Epidemiology"? Well, you could look it up on Wikipedia, of course, but, as long as you're going to trust some amateur hack to give you an encyclopedia entry, it might as well be me.

It started with shit. Seriously.

London in the 19th century was being overrun with the stuff. The flush toilet has been around in some form since ancient times. Sir Thomas Crapper contributed to its widespread development in the late 1800's, but did not actually invent the thing. What happened between ancient times and Crapper's fortune was the birth of Epidemiology.

In 1854, London was thrown into a panic as an epidemic of cholera hit the city. All sorts of things were blamed for it including the simple observation that the destitute were being hit the hardest, so it must just be that god doesn't like the lower class. Talk to Mr. Trump about Mexicans if you think we've progressed since then.

John Snow (no, not the dude from Game of Thrones) was a doctor in London who surmised that it might have something to do with the fact that people were literally drinking their own waste; the standard way to get rid of the stuff was to dump it out your window into the gutter where it could contaminate the water supply.

The powers that be were not impressed with his argument so he did what any self-respecting man of the Age of Reason would do: he gathered more data. He plotted where each victim lived on a map. He then looked for the nearest well. No luck; it appeared several wells must be contaminated and that would be quite a coincidence if the root cause was a breech allowing waste water into the fresh supply.

Then comes the brilliant part. Like most brilliancies, it's obvious in hindsight. People don't take helicopters to fetch water, they walk. So, he redrew his map, this time using minimum walking distance as his distance metric. He accounted for shortcuts like back doors and allyways. Blam. A single well was isolated. In the face of the newly presented evidence, officials issued an order to close the well in question. The epidemic ended as quickly as it had started and a new field of study was born. (And a new industry - London became the pioneering city for modern sanitation).

Like every field, Epidemiology has come a long way in the subsequent 150 years, but it's still basically about the same thing: extracting information embedded in disease and mortality data to suggest (not prove) a root cause. My master's work looked at methods to prioritize public health investigations based on the likelihood (or, unlikelihood) of the spatial distribution of disease cases. It was a minor result, for sure, but some very practical good came out of it. The New York State Department of Health used the results to streamline their investigation process and they did, in fact, turn up some nasty toxins.

Epidemiologists have been prying pearls of knowledge out of data for over a century, but the Data Science crowd seems to think they're inventing all this stuff. Read a typical Data Science paper and you'll see few, if any, references to journals like Biometrika (looks Russian, but it's actually an Oxford publication) or The American Journal of Epidemiology (the latter being where I was published). It's a gold mine. Just read all the significant results from the last 40 years in Epidemiology and you'll blow the doors off these newbies.

Please note that I'm not suggesting that the Data Science crowd is a bunch of idiots. They're just not using a lot of work that's already been done because it's been published in places they don't generally care about. In fairness, the vast majority of Epidemiology literature is very specific to particular applications so, unless you're in the trenches trying to sort out why diabetes rates are inverse correlated with income or some other such problem, you wouldn't bother reading it. However, if you want to get in on this Data Science thing (And, why not? It's a booming field), a review of the milestone articles is an easy way to get a leg up.

Sunday, February 28, 2016

What's old is new again

No great pearls of wisdom today. Spent most of the evening finishing off my homework.

I will note, however, that if one was to be seeking a PhD in Data Science, clustering would be the quickest route to that provided you knew about a field called Epidemiology. I'm reading all this stuff for my Data Mining class and the citations are all between 1990 and 2010 and I can't help thinking, "we knew all that in the 80's." It's almost enough to make me want to pivot my thesis topic but, for now, I'll carry on with the problem at hand.

Saturday, February 27, 2016

2005 Chicago Snowgaine

Plenty warm out today so this week's off-day throwback race report is to a run in more seasonable conditions. And, while she doesn't read this blog, I'll still give a "Happy Birthday" shout out to my wife Kate.

Run February 19, 2005

I'm not quite sure why it took David and I 8 years to try a team orienteering event together, but it obviously falls in the "really good idea" bucket. We've done navigation sections together as part of adventure races, but since you only get one set of maps, it wasn't so obvious how well we work together.

We arrived separately since I was staying in Chicago for the weekend. When I get there, David has already sketched out a preliminary route. It looks ugly - lots of doglegs - but I'm only able to find a few small improvements. The route ends with a road run of nearly two miles. To make that pay off, we'll have to hammer it pretty hard. One thing is for sure: the winning time will be well under the 6-hour limit.

Competition map: North, South.

Just before the start, meet director Kathy Bullard calls us all together for some instructions. As she counts down to the start I realize David is nowhere around. David always seems to be doing something right before the start of these things, even today when he was the first person to register. As she says "go" David pops out from behind his van and we head off.

On the way to our first control we are passed by Michael Eglin and Dimitry Perchonok. I expect that they will be our primary competition. Michael is as fast as we are and a better navigator on a true orienteering map. This map is a touched up version of USGS so our experience on rough maps evens the odds. I don't know whether Dimitry will help or slow him down.

Although the temperature is in the high 20's, the ground is frozen so solid it is like running on concrete. Make that concrete covered with ice, sticks, and leaves. Footing is bad in the valleys, worse on ridges, and nearly impossible on hillsides.

We get to the fourth control at the same time as Eglin/Perchonok and open up a slight lead heading up the hill. On the other side of the ridge, we drop into a steep reentrant. We can hardly go three steps without slipping and bruising something. Our adversaries take the power line cut down to the road and beat us to the next control by about 30 seconds. A few things are becoming clear:
  • We really should have worn spikes.
  • We are going to need to pick it up if we're going to get rid of Eglin/Perchonok (who are wearing spikes).
  • Spikes would be a good thing.
  • We need to give excessive preference to routes that follow trails and fields so we can run fast without spikes.
  • It sure would be nice to have some spikes.
Speaking of spikes (now in the navigation sense), another fact is becoming apparent: when David and I are navigating together on ridge and valley terrain (which is 90% of our terrain in St. Louis), we basically spike everything regardless of how fast we go. This is both good and bad. Good because it means we can go really fast. Bad because we now have no excuse not to go really fast.

After the fifth control (25), we head to 19 while the others go to 23. We crank up the pace on the route from 19 to 23 (mostly trail and field) and continue to push to 20. Just as we are getting there, we see Michael and Dimitry leaving the control. Encouraged by the fact that we are nearly a control ahead, we push even more.

Since the out and back strategy on trails and fields is working so well, we switch the order of 13 and 14, making 15-14-13 one big dogleg. 13-12-11 also becomes a dogleg as we stay in the valley rather than going over the ridge between 12 and 11.

We're running about as fast as I can go for this distance. Although we are covering a lot of extra ground to avoid the hillsides, we're still doing 8:00/Km against the straight line. We have one little boom at 8 as we both drift right of the control. We recover quickly and the error is less than two minutes.

Our most radical departure from the red line comes at control 4 where we decide to take the road all the way around to 3. This would have been great if the road hadn't been covered in ice. We have to take tiny steps as we half run/half slide down the hill.

When we get to 2 (our last control that requires any real navigation) we spend a few seconds taking in the view off the cliffs. Part of this is because it really is a nice view and part is because we know what comes next - the road run to the finish.

David scrambles down to the road ahead of me and I tell him to get moving since he has the punch card. I follow about 30 meters behind. If any further evidence was needed as to how evenly matched we are, we run the entire leg up to 10 spaced exactly like that. As David punches 10, I catch up and we run into the finish together.

Our strategy has paid off huge. We are the only team to use the long road run. Avoiding the hillsides near the end put us way out in front. Michael and Dimitry show up 43 minutes later and the third place team of Maricel Olaru and Natalia Babeti is nearly an hour back. Five teams and one solo competitor manage to sweep the course in under six hours.

Kathy is a little disappointed that we got done in just over 3 hours. I don't know that you could make the course much longer without adding a lot of controls. The area just isn't that big. I suggest that if she really wants to increase the distance, she could do a map exchange.

While I would have enjoyed a longer course, I also found it great fun to hammer full-on for three hours. We certainly wouldn't have been able to hold that pace for six. I felt like David could have gone a bit faster without me, although he might have had to slow down to navigate. At any rate, it was fun to put the Buckley-Frei Death Match aside for a day and work together. So much fun that I think we should do it more often than once every 8 years.

Friday, February 26, 2016

About that class thing...

Yes, I am taking classes this semester, not just writing algorithms of my own choosing. Spent a couple hours working the chapter 5 problems for Bayesian Stats this evening (homework is due Monday; this kind of crept up on me). A bit frustrating to spend two hours writing R programs to solve problems that I could do by hand in 15 minutes but I do get it. No point in programming stuff you don't understand until you're really sure you can validate the stuff you do understand.

Counting the cost

Missed yesterday's post due to yet another packed day. Straight to work upon waking so I could get out of there in time to meet with my adviser at 3PM. Then my Data Mining class at 4. Followed that with a short track workout (4x800m@3:00 plus a warm up and warm down). Then, straight to choir practice (I've pretty much quit choir, but a long-time associate pastor of our church passed away last week and I felt joining them for his funeral this Saturday was a fairly high calling). Finally, a short bit of research and off to bed. So, I'll enter a somewhat off-topic post this morning and may have some academic meat later on today.

It's only somewhat off topic because, as I've mentioned before, the big casualty in all this has been my running. That's a good thing, really. Something had to give and better my running performances than something that really mattered. That said, running matters to me.

I'll begin by acknowledging that there is absolutely no reason why somebody needs to be their absolute best at anything, much less a recreational pursuit. However, if one is trying to do something well, looking for shortcuts is not the way to go. My impending performance at The Woodlands serves as a illustrative case. Accepting the elite invite to The Woodlands Marathon last fall seemed like a good idea at the time. I was coming off a sub-3 and age group win at Milwaukee. With winter break providing six weeks of decent training time, it seemed like I could hold some decent form together.

We'll find out in 8 days just how far we need to stretch the definition of "decent form" to make that statement true but, for the purposes of argument, let's assume that all my pre-race predictors are right and I run somewhere around 3:03.

Taken in isolation, Runner's World would use that as a case study for their typical "Here's a super busy guy who does great on less training" story. Only 12 weeks of prep totaling barely 800 miles, but those two quality workouts per week had him beating his Boston Qualifying time by half an hour. They might even go further to suggest that I could qualify for Boston on even less if I just turned up the intense workouts a bit more.

Ah, but that's the rub. I probably could. But not because 40 miles a week and some track work is enough. It's because I have this huge base from when I had time to train right. My best marathons have been preceded by around 1700 miles in 18 weeks. The only reason I'll still be way under the BQ standard is that fitness doesn't leave you overnight (it just feels that way). I know, as does just about every credible coach, that running marathons well means running twice a day every day for years.

While it'a a surprising fact of biology that one can get reasonably close to optimal performance on much less, the running press has interpreted this as meaning that the extra effort doesn't really pay off. Nothing could be further from the truth. Five minutes may not sound like much, but the winning time for 50+ at nearly every mid-sized city marathon is between 2:55 and 3:00. That five minutes has demoted me from someone who had a legitimate shot of winning my age at 90% of the world's marathons to someone who doesn't have a chance. It's not coincidental that the winning times are so uniform. It's a diminishing returns game which compresses the results into a very small range. Those tight times are indicative that the top runners are engaged in some pretty serious marathon prep. You either do it or you lose to someone who did.

Again, I'm not super depressed about this. I knew that grad school would be the end of my competitive athletics. I'm merely making a point that so many people seem to miss. There simply are no shortcuts. If you want to be good at something, do so. But count the cost, because it is high.

Wednesday, February 24, 2016

Sampling harness complete

In a perfect world, I'd be further along, but I did get to a tangible milestone today and sleep is important, too. My sampling harness is complete. The architecture is pretty straightforward. As it's basically a command-line app, there's no presentation layer; just a model sitting on top of a data layer.

I've only implemented csv files in the data layer as that's all I'm working with right now. When I start using the full dataset, I'll swap that out for HDFS and host it in the cloud. Aside from basic read/write, the data layer also does the stratification of the data before presenting it to the model. This effectively re-blocks the data in the csv case. When I go to Hadoop, I'll do the stratification as part of the partitioning of the Parquet files so the data layer can leverage the parallelism inherent in the file system.

The model currently implements two sampling methods; a sequential scan of all blocks and a random sampling of blocks. Both run through the whole data set right now. I probably won't bother putting a stopping rule on either since I'm only using them as a baseline for convergence. If I do put a user interface on this thing, I'll factor the sampling methods out into a controller layer. For now, I'm just being careful not to build sampling rules into the true model classes. That separation will likely be important when I go to parallelize things.

It's just a start, but it works fine and the whole thing has unit test coverage. That's actually pretty important because once I start running large data sets through this thing, it's going to be very difficult to validate that each step is working right. Having those test cases to fall back on for regression when making changes will likely save a ton of time.

Next up, obviously, is to code my algorithm (need to come up with a better name than "my algorithm"). Most of it is pretty simple, the only difficulty is in computing the posteriors. I'll probably use a heuristic for starters just to get it working and then work through the math with my adviser when I get back from Texas.

Tuesday, February 23, 2016

Moving on

I'm pretty much ready to shelve the Metropolis algorithm. I'll concede that when very little is known about the distribution of the observations, it does have some merit, but that's not really the problem I'm trying to solve. I'm looking for fast convergence and this is clearly not it. Rather than demonstrate on my data which is too irregular to be a good illustration, I'll show what I mean on a the uniform distribution.


The blue bars show the count of actual observations from a discrete uniform (1,100) distribution. The red bars show the Metropolis count where the current observation, xi, is kept if it is larger than the next, xi+1, with probability xi+1/xi. If xi+1 is larger, Metropolis always moves to the next observation.

As one would expect, Metropolis correctly shifts the sampling roughly in proportion to the mass of the observations (I'm using the term mass to indicate density times the value). However, the convergence is terrible. This is a sample set of 10000 observations on the smoothest possible distribution. One can reasonably ask for pretty nice convergence here. Yet, look at the deviation from the expected value for the counts for each possible value of X:



The actual observations have a few spikes, but are in line with expectations (they'd better be, or I need a new random number generator). Metropolis does great in the area of the distribution where there is very little mass, but look at those giant spikes on the right. That's in the heavy part of the distribution where I can't afford to be off by much. This isn't a fluky set of observations, either. I ran it quite a few times with different seeds and the graphs looked pretty much the same each time. If Uniform is off by this much after 10000 observations, my data is going to need millions of hits to converge. Given the difficulties of random sampling a large database, there's no way that's going to outperform a full scan until the data size is way into the billions. My full data set is that large, but individual queries will often return "only" several million rows. This just isn't looking like a good answer.

So, here's my revised algorithm that simultaneously monitors convergence of the posterior distributions of the extreme value (to avoid missing big observations) and the sum (which is the value we really want):

As before, when loading the data, stratify it by log2 and group into blocks large enough to get some decent mixing of attributes. Might have to encourage such mixing by employing a hash to the block assignment, but I'm not worried about that right now.

Then, on query, start with a prior (ideally from actuals) P(|Xs|) which is the density of the sum of the absolute values of the relevant measures from a block from stratum s (note that, despite the somewhat misleading notation, this is the sum of absolute values, not the absolute value of the sum). Set P(Xs) to be some convenient distribution centered at zero, which I think will work fine as these are generally offsetting accounting entries, though I may need to generate a better prior for the actual sum as well. Then:
  • Pick a stratum s in proportion to E(|Xs|).
  • Sample a block from there to get both the sum of the absolute values |D| and the actual sum D.
  • Update the posteriors P(|Xs| |D) and P(Xs | D). May also need to update the nearby posteriors, but, for now, I'm treating them as independent even though they aren't.
  • Keep going until the (1-α) HDI on all the posteriors P(Xs | D) is within the specified tolerance.
  • Multiply E(Xs) by the number of blocks in each stratum and sum to get your answer.
This is essentially stratified Gibbs sampling with an update on the distribution at each step. That strikes me as obvious enough that I did a cursory review to see if there's already a paper out there somewhere for this and didn't find one (though I did turn up some related papers from Importance Sampling that I should probably read).

For a first cut, I'm not going to stress too much about computing posteriors; I'll just assume something mathematically convenient with a tail to the right. The goal is to get something working in the next few days so I can then refine it while I'm in Texas.