]> code.communitydata.science - stats_class_2020.git/blob - psets/pset3-worked_solution.rmd
e4b6d32985e421c6b90e7d5890953ffc4dddd8ce
[stats_class_2020.git] / psets / pset3-worked_solution.rmd
1 ---
2 title: "Problem set 3: Worked solutions"
3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
4 author: "Aaron Shaw"
5 date: "October 12, 2020"
6 output:
7   html_document:
8     toc: yes
9     toc_depth: 3
10     toc_float:
11       collapsed: true
12       smooth_scroll: true
13     theme: readable
14   pdf_document:
15     toc: yes
16     toc_depth: '3'
17 ---
18
19 ```{r setup, include=FALSE}
20 knitr::opts_chunk$set(echo = TRUE, message=FALSE, tidy='styler')
21 ```
22
23 # Programming challenges  
24 ## PC1. Learn about the data.
25
26 The documentation makes it sound like the Illinois data is particularly messy/incomplete. It explains that the data is limited to 2012-2017 (instead of extending back to 2004) "due to format issues and relevance." The notes point out that data related to searches vary substantially year-to-year, suggesting that the measurement was inconsistent and that caution should be exercised in the analysis of this measure in particular. The `search_conducted` variable seems to include searches of any kind (vehicles, driver, passengers). I also notice that the codebook does not provide much insight into *how* `subject_race` and `subject_sex` information were collected, measured, and/or recoded. I think it's best to presume that the records reflect the judgments/notes/documents accessed by the officer(s) and/or police forces involved rather than the perspectives of the subjects of the traffic stops. There are likely some inconsistencies in these measures as well.
27
28 ## PC2. Import, explore, clean
29
30 I'll use the tidyverse `read_csv()` function to do this. While I'm at it, I'll import the `lubridate` library because I already know this is time series data and there are a bunch of helpful functions in `lubridate` that I will probably want to use.
31
32 ```{r import}
33 library(tidyverse)
34 library(lubridate)
35 data.dir <- "https://communitydata.science/~ads/teaching/2020/stats/data/week_05"
36 filename <- "il_statewide_sample-2020_04_01.csv"
37 ilstops <- read_csv(url(paste(data.dir, filename, sep="/")))
38
39 ilstops
40 ```
41
42 ### Recode and explore key variables
43
44 Right off the bat, I see a lot of `NA` values and some things likely to need recoding. Let's look a little more closely at the variables we are going to work with.
45 ```{r class check}
46
47 key_variables <- c("date", "vehicle_year", "subject_race", "subject_sex", "search_conducted")
48
49 lapply(ilstops[key_variables], class)
50 ```
51 Indeed, I want to create some factors for the `subject_race` and `subject_sex` variables. 
52 ```{r recode factors}
53 ilstops$subject_race <- factor(ilstops$subject_race)
54 ilstops$subject_sex <- factor(ilstops$subject_sex)
55 ```
56 ## PC3 Summarize key variables
57 Now, let's look at summaries of the key variables. Keep in mind that these are all variables calculated over all traffic stops in the dataset.  
58 ```{r summarize}
59 lapply(ilstops[key_variables], summary)
60 ```
61 Note that the number of missing (`NA`) values is quite different depending on the variable.
62
63 ### Recode nonsense `vehicle_year` values
64
65 Review those summaries and you might notice something a little funny about the `vehicle_year` variable. Doesn't this data end in 2017? None of the `vehicle_year` values should be more recent than 2018 (if we're feeling generous since car manufacturers sometimes start selling "next year's car" at the end of a given year?). I'm also modestly skeptical that the police pulled over a car from 1900, but at least that's a plausible value! In the abence of additional information about the provenance of this variable or reasons why values larger than 2018 might be reasonable, I will go ahead and recode anything impossibly new as missing (`NA`).
66 ```{r recode vehicle_year}
67 ilstops$vehicle_year[ilstops$vehicle_year > 2018] <- NA
68
69 summary(ilstops$vehicle_year)
70 hist(ilstops$vehicle_year)
71 ```
72
73 ## PC4. Summarize bivariate relationships  
74 For bivariate summaries and comparisons, I'll walk through things variable by variable. In all of these cases except for the `date` variable, I'm considering the bivariate relationships in aggregate (rather than taking the "over time" character of the data into account). This is good to keep in mind since it might mean that the summaries mask temporal trends.
75
76 ### Searches by `vehicle year`  
77 I'll start with a continuous predictor, `vehicle_year`. Since the outcome variable is a dichotomous category, I can look at this by calculating summary statistics within the two outcomes or create a pair of plots. I love visual summaries, so I'll start with side-by-side boxplots:
78
79 ```{r pair boxplot vehicle year x searches}
80 ilstops %>%
81   filter(!is.na(vehicle_year) & !is.na(search_conducted)) %>%
82   ggplot(aes(x=search_conducted, y=vehicle_year)) +  geom_boxplot()
83
84 ```
85
86 It looks like stops of older vehicles are, on average, a little bit more likely to result in a search, but stops involving the very oldest vehicles tend to result in fewer searches.
87
88 With that in mind, I'll calculate summary statistics across the two groups. One way to do this is with a call to the `tapply` function nested within a call to the `with()` command. I don't think we've used `with` yet, but it can allow me to identify the object I want to use (the first argument to `with()`) to run some command on (the second argument, in this case the call to `tapply()`).
89
90 ```{r crosstab vehicle year x searches}
91 with(ilstops, 
92      tapply(vehicle_year, search_conducted, summary)
93      )
94 ```
95 Indeed, as I suspected based on the histograms, the center of the distribution for stops resulting in searches is a little lower (older vehicles) than among stops not resulting in searches. I'll calculate the standard deviations with a separate command:
96 ```{r}
97 with(ilstops, tapply(vehicle_year, search_conducted, sd, na.rm=TRUE))
98 ```
99
100 The stops not resulting in searches (`FALSE`) have a larger/wider spread of vehicle ages than those that do result in searches (`TRUE`). This difference is likely driven by the larger number of old vehicle outliers inovlved in stops that do not result in searches. This is a nice example of how the visual and numerical summaries might help convey a somewhat subtle feature of my data.
101
102 ### Searches by `subject_sex`
103
104 Next, we'll focus on the `subject_sex` variable. Both the outcome and this predictor variable are categorical measures that take only two (non-missing) values. That means that I can focus on creating contingency tables to report the raw numbers of traffic stops in each of the four possible combinations of the two variables. I've provided a couple of methods here:
105 ```{r crosstab searches and sex}
106 table(ilstops$subject_sex, ilstops$search_conducted)
107
108 ## the xtabs() command produces nicely formatted output that includes both variable names
109 xtabs(~ subject_sex + search_conducted, ilstops)
110
111 ```
112 One way to get the proportions in Base-R is to call the `proportions()` function on a contingency table. 
113 ```{r props searches}
114 proportions(
115   xtabs(~ subject_sex + search_conducted, ilstops)
116   )
117 ```
118
119 That's a lot of significant digits...especially for proportions. I can clean it up a bit by nesting everything within a call to the `round()` function and restricting the output to just two significant digits.
120
121 ```{r pretty prop searches}
122 round(
123   proportions(
124     xtabs(~ subject_sex + search_conducted, ilstops)
125     )
126   , digits=2)
127 ```
128 Much better! However, notice that these proportions are all calculated against the *total* number of non-missing values. I can specify row or column proportions with an additional argument to the `proportions` function to specify which margins to use. Following R conventions, `margin=1` calculates proportions row-wise and `margin=2` calculates them column-wise: 
129
130 ```{r marginal proportions}
131
132 ## Row-wise proportions
133 round(
134   proportions(
135     xtabs(~ subject_sex + search_conducted, ilstops),1
136     )
137   , digits=2)
138
139 ## Column-wise proportions
140 round(
141   proportions(
142     xtabs(~ subject_sex + search_conducted, ilstops),2
143     )
144   , digits=2)
145 ```
146 These different marginal proportions support distinct statements about the data. Here is an example for each:   
147
148 * **Within the rows (`subject_sex`):** *the proportion of searches in stops of male drivers is twice as large as the proportion in stops of female drivers.*  
149 * **Within the columns (`search_conducted`):** *stops that resulted in searches involved male drivers almost $75\%$ of the time.*  
150
151 ### Searches by `subject_race`
152
153 The `subject_race` variable is also categorical, but takes one of five non-missing values in this data. As a result, the basic structure/approach for my summaries is almost identical to that pursued with `subject_sex`.
154
155 First, tabular summaries of raw counts and proportions:
156 ```{r}
157 xtabs(~ subject_race + search_conducted, ilstops)
158 ```
159 In this case, I think the marginal proportions are far more interesting than the aggregate proportions. Here are proportions of searches within each category of `subject_race` (row-wise given the way I've structured the call to `xtabs`):
160 ```{r prop searches within group}
161 ## Proportions of searches among stops within groups
162 round(
163   proportions(
164     xtabs(~ subject_race + search_conducted, ilstops), 1)
165   , digits=2)
166 ```
167 Notice that the proportion of black and hispanic drivers with stops resulting in searches is higher than within other groups.
168
169 Now, I'll calculate the proportions of each `subject_race` category among stops resulting in searches and stops that do not result in searches:
170 ```{r prop group within searches}
171 ## Proportions of groups within those stops (not) resulting in searches
172 round(
173   proportions(
174     xtabs(~ subject_race + search_conducted, ilstops), 2)
175   , digits=2)
176 ```
177 Here the noteworthy comparisons again arise within the black and hispanic categories. Both groups account for a substantially larger proportion of stops resulting in searches vs. those that do not result in searches. 
178
179 For the sake of completeness/comparison, here's a way to do similar cross-tabulations in single a chunk of tidyverse code. I include conditional proportions of all stops to facilitate comparison with some of the tables I created earlier as well:  
180 ```{r tidyverse crosstabs}
181 ilstops %>%
182   group_by(subject_race) %>%
183   filter(!is.na(subject_race)) %>%
184   summarize(
185     n_stops = n(),
186     prop_total_stops = round(n() / nrow(ilstops), digits=3),
187     across(
188       search_conducted,
189       list(
190         sum = ~ sum(.x, na.rm=TRUE),
191         over_n_stops = ~ round(mean(.x, na.rm=TRUE), digits=3)  
192         )
193       )
194     ) %>%
195       arrange(desc(n_stops))
196 ```
197
198 Notice the use of `across()` within a call to `summarize()` provides one way to calculate conditional summariy info. I also use a nested call to `arrange()` and `desc()` at the end to sort my results in descending order by one of the columns.  
199
200 ### Searches by `date`
201
202 Last, but not least, I can inspect how the number of searches is distributed over time. I can do this a few ways, but once again will start with a visual summary because that is how I prefer to approach things. Remember that `search_conducted` is a dichotomous categorical measure, so a boxplot is a good starting point.
203 ```{r searches over time}
204 ilstops %>%
205   ggplot(aes(x=search_conducted, y=date)) + geom_boxplot()
206   
207 ```
208
209 How considerate of ggplot to include the `NA` category by default! Notice how seemingly all of the missing values come from 2012-2014? There's some of that missing data the SOPP folks warned us about in their documentation. Let's inspect things a bit more closely:
210
211 ```{r explore missing values in searches}
212 ilstops %>%
213   group_by(search_conducted) %>%
214   summarize(
215     min = min(date, na.rm=TRUE),
216     max = max(date, na.rm=TRUE)
217   )
218 ```
219
220 All the missing data for `search_conducted` comes within a two-year range. 
221
222 I can use `filter` to look at the `search_conducted` variable within that two year range alone:
223 ```{r missing search values pre 2014}
224 ilstops %>%
225   filter(date < as_date('2014-01-11')) %>%
226   group_by(search_conducted) %>%
227   summarize(
228     n = n()
229   )
230 ```
231 *Most* of the `search_conducted` data between 2012-2014 isn't missing, even if that's the only period within which missing values for `search_conducted` appear. There's not much we can do about this, but it's worth taking into consideration as we go.
232
233
234 ## PC5 Searches within `subject_race` groups over time
235
236 The goal here is to synthesize pieces of the work we did in PC3 and PC4 to analyze subgroups of the stops and the stops resulting in searches over time. I'll start by reorganizing/tidying the data to support the creation of my time series plots.
237
238 ### 5.1 Binning data within time periods
239
240 I want to plot this data as a time series of some numeric measure, but right now my outcome variable is a categorical indicator. How can I plot TRUE/FALSE values over time? One good way is to bin my data into meaningful time-periods (e.g., months or weeks) and count the number of TRUE values (and/or calculate proportions of TRUE values) within each bin. 
241
242
243 Let's do that using some of the handy features of `lubridate` for working with date objects. I have demonstrated a few approaches to binning and grouping in the R tutorial materials, but here's yet another that takes advantage of some `lubridate` features (specifically the call to `format` within the `mutate` step below) quite effectively:
244
245 ```{r bin to months}
246 ## library(lubridate)
247
248 ilstops_monthly <- ilstops %>%
249   filter(!is.na(subject_race)) %>%
250   mutate(
251     yearmonth = format(date, '%Y-%m')
252   ) %>%
253   group_by(yearmonth) %>%
254   summarize(
255     stops = n(),
256     searches = sum(as.numeric(search_conducted), na.rm=T),
257     prop_searched = round(searches / stops,3)
258     )
259
260 ilstops_monthly
261 ```
262 Looks great, except do you notice how my new `yearmonth` variable is no longer recognized as a date or date-time object? I think I can fix that with the handy `parse_date_time` function from `lubridate`:
263 ```{r clean binned dates}
264 ilstops_monthly$date <- parse_date_time(ilstops_monthly$yearmonth, "ym")
265 class(ilstops_monthly$date)
266 ```
267
268 Now I can plot it since `ggplot` knows what to do with objects in the `POSIXct` class. I'm also going to add a type of layer we haven't worked with yet (`geom_smooth`) which incorporates a smoothed conditional mean to the plot (appropriately enough in this dataset, in the form of a thin blue line).
269
270 ```{r plot monthly stops}
271 ilstops_monthly %>%
272   ggplot(aes(x=date, y=stops)) + geom_line() + geom_smooth() + labs(title="Illinois traffic stops", y="stops (per month)")
273          
274 ```
275
276 Turns out traffic stops are not distributed smoothly across the months of the year (I notice, in particular, the big dips and spikes around the end of each year). It also looks like there is a slight $\cup$-shaped trend overall, but the shifts do not look very large relative to either the month-to-month variations or the total number of stops per year.
277
278 Here's a similar plot for searches:
279 ```{r plot monthly searches}
280 ilstops_monthly %>%
281   ggplot(aes(x=date, y=searches)) + geom_line() + geom_smooth() + labs(title="Illinois traffic stops that result in searches",
282                                                                        y="searches (per month)")
283          
284 ```
285 Turns out searches are also not distributed smoothly across the months of the year (again, the big dips and spikes around the end of the year). In this case, there seems to be a very slight downward trend on average, but this decrease mostly appears in the first couple of years of the data and does not look very large relative to the month-to-month variations or the total number of searches per year (maybe ten fewer searches per month on average during the latter half of the dataset?).
286
287 Let's try looking at the proportion of stops that result in searches.
288
289 ```{r plot monthly prop searched}
290 ilstops_monthly %>%
291   ggplot(aes(x=date, y=prop_searched)) + geom_line() + geom_smooth() + labs(title="Proportion of Illinois traffic stops that result in searches",  y="proportion searched\n(per month)")
292 ```
293 Still quite spiky! Interestingly, we see an inverse-$\cap$ trend in the average here (contrast with the raw counts of traffic stops above), but these fluctuations are, again, very small relative to the underlying variation or the range of measure. Also, notice the narrow y-axis scale sort of hides the fact that almost all the values fall between $4-6\%$. In other words, all of the variation here is within quite a narrow range.
294
295 ### 5.2 Plot stops within `subject_race` categories
296
297 To do this, I'll need to go back to the code that created my `ilstops_monthly` tibble above and incorporate the conditional summary data. To do this, I'll "explode" the aggregated summaries from the old tibble by inserting an extra variable to `group_by()` before making my call to `summarize`. Notice that I go ahead and cleanup the dates again too (since I know to anticipate that issue)
298
299 ```{r bin to subgroup months}
300 ilstops_monthly_grouped <- ilstops %>%
301   filter(!is.na(subject_race)) %>%
302   mutate(
303     yearmonth = format(date, '%Y-%m')
304   ) %>%
305   group_by(yearmonth, subject_race) %>%
306   summarize( 
307     stops = n(),
308     searches = sum(as.numeric(search_conducted), na.rm=T),
309     prop_searched = round(searches / stops,3)
310     )
311
312 ilstops_monthly_grouped$date <- parse_date_time(ilstops_monthly_grouped$yearmonth, "ym")
313
314 ilstops_monthly_grouped
315
316 ```
317 Great! Now, I need to add the values for total stops and total searches each month (note that this will be the same for all five groups/rows within any given value of `date` or `yearmonth`) so that I can calculate the proportion of total stops and total searches accounted for by each group. There are a bunch of ways we could go about this. Since I've made it this far with tidyverse code, I'll keep going in that direction and use the `mutate` verb to keep all my existing rows while also creating some new ones. Notice that I have to do this twice because my first call to `group_by` is only by `yearmonth` so that I can create the monthly totals, whereas my second call to `group_by` once again looks wihtin `yearmonth` and `subject_race` to calculate the proportions. In order to make sure that nothing weird happens along the way I insert a call to `ungroup` as well (which does exactly what you might hope following an operation performed on some data that's already grouped in some way).
318
319 ```{r create total stops, total searches and prop for each within groups}
320 ilstops_monthly_grouped <- ilstops_monthly_grouped %>%
321   group_by(yearmonth) %>%
322   mutate(
323     total_stops = sum(stops),
324     total_searches = sum(searches)
325   ) %>%
326   ungroup() %>%
327   group_by(yearmonth, subject_race) %>%
328   mutate(
329     prop_total_stops = stops / total_stops,
330     prop_total_searches = searches / total_searches
331   )
332
333 ilstops_monthly_grouped
334
335 ## I think the html drops a few of our new columns, so I'll inspect them manually here
336 head(ilstops_monthly_grouped[,c("yearmonth", "subject_race", "total_searches", "prop_total_stops", "prop_total_searches")])
337 ```
338
339 Now my plotting code will look pretty similar to the univariate plots above, just with the addition of another aesthetic...
340
341 ```{r plot monthly subgroup stops}
342 ilstops_monthly_grouped %>%
343   ggplot(aes(x=date, y=stops, color=subject_race)) + geom_line() + geom_smooth() + labs(title="Illinois traffic stops",  y="stops (per month)")
344          
345 ```
346
347 ### 5.3 Plot searches within `subject_race` categories
348 ```{r plot monthly subgroup searches}
349 ilstops_monthly_grouped %>%
350   ggplot(aes(x=date, y=searches, color=subject_race)) + geom_line() + geom_smooth() + labs(title="Illinois traffic stops resulting in searches",  y="searches (per month)")
351
352 ```
353
354
355 ### 5.4 Plot proportion of searches within `subject_race` categories
356
357 The question for this one only asked about proportions of total searches accounted for by each group of `subject_race`, when several other interesting baselines/comparisons are available. For example, the proportion of total stops across groups as well as stops resulting in searches *within* each group also strike me as potentially illustrative. I'll plot each of these proportions just for the sake of completeness. Let's start with the proportion of total stops plotted across groups.
358 ```{r plot monthly subgroup prop total stops}
359 ilstops_monthly_grouped %>%
360   ggplot(aes(x=date, y=prop_total_stops, color=subject_race)) + geom_line() + geom_smooth() + labs(title="Proportion of total Illinois traffic stops by race/ethnicity",  y="proportion of total stops\n(per month)\n")
361 ```
362
363 Onwards to the proportion of stops resulting in searches calculated *within* each group of `subject_race`:
364
365 ```{r plot monthly subgroup prop searched}
366 ilstops_monthly_grouped %>%
367   ggplot(aes(x=date, y=prop_searched, color=subject_race)) + geom_line() + geom_smooth() + labs(title="Proportion of Illinois traffic stops that result in searches",  y="proportion searched\n(calculated per month)")
368 ```
369
370 Hmm, that looks pretty messy. The proportions compress everything into a narrow region of the y-axis. Also, the data for some of these groups (especially "other") fluctuates all over the place, which is a distration since that's the smallest category. 
371
372 Instead of plotting all the proportions on the same figure, it might be easier to see/compare variations if I use facets and don't include the "other" category (because the y-axis scale is just so different).
373
374 ```{r plot monthly subgroup prop searched facets}
375 ilstops_monthly_grouped %>%
376   filter(subject_race != "other") %>%
377   ggplot(aes(x=date, y=prop_searched, color=subject_race)) + geom_line() + geom_smooth() + facet_grid(rows=vars(subject_race)) + labs(title="Proportion of Illinois traffic stops that result in searches by race/ethnicity",  y="proportion searched\n(per month)\n")
378 ```
379
380 The title runs a little long and facet label for "asian/pacific islander" gets truncated. If I were preparing this for publication I'd fix those issues, but for now, I'm going to acknowledge them and move on.
381
382 Now, let's look at the proportion of total searches.
383
384 ```{r plot monthly subgroup prop total searches}
385 ilstops_monthly_grouped %>%
386   ggplot(aes(x=date, y=prop_total_searches, color=subject_race)) + geom_line() + geom_smooth() + labs(title="Proportion of total Illinois traffic stops resulting in searches by race/ethnicity",  y="proportion of total searches\n(per month)\n")
387 ```
388
389 ## PC6. Calculate baseline population proportions
390
391 As recommended, I'll calculate this using the [`county_complete`](https://www.openintro.org/data/index.php?data=county_complete) provided as part of the `openintro` package to build the relevant variables for Illinois in 2010. You'll want to review the [list of variables](https://www.openintro.org/data/index.php?data=county_complete) to figure out what you need to use here. 
392
393 Note that what we're looking for are *statewide-proportions* of race/ethnic groups in 2010 within the same categories in the SOPP traffic stop data for Illinois. The `county_complete` dataset does not provide these values directly, but instead has county-level populations and county-level proportions of slightly different race/ethnic groups over multiple different years across all the counties in the United States. That means that we'll need to take some steps (like the following) to generate state-level proportion estimates:  
394
395 * Restrict the dataset to Illinois counties only.  
396 * Multiply county-level populations by county-level proportions for salient groups to create county-level population counts within groups.  
397 * Sum county-level populations within groups to create state-level group populations.  
398 * Divide state-level group populations by total state population (sum of all county populations) to calculate state-level proportions.  
399
400 In the code chunk below, I'll tackle these steps in a very "verbose" and redundant way:
401 ```{r verbose county baselines}
402 library(openintro)
403 data(county_complete)
404
405 il_county <- county_complete[county_complete$state == "Illinois",]
406
407 ## total state population is easy
408 il_pop2010 <- sum(il_county$pop2010)
409
410 ## Turn IL county-level proportions into counts
411 il_county$white2010 <- il_county$pop2010*il_county$white_2010
412 il_county$black2010 <- il_county$pop2010*il_county$black_2010
413 il_county$hispanic2010 <- il_county$pop2010*il_county$hispanic_2010
414
415 ## We'll have to sum these two groups to approximate "Asian / Pacific Islander" from SOPP data
416 il_county$asian2010 <- (il_county$pop2010*il_county$asian_2010)
417 il_county$pi2010 <- (il_county$pop2010*il_county$pac_isl_2010)
418
419 ## Then go through summing and dividing by the total population
420 ##
421 ##prop.white
422 sum(il_county$white2010) / il_pop2010
423 ##prop.black 
424 sum(il_county$black2010) / il_pop2010
425 ##prop.hispanic 
426 sum(il_county$hispanic2010) / il_pop2010
427 ##prop.api
428 (sum(il_county$asian2010, na.rm=TRUE)+ sum(il_county$pi2010, na.rm=TRUE)) / il_pop2010
429 ```
430
431 That is all perfectly accurate; however, the code is verbose and redundant because it calculates each (intermediate and group) value with separate lines of code that are very repetitive. There are many more concise ways to go about this, all of which make it faster, less prone to typing mistakes, and easier to extend. How to proceed? Anytime I find myself repeating the same steps of analysis more than once, I consider creating a function and then calling that function within a call to an `*apply()` function. Here's what that might look like here:
432 ```{r county baselines via function}
433 gen_subgroup_prop <- function(group_var, d=il_county){
434   total_pop2010 <- sum(d["pop2010"]) # state population total
435   county_ests <- d["pop2010"]*d[group_var] # county-level counts within groups
436   group_prop <- sum(county_ests, na.rm=TRUE) / total_pop2010 # state-level proportions within groups
437   return(round(group_prop, 4)) # prettify output
438 }
439
440 ## give it a try:
441 gen_subgroup_prop("hispanic_2010")
442
443 ## Now let's tie everything together
444 groups <- c("white_2010", "black_2010", "hispanic_2010", "asian_2010", "pac_isl_2010")
445
446 sapply(groups, gen_subgroup_prop)
447
448 ```
449
450 I need to add those last two values together to reproduce the "asian/pacific islander" category from the SOPP data. Note also that there are slight differences between my results here and my results above. Any differences **should** be due entirely to rounding and, in a scenario where this analysis was being prepared for publication and/or informing a policy decision, I'd ansp;ite;y want to check that and think carefully how many significant digits to include in my reported estimates. The Illinois population is sufficiently large that even a few hundredths of a percent turn into very meaningful numbers of people!
451
452 # Statistical questions
453
454 Note that my "solutions" below reflect possible interpretations, but are not intended to be exhaustive or exclusive of other possible interpretations. We'll discuss these in our class meeting.
455
456 ## SQ1. Interpret analysis from PCs3-5  
457
458 Some interpretation appeared above alongside specific results. Here is a brief summary of several striking points.  
459
460 * Overall, traffic stops and searches in Illinois between 2012-2018 were distributed unequally across different racial/ethnic groups.  
461     + The proportion of stops resulting in searches is higher among stops of black and hispanic drivers ($7\%$ and $6\%$ respectively vs. $4\%$ or less for all other groups).  
462     + Black and hispanic drivers account for a higher proportion of stops resulting in searches than stops not resulting in searches ($30\%$ and $18\%$ vs. $20\%$ and $13\%$ respectively).  
463
464 * Over time, the tpyical number of traffic stops, stops resulting in searches, and the proportion of stops resulting in searches fluctuate widely, but have fairly stationary central trends.  
465     + The typical number of monthly stops dipped a bit between 2014-2015.  
466     + The typical number of monthly stops resulting in searches may have decreased 2012-2015, but held flat 2015-2018.  
467     + The typical monthly proportion of stops resulting in searches increased a little 2012-2014, but then declined a little bit 2014-2017.  
468     
469 * Comparing across the racial/ethnic groups identified in the data, the *typical monthly numbers of stops and searches* vary substantially and exhibit distinct trends.  
470     + The typical number of black and hispanic drivers stopped per month has increased.  
471     + The typical number of white drivers stopped per month has decreased.  
472     + The typical number of asian/pacific islander and "other" category drivers per month has remained nearly flat.  
473     + The typical number of stops resulting in searches among black, white, and hispanic drivers per month have tended to be quite close to each other, and sometimes overlap,
474     + The typical number of stops resulting in searches among white drivers per month has increased slightly.
475     + The typical number of stops resulting in searches among black and hispanic drivers per month has decreased slightly.
476     
477 * The typical *proportion of stops resulting in searches within racial/ethnic groups* have shifted somewhat over time.  
478     + The proportion has increased slightly among white and asian/pacific islander drivers.  
479     + The proportion has decreased slightly among black and hispanice drivers.  
480     
481 * The typical *proportion of total stops accounted for by each racial/ethnic group* has shifted somewhat over time.  
482     + The proportion has decreased substantially among white drivers.  
483     + The proportion has increased among black drivers (substantially) and among hispanic drivers (slightly).  
484     
485 * The typical *proportion of total searches accounted for by each racial/ethnic group* has also shifted over time.  
486     + As of 2012, the proportion accounted for by black and white drivers was nearly identical (close to $40\%$ each)!
487     + The proprtion has increased substantially among white drivers (up to .  
488     + The proportion has decreased substantially among black and hispanic drivers. 
489
490 ## SQ2. Contextualize SQ1 interpretation in relation to PC6  
491
492 Several noteworthy comparisons come looking across the different proportions for traffic stops and searches and comparing those against the baseline population proportions accounted for by each racial/ethnic group category. Here are several specific observations:  
493
494 * The white proportion of the state population is consistently larger than either the proportion of white drivers involved in all traffic stops or the proportion of white drivers involved in all stops resulting in searches.  
495 * The hispanic proportion of the state population is about the same as the proportion of hispanic drivers involved in all traffic stops and usually a little smaller than the proportion of hispanic drivers involved in all stops resulting in searches (although this latter difference shrank over the period of data collection).  
496 * The black proportion of the state population is consistently and substantially smaller than the proportion of black drivers involved in all traffic stops and consistently and substantially smaller than the proportion of black drivers involved in all stops resulting in searches. This latter difference in particular stands out as it is almost more than double the baseline population proportion.  
497 * Overall, these results suggest that race/ethnicity and traffic stops and searches are not indepdendent in a statistical sense. The nature of their relationship is complex and could be due to a number of factors including biased policing practices, socioeconomic inequalities, different levels of behaviors linked to traffic stops within specific sub-groups, regional variations masked by these state-level analyses, political/cultural factors shaping police/driver interactions, and more.  
498
499 ## SQ3. Limitations and possible extensions  
500
501 Again, many possible things worth mentioning here, so I'll provide a few that stand out to me.  
502
503 * The generalizability of analysis focused on one state during one 6 year period is limited.
504 * Working with a random $1\%$ sample of the full dataset means that our results here could diverge from those we would find in an analysis of the full population of traffic stops in unpredictable ways. That said, even the very small sample is quite big and once you've read *OpenIntro* ยง5 you'll have some tools to estimate standard errors and confidence intervals around the various results from this analysis.   
505 * The data seem very prone to measurement errors of various kinds. In particular, I suspect the race/ethnicity classifications provided by officers are subject to some biases that are hard to identify and might also shift over time/region. The prevalence of missing values during the first two years of the dataset illustrate one aspect of this and may impact estimates of raw counts and proportions.  
506 * While the comparisons across racial/ethnic groups and between the traffic stops/searches and baseline population proportions illustrates a number of suggestive patterns, conclusive interpretation or attribution of those patterns to any specific cause or causes is quite difficult in the absence of additional information or assumptions. For one example, see my comments regarding statistical independence and the possible explanations in SQ2 above. 
507 * Extensions of this analysis might seek to investigate how some of the patterns identified in the aggregate sate-level data vary across sub-regions (e.g., counties or police districts) or even in comparison to other states.  
508 * Another extension might investigate how specific policy and/or socioeconomic changes may or may not relate to shifting patterns in the data.
509 * Yet another extension could try to identify some factor that introduces some sudden, uncontrolled (by the drivers or the police) change into the process of traffic stops that could make it possible to identify and estimate the effects of specific explanatory factors. The original SOPP paper is an extraordinary example of this that we can discuss! 

Community Data Science Collective || Want to submit a patch?