Juliaで統計解析：手法によって結論が反対になる例

データ解析の基本は、まずよく目で見ることである。[Rをはじめよう 生命科学のためのRStudio入門](https://www.yodosha.co.jp/yodobook/book/9784758120951/index.html)のデータセット（[ここ](http://r4all.org/books/datasets/)にある）に含まれる compensation.csv を例に、どう見るかで結論が反対になる例を見てみる。

このデータは、二つの条件でリンゴの幹の直径と収穫量を調べたものである。つまり、どちらの条件で収穫量が多くなるかを知るためのデータである。直径も収穫量も実数なので、まずは観測値の分布を見てみる。

```Julia
julia> using CSV, Gadfly, DataFrames, DataFramesMeta, Statistics, HypothesisTests
julia> dat = CSV.read("compensation.csv", allowmissing=:none)
40×3 DataFrame
│ Row │ Root     │ Fruit    │ Grazing  │
│     │ Float64⍰ │ Float64⍰ │ String⍰  │
├─────┼──────────┼──────────┼──────────┤
│ 1   │ 6.225    │ 59.77    │ Ungrazed │
⋮
│ 40  │ 8.53     │ 83.03    │ Grazed   │

julia> describe(dat)
3×8 DataFrame. Omitted printing of 1 columns
│ Row │ variable │ mean    │ min    │ median │ max      │ nunique │ nmissing │
│     │ Symbol   │ Union…  │ Any    │ Union… │ Any      │ Union…  │ Int64    │
├─────┼──────────┼─────────┼────────┼────────┼──────────┼─────────┼──────────┤
│ 1   │ Root     │ 7.18115 │ 4.426  │ 7.1235 │ 10.253   │         │ 0        │
│ 2   │ Fruit    │ 59.4105 │ 14.73  │ 60.875 │ 116.05   │         │ 0        │
│ 3   │ Grazing  │         │ Grazed │        │ Ungrazed │ 2       │ 0        │

julia> plot(dat, x = :Grazing, y = :Fruit, Geom.boxplot)
```
![compensation_1.png](https://qiita-image-store.s3.amazonaws.com/0/377279/25e0be25-74c0-c2c6-ce42-3e7ed6a46a70.png)
収穫量（データ中のFruit列）は、Grazing列の水準がGrazedになっている方が多い。四分位数のうち最低値以外はGrazedの方が大きいので、平均値もGrazedの方が大きいだろう。観測値はRootとFruitの二つのなので、次はそれで散布図を見てみる。

```Julia
julia> plot(dat, x="Root", y="Fruit", color="Grazing", Geom.point)
```
![compensation_2.png](https://qiita-image-store.s3.amazonaws.com/0/377279/a8f59b94-dca3-115c-463f-f51b8e48d5a4.png)
GrazingがGrazedでもUngrazedのどちらでも、それぞれ直線を当てはめられそうな感じだが、どっちの直線が上かというと、Ungrazedの方である。つまり、同じRootで比較すれば収穫量が多いのはUngrazedなリンゴの木である。箱ヒゲ図か散布図かで、判断は逆になる。

もちろん、統計検定でも同じことがおこる。箱ヒゲ図の比較は、ｔ検定である。

```Julia
julia> datG = @linq dat |> where(:Grazing .== "Grazed")   |> select(:Fruit, :Grazing)
(出力略)
julia> datU = @linq dat |> where(:Grazing .== "Ungrazed") |> select(:Fruit, :Grazing)
(出力略)
julia> pvalue(UnequalVarianceTTest(datG[:, :Fruit], datU[:, :Fruit]))
0.026888248295954684
```

このp値は両側検定のものである。「Grazedの方が平均値が大きいか大きくないか」を判断するならp値はこの半分であり、いずれにせよ有意水準が５％なら十分に小さい。つまり、平均値はGrazedの方が有意に大きい。

散布図では直線をあてはめるが、二本の直線を比較するのは共分散分析である。

```Julia
julia> using GLM
julia> model = lm(@formula(Fruit ~ Root * Grazing), dat)
(出力一部略)
Coefficients:
                          Estimate Std.Error  t value Pr(>|t|)
(Intercept)               -125.173   12.8112 -9.77062   <1e-10
Root                       23.2404   1.53077  15.1821   <1e-16
Grazing: Ungrazed          30.8057   16.8418  1.82912   0.0757
Root & Grazing: Ungrazed  0.756034   2.35411 0.321155   0.7500
```

二つの直線の係数の違いは、一番下の二行に出ている。Grazing: Ungrazed の横が切片の違い、Root & Grazing: Ungrazed の横が傾きの違いで、それぞれp値は 0.0757 と 0.7500 である。前者は「差はあるかもしれないがハッキリしない」、後者は「差はなさそう」という感じである。つまり二本の直線は平行だろう、切片はUngrazedの方が大きいかも、というところで、数学に全然詳しくない人向けにあえて用語の不正確さを飲み込んで砕いて言うと、直径と収穫量は比例していて、Grazingがどっちでも収穫量の増え方は同じ、でも同じ太さの木ならUngrazedの方が収穫が多いかもしれない、となろう。だったら、ｔ検定の結論と逆である。

箱ヒゲ図と散布図と、どっちをベースに判断を下すかで結論は逆になる。実際にはもちろん、どっちの図も見て、また他にどんな図が描けるか考えて判断を進めるべきである。

この例の場合は、Rootの値の分布範囲がGrazedとUngrazedであっていないのがよくない。Grazedの方がRootが大きいので、Fruitも大きい。比較するなら条件を揃える、というのは鉄則である。しかし、条件を揃えられない、たとえばもしUngrazedはRootが大きくならないからデータがない、とでも言われれば、そのメカニズムが次の研究対象になる。