Load all of the packages and the data
Loading Libraries
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.4.4 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.0
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
── Attaching packages ────────────────────────────────────── tidymodels 1.1.1 ──
✔ broom 1.0.5 ✔ rsample 1.2.0
✔ dials 1.2.0 ✔ tune 1.1.2
✔ infer 1.0.6 ✔ workflows 1.1.3
✔ modeldata 1.3.0 ✔ workflowsets 1.0.1
✔ parsnip 1.1.1 ✔ yardstick 1.3.0
✔ recipes 1.0.9
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ scales::discard() masks purrr::discard()
✖ dplyr::filter() masks stats::filter()
✖ recipes::fixed() masks stringr::fixed()
✖ dplyr::lag() masks stats::lag()
✖ yardstick::spec() masks readr::spec()
✖ recipes::step() masks stats::step()
• Learn how to get started at https://www.tidymodels.org/start/
Load the Data
aqidata <- read.csv ('aqidataset.csv' )
Restructuring the data, creating additional variables if needed
Turn each column of population estimate into a variable
pop2020 <- c (na.omit (aqidata$ POPESTIMATE2020))
pop2021 <- c (na.omit (aqidata$ POPESTIMATE2021))
pop2022 <- c (na.omit (aqidata$ POPESTIMATE2022))
pop2023 <- c (na.omit (aqidata$ POPESTIMATE2023))
Concatenate all of them into a single list
popestimate <- c (pop2020, pop2021, pop2022, pop2023)
Add the column into the dataframe and remove the separate popestiamte columns
aqi1 <- aqidata |>
select (! c (POPESTIMATE2020,POPESTIMATE2021,POPESTIMATE2022,POPESTIMATE2023)) |>
mutate (popestimate)
Do the same thing for popchanges
# Create npopchange variables
pchange2020 <- c (na.omit (aqidata$ NPOPCHG2020))
pchange2021 <- c (na.omit (aqidata$ NPOPCHG2021))
pchange2022 <- c (na.omit (aqidata$ NPOPCHG2022))
pchange2023 <- c (na.omit (aqidata$ NPOPCHG2023))
# concatenate
popchange <- c (pchange2020,pchange2021,pchange2022,pchange2023)
aqi2 <- aqi1 |>
select (! c (NPOPCHG2020,NPOPCHG2021,NPOPCHG2022,NPOPCHG2023)) |>
mutate (popchange)
Create a population density variable
aqi3 <- aqi2 |>
mutate (pop_per_sqmi = popestimate / ALAND_SQMI)
Now lets make some cool charts!!!
First, lets make a clarification. AQI data has not been taken for every single day for every single county. So our proportion of days with good AQI will be a quotient of Days with good AQI and days with AQI, NOT days with good AQI divided by the total number of days in the year. Because of this, we need to check out all of the counties and see how many days they had with AQI recorded.
ggplot (data = aqi3, aes (x = Days.with.AQI)) +
geom_boxplot ()
We can see from this boxplot that a vast majority of the observations (county for each year 2020-2023) have more than 100 Days with AQI on record. 100 is an arbitrary number and represents a minorty of all days throughout the year but it is the number I will be chosing for the immediate preceeding part of this study.
The reason I am choosing to explore this is to make sure that there will not be any observations that massively skew this study.
# First, we'll see how many counties have less than 100 days of AQI recorded
aqi_by_year <- aqi3 |>
group_by (Year) |>
count (aqidays100 = c (Days.with.AQI > 100 )) |>
filter (aqidays100 == FALSE )
aqi_by_year
# A tibble: 4 × 3
# Groups: Year [4]
Year aqidays100 n
<int> <lgl> <int>
1 2020 FALSE 7
2 2021 FALSE 5
3 2022 FALSE 3
4 2023 FALSE 43
Based of the information in this table, we can see that there are not very many observations in the years 2020, 2021, and 2022 that fall below the chosen cutoff of 100 days with AQI, but 2023 has far more counties which do not report enough days with AQI to meet this cutoff.
Lets filter the main dataframe by this cutoff to get a better idea of the data we’d be discarding for our plots to see if we should adjust this level, or if we should scrap.
aqi_parsed <- aqi3 |>
mutate (aqidays100 = c (Days.with.AQI > 100 )) |>
filter (aqidays100 == FALSE )
aqi_parsed
X State County Year Days.with.AQI Good.Days Moderate.Days
1 96 Colorado Fremont 2020 17 17 0
2 107 Colorado Pueblo 2020 96 91 5
3 359 Michigan Monroe 2020 61 58 3
4 547 Ohio Athens 2020 86 83 3
5 553 Ohio Columbiana 2020 48 48 0
6 766 Virginia Carroll 2020 56 56 0
7 819 West Virginia Marion 2020 74 69 5
8 96 Colorado Fremont 2021 86 85 1
9 359 Michigan Monroe 2021 61 57 4
10 553 Ohio Columbiana 2021 61 61 0
11 766 Virginia Carroll 2021 60 60 0
12 819 West Virginia Marion 2021 91 78 13
13 359 Michigan Monroe 2022 61 51 8
14 553 Ohio Columbiana 2022 61 61 0
15 766 Virginia Carroll 2022 56 56 0
16 27 Arkansas Arkansas 2023 90 71 19
17 28 Arkansas Ashley 2023 89 66 22
18 31 Arkansas Garland 2023 76 59 17
19 32 Arkansas Jackson 2023 91 70 21
20 36 Arkansas Union 2023 90 65 24
21 57 California Mendocino 2023 90 47 43
22 77 California Siskiyou 2023 90 82 8
23 96 Colorado Fremont 2023 80 79 1
24 97 Colorado Garfield 2023 90 85 5
25 105 Colorado Pitkin 2023 68 68 0
26 110 Colorado San Miguel 2023 76 75 1
27 136 Florida Okaloosa 2023 57 48 9
28 154 Georgia Clayton 2023 70 53 17
29 156 Georgia Coffee 2023 66 55 11
30 246 Indiana Spencer 2023 69 57 12
31 256 Iowa Cerro Gordo 2023 60 58 2
32 260 Iowa Lee 2023 60 47 13
33 266 Iowa Pottawattamie 2023 82 58 24
34 269 Iowa Woodbury 2023 78 54 23
35 359 Michigan Monroe 2023 40 37 3
36 398 Missouri Andrew 2023 91 39 46
37 399 Missouri Boone 2023 91 50 36
38 401 Missouri Callaway 2023 91 52 34
39 402 Missouri Cass 2023 91 56 31
40 405 Missouri Clinton 2023 91 38 43
41 410 Missouri Lincoln 2023 91 41 35
42 412 Missouri Perry 2023 91 43 37
43 553 Ohio Columbiana 2023 41 40 1
44 640 Pennsylvania Lawrence 2023 96 85 10
45 766 Virginia Carroll 2023 29 29 0
46 817 West Virginia Harrison 2023 53 43 9
47 819 West Virginia Marion 2023 59 50 8
48 827 Wisconsin Columbia 2023 96 50 28
49 830 Wisconsin Door 2023 96 64 25
50 832 Wisconsin Fond du Lac 2023 96 64 28
51 835 Wisconsin Jefferson 2023 96 52 30
52 837 Wisconsin Kewaunee 2023 96 57 33
53 839 Wisconsin Manitowoc 2023 96 69 26
54 840 Wisconsin Marathon 2023 96 59 28
55 844 Wisconsin Racine 2023 96 56 28
56 845 Wisconsin Rock 2023 96 43 32
57 847 Wisconsin Sheboygan 2023 96 56 27
58 850 Wisconsin Walworth 2023 96 51 29
Unhealthy.for.Sensitive.Groups.Days Unhealthy.Days Very.Unhealthy.Days
1 0 0 0
2 0 0 0
3 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0
7 0 0 0
8 0 0 0
9 0 0 0
10 0 0 0
11 0 0 0
12 0 0 0
13 1 0 0
14 0 0 0
15 0 0 0
16 0 0 0
17 1 0 0
18 0 0 0
19 0 0 0
20 1 0 0
21 0 0 0
22 0 0 0
23 0 0 0
24 0 0 0
25 0 0 0
26 0 0 0
27 0 0 0
28 0 0 0
29 0 0 0
30 0 0 0
31 0 0 0
32 0 0 0
33 0 0 0
34 0 1 0
35 0 0 0
36 6 0 0
37 5 0 0
38 5 0 0
39 4 0 0
40 10 0 0
41 15 0 0
42 11 0 0
43 0 0 0
44 1 0 0
45 0 0 0
46 0 1 0
47 0 1 0
48 16 2 0
49 6 1 0
50 4 0 0
51 14 0 0
52 5 1 0
53 1 0 0
54 8 1 0
55 11 1 0
56 19 2 0
57 12 1 0
58 16 0 0
Hazardous.Days Max.AQI X90th.Percentile.AQI Median.AQI Days.CO Days.NO2
1 0 27 20 11 0 0
2 0 63 45 19 0 0
3 0 66 40 17 0 0
4 0 61 41 24 0 0
5 0 31 20 10 0 0
6 0 15 12 7 0 0
7 0 67 47 26 0 0
8 0 62 28 15 0 0
9 0 77 45 19 0 0
10 0 35 24 16 0 0
11 0 33 16 9 0 0
12 0 80 53 30 0 0
13 1 596 59 26 0 0
14 0 40 23 12 0 0
15 0 23 13 7 0 0
16 0 78 59 35 0 0
17 0 107 58 34 0 0
18 0 63 56 34 0 0
19 0 73 61 35 0 0
20 0 102 59 39 0 0
21 0 87 72 49 0 0
22 0 60 50 38 0 0
23 0 71 23 13 0 0
24 0 71 50 43 0 0
25 0 46 28 13 0 0
26 0 54 20 12 0 0
27 0 84 64 42 0 0
28 0 97 60 38 0 0
29 0 87 56 30 0 0
30 0 84 55 34 0 0
31 0 56 41 19 0 0
32 0 97 60 35 0 0
33 0 97 61 38 0 0
34 0 164 59 35 0 0
35 0 82 45 26 0 0
36 0 122 97 54 0 0
37 0 143 90 49 0 0
38 0 133 90 49 0 0
39 0 143 80 47 0 0
40 0 133 101 54 0 0
41 0 140 108 51 0 0
42 0 150 101 51 0 0
43 0 74 26 11 0 0
44 0 119 54 41 0 0
45 0 45 19 7 0 0
46 0 162 60 29 0 0
47 0 165 58 31 0 0
48 0 159 115 50 0 0
49 0 159 100 45 0 0
50 0 122 84 44 0 0
51 0 140 105 49 0 0
52 0 172 97 46 0 0
53 0 133 80 44 0 0
54 0 166 100 46 0 0
55 0 151 108 47 0 0
56 0 177 126 54 0 0
57 0 182 112 49 0 0
58 0 143 122 49 0 0
Days.Ozone Days.PM2.5 Days.PM10 County.Specific ALAND_SQMI
1 0 0 17 Fremont County Colorado 1533.835
2 0 49 47 Pueblo County Colorado 2386.591
3 0 0 61 Monroe County Michigan 549.353
4 0 86 0 Athens County Ohio 503.643
5 0 0 48 Columbiana County Ohio 531.917
6 0 0 56 Carroll County Virginia 474.744
7 0 74 0 Marion County West Virginia 308.785
8 0 0 86 Fremont County Colorado 1533.835
9 0 0 61 Monroe County Michigan 549.353
10 0 0 61 Columbiana County Ohio 531.917
11 0 0 60 Carroll County Virginia 474.744
12 0 91 0 Marion County West Virginia 308.785
13 0 0 61 Monroe County Michigan 549.353
14 0 0 61 Columbiana County Ohio 531.917
15 0 0 56 Carroll County Virginia 474.744
16 0 90 0 Arkansas County Arkansas 992.140
17 0 89 0 Ashley County Arkansas 925.492
18 0 76 0 Garland County Arkansas 677.615
19 0 91 0 Jackson County Arkansas 633.611
20 0 90 0 Union County Arkansas 1039.067
21 21 69 0 Mendocino County California 3506.807
22 58 32 0 Siskiyou County California 6278.767
23 0 0 80 Fremont County Colorado 1533.835
24 90 0 0 Garfield County Colorado 2947.425
25 0 0 68 Pitkin County Colorado 970.701
26 0 0 76 San Miguel County Colorado 1286.725
27 52 0 5 Okaloosa County Florida 929.862
28 0 70 0 Clayton County Georgia 141.756
29 0 66 0 Coffee County Georgia 592.247
30 0 69 0 Spencer County Indiana 411.136
31 0 0 60 Cerro Gordo County Iowa 568.300
32 0 60 0 Lee County Iowa 517.550
33 0 78 4 Pottawattamie County Iowa 951.272
34 0 78 0 Woodbury County Iowa 872.896
35 0 0 40 Monroe County Michigan 549.353
36 91 0 0 Andrew County Missouri 432.632
37 91 0 0 Boone County Missouri 685.551
38 91 0 0 Callaway County Missouri 834.568
39 89 2 0 Cass County Missouri 696.615
40 91 0 0 Clinton County Missouri 418.940
41 91 0 0 Lincoln County Missouri 626.562
42 91 0 0 Perry County Missouri 474.355
43 0 0 41 Columbiana County Ohio 531.917
44 96 0 0 Lawrence County Pennsylvania 357.390
45 0 0 29 Carroll County Virginia 474.744
46 0 53 0 Harrison County West Virginia 416.008
47 0 59 0 Marion County West Virginia 308.785
48 96 0 0 Columbia County Wisconsin 765.600
49 96 0 0 Door County Wisconsin 481.936
50 96 0 0 Fond du Lac County Wisconsin 719.566
51 96 0 0 Jefferson County Wisconsin 556.492
52 96 0 0 Kewaunee County Wisconsin 342.454
53 96 0 0 Manitowoc County Wisconsin 589.302
54 96 0 0 Marathon County Wisconsin 1545.194
55 96 0 0 Racine County Wisconsin 332.611
56 96 0 0 Rock County Wisconsin 718.167
57 96 0 0 Sheboygan County Wisconsin 511.519
58 96 0 0 Walworth County Wisconsin 555.433
popestimate popchange pop_per_sqmi aqidays100
1 48878 -56 31.866531 FALSE
2 168337 172 70.534499 FALSE
3 154896 94 281.960779 FALSE
4 62393 -44 123.883386 FALSE
5 101642 -237 191.086203 FALSE
6 29148 -10 61.397300 FALSE
7 56240 24 182.133199 FALSE
8 49237 359 32.100584 FALSE
9 155534 638 283.122146 FALSE
10 101352 -290 190.541005 FALSE
11 29098 -50 61.291981 FALSE
12 56053 -187 181.527600 FALSE
13 155235 -299 282.577869 FALSE
14 100569 -783 189.068971 FALSE
15 29118 20 61.334108 FALSE
16 16307 -191 16.436188 FALSE
17 18262 -86 19.732207 FALSE
18 99784 -198 147.257661 FALSE
19 16784 130 26.489439 FALSE
20 37397 -416 35.990942 FALSE
21 89108 -899 25.410010 FALSE
22 42905 -881 6.833348 FALSE
23 50318 702 32.805354 FALSE
24 62707 451 21.275181 FALSE
25 16640 -207 17.142251 FALSE
26 7868 -140 6.114749 FALSE
27 218464 2000 234.942389 FALSE
28 298300 -99 2104.320099 FALSE
29 43317 70 73.140092 FALSE
30 45231 42 110.014691 FALSE
31 42406 -20 74.619039 FALSE
32 32565 -264 62.921457 FALSE
33 93179 114 97.952005 FALSE
34 105951 393 121.378721 FALSE
35 155045 -190 282.232007 FALSE
36 18127 111 41.899351 FALSE
37 189463 1742 276.366018 FALSE
38 44731 8 53.597790 FALSE
39 111732 1375 160.392756 FALSE
40 21548 245 51.434573 FALSE
41 64699 1517 103.260332 FALSE
42 18950 72 39.948983 FALSE
43 100182 -387 188.341414 FALSE
44 84472 -340 236.358040 FALSE
45 29239 121 61.588983 FALSE
46 64639 -113 155.379223 FALSE
47 55807 23 180.730929 FALSE
48 58091 -137 75.876437 FALSE
49 30562 11 63.415059 FALSE
50 103948 136 144.459299 FALSE
51 85743 18 154.077687 FALSE
52 20690 37 60.416874 FALSE
53 81331 133 138.012428 FALSE
54 138612 657 89.705241 FALSE
55 196613 652 591.119957 FALSE
56 164278 437 228.746239 FALSE
57 117752 91 230.200638 FALSE
58 105822 158 190.521629 FALSE
I need to review with Prof. Saidi on what I can/should do with this information. But this is a list of observations that should be datapoints that are taken with a grain of salt.
Here, we will additionally explore possible issues with this dataset and display a dot plot which shows the correlation (or lack thereof) between days with AQI and the proportion of good days.
ggplot (data = aqi3, aes (x = Days.with.AQI, y = Good.Days/ Days.with.AQI)) +
geom_point ()
This is a very strange looking plot with many noticeable peaks. But we can see that there is a fan shaped negative correlation showcasing the higher proportion of days with good AQI as the amount of days with good AQI decreases.
This plot reiterates the concern that counties with lower days reported could skew our results in studying.
ggplot (data = aqi3, aes (x = log (pop_per_sqmi), y = Median.AQI)) +
geom_point () +
geom_smooth (method = "lm" )
`geom_smooth()` using formula = 'y ~ x'