/
linear-interpolation.Rmd
247 lines (202 loc) · 10.2 KB
/
linear-interpolation.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
---
title: "Linearly Interpolate Farm(er) Characteristics & Build True-Test/True-Train Datasets"
author: "Britta Schumacher"
date: "Last updated: `r format(Sys.time(), '%d %B, %Y')`"
output:
html_document:
theme: yeti
toc: yes
toc_float: true
code_folding: hide
---
```{r setup, include = FALSE, cache = FALSE}
knitr::opts_chunk$set(
warning = FALSE,
message = FALSE
)
```
## Load packages & data
```{r}
library(tidyverse)
library(zoo)
yield <- readRDS("./data/data-raw/yield-modeling-panel.RDS")
yield <- yield %>% mutate(perc_corn = CORN_SQKM/AG_SQKM*100)
yr_df <- yield %>% select(GEOID, YEAR)
cty <- readRDS("./data/county.RDS")
cty <- as.data.frame(cty)
coterm <- yield %>% filter(GEOID %in% cty$GEOID) # check!
# bring in extra coa variables of interest
coa <- readRDS("./data/data-raw/BSDS_standatt_operated.RDS")
coa2 <- readRDS("./data/data-raw/COA-vars.RDS")
human <- coa2 %>% select(GEOID, county, year, chem, fert, labor_expense, machinery, insur_acres, gvt_prog, exp, occup, tenant, acres_per_op) %>%
filter(year == 2017 | year == 2012 | year == 2007) %>% rename(YEAR = year)
human <- merge(human, yr_df, by = c("GEOID","YEAR"), all = TRUE)
```
## Interpolate
### 2008 - 2011
```{r}
# interpolate 2008-2011 using the linear interpolation of 2007 and 2012 data for each GEOID
# NAs produced when no data is reported in 2007 or 2012
test <- human %>%
filter(YEAR == c(2007:2012)) %>%
group_by(GEOID) %>%
mutate(chem=ifelse(is.na(chem),na.approx(chem,na.rm=TRUE),chem),
fert=ifelse(is.na(fert),na.approx(fert,na.rm=TRUE),fert),
labor_expense=ifelse(is.na(labor_expense),na.approx(labor_expense,na.rm=TRUE),labor_expense),
machinery=ifelse(is.na(machinery),na.approx(machinery,na.rm=TRUE),machinery),
insur_acres=ifelse(is.na(insur_acres),na.approx(insur_acres,na.rm=TRUE),insur_acres),
gvt_prog=ifelse(is.na(gvt_prog),na.approx(gvt_prog,na.rm=TRUE),gvt_prog),
exp=ifelse(is.na(exp),na.approx(exp,na.rm=TRUE),exp),
occup=ifelse(is.na(occup),na.approx(occup,na.rm=TRUE),occup),
tenant=ifelse(is.na(tenant),na.approx(tenant,na.rm=TRUE),tenant),
acres_per_op=ifelse(is.na(acres_per_op),na.approx(acres_per_op,na.rm=TRUE),acres_per_op)) %>%
fill(county, .direction = "updown") %>%
filter(YEAR != 2012)
```
### 2013 - 2016
```{r}
# interpolate 2013-2016 using the average of 2012 and 2017 data for each GEOID
# NAs produced when no data is reported in 2012 or 2017
test2 <- human %>%
filter(YEAR == 2012 | YEAR == 2013 | YEAR == 2014 | YEAR == 2015 | YEAR == 2016 | YEAR == 2017) %>%
group_by(GEOID) %>%
mutate(chem=ifelse(is.na(chem),na.approx(chem,na.rm=TRUE),chem),
fert=ifelse(is.na(fert),na.approx(fert,na.rm=TRUE),fert),
labor_expense=ifelse(is.na(labor_expense),na.approx(labor_expense,na.rm=TRUE),labor_expense),
machinery=ifelse(is.na(machinery),na.approx(machinery,na.rm=TRUE),machinery),
insur_acres=ifelse(is.na(insur_acres),na.approx(insur_acres,na.rm=TRUE),insur_acres),
gvt_prog=ifelse(is.na(gvt_prog),na.approx(gvt_prog,na.rm=TRUE),gvt_prog),
exp=ifelse(is.na(exp),na.approx(exp,na.rm=TRUE),exp),
occup=ifelse(is.na(occup),na.approx(occup,na.rm=TRUE),occup),
tenant=ifelse(is.na(tenant),na.approx(tenant,na.rm=TRUE),tenant),
acres_per_op=ifelse(is.na(acres_per_op),na.approx(acres_per_op,na.rm=TRUE),acres_per_op)) %>%
fill(county, .direction = "updown")
```
### 2018
```{r}
# interpolate 2018 using 2017 data for each GEOID
# NAs produced when no data is reported in 2017
test3 <- human %>%
filter(YEAR == 2017 | YEAR == 2018) %>%
group_by(GEOID) %>%
mutate(chem=ifelse(is.na(chem),mean(chem,na.rm=TRUE),chem),
fert=ifelse(is.na(fert),mean(fert,na.rm=TRUE),fert),
labor_expense=ifelse(is.na(labor_expense),mean(labor_expense,na.rm=TRUE),labor_expense),
machinery=ifelse(is.na(machinery),mean(machinery,na.rm=TRUE),machinery),
insur_acres=ifelse(is.na(insur_acres),mean(insur_acres,na.rm=TRUE),insur_acres),
gvt_prog=ifelse(is.na(gvt_prog),mean(gvt_prog,na.rm=TRUE),gvt_prog),
exp=ifelse(is.na(exp),mean(exp,na.rm=TRUE),exp),
occup=ifelse(is.na(occup),mean(occup,na.rm=TRUE),occup),
tenant=ifelse(is.na(tenant),mean(tenant,na.rm=TRUE),tenant),
acres_per_op=ifelse(is.na(acres_per_op),mean(acres_per_op,na.rm=TRUE),acres_per_op)) %>%
fill(county, .direction = "updown") %>%
filter(YEAR == 2018)
h <- rbind(test, test2, test3)
is.nan.data.frame <- function(x)
do.call(cbind, lapply(x, is.nan))
h[is.nan(h)] <- NA
```
### Which counties still have NA values?
```{r}
# Which counties have NAs after taking the average between census years?
na_vals <- filter_all(h, any_vars(is.na(.)))
na_count <- sapply(h, function(x) sum(is.na(x)))
# find counties with 12 NAs/column and remove; these counties have no data for these variables across time (2007 - 2012 - or 2017)
remove <- h %>% group_by(GEOID) %>% summarise_all(funs(sum(is.na(.)))) %>% filter_all(any_vars(. == 12)) # remove n = 234 counties
remove <- merge(remove, coa, by = "GEOID")
remove <- remove %>% select(GEOID, county.y) %>% unique()
write.csv(remove, "./results/impute-remove.csv")
h <- h %>% filter(!GEOID %in% remove$GEOID) # n = 2874 remaining counties with at least some data!
```
### 2018 using interpolated values from 2012-2017
```{r}
# interpolate 2018 using interpolated values from 2012-2017 for each GEOID
# NAs remain when no data is reported in 2012 or 2017
test4 <- h %>%
filter(YEAR == 2017 | YEAR == 2018) %>%
group_by(GEOID) %>%
mutate(chem=ifelse(is.na(chem),na.approx(chem,na.rm=TRUE),chem),
fert=ifelse(is.na(fert),na.approx(fert,na.rm=TRUE),fert),
labor_expense=ifelse(is.na(labor_expense),na.approx(labor_expense,na.rm=TRUE),labor_expense),
machinery=ifelse(is.na(machinery),na.approx(machinery,na.rm=TRUE),machinery),
insur_acres=ifelse(is.na(insur_acres),na.approx(insur_acres,na.rm=TRUE),insur_acres),
gvt_prog=ifelse(is.na(gvt_prog),na.approx(gvt_prog,na.rm=TRUE),gvt_prog),
exp=ifelse(is.na(exp),na.approx(exp,na.rm=TRUE),exp),
occup=ifelse(is.na(occup),na.approx(occup,na.rm=TRUE),occup),
tenant=ifelse(is.na(tenant),na.approx(tenant,na.rm=TRUE),tenant),
acres_per_op=ifelse(is.na(acres_per_op),na.approx(acres_per_op,na.rm=TRUE),acres_per_op)) %>%
fill(county, .direction = "updown") %>%
filter(YEAR == 2018)
h <- h %>% filter(YEAR != 2018)
h <- rbind(h, test4)
# Which counties STILL have NAs?
na_vals <- filter_all(h, any_vars(is.na(.)))
na_count <- sapply(h, function(x) sum(is.na(x)))
```
### 2007 - 2012 using interpolated values from 2012-2017
```{r}
# interpolate 2007 to 2012 using interpolated 2012 data
# NAs remain when no data is reported in 2012 or 2017
test5 <- h %>%
filter(YEAR == c(2007:2012)) %>%
group_by(GEOID) %>%
mutate(chem=ifelse(is.na(chem),na.approx(chem,na.rm=TRUE),chem),
fert=ifelse(is.na(fert),na.approx(fert,na.rm=TRUE),fert),
labor_expense=ifelse(is.na(labor_expense),na.approx(labor_expense,na.rm=TRUE),labor_expense),
machinery=ifelse(is.na(machinery),na.approx(machinery,na.rm=TRUE),machinery),
insur_acres=ifelse(is.na(insur_acres),na.approx(insur_acres,na.rm=TRUE),insur_acres),
gvt_prog=ifelse(is.na(gvt_prog),na.approx(gvt_prog,na.rm=TRUE),gvt_prog),
exp=ifelse(is.na(exp),na.approx(exp,na.rm=TRUE),exp),
occup=ifelse(is.na(occup),na.approx(occup,na.rm=TRUE),occup),
tenant=ifelse(is.na(tenant),na.approx(tenant,na.rm=TRUE),tenant),
acres_per_op=ifelse(is.na(acres_per_op),na.approx(acres_per_op,na.rm=TRUE),acres_per_op)) %>%
fill(county, .direction = "updown") %>%
filter(YEAR != 2012)
h <- h %>% filter(YEAR != 2007 & YEAR != 2008 & YEAR != 2009 & YEAR != 2010 & YEAR != 2011)
h <- rbind(h,test5)
h[is.nan(h)] <- NA
na_vals <- filter_all(h, any_vars(is.na(.)))
na_count <- sapply(h, function(x) sum(is.na(x)))
```
### 2012 - 2018 using interpolated values from 2007-2011
```{r}
# fill in 2012-2018 with 2011 data (interp from 2007) in counties that have no data reported in 2012 or 2017
test6 <- h %>%
filter(YEAR == 2011 | YEAR == 2012 | YEAR == 2013 | YEAR == 2014 | YEAR == 2015 | YEAR == 2016 | YEAR == 2017 | YEAR == 2018) %>%
group_by(GEOID) %>%
mutate(chem=ifelse(is.na(chem),na.approx(chem,na.rm=TRUE),chem),
fert=ifelse(is.na(fert),na.approx(fert,na.rm=TRUE),fert),
labor_expense=ifelse(is.na(labor_expense),na.approx(labor_expense,na.rm=TRUE),labor_expense),
machinery=ifelse(is.na(machinery),na.approx(machinery,na.rm=TRUE),machinery),
insur_acres=ifelse(is.na(insur_acres),na.approx(insur_acres,na.rm=TRUE),insur_acres),
gvt_prog=ifelse(is.na(gvt_prog),na.approx(gvt_prog,na.rm=TRUE),gvt_prog),
exp=ifelse(is.na(exp),na.approx(exp,na.rm=TRUE),exp),
occup=ifelse(is.na(occup),na.approx(occup,na.rm=TRUE),occup),
tenant=ifelse(is.na(tenant),na.approx(tenant,na.rm=TRUE),tenant),
acres_per_op=ifelse(is.na(acres_per_op),na.approx(acres_per_op,na.rm=TRUE),acres_per_op)) %>%
fill(county, .direction = "down") %>%
filter(YEAR != 2011)
h <- h %>% filter(YEAR != 2012 & YEAR != 2013 & YEAR != 2014 & YEAR != 2015 & YEAR != 2016 & YEAR != 2017 & YEAR != 2018)
h <- rbind(h,test6)
h[is.nan(h)] <- NA
na_vals <- filter_all(h, any_vars(is.na(.)))
na_count <- sapply(h, function(x) sum(is.na(x)))
# SAVE INTERPOLATED DATA
saveRDS(h, "./data/data-out/farm-er-linearinterp.RDS")
h <- readRDS("./data/data-out/farm-er-linearinterp.RDS")
```
## Yield data
```{r}
# how many counties report yield?
y_present <- coterm %>% filter(GEOID %in% h$GEOID,
!is.na(YIELD)) # reduces to 16983 unique county-year observations
bio_farm_train <- merge(y_present, h , by = c("GEOID", "YEAR"))
saveRDS(bio_farm_train, "./data/data-out/true-train.RDS")
# how many counties are missing yield?
y_missing <- coterm %>% filter(GEOID %in% h$GEOID,
is.na(YIELD)) # 14631
bio_farm_infill <- merge(y_missing, h, by = c("GEOID", "YEAR"))
saveRDS(bio_farm_infill, "./data/data-out/true-test-infillobs.RDS")
# for counties reporting yield, the number of years with reports (range = 1 - 11)
y_count <- coterm %>% group_by(GEOID) %>% filter(!is.na(YIELD)) %>% summarise(sum = sum(length(YIELD)))
```