Messing Around with DWI Arrests, TABC Mixed Beverage Taxes, and Ridesharing, in Austin, TX

(Playing along at home? Download the .Rmd for a list of libraries you’ll need!)

0.1 Background

So, as you may or may not know, Uber and Lyft, two major ridesharing services, stopped offering rides in Austin, TX as of May 9, 2016 after a contentious public vote which you can read about here: http://www.bizjournals.com/austin/news/2016/05/07/uber-lyft-defeated-in-prop-1-referendum.html

A common argument is that the presence of ride-sharing services such as Uber or Lyft decreases the incidence of drunken driving in a city in which they operate; it sounds reasonable enough to say that if there are more options for getting home after drinking, at least one of them might be more attractive than driving home in a potentially compromised state.

There have been a few articles and a couple papers studying the availability of ridesharing services and drunk driving (see references below), but in my opinion they fall short in one or more ways:

• They fail to count all available records of drunk driving (like if, say, a study only counted drunk driving fatalities)
• They fail to account for trends in drunk driving frequency (like if, say, more people than average drove drunk on certain holidays), or aggregate time-series data too broadly
• They fail to control for exogenous effects on drunk driving frequency (like if, say, the population of a core drinking demographic increased)
• They fail to properly treat the regression discontinuity between data before and after ridesharing was available in Austin (more on this later)

I’ve managed to get a hold of three datasets that I think will be helpful:

1. APD Incident Reports for 2008-2011, Jan-July 2015, Jan-July 2016 (including a number of intoxication-related crimes)
2. TABC Mixed Beverage Tax Receipts for Austin, 2008-2016.
3. APD Chief’s Monthly Reports, 2010-2016, with numbers of DWI and Narcotics arrests.

The APD (Austin Police Department) Incident Reports are publicly available for years 2008-2011, and for ‘Year to Date’ which supposedly includes the last 18 months of reported crime in Austin but actually only includes 8 months of 2016 and a small number of reports from 2015. I was able to obtain an archived version of the data from late 2015, which has more complete data for Jan-Jul 2015.

The APD Chief’s Monthly Report is a high-level summary of crime in Austin for the month, including counts for high-profile crimes including DWI. It has no time-of-day or location information, and no info other than just that a DWI occurred, but is another source of data to compare to the incident reports, which may be spotty. It additionally contains data on narcotics arrests, which we can check to see if we can use it as a control variable - narcotics crimes may follow the same broad trends as DWI, but are potentially unaffected by ridesharing.

It’s worth noting that the study of DWI arrest numbers as a proxy for the incidence of drunk driving is an open area of research - some studies prefer drunk driving crashes, although Austin has not released enough monthly statistics on drunk driving crashes to make it a useful for what I want to do. We’ll assume DWI arrests are a viable proxy for drunk driving here. For a more detailed treatment of this, see the appendix.

The TABC (Texas Alcoholic Beverage Commission) Mixed Beverage Tax is a fixed tax charged on alcoholic beverages sold by a holder of a Mixed Beverage License in the state of Texas. I suspect that Mixed Beverage Tax revenues are associated with Units of Alcohol Consumed in Austin, and we’ll check to see if any relationship exists between these tax revenues and DWI counts.

0.2 Dates of Interest

