Opioids in MA, 2017

I am surprised that the death rate from drug overdoses is so high on Cape Cod (Barnstable County). The simplistic reading of the opioid crisis usually associates it with deprivation and poverty; the typical narrative does not particularly associate the crisis with affluent or vacation communities. But Barnstable's age-adjusted opioid death rate of 53.8/100,000 puts it at about 2.5x the national average
A couple of technical points:
  • To protect individuals' privacy, the CDC does not report counts in counties with low numbers of deaths. The number of deaths in Franklin County and Nantucket County fall below this threshold, hence they are absent from the above map.
  • The rates of almost all causes of deaths vary with age. Age adjustment is a technique that removes the effects of age from crude rates and allows meaningful comparisons between populations of different age profiles. See the CDC's description of the process here.
  • The map shows rates not raw counts. See the wikipedia entry on choropleths for the justification.
So what is causing these high rates of drug overdoses? It is impossible to tell from the CDC data alone. A good first stop for other data that could contribute to a model is the US census. More on that later. The R-code that produces the above map is available in my git repository here.

Opioid deaths: 2017 figures

The day after my earlier post I found that the CDC had released the 2017 files. Here is the updated version of the plot that shows the histories of death rates in the 12 counties with the highest rates. Click to enlarge.

The numbers just keep going up. I know this is a highly selective snapshot (by definition I am looking at the worst hit counties) but the full national picture is sobering. From the CDC's 2017 report:

  • Average rate in 2017 was 10% higher in 2017 than in 2016
  • Amongst the synthetic opioids (e.g. fentanyl and tramadol) the increase was 45%

Opioid deaths: rates or counts?

The CDC have just updated their annual report "Drug Overdose Deaths in the United States, 1999-2017". The average age-adjusted rate in 2017 is almost 22 per 100,000; up 10% on the 2016 figures. Some individual states are well above these numbers: the highest death rates are in West Virginia (58 per 100,000), Ohio (46), Pennsylvania (44) and District of Columbia (44).

What does this look like at the county level? The CDC make the Detailed Mortality Files—on which the above report is based—available through WONDER, the Wide-ranging OnLine Data for Epidemiologic Research. Unfortunately this has not been updated with the 2017 numbers, so a county-level analysis can only go through to 2016.

The counties with the highest death rates are shown below (click to enlarge).

As a New Mexico resident I am not surprised that Rio Arriba is at the top. Its largest town, Española, has been identified as a "drug capital of America" for a decade or more (for example, see this Forbes article from 2009). The West Virginia counties are not surprising either, given the state figures.

But the figures that really shock me are those for the counties in Maryland and Ohio. Yes, the rates are somewhat lower (though still north of 50 per 100,000) but the populations of these counties are much higher than you see in counties in NM and WV. The numbers of deaths in Baltimore, Montgomery, Butler, and Clermont Counties are—to my eyes—startlingly high.

