깔끔한 데이터 원칙으로 작물 수확량 🌽🍚🌾 모델링

최근에 저는 screencasts 프레임워크를 시작하는 것부터 더 복잡한 모델 튜닝에 이르기까지 사용 방법을 보여주는 tidymodels을 게시하고 있습니다. 오늘의 스크린캐스트에서는 작물 수확량에 대한 이번 주 #TidyTuesday dataset를 사용하여 많은 모델을 구축하는 작업에 깔끔한 데이터 원칙을 유창하게 적용하는 방법을 살펴봅니다.



다음은 비디오 대신 또는 비디오에 추가하여 읽기를 선호하는 사람들을 위해 비디오에서 사용한 코드입니다.

데이터 탐색

우리의 모델링 목표는 crops yields are changing around the world using this week’s #TidyTuesday dataset 방법을 추정하는 것입니다. 관심 있는 국가-작물 조합에 대해 많은 모델을 구축할 수 있습니다.

먼저, 이번 주의 두 가지 데이터 세트를 읽어 보겠습니다.

library(tidyverse)

key_crop_yields <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-09-01/key_crop_yields.csv")
land_use <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-09-01/land_use_vs_yield_change_in_cereal_production.csv")

land_use 데이터 세트를 사용하여 인구가 가장 많은 국가를 찾을 것입니다. 이름 벡터를 만들어 봅시다.

top_countries <- land_use %>%
  janitor::clean_names() %>%
  filter(!is.na(code), entity != "World") %>%
  group_by(entity) %>%
  filter(year == max(year)) %>%
  ungroup() %>%
  slice_max(total_population_gapminder, n = 30) %>%
  pull(entity)

top_countries


## [1] "China" "India"                       
## [3] "United States" "Indonesia"                   
## [5] "Pakistan" "Brazil"                      
## [7] "Nigeria" "Bangladesh"                  
## [9] "Russia" "Mexico"                      
## [11] "Japan" "Ethiopia"                    
## [13] "Philippines" "Egypt"                       
## [15] "Vietnam" "Democratic Republic of Congo"
## [17] "Germany" "Turkey"                      
## [19] "Iran" "Thailand"                    
## [21] "United Kingdom" "France"                      
## [23] "Italy" "South Africa"                
## [25] "Tanzania" "Myanmar"                     
## [27] "Kenya" "South Korea"                 
## [29] "Colombia" "Spain"

이제 관심 있는 국가 및 작물에 대한 깔끔한 버전의 작물 수확량 데이터를 만들어 보겠습니다.

tidy_yields <- key_crop_yields %>%
  janitor::clean_names() %>%
  pivot_longer(wheat_tonnes_per_hectare:bananas_tonnes_per_hectare,
    names_to = "crop", values_to = "yield"
  ) %>%
  mutate(crop = str_remove(crop, "_tonnes_per_hectare")) %>%
  filter(
    crop %in% c("wheat", "rice", "maize", "barley"),
    entity %in% top_countries,
    !is.na(yield)
  )

tidy_yields


## # A tibble: 6,032 x 5
## entity code year crop yield
## <chr> <chr> <dbl> <chr> <dbl>
## 1 Bangladesh BGD 1961 wheat 0.574
## 2 Bangladesh BGD 1961 rice 1.70 
## 3 Bangladesh BGD 1961 maize 0.799
## 4 Bangladesh BGD 1961 barley 0.577
## 5 Bangladesh BGD 1962 wheat 0.675
## 6 Bangladesh BGD 1962 rice 1.53 
## 7 Bangladesh BGD 1962 maize 0.738
## 8 Bangladesh BGD 1962 barley 0.544
## 9 Bangladesh BGD 1963 wheat 0.607
## 10 Bangladesh BGD 1963 rice 1.77 
## # … with 6,022 more rows

이 데이터 구조는 시간 경과에 따른 작물 수확량을 표시하는 데 적합합니다!

tidy_yields %>%
  ggplot(aes(year, yield, color = crop)) +
  geom_line(alpha = 0.7, size = 1.5) +
  geom_point() +
  facet_wrap(~entity, ncol = 5) +
  scale_x_continuous(guide = guide_axis(angle = 90)) +
  labs(x = NULL, y = "yield (tons per hectare)")

모든 국가가 모든 작물을 생산하는 것은 아니지만 전체 작물 수확량이 증가하고 있음을 알 수 있습니다.

많은 모델

이제 각 국가-작물 조합에 선형 모델을 적용해 보겠습니다.

library(tidymodels)

tidy_lm <- tidy_yields %>%
  nest(yields = c(year, yield)) %>%
  mutate(model = map(yields, ~ lm(yield ~ year, data = .x)))