• June 3, 2014: Lyft launches in Austin. (https://twitter.com/AustinLyft/status/473857057720766464)

• June 4, 2014 : Uber launches in Austin. (https://twitter.com/Uber_ATX/status/474199740217716736)

• May 9, 2016: Uber and Lyft stop giving rides in Austin.

We’ll call June 4 and May 9 treatment dates - the inflection points we’ll investigate to see if the data on either side is substantially different. When dealing with monthly data we’ll have to say that June 2014 is the first month of the treatment and April 2016 is the last month.

You can find the same data I used at these sites:

• https://data.austintexas.gov/browse?q=apd&sortBy=relevance
• https://www.comptroller.texas.gov/transparency/open-data/search-datasets/
• https://austintexas.gov/page/chiefs-monthly-reports
• https://durs.carto.com/tables/apd_incident_extract_2015csv/public

Or you can pull it from my Github, http://www.github.com/ianwells/atx-dwi

Alright, Let’s get to it.

1 The Data

1.1 DWI Arrests

It’s always a good idea to just poke around in your data and see what’s going on. Here’s a quick glimpse of what the raw APD data looks like:

sample_n(dwi.raw,10)

Cool. We have a bunch of crimes, pre-filtered to intoxication-related ones and the time they approximately occurred. I’ve verified that time and location are unique for each crime (no double-bookings, although there are a few rare cases where two drunk drivers were arrested at the same time), so we can be confident that one incident report = one confirmed instance of drunk driving.

You can check this against the full list of crimes available to commit in Austin, which is available on that same Austin Data Portal site, I’m not including it here because it’s quite long, please trust that I’ve got all of them.

levels(dwi.raw$crime)

[1] "CRASH/INTOX MANSLAUGHTER"     "CRASH/INTOXICATION ASSAULT"  
[3] "DRIVING WHILE INTOX / FELONY" "DUI - AGE 17 TO 20"          
[5] "DWI"                          "DWI  .15 BAC OR ABOVE"       
[7] "DWI 2ND"                     

CRASH/INTOX MANSLAUGHTER

CRASH/INTOXICATION ASSAULT

DRIVING WHILE INTOX / FELONY

DUI - AGE 17 TO 20

DWI

DWI  .15 BAC OR ABOVE

DWI 2ND

If you’re familiar with this data set you’ll notice I’ve removed a few crimes:

1. BOATING WHILE INTOXICATED - This rare crime doesn’t seem like it could be avoided by hailing an Uber.
2. DWI - DRUG RECOGNITION EXPERT - You don’t have to have any alcohol in your system to be charged with a DWI: in this case, an officer certified in recognizing the effects of some other substance has decided you shouldn’t be driving. While I imagine you could avoid trouble by hailing a ride in this case, I don’t think I could relate it to bar sales, so this stat is going to be relegated to some other study.
3. CRASH/LEAVING THE SCENE - This is a tricky one. Certainly a non-zero percentage of hit-and-run offenses are committed by drunk drivers, but this charge doesn’t offer any more details, and worse, there are more of this kind of crime than all the others put together. How should we count them, if at all? For now, let’s leave them all out, as they’re not mentioned specifically in the chief’s monthly report. We might be able to project how many of these crimes involve alcohol at a later date, based on when and where they occurred, but we’ll leave that for another analysis. See the appendix for further discussion.

Let’s take a look at the relative frequency of these offenses. Periodicity can be exploited later on.

qplot(hour(date), data=dwi.raw, geom="histogram", bins = 24)

Perhaps not surprisingly, nearly all DWI offenses occur between 8 PM and 5 AM.

qplot(wday(date), data=dwi.raw, geom="histogram", bins = 7)

Again, looks like Fridays and Saturdays are popular days to get caught drunk driving.

We’re going to look at the whole series, and then a histogram of monthly data to see if there’s a seasonal trend (you might expect more DWI’s to occur in December, January, or July months with holiday weekends that involve a lot of drinking)

Note: we have incomplete data for July 2015 and August 2016, as well as an old offense in 2014, so we’ll remove those months now. We’ll leave 2015 and 2016 out of the monthly histogram because we only have data from a few months for those years.

dwi.raw$year <- year(dwi.raw$date)
dwi.raw$month <- month(dwi.raw$date)
  
dwi.monthly <- dwi.raw %>% group_by(year,month)
dwi.monthly <- summarize(dwi.monthly,count=n())
dwi.monthly$date <- as.Date(paste0(dwi.monthly$year,'-',dwi.monthly$month,'-01'))

dwi.monthly <- dwi.monthly[-56,] #remove partial July 2015 data
dwi.monthly <- dwi.monthly[-63,] #remove partial August 2016 data
dwi.monthly <- dwi.monthly[-49,] #remove lone 2014 arrest

label.months <- c('J','F','M','A','M','J','J','A','S','O','N','D')

ggplot(data=dwi.monthly, aes(x = date, y = count, group = year(date), color = factor(year(date)))) + geom_point()

dummy <- filter(dwi.monthly,date < '2015-01-01') %>% group_by(month(date)) %>% summarize(count = sum(count)) %>% select(count) 
dummy1 <- mean(dummy$count) #1824
dummy2 <- sd(dummy$count) #99

ggplot(data = filter(dwi.monthly,date < '2015-01-01'), aes(factor(month(date)),weight = count) ) + geom_bar() + geom_hline(yintercept=1824, color = 'red') +
geom_hline(yintercept=1724, color = 'gray') +
geom_hline(yintercept=1924, color = 'gray')

So we can get a sense that, for the years 2008-2011, monthly DWI’s seem to be decreasing (good job APD), and as for a monthly trend we see an annual bump in December-January, something in March, as well as something going on in August-October, although all seem to be within one standard deviation of our monthly mean.

Now, I’m curious to see how the Chief’s Monthly Report compares.

apd.raw$date <- as.Date(apd.raw$date)
apd.raw$year <- year(apd.raw$date)
apd.raw$month <- month(apd.raw$date)

apd.monthly <- apd.raw
  
ggplot(data=apd.monthly, aes(x = date, y = count, color = factor(year))) + geom_point() 

dummy <- apd.monthly %>% group_by(month(date)) %>% summarize(count = sum(count)) %>% select(count) 
dummy1 <- mean(dummy$count) #3290
dummy2 <- sd(dummy$count) #145

ggplot(data = apd.monthly, aes(factor(month(date)),weight = count) ) + geom_bar() + geom_hline(yintercept=3290, color = 'red') +
geom_hline(yintercept=3145, color = 'gray') +
geom_hline(yintercept=3435, color = 'gray')

This has similar bumps, with a weaker September. Let’s plot both at the same time.

both.monthly <- merge(apd.monthly,dwi.monthly, by = 'date', all = TRUE);
both.monthly$apd <- both.monthly$count.x
both.monthly$dwi <- both.monthly$count.y

both.monthly <- both.monthly %>% select(date,apd,dwi)
both.melted <-melt(both.monthly, id.vars="date")

ggplot(data=both.melted, mapping = aes(x=date,y=value, color = variable)) + geom_point()

Okay, so it looks like the chief’s report is a little different than the incident report database. That’s annoying, but based on the disclaimers on both websites, and that the chief’s report is updated more often, I’m going to call it and say the chief’s data is what I’m using moving forward. I’d like those extra two years of data that the database has though - and since the data is so historic, and the difference is so consistent, what I’m going to do is nudge the data up slightly so that it lines up.

both.monthly$diff <- both.monthly$apd - both.monthly$dwi
mean(filter(both.monthly,date < '2012-01-01', date > '2010-01-01')$diff)

[1] 14.86957

both.monthly.nudge <- both.monthly
both.monthly.nudge$dwi = both.monthly.nudge$dwi + 15 #approximate mean difference

both.melted.nudge <-melt(both.monthly.nudge, id.vars="date")
ggplot(data=both.melted.nudge, mapping = aes(x=date,y=value, color = variable)) + geom_point()

Alright, that’s a little cleaner. We’ll move forward with this set of monthly data until I can get more comprehensive arrest records, being able to use weekly periodicity might be helpful.

1.2 Narcotics Arrests

The Chief’s Monthly Report lists all Part 1 Offenses (according to the Uniform Crime Reporting standards established by the FBI, mostly violent felonies) and a few key Part 2 Offenses, namely DWI, narcotics, prostitution, and weapons-law related crimes. Prostitution and weapons crimes occur sparsely compared to DWI and narcotics, so we’ll import the narcotics arrest data in hopes that it will allow us to control for sweeping effects across all crime categories - things like the number of officers working, motivation of the police force, population of the city, etc.

ggplot(data=apd.monthly, aes(x = date, y = narc, color = factor(year))) + geom_point() 

qplot(factor(month(date)), data=apd.monthly, geom="bar", weight = narc)

Not a lot to say here - there’s a strong January but not a lot else that lines up with the DWI data.

I’ve saved out our final data set here to save time.

dwi <- read.csv('data/dwi-final.csv')
#sapply(dwi,class)
dwi$date <- as.Date(dwi$date)
dwi$count <- as.numeric(as.character(dwi$count)) #cast strictly
dwi$narc <- as.numeric(as.character(dwi$narc))
dwi$treatment <- as.factor(dwi$treatment) #0: pre uber, 1:uber, 2:post uber

1.3 TABC Mixed Beverage Tax

Let’s have a look at what’s going on with the TABC data. Now, I’ve already done a bit of work with this data set in the past (see http://tabcmap.s3-website-us-east-1.amazonaws.com), so I’ve just pulled values from my database of monthly revenues calculated from monthly tax receipts.

Here’s what a few rows look like:

sample_n(tabc.raw,10)

Let’s do the same series and histogram plots to get a sense of what’s in there.

There’s a lone data point from 1999, let’s clean that up too, and let’s scale revenues into the millions of dollars.

tabc.monthly <- tabc.raw[-40,] #remove lone 1999 point

tabc.monthly$date <- as.Date(tabc.monthly$date)
tabc.monthly$year <- year(tabc.monthly$date)
tabc.monthly$month <- month(tabc.monthly$date)
tabc.monthly$rev <- round(tabc.monthly$rev/1000000,2)

ggplot(data=tabc.monthly, aes(x = date, y = rev, group = year, color = factor(year))) + geom_point()

qplot(factor(month(date)), data=tabc.monthly, geom="bar", weight = rev)

Look at that growth! Notice that spike in March? That’s probably SXSW, a major annual music, technology, and film conference. Saint Patrick’s day is in March, too. See that spike in October? I bet that’s Austin City Limits, a major music festival, and also Halloween. This data is already good to go, so let’s rename it to keep things clean and keep going.

tabc <- select(tabc.monthly,date,rev)

1.4 DWI per $MM

Let’s see about a new variable - DWI per Millions of Dollars Bar Revenue. I think that this variable may be useful as it incorporates things you would expect to vary with bar revenue - drinking-age population and popularity of drinking - into the DWI count. All things the same, a town with a lower DWI-per-drink-consumed number is going better than one with a higher ratio. We may be able to link ridesharing to an improvement in this ratio more convincingly than to just raw DWI’s.

dwi <- merge(dwi,tabc, by = 'date')

dwi$count_per_mm <- dwi$count/dwi$rev

ggplot(data=dwi, aes(x = date, y = count_per_mm, group = year(date), color = factor(year(date)))) + geom_point()

ggplot(data=dwi, aes(x = rev, y = count, color = year(date))) + geom_point() + geom_smooth(method = "lm")

summary(lm(data = dwi, count ~ rev))

Call:
lm(formula = count ~ rev, data = dwi)

Residuals:
     Min       1Q   Median       3Q      Max 
-124.302  -36.895    2.742   37.032  112.456 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 471.91537   21.21465  22.245   <2e-16 ***
rev           0.01024    0.47007   0.022    0.983    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 51.7 on 105 degrees of freedom
Multiple R-squared:  4.519e-06, Adjusted R-squared:  -0.009519 
F-statistic: 0.0004745 on 1 and 105 DF,  p-value: 0.9827

So while DWI/$MM is decreasing over time, it definitely has periods where it’s still dominated by a rising DWI rate, and it would appear that there’s no direct correlation between DWI counts and bar revenue - which makes sense, as they’re both clearly non-stationary.

qplot(factor(month(date)), data=dwi, geom="bar", weight = count_per_mm)

This variable does however have a much more obvious yearly peak in January - we may be able to exploit this entangled variable more easily than just the raw DWI, if we can isolate the effect of ridesharing on TABC revenue (which we’ll talk about later).

2 The Analysis

Alright, so now that we’ve got some data, what can we do with it? How can we establish one way or the other that the advent of ridesharing affected the number of DWI arrests in Austin, or the number of DWI’s per dollar of alcohol sold? There are a few things I’m going to look at.

2.1 Regression Discontinuity

Typically, when you set out to prove something with science, you sit down and design an experiment. We don’t have the resources to do a city-scale double-blind ridesharing/drunk-driving experiment, so we’ll have to use the next best things: math and historical data.

Regression Discontinuity Design is a method of establishing causality without the benefit of experimental design. The key assumption is that, with regards to a time series, the data on either side of a treatment inflection are otherwise similar - nothing else special has happened. If you can get behind that, you could claim that if there’s a kink in your data, it’s because of the treatment you’re investigating.

Using this technique on time series data can be tricky, especially when you’re working with 8 years of data, because surely something else besides your treatment happened over that time frame. Let’s remind ourselves what the DWI count looks like:

ggplot(data = dwi, mapping = aes(x=date,y=count, color = factor(year(date)))) + geom_point()

You can probably see a lot of jumps in the data - from 2009-2011 there’s a rapid rise in DWI, followed by a swift drop and then a leap up again in 2013, and a gradual drop into 2016. We can’t honestly include all that data, so we’re going to limit things to about two years on either side of the first treatment (Uber and Lyft open for business in Austin), and the 8 months we have after the second treatment (Uber and Lyft cease operations).

First, we’ll plot regressions on either side of the treatment and see how they look - with just lm to start and then with a specialized library.

ggplot(data = filter(dwi,date > '2012-05-01'),mapping = aes(date, count, color = treatment) )+ geom_point() + geom_smooth(method = 'lm')

Visually, you might be tempted to get a sense of the effect of ridesharing over the period marked in green, but let’s keep going. We’re going to use the rddtools library and see if there is a statistically significant change in the regressions before and after each discontinuity, we’ll also compare rddtools to the output of lm .

dwi$index <- seq.int(nrow(dwi))

dwi.w <- filter(dwi, date < '2016-05-01', date > '2012-05-1')
dwi.rdd.w <- rddtools::rdd_data(y = dwi.w$count,x = dwi.w$index,cutpoint = 77)

dwi.wo <- filter(dwi, date < '2016-12-31', date > '2014-05-30')
dwi.rdd.wo <- rddtools::rdd_data(y = dwi.wo$count,x = dwi.wo$index,cutpoint = 100)

lm.w<-lm(count~I(index-77)*treatment,dwi.w) 
reg_para.w <- rdd_reg_lm(rdd_object = dwi.rdd.w, order = 1)

lm.wo<-lm(count~I(index-8)*treatment,dwi.wo) 
reg_para.wo <- rdd_reg_lm(rdd_object = dwi.rdd.wo, order = 1)


print(reg_para.w)

### RDD regression: parametric ###
    Polynomial order:  1 
    Slopes:  separate 
    Number of obs: 47 (left: 23, right: 24)

    Coefficient:
  Estimate Std. Error t value Pr(>|t|)
D  -6.3097    18.7822 -0.3359   0.7386

plot(reg_para.w)

#summary(lm.w) #for a sanity check

For the first treatment, RDD tells us (through the P-value of the regression, 0.7386) the that the drop is small enough relative to the variance in the data that any bump and change in the slope of the regression is too small to be attributed to anything other than randomness - at least when modeled with a straight line. We’ll consider that a straight line isn’t the best way to fit this data later on - remember we spent some time investigating periodic patterns.

It’s also worth noting that our primary assumption here - that everything is the same on either side of the treatment - might not hold, as the popularity and capacity of ridesharing probably grew over time. RDD may not be a rigorous way to inspect this side of the treatment.

print(reg_para.wo)

### RDD regression: parametric ###
    Polynomial order:  1 
    Slopes:  separate 
    Number of obs: 30 (left: 22, right: 8)

    Coefficient:
  Estimate Std. Error t value Pr(>|t|)   
D  -64.579     23.084 -2.7975 0.009565 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(reg_para.wo)

#summary(lm.wo) #for a sanity check

For the second treatment, the jump in the regression is more pronounced and the slope changes wildly - our P values are low enough (< 1% chance that this result is due to randomness) to conclude that the change may be due to the discontinuity in May 2016 and not random chance. I think this might be a more valid case for RDD than the first side of the treatment - Uber and Lyft quite literally stopped offering ridesharing the Monday after the vote came down. This is a point in favor of ridesharing, but we’re not out of techniques yet.

2.2 Exponential Smoothing

Our regression discontinuity design fit 1st degree polynomials to the DWI data, but what if the behavior is more complicated than that? We observed earlier that there may be a periodic component to the data, what if ridesharing disrupts busy DWI months but not slower months, would a simple linear regression uncover that? As we mentioned, Uber and Lyft’s ridesharing capacity grew as they added drivers (certainly a metric I’d like to get my hands on), how can we investigate a higher-degree change in DWI?

Exponential Smoothing is a technique used in a whole class of methods for predicting future values in a time series based on manipulations of past values. Holt-Winters Forecasting is a method of predicting time-series data that has strong periodic components by separately smoothing long-term trend, short-term seasonality, and local components to create a forecast that takes into account predictable variations in the data. We saw above that our DWI data probably has some yearly trends in it - the data shows it a little, and you can Google around for articles about DWI in Austin and see for yourself that January is a special month for the police (look for ‘no-refusal weekend’). We also got a sense that a few years saw a continuous reduction in DWI arrests, perhaps as broad anti-DWI measures kicked in. We’ll decompose the DWI and DWI/$ really quick and you can see for yourself:

plot(stl(ts(dwi$count, start = c(2008, 1), frequency = 12), s.window="periodic"))

plot(stl(ts(dwi$count_per_mm, start = c(2008, 1), frequency = 12), s.window="periodic"))

You can see seasonality contributes to a swing of 70 DWI arrests, or 2.5 arrests per $MM. Additionally, you can see a highly smoothed trend for all the data over our entire time fame, and the remainder that this decomposition doesn’t think is attributable to a pattern.

2.2.1 Pre-Uber, and Validating Holt-Winters

What I’m going to do is feed DWI data from January 2008 - May 2014 (the last month before Uber and Lyft were widely available in Austin) into a Holt-Winters model (which will use this same stl seasonal decomposition) and forecast DWI counts for the next few months (we can roll the data back to 2008 again because we’re not concerned with local effects so much as trends inherent in DWI counts).

I’m then going to compare those to the actual data, with the reasoning that since the model takes into account yearly cycles and broad trends, it will come up with a prediction for DWI counts as if ridesharing never appeared. We can call the difference between the prediction and the real data the number of DWI’s ridesharing prevented - or maybe caused.

Let’s get confident in this method - first I’m going to take some untreated months and see if we can predict them reliably.

dwi.pre <- filter(dwi,date < '2013-10-31')
dwi.post <- filter(dwi,date > '2013-10-31')

dwi.pre.count.ts <- ts(dwi.pre$count, start = c(2008, 1), frequency = 12)
dwi.post.count.ts <- ts(dwi.post$count, start = c(2013, 6), frequency = 12)
#dwi.post.count.ts <- ts(dwi.post$count, start = c(2013, 11), frequency = 12) futher back forecast

dwi.pre.count.hw <- HoltWinters(dwi.pre.count.ts)

forecast <- predict(dwi.pre.count.hw, n.ahead = 6, prediction.interval = T, level = 0.9)

#acf((dwi.pre.count.hw$fitted - dwi.pre.count.ts)[,1]) #confirm no autocorrelations in residual
#hist((dwi.pre.count.hw$fitted - dwi.pre.count.ts)[,1],30) #confirm residual looks normal

plot(dwi.pre.count.hw, forecast); lines(dwi.post.count.ts,col='green')

effect <- forecast - dwi.post.count.ts

#ttest for 12month
#pre.residual <- as.numeric((dwi.pre.count.hw$fitted - dwi.pre.count.ts)[,1])
#post.residual <- as.numeric(effect)
#t.test(pre.residual,post.residual) #.59 too high

pct.effect <- cumsum(effect)/cumsum(forecast)
pct.effect[6]

[1] 0.0474804

What we’re looking at here is real data before our treatment (in black), and a fitted model with a forecast extending into our treatment (in red), with two 90% confidence intervals (these are called prediction intervals in this context, they’re in blue), and real treatment data (in green). We’ll plot all our Holt-Winters models like this.

Not bad. The model is a little rough - it fits it approximately but fails to catch all the peaks - the data may not have been periodic enough for this model to be effective, or there may be too much non-seasonal noise.

A failure to fully model the seasonality would have shown up as an autocorrelation or non-normality in the residuals - the difference between the model’s fit and the real data, over the set of data we used to train it. If you’re curious you can uncomment the acf and hist calls and see for yourself - otherwise trust me that I’ve verified that this model fully captures the seasonality available in the training set.

Using the data from January 2008 to October 2013, it looks like we can predict cumulative DWI’s for the next 6 months within 1.5%. That’s pretty close, and that should make us feel better about the numbers coming out, with the caveat that cumulative DWI predictions over longer time frames may be more accurate than picking individual months. (If you’re reading the commented-out bits of code, you’ll note that I also validated predictions 12 months out within 1%). Let’s move the end of the training data to May 2014 (the last month without Uber and Lyft) and see how things look (we’ll improve the model a little more by adding those extra 6 months to the training set as well).

dwi.pre <- filter(dwi,date < '2014-05-31')
dwi.post <- filter(dwi,date > '2014-05-31')

dwi.pre.count.ts <- ts(dwi.pre$count, start = c(2008, 1), frequency = 12)
dwi.post.count.ts <- ts(dwi.post$count, start = c(2014, 6), frequency = 12)

dwi.pre.count.hw <- HoltWinters(dwi.pre.count.ts)

forecast <- predict(dwi.pre.count.hw, n.ahead = 24, prediction.interval = T, level = 0.9)

#acf((dwi.pre.count.hw$fitted - dwi.pre.count.ts)[,1]) #confirm no autocorrelations in residual
#hist((dwi.pre.count.hw$fitted - dwi.pre.count.ts)[,1],30) #confirm residual looks normal

plot(dwi.pre.count.hw, forecast); lines(dwi.post.count.ts,col='green')

effect <- forecast - dwi.post.count.ts

pct.effect <- cumsum(effect)/cumsum(forecast)
pct.effect[12]

[1] 0.07247797

pct.effect[3]

[1] 0.01613347

For the first 12 months Uber and Lyft were available in Austin, the real cumulative DWI count was 440 arrests less than the predicted count - an effect of about 7% (over the whole 24-month period, we see 8%). For the first 3 months, where our error bars are much tighter (and perhaps while Uber and Lyft were still growing in capacity), we’re only 24 DWI’s less than predicted - a decrease of 1.6%.

Now you’re probably wondering - the real values are still well within these prediction bands, what does that mean? If our real data was a few points, it would mean we probably couldn’t count on it to tell us anything - it could be well within the model’s tolerance for random chance, and we wouldn’t be seeing a treatment effect, just a coincidence. But we have a lot of points, and all of them are below the forecast value - can we tell if these values are meaningfully different than the forecast?

We can. Since we’ve validated that our residuals are normally distributed, we can use a t-test to check if the fit-vs-training residuals are different from the forecast-vs-real residuals, which would tell us that the treatment effect we’re seeing is real, with some degree of confidence.

pre.residual <- as.numeric((dwi.pre.count.hw$fitted - dwi.pre.count.ts)[,1])
post.residual <- effect[,1]
t.test(pre.residual,post.residual)

    Welch Two Sample t-test

data:  pre.residual and post.residual
t = -4.4255, df = 44.764, p-value = 6.102e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -69.96023 -26.19307
sample estimates:
mean of x mean of y 
-6.507904 41.568743 

#qplot(as.numeric(pre.residual), geom="density")
#qplot(as.numeric(post.residual), geom="density")

The p-value of 0.00006 tells us these results are significant. We’ll keep an eye on these behind the scenes.

We can try a potentially better-fitting model (really, optimized for a different criterion) from the exponential smoothing family of models with ets.

dwi.pre.count.ets <- ets(dwi.pre.count.ts, model = 'ZZZ')
forecast.ets <- predict(dwi.pre.count.ets, n.ahead = 24, prediction.interval = T, level = 0.9)

plot(forecast.ets); lines(dwi.post.count.ts,col='green')

effect.ets <- forecast.ets$mean - dwi.post.count.ts

#pre.residual <- as.numeric((dwi.pre.count.ets$fitted - dwi.pre.count.ets)[,1])
#post.residual <- as.numeric(effect.ets)
#t.test(pre.residual,post.residual) #e-7

pct.effect.ets <- cumsum(effect.ets)/cumsum(forecast.ets$mean)
pct.effect.ets[12]

[1] 0.08457409

pct.effect.ets[3]

[1] 0.05609513

This plot looks like our Holt-Winters plots, but omits the fitted curve before the treatment, and the forecast is in blue.

With much tighter prediction intervals, ets suggests a drop of 8.5% over 12 months (9.7 over the full 24) and 5.6% over 3 months.

I’d like to look at our DWI per Revenue variable, but first we have to answer, how much did ridesharing affect bar revenues? If it’s a lot, we’ll have a hard time attributing a change in DWI/$ to a change in DWI’s.

dwi.pre.rev.ts <- ts(dwi.pre$rev, start = c(2008, 1), frequency = 12)
dwi.post.rev.ts <- ts(dwi.post$rev, start = c(2014, 6), frequency = 12)

dwi.pre.rev.ets <- ets(dwi.pre.rev.ts, model = 'ZZZ')
forecast.ets <- predict(dwi.pre.rev.ets, n.ahead = 24, prediction.interval = T, level = 0.9)

plot(forecast.ets); lines(dwi.post.rev.ts,col='green')

effect.ets <- forecast.ets$mean - dwi.post.rev.ts

#pre.residual <- as.numeric((dwi.pre.rev.ets$fitted - dwi.pre.rev.ts))
#post.residual <- as.numeric(effect.ets)
#t.test(pre.residual,post.residual) #.001

pct.effect.ets <- cumsum(effect.ets)/cumsum(forecast.ets$mean)
pct.effect.ets[3]

[1] -0.02926871

pct.effect.ets[12]

[1] -0.03358013

Looks like we see a 3% rise in TABC sales due to ridesharing - which seems easy to reason around, if you don’t have to worry about driving you could be drinking more. I don’t want to spend too much time focusing on TABC sales and ridesharing, but I’ll mention that according to this, that’s an increase of about $68,000,000 in revenue to local businesses over 24 months, or roughly $9,500,000 in extra taxes.

dwi.pre.count_per_mm.ts <- ts(dwi.pre$count_per_mm, start = c(2008, 1), frequency = 12)
dwi.post.count_per_mm.ts <- ts(dwi.post$count_per_mm, start = c(2014, 6), frequency = 12)

dwi.pre.count_per_mm.ets <- ets(dwi.pre.count_per_mm.ts, model = 'ZZZ')
forecast.ets <- predict(dwi.pre.count_per_mm.ets, n.ahead = 24, prediction.interval = T, level = 0.9)

plot(forecast.ets); lines(dwi.post.count_per_mm.ts,col='green')

effect.ets <- forecast.ets$mean - dwi.post.count_per_mm.ts

#pre.residual <- as.numeric((dwi.pre.count_per_mm.ets$fitted - dwi.pre.count_per_mm.ts))
#post.residual <- as.numeric(effect.ets)
#t.test(pre.residual,post.residual) #.0007


pct.effect.ets <- cumsum(effect.ets)/cumsum(forecast.ets$mean)
pct.effect.ets[12]

[1] 0.1334951

pct.effect.ets[3]

[1] 0.1029293

ETS tells us a whopping 13% drop over 12 months (10% over 3, 17% over 24) in DWI/$MM is attributable to ridesharing - well beyond the change we saw from just the rise in $MM.

2.2.2 Post-Uber

Let’s try the other side - fit a model from data starting from Uber’s first month, and then stopping on Uber’s last month. If the prediction is lower than the real data, we might be able to say that ridesharing was having an effect, but there was some other factor contaminating the data from 2008-2014 (say an organized campaign against drunk driving or something).

We’re going to cheat a tiny bit because we need another sample of data (Holt-Winters requires 2 full years of data here) - so I’m going to say Uber actually opened for business a month earlier than it did. Hopefully the first month was relatively slow as people started adopting to ridesharing, and it won’t make much of a difference to stretch the effect out over one more month in addition the 23 months Uber and Lyft were in Austin.

dwi.during <- dwi %>% filter(date >= '2014-5-01') %>% filter(date < '2016-05-01')
dwi.exit <- filter(dwi,date >= '2016-05-01')

dwi.during.count.ts <- ts(dwi.during$count, start = c(2014, 5), frequency = 12)
dwi.exit.count.ts <- ts(dwi.exit$count, start = c(2016, 5), frequency = 12)

dwi.during.count.hw <- HoltWinters(dwi.during.count.ts)

forecast <- predict(dwi.during.count.hw, n.ahead = 12, prediction.interval = T, level = 0.9)
plot(dwi.during.count.hw, forecast); lines(dwi.exit.count.ts,col='green')

effect <- forecast - dwi.exit.count.ts

#pre.residual <- as.numeric((dwi.during.count.hw$fitted - dwi.during.count.ts))
#post.residual <- as.numeric(effect)
#t.test(pre.residual,post.residual) #e-9

pct.effect <- cumsum(effect)/cumsum(forecast)
pct.effect[7]

[1] -0.04835785

pct.effect[3]

[1] -0.06414393

This looks better - the model fits better (we’ll stick with this over ets), the confidence interval is tighter, and the treatment effect is 80 DWI’s over the first 3 months, a decrease of about 6%. In other words, not having Uber and Lyft may have resulted in an additional 80 arrests in the first 3 months. After 7 months (taking us to the end of our data), we observe a smaller 5% decrease - perhaps we can attribute this to the smaller ridesharing services (like Fare) that sprung up in Austin after Uber and Lyft stopped offering rides.

This gives us a number to associate with the qualitative change we saw with rdd earlier. Cool!

Let’s do the same with DWI/$MM:

dwi.during.rev.ts <- ts(dwi.during$rev, start = c(2014, 5), frequency = 12)
dwi.exit.rev.ts <- ts(dwi.exit$rev, start = c(2016, 5), frequency = 12)

dwi.during.rev.hw <- HoltWinters(dwi.during.rev.ts)

forecast <- predict(dwi.during.rev.hw, n.ahead = 12, prediction.interval = T, level = 0.9)
plot(dwi.during.rev.hw, forecast); lines(dwi.exit.rev.ts,col='green')

effect <- forecast - dwi.exit.rev.ts

#pre.residual <- as.numeric((dwi.during.rev.hw$fitted - dwi.during.rev.ts))
#post.residual <- as.numeric(effect)
#t.test(pre.residual,post.residual) #0.0004

pct.effect <- cumsum(effect)/cumsum(forecast)
pct.effect[7]

[1] 0.07915371

pct.effect[3]

[1] 0.06589074

Bad news for bars - Holt-Winters suggests that revenues were depressed by 6.5% over the first 3 months without Uber and Lyft. That’s $12,000,000 - and nearly $1,700,000 in lost tax revenue. It’s $34,000,000 over the whole 7 months we have data for.

dwi.during.count_per_mm.ts <- ts(dwi.during$count_per_mm, start = c(2014, 5), frequency = 12)
dwi.exit.count_per_mm.ts <- ts(dwi.exit$count_per_mm, start = c(2016, 5), frequency = 12)

dwi.during.count_per_mm.hw <- HoltWinters(dwi.during.count_per_mm.ts)

forecast <- predict(dwi.during.count_per_mm.hw, n.ahead = 12, prediction.interval = T, level = 0.9)
plot(dwi.during.count_per_mm.hw, forecast); lines(dwi.exit.count_per_mm.ts,col='green')

effect <- forecast - dwi.exit.count_per_mm.ts

#pre.residual <- as.numeric((dwi.during.count_per_mm.hw$fitted - dwi.during.count_per_mm.ts))
#post.residual <- as.numeric(effect)
#t.test(pre.residual,post.residual) #0.018 <- this one is a little high but still acceptable


pct.effect <- cumsum(effect)/cumsum(forecast)
pct.effect[7]

[1] -0.1268825

pct.effect[3]

[1] -0.1481808

14% higher DWI/$MM over the first 3 months without Uber, or 12% through 2016, again exceeding the effect of not having ridesharing on revenue alone.

We’ve got one more trick up our sleeves.

2.3 Bayesian Structural Time-Series

2.3.1 With Google CausalImpact

Google has a great library for testing the causal effect of a given treatment on a time series. Designed for testing the effectiveness of advertising, I think this library will help us out - it’s designed for inference on non-experimental data, which is what we’ve got.

For this, in contrast to regression discontinuity or exponential smoothing methods, we require a response and predictor - one thing that we think a treatment affected and one thing that we are confident it didn’t. Let’s use DWI count as a response, and TABC revenue as a control, as it appears TABC sales growth is less affected than DWI, per our experiments in the previous section, and also narcotic arrests, which we don’t think will be shifted at all.

We need a well-behaved stretch of data (we’re just studying our treatment, not other inflections) to use as a control - so let’s re-use the same cutoff we used for our regression discontinuity design (2012-05-01). Let’s plot it to see if it it looks like anything strange happened to our control data over our period of interest.

uber <- geom_rect(mapping = aes(xmin=as.Date('2014-06-01'), xmax = as.Date('2016-05-01')), ymin=-Inf, ymax=+Inf, color="grey20", alpha=0.3, inherit.aes = FALSE)

ggplot(data=melt(select(filter(dwi,date>="2012-04-01"),narc,count,rev,date),id = 'date'), mapping = aes(x=date, y=value, color = variable)) + uber + geom_point() + geom_smooth(method = 'loess', se = FALSE)

Not the best, but it’s what we’ve got. I’m especially concerned by the trend in TABC revenue and the little uptick in narcotics arrests in late 2016, and generally concerned about the jumpiness of the data. Major drug busts can change the availability of drugs in a city overnight, and could reduce the number of drug arrests in subsequent months. If CausalImpact doesn’t like it it should tell us, we’ll see how it goes.

(If you’re running this at home note that you’ll need the development version of devtools to install CausalImpact due to a bug in the CRAN version of devtools. Get it with devtools::install_github("hadley/devtools") )

2.3.2 TABC as a Control

dwi.ci <- select(filter(dwi,date>="2012-05-01"),date,count,rev)

pre.ci <- as.Date(c("2012-05-01", "2014-05-31"))
post.ci <- as.Date(c("2014-06-01", "2016-05-01"))

impact <- CausalImpact(dwi.ci, pre.ci, post.ci,  model.args = list(prior.level.sd = 0.1,nseasons = 12, season.duration = 1)) 
#i've turned up the sd since this data is more scattered than web traffic data
#notice we've set it to look for a trend with a period of one year

plot(impact)

impact

Posterior inference {CausalImpact}

                         Average        Cumulative    
Actual                   463            11114         
Prediction (s.d.)        507 (15)       12172 (355)   
95% CI                   [479, 539]     [11485, 12926]
                                                      
Absolute effect (s.d.)   -44 (15)       -1058 (355)   
95% CI                   [-75, -15]     [-1812, -371] 
                                                      
Relative effect (s.d.)   -8.7% (2.9%)   -8.7% (2.9%)  
95% CI                   [-15%, -3%]    [-15%, -3%]   

Posterior tail-area probability p:   0.002
Posterior prob. of a causal effect:  99.8%

For more details, type: summary(impact, "report")

CausalImpact attributes a 8.7% drop in DWI’s to ridesharing (with 99.9% certainty, it says), pretty close to what the last model said, but can we trust it? We validated Holt-Winters, let’s validate this as well.

Our control data might have trends and other unrelated treatments in it, or this model might not be as good at handling seasonality. If we give CasualImpact a ‘fake’ treatment date when we know nothing happened - say, in the end of 2013 - it shouldn’t report any change in DWI due to this made-up treatment.

pre.ci <- as.Date(c("2012-05-01", "2013-11-30"))
post.ci <- as.Date(c("2013-12-01", "2014-05-01"))

impact <- CausalImpact(dwi.ci, pre.ci, post.ci,  model.args = list( prior.level.sd = 0.1, nseasons = 12, season.duration = 1))

plot(impact)

impact

Posterior inference {CausalImpact}

                         Average        Cumulative  
Actual                   504            3022        
Prediction (s.d.)        516 (14)       3097 (82)   
95% CI                   [488, 544]     [2929, 3262]
                                                    
Absolute effect (s.d.)   -12 (14)       -75 (82)    
95% CI                   [-40, 16]      [-240, 93]  
                                                    
Relative effect (s.d.)   -2.4% (2.7%)   -2.4% (2.7%)
95% CI                   [-7.8%, 3%]    [-7.8%, 3%] 

Posterior tail-area probability p:   0.18
Posterior prob. of a causal effect:  82%

For more details, type: summary(impact, "report")

Alright, a 2.4% drop (80% sure, it’ll tell you if you check, which isn’t the best). That’s not a lot, but it’s more than the 0% we’d expect. We may not be able to use these results as they are, at least with TABC revenues as a control. The upward trend in revenues through our treatment period and the downward trend in DWI may be distorting the model and may always tell us there’s a drop in DWI arrests.

I’ll note here that other randomly-picked dates for the fake treatment might show different results - we’ve compensated for the higher-than-normal variance in our control data, but if you land on a peak CausalImpact can tell you there’s more of a drop - up to the scale of the real treatment, but with less probability of a causal effect - or no drop at all, or even a slight rise, all with low posterior probabilities. If you’re inclined, I encourage you to play around with it. I’ve picked this one because it’s close to a typical value for the pre-treated range - but a lack of a consistent result from this test should make us feel a little weird about how well suited our data is to this model.

2.3.3 Narcotics Arrests as a Control

Since we were worried TABC revenues are a bad control variable, let’s see if narcotics arrests are any better.

dwi.ci.n <- filter(select(dwi,date,count,narc),date >= '2012-05-01')

pre.ci <- as.Date(c("2012-05-01", "2014-05-31"))
post.ci <- as.Date(c("2014-06-01", "2016-05-09"))

impact.n <- CausalImpact(dwi.ci.n, pre.ci, post.ci,model.args = list(prior.level.sd = 0.1,nseasons = 12, season.duration = 1))
plot(impact.n)

impact.n

Posterior inference {CausalImpact}

                         Average         Cumulative     
Actual                   463             11114          
Prediction (s.d.)        498 (14)        11964 (334)    
95% CI                   [471, 527]      [11314, 12644] 
                                                        
Absolute effect (s.d.)   -35 (14)        -850 (334)     
95% CI                   [-64, -8.3]     [-1530, -200.3]
                                                        
Relative effect (s.d.)   -7.1% (2.8%)    -7.1% (2.8%)   
95% CI                   [-13%, -1.7%]   [-13%, -1.7%]  

Posterior tail-area probability p:   0.005
Posterior prob. of a causal effect:  99.5%

For more details, type: summary(impact, "report")

Using narcotics arrests as a control, ridesharing can account for a 7% drop in DWI arrests, with about the same certainty. Alright, that’s almost the difference we saw earlier - let’s double check on a fake treatment date:

pre.ci <- as.Date(c("2012-05-01", "2013-11-30"))
post.ci <- as.Date(c("2013-12-01", "2014-05-01"))

impact.n <- CausalImpact(dwi.ci.n, pre.ci, post.ci,model.args = list(prior.level.sd = 0.1, nseasons = 12, season.duration = 1))
plot(impact.n)

impact.n

Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   504             3022         
Prediction (s.d.)        520 (14)        3118 (85)    
95% CI                   [493, 548]      [2956, 3286] 
                                                      
Absolute effect (s.d.)   -16 (14)        -96 (85)     
95% CI                   [-44, 11]       [-264, 66]   
                                                      
Relative effect (s.d.)   -3.1% (2.7%)    -3.1% (2.7%) 
95% CI                   [-8.5%, 2.1%]   [-8.5%, 2.1%]

Posterior tail-area probability p:   0.136
Posterior prob. of a causal effect:  86%

For more details, type: summary(impact, "report")

Less impact (-3%). 87% sure though - which most statisticians will throw out. This control looks like it’s better, but I’m not super into it. What else could we do?

2.3.4 Both

dwi.ci.b <- filter(select(dwi,date,count,rev,narc),date >= '2012-05-01')

pre.ci <- as.Date(c("2012-05-01", "2014-05-31"))
post.ci <- as.Date(c("2014-06-01", "2016-05-09"))

impact.b <- CausalImpact(dwi.ci.b, pre.ci, post.ci,model.args = list(prior.level.sd = 0.1,nseasons = 12, season.duration = 1))
plot(impact.b)

impact.b

Posterior inference {CausalImpact}

                         Average        Cumulative    
Actual                   463            11114         
Prediction (s.d.)        504 (16)       12103 (377)   
95% CI                   [473, 535]     [11361, 12833]
                                                      
Absolute effect (s.d.)   -41 (16)       -989 (377)    
95% CI                   [-72, -10]     [-1719, -247] 
                                                      
Relative effect (s.d.)   -8.2% (3.1%)   -8.2% (3.1%)  
95% CI                   [-14%, -2%]    [-14%, -2%]   

Posterior tail-area probability p:   0.005
Posterior prob. of a causal effect:  99.5%

For more details, type: summary(impact, "report")

A 8.2% effect. Some validation:

pre.ci <- as.Date(c("2012-05-01", "2013-11-30"))
post.ci <- as.Date(c("2013-12-01", "2014-05-01"))

impact.b <- CausalImpact(dwi.ci.b, pre.ci, post.ci,model.args = list(prior.level.sd = 0.1,nseasons = 12, season.duration = 1))
plot(impact.b)

impact.b

Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   504             3022         
Prediction (s.d.)        521 (15)        3126 (87)    
95% CI                   [491, 549]      [2948, 3295] 
                                                      
Absolute effect (s.d.)   -17 (15)        -104 (87)    
95% CI                   [-46, 12]       [-273, 74]   
                                                      
Relative effect (s.d.)   -3.3% (2.8%)    -3.3% (2.8%) 
95% CI                   [-8.7%, 2.4%]   [-8.7%, 2.4%]

Posterior tail-area probability p:   0.112
Posterior prob. of a causal effect:  89%

For more details, type: summary(impact, "report")

A 3.3% effect at 88%. We’re getting semi-workable results - the real inflection date with controls just about matches our exponential smoothing models, and we can sort of discount the fake inflection we’re testing with, but not very confidently. We could mess around with a custom multi-seasonal BSTS model and feed it to CausalImpact, I’ll call that an area for future testing.

We’ll have to revisit this model if we get more data.

3 Conclusions

Our conclusions were a bit mixed, but promising. Investigating with RDD suggested that there might be an intervention effect due to Uber and Lyft ceasing ridesharing in Austin, and exponential smoothing techniques suggested measurable drops in DWI counts and DWI counts per MM bar revenue due to the introduction of ridesharing (as well as a small uptick in bar revenue), as well as a subsequent increase in DWI counts and DWI per MM after Uber and Lyft stopped offering rides, once you took into account broad trends and periodic cycles in both data sets.

CausalImpact analysis didn’t work out to be something we could use - both TABC revenue and narcotics arrests data weren’t suitable controls. They did sort of corroborate the results from the other models, though.

Taking into account our results, the models suggest attributing about a 8% drop in DWI arrests in Austin from June 2014 - May 2016 to ridesharing, and about a 5% increase in DWI arrest from May 2016 - December 2016 to an abrupt cessation of Uber and Lyft’s ridesharing services. An 8% drop in DWI over the two years Uber and Lyft were offering ridesharing services accounts for nearly 1000 DWI’s.

Similarly, our new DWI/TABC Mixed Beverage Revenue (Millions) metric is decreased by about 17% by ridesharing, and increased by about 14% without it.

As a little bonus, modelling just bar sales showed that Austin bar revenues may have been improved by as much as $68,000,000 due to the availablity ridesharing over two years, adding to the TABC’s tax coffers by nearly $9,500,000 over that time. Without Uber and Lyft, bars may have seen as much as a $34,000,000 depression in sales over 7 months, shorting the TABC $4,700,000 in taxes. We didn’t study this as much as DWI’s, though, and unobserved interventions could account for some of these amounts without changing our conclusions about DWI counts directly.

4 Future Work and Caveats

With more fine-grained time series data, or some certified-extra-accurate data from the APD, I think we could improve on this analysis. My original thoughts were to see if ridesharing had a more specific effect on certain kinds of DWI’s - say specifically 1st offenses - and in certain places, but the data was not complete enough for that. Additionally, Austin’s crime data in general may be unreliable (that’s certainly a caveat - see http://www.mystatesman.com/news/local/why-the-austin-police-department-dwi-crash-statistics-keep-changing/GNv8LTrbEgTd7KEObnPwgL/), although we did compare two sources for crime data. Updated data from after the treatment period could alter the models’ results substantially. I would love to get something like number of rideshare rides taken per month, too.

If you’ve got some comments or suggestions regarding this work, contact me! Pull the data and models down and mix them up yourself if you’d like, I’d love to see any new results or conflicting conclusions.

5 Appendix

5.1 DWI Arrests as a Proxy for Drunk Driving

We could be concerned that, rather than being a reasonable proxy for drunk driving, the number of DWI arrests could be constrained by the resources available to the APD - what if drunk driving is out of control, and the APD has hit some sort of arrest ceiling? If this was the case, we would expect to see an unchanging average amount of monthly arrests with some noise, possibly varying with the size of the police force at their maximum capacity for making DWI arrests.

After some quick Googling, we have a few sworn-officer counts for the APD: about 1200 in 2002, 1600 in 2010, and 1800 in 2016, so we can say that the police force is growing. Reading up on the APD also indicates that the number of officers familiar with handling a DWI arrest increased over that time frame as well. Let’s remind ourselves what the DWI incident count looks like:

ggplot(data = dwi, aes(x=date,y=count)) + geom_point() + geom_smooth(method = 'loess')

Even as the APD grows, the amount of arrests varies wildly (a range of nearly 200 arrests per month), in addition to holding several increasing and decreasing trend periods we already studied. From this, I think it’s safe to say that the number of arrests does track drunk driving without saturation.

5.2 Crash/Leaving The Scene

We can assume, I think, that some nonzero percentage of hit-and-run incidents, especially those happening during peak drunk driving hours, involve a drunk driver who was not caught and charged with a DWI. If this number is high enough, and/or fluctuating enough, it could cast doubt on the validity of DWI arrests as a proxy for drunk driving by suggesting that the arrest count is meaninglessly small and invariant with respect to some actual massive population of drunk drivers.

Additionally, if we can confidently extract a portion of these incidents and label them as probable drunk-driving incidents, and this portion is significantly large and not in line with larger trends across our treatment area, it may change the output of our models if we include it. This second one is a stretch and relies on another big assumption (that we can identify a drunk driving hit-and-run), but let’s take a crack at it.

har.raw <- read.csv('data/hitandrun-raw.csv')
har.raw$hour <- substr(har.raw$date,12,13)
#har.raw$longdate <- parse_date_time(har.raw$date,'%Y-%m-%d %H:%M:%S',tz ='America/Chicago')
har.raw$date <- as.Date(har.raw$date)

#qplot(hour(longdate), data=har.raw, geom="histogram", bins = 24)
#qplot(wday(longdate), data=har.raw, geom="histogram", bins = 24)

har.raw$year <- year(har.raw$date)
har.raw$month <- month(har.raw$date)

har.monthly <- har.raw %>% group_by(year,month)
har <- summarize(har.monthly,count=n())
har$date <- as.Date(paste0(har$year,'-',har$month,'-01'))

har <- filter(har,count > 50) #remove a few months outside of the 1-6/2015 and 1-7/2016 range without complete data
har <- har[-63,] 
har <- har[-87,] 

qplot(x=date,y=count,data = har)

There are a lot of incidents here, but not enough to totally swamp the data, so our first concern is ruled out, I think. I’m a little concerned about the jump in hit and runs in the later part of the data though, even though it’s noisy and the incident report is known to disagree with the chief’s report towards the end of the time frame - let’s take a stab at picking the incidents most likely to be alcohol related, by filtering just the ones that occur when DWI’s are most likely to be reported.

har.raw <- filter(filter(har.raw, hour > 19),hour < 6) #only hit-and-runs during prime drinking and driving hours

har.monthly <- har.raw %>% group_by(year,month)
har <- summarize(har.monthly,count=n())
har$date <- as.Date(paste0(har$year,'-',har$month,'-01'))

har <- filter(har,count > 50) #remove a few months outside of the 1-6/2015 and 1-7/2016 range without
har <- har[-63,] 
har <- har[-87,] 

qplot(x=date,y=count,data = har)

#ggplot(data = filter(dwi.monthly,date < '2015-01-01'), aes(factor(month(date)),weight = count) ) + geom_bar()

ggplot(data = filter(har,date < '2015-01-01'), aes(factor(month(date)),weight = count) ) + geom_bar()

That brings the mean down considerably and levels out the trend - although the histogram doesn’t match our incident report DWI histogram.

If we were to further sample this data for likely drunk driving hit-and-runs, we’d have to weight it to match the DWI histogram, and then apply some flat ‘this many hit-and-runs in prime hit-and-run timeframes involve alcohol’ percentage. I don’t have a number for that, we’ll just call it 20%.

filter<-data.frame(month=c(1,2,3,4,5,6,7,8,9,10,11,12),trim=c(1,.70,.7,.60,.64,.56,.58,.62,.76,.72,.73,.81))

har <- merge(har,filter,id ='month')
har$count <- har$count * har$trim
ggplot(data = filter(har,date < '2015-01-01'), aes(factor(month(date)),weight = count) ) + geom_bar()

har$trim <- NULL

har$count = har$count * .20

qplot(x=date,y=count,data = har)

And it’s less than or about the same as the difference between the APD Incident Report and the Chief’s Monthly Report. Maybe a more official estimate of hit-and-run drunk driving may have made it into the Chief’s report, accounting for the discrepancy.

For the range 2008-2012, the difference is small, in line with the prevailing decreasing trend, and as it contributes only to seasonal models I’m comfortable discounting it entirely. For the months in 2015 and 2016, I’m not sure what to do with it - I don’t think it’s a good idea to add it to the Chief’s report, and if it only makes up the difference between the APD Incident Report and the Chief’s report, and since I don’t even have it for the rest of the time frame - I think it’s safe to ignore it and if more accurate data presents itself, test our conclusions again.

5.3 Downtown and All Other Arrest Data

Late in the development of this writeup I noticed that the Chief’s Monthly Report was broken into broad geographic sections, and also considered that the ‘all other Part 2 crimes’ count available on the report would be a better control than narcotics arrests, as it might have a lower variance, and be a better representation of ‘background’ crime. I’ll quickly look at some of this data.

bonus <- read.csv('data/other.csv')
head(bonus)

bonus$date <- as.Date(bonus$date)
bonus$dt.count <- as.numeric(as.character(bonus$dt.count)) #cast strictly
bonus$count <- as.numeric(as.character(bonus$count)) 
bonus$other <- as.numeric(as.character(bonus$other)) 

ggplot(data = bonus, aes(factor(month(date)),weight = dt.count) ) + geom_bar()

ggplot(data = bonus, aes(y=dt.count,x=date)) + geom_point()

ggplot(data = bonus, aes(y=dt.count/count,x=date)) + geom_point()

bonus.pre <- filter(bonus,date < '2014-05-31')
bonus.post <- filter(bonus,date > '2014-05-31')

bonus.pre.ts <- ts(bonus.pre$dt.count, start = c(2010, 1), frequency = 12)
bonus.post.ts <- ts(bonus.post$dt.count, start = c(2014, 6), frequency = 12)

bonus.hw <- HoltWinters(bonus.pre.ts)
forecast <- predict(bonus.hw, n.ahead = 24, prediction.interval = T, level = 0.9)

plot(bonus.hw,forecast); lines(bonus.post.ts,col='green')

effect <- forecast - bonus.post.ts

#pre.residual <- as.numeric((bonus.hw$fitted - bonus.pre.ts))
#post.residual <- as.numeric(effect)
#t.test(pre.residual,post.residual) #e-16

pct.effect <- cumsum(effect)/cumsum(forecast)
pct.effect[24]

[1] 0.1412853

Cool. Some seasonality, a plot that looks like a scaled version of the full DWI plot, and a sense that downtown DWI’s generally account for 10-14% of all DWI’s in Austin.

A quick glance at just the downtown data suggests that ridesharing had the effect of lowering DWI counts downtown by 14% over 24 months (9% over 12.)

ggplot(data = bonus, aes(y=other,x=date)) + geom_point() + geom_smooth()

`geom_smooth()` using method = 'loess'

ggplot(data = bonus, aes(factor(month(date)),weight = other) ) + geom_bar()

bonus.b <- filter(select(bonus,date,count,other),date >= '2012-05-01')

pre.ci <- as.Date(c("2012-05-01", "2014-05-31"))
post.ci <- as.Date(c("2014-06-01", "2016-05-09"))

impact.b <- CausalImpact(bonus.b, pre.ci, post.ci,model.args = list(prior.level.sd = 0.1,nseasons = 12, season.duration = 1))
plot(impact.b)

impact.b

Posterior inference {CausalImpact}

                         Average         Cumulative     
Actual                   463             11114          
Prediction (s.d.)        496 (13)        11900 (316)    
95% CI                   [471, 522]      [11295, 12534] 
                                                        
Absolute effect (s.d.)   -33 (13)        -786 (316)     
95% CI                   [-59, -7.5]     [-1420, -180.9]
                                                        
Relative effect (s.d.)   -6.6% (2.7%)    -6.6% (2.7%)   
95% CI                   [-12%, -1.5%]   [-12%, -1.5%]  

Posterior tail-area probability p:   0.005
Posterior prob. of a causal effect:  99.5%

For more details, type: summary(impact, "report")

pre.ci <- as.Date(c("2012-05-01", "2013-11-30"))
post.ci <- as.Date(c("2013-12-01", "2014-05-01"))

impact.b <- CausalImpact(bonus.b, pre.ci, post.ci,model.args = list(prior.level.sd = 0.1,nseasons = 12, season.duration = 1))
plot(impact.b)

impact.b

Posterior inference {CausalImpact}

                         Average         Cumulative   
Actual                   504             3022         
Prediction (s.d.)        491 (17)        2943 (101)   
95% CI                   [456, 522]      [2737, 3130] 
                                                      
Absolute effect (s.d.)   13 (17)         79 (101)     
95% CI                   [-18, 47]       [-108, 285]  
                                                      
Relative effect (s.d.)   2.7% (3.4%)     2.7% (3.4%)  
95% CI                   [-3.7%, 9.7%]   [-3.7%, 9.7%]

Posterior tail-area probability p:   0.205
Posterior prob. of a causal effect:  80%

For more details, type: summary(impact, "report")

Other Part 2 crimes seem to be on the decline, although there looks like there might be an inflection in late 2013 - early 2014. We see a little seasonality, possibly consistent with studies that have shown a link between warm weather and crime. CausalImpact tells us that we see about a 6% drop in DWI, and it validates much better than our other controls for our selected fake intervention date, but not much better for other dates. I’ll leave it here for future research.

6 References

Deterring Drunk Driving Fatalities: An Economics of Crime Perspective Bruce L. Benson and David W. Rasmussen

Inferring Causal Impact Using Bayesian Structural Time-series Models Kay H. Brodersen, Fabian Gallusser, Jim Koehler, Nicolas Remy, Steven L. Scott

Regression Discontinuity Designs in Economics David S. Lee and Thomas Lemieux

Program Description: APD DWI Enforcement https://one.nhtsa.gov/people/injury/research/AustinPoliceDWI/pages/ProgDesc.html

7 Acknowledgements

This work would have been much worse without the intensely appreciated input of its peer reviewers, who know who they are :)

8 Disclaimers

8.1 APD

The Austin Police Department does not assume any liability for any decision made or action taken or not taken by the recipient in reliance upon any information or data provided.

8.2 TABC

The Comptroller of Public Accounts cannot vouch for the data or analysis derived from data after it has been retrieved from the Comptroller’s Web site.

Some numbers in the R may vary slightly from the numbers in the text, as some models are rolled at render-time and may involve some randomness.

9 License

The MIT License (MIT)

Copyright (c) 2017 Ian Wells

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.