This raises the question of how to measure "worst". Is it death rate (i.e. deaths per 100,000 population), deaths, or some combination? Although the use of rate has a mathematical appeal, I think it carries an inherent assumption that may not be correct. Specifically, that if county A has 10 times the population of county B, the threshold at which local services get overloaded is also 10 times higher in A than it is in B. I suspect (though I don't have evidence for this) that this is not the case and that the threshold in A would be substantially less than 10 times the threshold in B.

Note: The R code for the above plots and a discussion of assumptions is available here.

Python's magic numbers

I've been troubleshooting a problem that arose while running a pyc file: RuntimeError: Bad magic number in .pyc file

OK, so how do I find the magic number that is in the myfile.pyc and how do I find the magic number that was expected?

The magic number is the first two bytes of the file. I can get them using xxd:

robert@dante ~/MathScripts
$ xxd myfile.pyc | head -n 1
0000000: 330d 0d0a 33c7 925b 4c29 0000 e300 0000  3...3..[L)......

In other words, the magic number of the python36 version of myfile.pyc is 330d 0d0a.

The python27 version has a different magic number:

$ xxd ../MathScripts_py27/myfile.pyc | head -n 1
0000000: 03f3 0d0a 24d4 f75a 6300 0000 0000 0000  ....$..Zc.......

…namely 03f3 0d0a.

Having got the magic number of my byte-code file, how do I discover the magic number that my current python is expecting? The commands are different in python2 and 3. In python2:

robert@dante ~
$ source activate py27
(py27)
robert@dante ~
$ python
Python 2.7.15 |Anaconda, Inc.| (default, May  1 2018, 23:32:55)
[GCC 7.2.0] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import imp
>>> imp.get_magic().encode('hex')
'03f30d0a'

And in python3:

robert@dante ~
$ source activate py36
(py36)
robert@dante ~
$ python
Python 3.6.5 |Anaconda, Inc.| (default, Apr 26 2018, 08:42:37)
[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import imp
>>> imp.get_magic().hex()
'330d0d0a'

There's a full list of magic numbers here.

Circular network, fully connected

Here's my version of a design taken from the Jared Tarbell talk that I referenced earlier.

31 points, radius = [1., 1.8, 2.2, 3.1, 3.5, 3.7, 4.]

The recipe is simple:
* Fully connect points that are equispaced around the circumference of a circle
* Repeat for circles of various radii

This exercise posed an interesting challenge. Full details in my repository here, but the short version is that I thought I could take a clever approach using Processing's scale transform and—well—it didn't pan out.

So it was pretty obvious from the symmetry of the pattern that I would want to do my math in polar coordinates and then convert to cartesian at the last moment. And rather than recalculating the points around the edge of the circle whenever I change the radius of the circle, I can leave the points' r and theta (and their equivalent x and y) untouched and use Processing's scale() transform to change the underlying basis vectors and hence the points' position on the screen. Neat! No recalculations!

Unfortunately, scale() scales everything. It's not just that all the basis vectors are scaled, the actual screen representation of a pixel is scaled too. So if I call scale(8) all of a sudden any line (say) that used to be one pixel wide is now 8 pixels wide. And although I can downscale the resulting line widths, I cannot do so below the new representation of 1 pixel. For example, the minimum value that I can pass to strokeWeight() is 1. This is 1 transformed pixel, not 1 actual pixel. So when I use a large value in scale(), my lines suddenly start to get thick and they stay thick. Here's a detail from an early iteration of the design:

The solution is not to use scale(): instead simply recalculate the points

More curves; or "how do I google something whose name I don't know?"

I have a vague childhood memory of visiting a museum and seeing some Victorian contraption of pendulums and pulleys that would produce lovely swirling, swooshing designs from the pen clamped at its center. The designs weren't Lissajous figures, they were far more complex than that. What were they? And what were the devices that produced them?

Sand pendulums? No, I think that was from "Vision On" with Tony Hart. But searching for "sand pendulums" did take me to the archive of the Bridges Organization, which "oversees the annual Bridges conference on mathematical connections in art, music, architecture, education, and culture". OK, I feel I should have known about that and am ashamed that I did not. Moving on…

My point of entry into the archive was a paper by Douglas McKenna, "From Lissajous to Pas de Deux to Tattoo: The Graphic Life of a Beautiful Loop". Victory! It mentions the device I was trying to recall. The harmonograph. Yes, those are the patterns I remember, those are the devices.

Googling for examples is much easier now I know what the things are called. It turns out that several artists are working in this space and producing delightful work. Here's an example from Rick Hall (it's also the source of the image at the top of this post):

Source: http://rickhall.co.uk/harmonograph-2/

Now I just need to work out how to build one of these creatures.

Spirograph math

Yesterday's post about Mary Wagner's giant spirographs prompted two lines of thought:

  1. Wouldn't it be cool to actually make some of those cogs?
  2. Wouldn't it be cool to prototype them in code?

It didn't take me long to put the first of these on ice. My local maker space has all the necessary kit to convert a sketch on my laptop into something solid but oh my, gears are hard. There seem to be all sorts of ways of getting a pair of gears to grind immovably into one another. I suspect I'd get to explore that space thoroughly before actually getting something to work. I don't have the patience for that. (If you are more patient than me, check out this guide to gear constrution from Make.)

Let's turn to the code then.

There are plenty of spirograph examples out in the wild (the one by Nathan Friend is probably most well known). But if I want to build one of my own, that's going to need some math.

Wikipedia didn't let me down. Apparently the pattern produced by a spirograph is an example of a roulette, the curve that is traced out by a point on a given curve as that curve rolls without slipping on a second fixed curve. In the case of a spirograph we have two roulettes:

  1. A circle rolling without slipping inside a circle (this is a hypotrochoid); and
  2. A circle rolling without slipping outside a circle (this is an epicycloid)

Here's a hypotrochoid:

Source: https://commons.wikimedia.org/wiki/File:HypotrochoidOutThreeFifths.gif

And here's an epicycloid:

Source: https://commons.wikimedia.org/wiki/File:EpitrochoidOn3-generation.gif

The math is straightforward (here and here) so no real barrier to coding this up, right?

Hmm. As usual one line of enquiry opens up several others. Sure, spirographs look nice but they are rather simple, no? What are our options if we want something more complex but still aesthetically satisfying? Wouldn't it be more fun to play with that? I mean the math couldn't be much harder, right?

More on that tomorrow…

Belemnite battlefields

I have been tinkering around with the code in Pearson's "Generative Art" and came up with this and the image above.

The author's intention was to show how complex patterns could spontaneously emerge from the interaction of multiple components that individually exhibit simple behavior. In this case, a couple of dozen invisible disks track across the display; whenever the disks intersect they cause a circle to suddenly appear at the point of intersection. So a little like the bubble chamber that I'd used in physics lab at university. You don't observe the particle directly, you observe the evidence of its passage through a medium. (My code is here if you want to play around with it.)

But what the images really remind me of are the belemnites that Cynthia and I had tried (and failed) to find in the fossil beds of Dorset and Yorkshire.

Belemnites are the fossilized remains of little squid-like creatures. Most of the times you find individuals—a small one looks a little like a shark's tooth—but sometimes you find a mass of them in a single slab, a "belemnite battlefield".

Source: thefossilforum.com

How do these battlefields occur? I found a paper by Doyle and MacDonald from 1993 in which the authors identify five pathways: post-spawning mortality, catastrophic mass mortality, predation concentration, stratigraphical condensation, and resedimentation. I don't know which occurred in the slab shown above but hey, I now know more potential pathways than I did before I read that paper.

And that's why art is fun. One moment you are messing around with RGB codes to get just the right sandy color for the background of a picture, the next you are researching natural phenomena that until an hour ago you never knew even had a name.