散布図と回帰分析によるワイン生産地の相対評価

## 概要
[J.S.A.ソムリエ](https://www.sommelier.jp/exam/exam_guidance.html)の学科試験では世界各地の「ブドウ栽培面積」や「ワイン生産量」が出題されますが、ソムリエ教本上では生産地ごとの比較が無いため、この記事では散布図と回帰分析を用いて各生産地を相対的に評価してみます。

## データソース
GitHubのリポジトリに、2020年のブドウ栽培面積・ワイン生産量をCSVファイルとして配置しました。

https://github.com/29e/statistic/blob/main/source/wine_production_2020.csv

## データの読み込み

GitHub上のCSVファイルををread_csv()で直接読み込みます。

```R
library(tidyverse)
library(ggrepel)
library(gridExtra)

df_raw <- read_csv("https://raw.githubusercontent.com/29e/statistic/main/source/wine_production_2020.csv")

df <- df_raw %>% 
  filter(!is.na(wine_production) == TRUE) # 一部の国で欠損があるためfilterで除外
```

## すべての生産地の散布図

ggplotのgeom_point()を使ってすべての生産地をプロットしてみます。

```R
ggplot(df, aes(x = vineyard_surface_area, y = wine_production, color = country)) +
  geom_text_repel(aes(label = region), size = 3) +
  geom_point()
```

![Rplot1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/147689/6230ac53-df13-08ee-c9a8-56cc87744309.png)

X軸をブドウ栽培面積、Y軸をワイン生産量とした散布図がプロットできました。

右上にイタリア、フランス、スペイン、アメリカが位置しており、これらの国が突出して生産規模が大きいようです。

単位面積あたりの収量（=生産量÷栽培面積）をフランスと比較すると、イタリアは収量が多く、スペインは収量が低いということが分かります。

上位の4生産地の影響でそれ以外の生産地の分布がわかりづらくなってしまっているので、少し調整してみます。

## 主要国を除外する

イタリア、フランス、スペイン、アメリカを除外します。

```R
df_fix <- df %>% 
  filter(region != "France" & region != "Italia" & region != "Spain" & region != "United States")

ggplot(df_fix, aes(x = vineyard_surface_area, y = wine_production, color = country)) +
  geom_text_repel(aes(label = region), size = 3) +
  geom_point()
```

![Rplot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/147689/79da23ae-236a-4ba4-1d48-51254b14e673.png)

さきほどの散布図では見えなかった生産地が可視化されました。

ラングドック・ルーションのワイン生産量はフランス全体の40%を占めており、この図でも南アフリカ・ドイツ・ポルトガル・ルーマニアなどの一国よりも生産規模が大きくなっていることが分かります。

散布図からは一定の線形の関係が読み取れるので、回帰分析での当てはまりの良さを確認してみます。
```R
lm(wine_production ~ vineyard_surface_area, df_france) %>% 
  summary()

## 
## Call:
## lm(formula = wine_production ~ vineyard_surface_area, data = df_fix)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3884772  -418758  -297348    28915  5742747 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.465e+05  3.322e+05   1.043    0.302    
## vineyard_surface_area 5.042e+01  5.164e+00   9.762 5.53e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1734000 on 48 degrees of freedom
## Multiple R-squared:  0.6651, Adjusted R-squared:  0.6581 
## F-statistic: 95.31 on 1 and 48 DF,  p-value: 5.531e-13
```
決定係数は0.66となり、少し弱い数値になりました。

フランス・イタリア・その他の国 で分けた形でも可視化してみましょう。

## フランスの散布図と回帰直線

フランスのデータで絞り、散布図と合わせて回帰線も引いてみます。

ggplot2では stat_smooth()で回帰線が設定可能で、method引数では回帰モデルが指定できます。ここでは線形回帰としてlmを指定しました。

```R
df_france <- df_fix %>% 
  filter(country == "France") 

plot_france <- ggplot(df_france, aes(x = vineyard_surface_area, y = wine_production)) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  ggtitle("France") +
  geom_text_repel(aes(label = region), size = 3) +
  geom_point() +
  stat_smooth(method = "lm", aes(x = vineyard_surface_area, y = wine_production))

plot_france
```

![Rplot3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/147689/7f7e3cdb-47ed-7954-9eae-00f6426d5e1e.png)

さきほどの散布図と比較するとばらつきが減り、直線的になったように見えます。

ラングドック・ルーション地方が回帰直線よりも上にあるのは（単位面積あたりの収量が多いのは）、I.G.P.ワインの生産量の多さが起因しているのかもしれません。
この地方はフランス全体のI.G.P.ワインの80%を生産しているためI.G.P.比率が高く（A.O.C.比率が低く）、他の地域よりも生産に関する制約が緩いものと思われます。

逆にボルドー地方は回帰直線よりも下にあり、これは品質向上のための収量制限が影響しているのかもしれません。
※ボルドーはフランスを代表する銘醸地で、A.O.C.ワインの生産量が国内第1位です。

ここでも決定係数を確認してみます。

```R
lm(wine_production ~ vineyard_surface_area, df_france) %>% 
  summary()

## 
## Call:
## lm(formula = wine_production ~ vineyard_surface_area, data = df_france)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1328673  -111553    78275   188482   808393 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -71260.421 224404.460  -0.318    0.757    
## vineyard_surface_area     54.719      3.123  17.521  2.2e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 593300 on 11 degrees of freedom
## Multiple R-squared:  0.9654, Adjusted R-squared:  0.9623 
## F-statistic:   307 on 1 and 11 DF,  p-value: 2.199e-09
```

フランスのモデルでは決定系数が0.96となり非常に当てはまりがよくなりました。
原産地呼称保護の下で、どの生産地も収量がコントロールされているのかもしれません。

vineyard_surface_areaのEstimate（係数）は54となっており、「面積が1増えるごとにワイン生産量が54増える（つまり単位面積あたりの収量が54）」と解釈できます。

ソムリエ・ワインエキスパートの試験ではブドウ栽培面積を「ボルドーが11万ha、ブルゴーニュが3万ha、、」といった形で覚えると思いますが、実はそれぞれの数値に54を掛ければワイン生産量となるわけです。
なのでフランスにおいてはワイン生産量を覚える必要が無くなりますね。

## イタリアの散布図と回帰直線

次にイタリアで見てみましょう。

```R
df_italia <- df_fix %>% 
  filter(country == "Italia") 

plot_italia <- ggplot(df_italia, aes(x = vineyard_surface_area, y = wine_production)) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  ggtitle("Italia") +
  geom_text_repel(aes(label = region), size = 3) +
  geom_point() +
  stat_smooth(method = "lm", aes(x = vineyard_surface_area, y = wine_production))

plot_italia
```

![Rplot4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/147689/d14b334e-4d88-4cfb-f85e-bcca42f7a1aa.png)

イタリアではバラつきが多くなっているようです。

ヴェネト、プーリア、エミリア・ロマーニャは回帰直線の上側に位置しています。

ヴェネトは国内第1位の生産量を誇り、プーリアはヨーロッパでも有数なバルクワイン供給地、エミリアではランブルスコ（赤の微発泡）を大量生産していた過去があり、これらの生産地は面積あたりの収量が多いのかもしれません。

フランスのラングドック・ルーション地方もそうですが、生産規模が大きくなると生産効率が重視され、結果として面積あたりの収量が増える傾向にありそうです。

逆にシチリア（島）は回帰直線の下にあるようで、これは上質なワインの多さを意味しているのでしょうか。気になります。

決定係数は0.70となり、フランスよりも低い数値となりました。

```R
lm(wine_production ~ vineyard_surface_area, df_italia) %>% 
  summary()

## 
## Call:
## lm(formula = wine_production ~ vineyard_surface_area, data = df_italia)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4609344  -395406   -97138   127908  3327058 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -157631.7   561180.1  -0.281    0.782    
## vineyard_surface_area      81.5       12.4   6.573 3.56e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1744000 on 18 degrees of freedom
## Multiple R-squared:  0.7059, Adjusted R-squared:  0.6895 
## F-statistic:  43.2 on 1 and 18 DF,  p-value: 3.564e-06
```


## その他の国の散布図と回帰直線

最後に、フランス・イタリア以外の国を確認します。

```R
df_other <- df_fix %>% 
  filter(country == "Other") 

plot_other <- ggplot(df_other, aes(x = vineyard_surface_area, y = wine_production)) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  ggtitle("Other Countries") +
  geom_text_repel(aes(label = region), size = 3) +
  geom_point()+
  stat_smooth(method = "lm", aes(x = vineyard_surface_area, y = wine_production))

plot_other
```

![Rplot5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/147689/8a1d4264-a6dc-f79b-3065-ebbb29392761.png)

南アフリカとドイツが回帰直線よりも上に位置していました。

南アフリカは生産量の半数が輸出、そのうち6割がバルクで、国際市場を意識した大量生産のワインが多いようです。
ナチュラル・ワインのカテゴリにも南アフリカのワインが多くみられますが、国全体としてはまだ「質より量」の状態から抜け出せていないのかもしれません。

ドイツはよくわかりませんが、ミュラー・トゥルガウのような収量の高いブドウ品種の比率が多いのでしょうか。


```R
lm(wine_production ~ vineyard_surface_area, df_other) %>% 
  summary()

## 
## Call:
## lm(formula = wine_production ~ vineyard_surface_area, data = df_other)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2138131  -532212  -382974    69236  4873095 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           4.784e+05  6.090e+05   0.786 0.444296    
## vineyard_surface_area 3.944e+01  7.987e+00   4.938 0.000179 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1838000 on 15 degrees of freedom
## Multiple R-squared:  0.6191, Adjusted R-squared:  0.5937 
## F-statistic: 24.38 on 1 and 15 DF,  p-value: 0.0001788
```

決定係数は0.61となりました。

## 散布図の比較
それぞれの散布図を並べて比較してみます。
gridExtraのgrid.arrange()では ggplotオブジェクトを並べて描画することができます。

```R
grid.arrange(plot_france, plot_italia, plot_other, ncol = 2)
```
![Rplot6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/147689/bd224f34-de81-cb65-dab1-e20708b69485.png)

フランスとイタリアを比較するとフランスの方が直線関係が強くなっており、国全体としてフランスの方が収量がコントロールされているのかもしれません。

ただの2変数での散布図ですが、何か示唆が得られそうな興味深い結果となりました。

## おわりに
散布図と回帰線を用いて各生産地を相対的に評価してみました。

ソムリエ・ワインエキスパートの勉強量はあまりにも多いので各生産地の横断的な比較がしづらいのですが、自分なりの軸でビジュアライズしてみるのは面白いかもしれません。

## 関連URL
https://www.sommelier.jp/exam/exam_guidance_detail2.html