tidy_lm


## # A tibble: 111 x 5
## entity code crop yields model 
## <chr> <chr> <chr> <list> <list>
## 1 Bangladesh BGD wheat <tibble [58 × 2]> <lm>  
## 2 Bangladesh BGD rice <tibble [58 × 2]> <lm>  
## 3 Bangladesh BGD maize <tibble [58 × 2]> <lm>  
## 4 Bangladesh BGD barley <tibble [58 × 2]> <lm>  
## 5 Brazil BRA wheat <tibble [58 × 2]> <lm>  
## 6 Brazil BRA rice <tibble [58 × 2]> <lm>  
## 7 Brazil BRA maize <tibble [58 × 2]> <lm>  
## 8 Brazil BRA barley <tibble [58 × 2]> <lm>  
## 9 China CHN wheat <tibble [58 × 2]> <lm>  
## 10 China CHN rice <tibble [58 × 2]> <lm>  
## # … with 101 more rows

다음으로 tidy() 해당 모델을 사용하여 계수를 가져오고 그 동안 다중 비교를 위해 p-값을 조정합니다.

slopes <- tidy_lm %>%
  mutate(coefs = map(model, tidy)) %>%
  unnest(coefs) %>%
  filter(term == "year") %>%
  mutate(p.value = p.adjust(p.value))

slopes


## # A tibble: 111 x 10
## entity code crop yields model term estimate std.error statistic p.value
## <chr> <chr> <chr> <list> <lis> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Bangla… BGD wheat <tibbl… <lm> year 0.0389 0.00253 15.4 5.11e-20
## 2 Bangla… BGD rice <tibbl… <lm> year 0.0600 0.00231 26.0 6.05e-31
## 3 Bangla… BGD maize <tibbl… <lm> year 0.122 0.0107 11.3 1.82e-14
## 4 Bangla… BGD barl… <tibbl… <lm> year 0.00505 0.000596 8.47 4.34e-10
## 5 Brazil BRA wheat <tibbl… <lm> year 0.0366 0.00222 16.5 2.55e-21
## 6 Brazil BRA rice <tibbl… <lm> year 0.0755 0.00490 15.4 4.96e-20
## 7 Brazil BRA maize <tibbl… <lm> year 0.0709 0.00395 18.0 4.37e-23
## 8 Brazil BRA barl… <tibbl… <lm> year 0.0466 0.00319 14.6 5.05e-19
## 9 China CHN wheat <tibbl… <lm> year 0.0880 0.00141 62.6 1.72e-51
## 10 China CHN rice <tibbl… <lm> year 0.0843 0.00289 29.2 1.47e-33
## # … with 101 more rows

결과 살펴보기

이제 우리는 이 모델링의 결과를 시각화할 수 있습니다. 이 모델링은 전 세계적으로 작물 수확량이 어떻게 변하고 있는지 추정합니다.

library(ggrepel)
slopes %>%
  ggplot(aes(estimate, p.value, label = entity)) +
  geom_vline(
    xintercept = 0, lty = 2,
    size = 1.5, alpha = 0.7, color = "gray50"
  ) +
  geom_point(aes(color = crop), alpha = 0.8, size = 2.5, show.legend = FALSE) +
  scale_y_log10() +
  facet_wrap(~crop) +
  geom_text_repel(size = 3, family = "IBMPlexSans") +
  theme_light(base_family = "IBMPlexSans") +
  theme(strip.text = element_text(family = "IBMPlexSans-Bold", size = 12)) +
  labs(x = "increase in tons per hectare per year")

x축에는 이러한 모델의 기울기가 있습니다. 대부분의 국가는 작물 수확량이 증가함에 따라 긍정적인 측면에 있습니다. 국가가 오른쪽으로 갈수록 이 기간 동안 작물 수확량의 증가가 더 큽니다. 옥수수 수확량이 가장 많이 증가했습니다.

y축에는 p-값이 있습니다. 이 값은 관계가 없다는 가정(시간에 따른 변화 없음)에서 우리가 보는 효과가 얼마나 놀라운지를 측정합니다. 도표에서 아래쪽에 있는 국가는 p-값이 더 작습니다. 우리는 그것들이 실제 관계임을 더 확신합니다.

이를 확장하여 이러한 모델이 glance()의 데이터에 얼마나 잘 맞는지 확인할 수 있습니다. 많은 하위 그룹의 변화를 한 번에 추정하기 위해 통계 모델을 사용하는 이 접근 방식은 많은 상황에서 저에게 큰 도움이 되었습니다!