Julia 入門05 -- R と比較して（線形モデルによる回帰分析）

以下は，R 入門のために書かれた

https://ai-trend.jp/programming/r-beginner/r-beginner/

を，Julia ではこんな風になるよということで書いたものである。

R の入門というページは多く存在する。

R に入門しようとしている人 ... Julia にしませんか？

---


# 線形モデルによる回帰分析

## `llsq()` 関数を用いて単回帰分析

説明変数 `x` を用いて目的変数 `y`を予測する。

以下の例で，`x`, `y` は実数とするので，簡単のため 最初の要素だけ小数点を付けて定義する。<br>
なお，Julia では行末に `;` は不要であるが，REPL の場合に評価結果を表示させないために付けることがある。


```julia
x = [10., 20, 10, 40, 50, 10, 10, 20, 70];
y = [1., 3, 4, 10, 5, 1, 3, 14, 21];
```

Julia でも何通りかあるが，`MultivariateStats` にある `llsq()` を用いる。[^1]

[^1]: Linear Least Square Regression の略


```Julia
using MultivariateStats
```

また，基本的に `llsq()` は重回帰分析も行えるので，説明変数はたとえ 1 個（つまり単回帰分析）の場合であっても，n × 1 行列にしておかねばならない。<br>
ベクトル `x` を `reshape()` によって行列に変換する。


```julia
X = reshape(x, length(y), :);
```

`llsq()` はデフォルトで定数項を含めるモデル（`y = ax + b`）を採用する。[^2]

[^2]: 原点を通るモデルの場合は `bias=false` を指定する。


```julia
slope, intercept = llsq(X, y);
```

戻り値は，傾き(slope) と切片(intercept) である。<br>
上の例に示すように，左辺に `slope, intercept` と書いておけば，結果がそれぞれの変数に代入される。


```julia
println("slope = $slope, intercept = $intercept")
```

    slope = 0.23596491228070177, intercept = 0.5964912280701752


`llsq()` は傾きと切片の結果しか返さないので，偏回帰係数の有意性検定などをしたければ，自分でプログラムしなければならない。

## `GLM` による分析

`GLM` によれば，R と同じような感じで回帰分析を行える。


```julia
using GLM
```

ただし，データはデータフレームにしておく必要がある。


```julia
using DataFrames
df = DataFrame(x=x, y=y)
```




<table class="data-frame"><thead><tr><th></th><th>x</th><th>y</th></tr><tr><th></th><th>Float64</th><th>Float64</th></tr></thead><tbody><p>9 rows × 2 columns</p><tr><th>1</th><td>10.0</td><td>1.0</td></tr><tr><th>2</th><td>20.0</td><td>3.0</td></tr><tr><th>3</th><td>10.0</td><td>4.0</td></tr><tr><th>4</th><td>40.0</td><td>10.0</td></tr><tr><th>5</th><td>50.0</td><td>5.0</td></tr><tr><th>6</th><td>10.0</td><td>1.0</td></tr><tr><th>7</th><td>10.0</td><td>3.0</td></tr><tr><th>8</th><td>20.0</td><td>14.0</td></tr><tr><th>9</th><td>70.0</td><td>21.0</td></tr></tbody></table>



R の `lm()` 関数と同じような書き方で分析できる。


```julia
results = lm(@formula(y ~ x), df);
println(results)
```

    StatsModels.DataFrameRegressionModel{LinearModel{LmResp{Vector{Float64}}, DensePredChol{Float64, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}}}, Matrix{Float64}}
    
    Formula: y ~ 1 + x
    
    Coefficients:
    ────────────────────────────────────────────────────
                 Estimate  Std.Error   t value  Pr(>|t|)
    ────────────────────────────────────────────────────
    (Intercept)  0.596491  2.60529    0.228954    0.8254
    x            0.235965  0.0773885  3.04909     0.0186
    ────────────────────────────────────────────────────


結果も同じように表示される。<br>
多分，R の関数を呼んでいるのだろう。
偏回帰係数（定数項），それらの標準誤差，有意性検定の t 値および P 値が表示される。

個別の統計量は以下のようにして求めることができる。<br>
偏回帰係数


```julia
coef(results)
```




    2-element Vector{Float64}:
     0.5964912280701711
     0.23596491228070188



モデルの逸脱度（deviance)


```julia
deviance(results)
```




    159.30701754385967



残差の自由度


```julia
dof_residual(results)
```




    7.0



モデルの各係数の標準誤差


```julia
stderror(results)
```




    2-element Vector{Float64}:
     2.6052863639549413
     0.07738853676061858



係数の共分散行列


```julia
vcov(results)
```




    2×2 Matrix{Float64}:
      6.78752   -0.159706
     -0.159706   0.00598899



予測値


```julia
predict(results)
```




    9-element Vector{Float64}:
      2.95614035087719
      5.315789473684209
      2.95614035087719
     10.035087719298247
     12.394736842105265
      2.95614035087719
      2.95614035087719
      5.315789473684209
     17.1140350877193



実測値と予測値の図を描いてみよう


```julia
using Plots
scatter(x, y, alpha=0.2, label="observed")
scatter!(x, predict(results), alpha=0.2, label="predicted")
x2 = minimum(x):0.1:maximum(x)
predicted = coef(results)[2] .* x2 .+ coef(results)[1]
plot!(x2, predicted, label="predicted", legend=:bottomright)
```

![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/e243ea16-fa9e-147f-4bac-bc6e246722e9.png)




## 重回帰分析

ここでは，例示のためのデータセットが使われているので，まずは下記の URL から，データをダウンロードせよと言われている。
https://ai-trend.jp/wp-content/uploads/2016/07/sample-data.csv
Mac OS の場合には `/Users/foo/Dwnloads/sample-data.csv` にダウンロードされる。

このファイルは邪悪な Window ファイル(cp932 エンコーディング)なので（Mac OS では）普通の `CSV.read()` では読めない。<br>
これを読み取るために以下のように指定する。


```julia
using DataFrames, CSV
using StringEncodings
df = CSV.File(open(read, "/Users/foo/Downloads/sample-data.csv", enc"cp932")) |> DataFrame
```




<table class="data-frame"><thead><tr><th></th><th>番号</th><th>年齢</th><th>血圧</th><th>肺活量</th><th>性別</th><th>病気</th><th>体重</th></tr><tr><th></th><th>Int64</th><th>Int64</th><th>Int64</th><th>Int64</th><th>String</th><th>Int64</th><th>Int64</th></tr></thead><tbody><p>30 rows × 7 columns</p><tr><th>1</th><td>1</td><td>22</td><td>110</td><td>4300</td><td>M</td><td>1</td><td>79</td></tr><tr><th>2</th><td>2</td><td>23</td><td>128</td><td>4500</td><td>M</td><td>1</td><td>65</td></tr><tr><th>3</th><td>3</td><td>24</td><td>104</td><td>3900</td><td>F</td><td>0</td><td>53</td></tr><tr><th>4</th><td>4</td><td>25</td><td>112</td><td>3000</td><td>F</td><td>0</td><td>45</td></tr><tr><th>5</th><td>5</td><td>27</td><td>108</td><td>4800</td><td>M</td><td>0</td><td>80</td></tr><tr><th>6</th><td>6</td><td>28</td><td>126</td><td>3800</td><td>F</td><td>0</td><td>50</td></tr><tr><th>7</th><td>7</td><td>28</td><td>126</td><td>3800</td><td>F</td><td>1</td><td>43</td></tr><tr><th>8</th><td>8</td><td>29</td><td>104</td><td>4000</td><td>F</td><td>1</td><td>55</td></tr><tr><th>9</th><td>9</td><td>30</td><td>125</td><td>3600</td><td>F</td><td>1</td><td>47</td></tr><tr><th>10</th><td>10</td><td>31</td><td>120</td><td>3400</td><td>F</td><td>1</td><td>49</td></tr><tr><th>11</th><td>11</td><td>32</td><td>116</td><td>3600</td><td>M</td><td>1</td><td>64</td></tr><tr><th>12</th><td>12</td><td>32</td><td>124</td><td>3900</td><td>M</td><td>0</td><td>61</td></tr><tr><th>13</th><td>13</td><td>33</td><td>106</td><td>3100</td><td>F</td><td>0</td><td>48</td></tr><tr><th>14</th><td>14</td><td>33</td><td>134</td><td>2900</td><td>F</td><td>0</td><td>41</td></tr><tr><th>15</th><td>15</td><td>34</td><td>128</td><td>4100</td><td>M</td><td>1</td><td>70</td></tr><tr><th>16</th><td>16</td><td>36</td><td>128</td><td>3420</td><td>M</td><td>1</td><td>55</td></tr><tr><th>17</th><td>17</td><td>37</td><td>116</td><td>3800</td><td>M</td><td>1</td><td>70</td></tr><tr><th>18</th><td>18</td><td>37</td><td>132</td><td>4150</td><td>M</td><td>1</td><td>90</td></tr><tr><th>19</th><td>19</td><td>38</td><td>134</td><td>2700</td><td>F</td><td>0</td><td>39</td></tr><tr><th>20</th><td>20</td><td>39</td><td>116</td><td>4550</td><td>M</td><td>1</td><td>86</td></tr><tr><th>21</th><td>21</td><td>40</td><td>120</td><td>2900</td><td>F</td><td>1</td><td>50</td></tr><tr><th>22</th><td>22</td><td>42</td><td>130</td><td>3950</td><td>F</td><td>1</td><td>65</td></tr><tr><th>23</th><td>23</td><td>46</td><td>126</td><td>3100</td><td>M</td><td>0</td><td>58</td></tr><tr><th>24</th><td>24</td><td>49</td><td>140</td><td>3000</td><td>F</td><td>0</td><td>45</td></tr><tr><th>25</th><td>25</td><td>50</td><td>156</td><td>3400</td><td>M</td><td>1</td><td>60</td></tr><tr><th>26</th><td>26</td><td>53</td><td>124</td><td>3400</td><td>M</td><td>1</td><td>71</td></tr><tr><th>27</th><td>27</td><td>56</td><td>118</td><td>3470</td><td>M</td><td>1</td><td>62</td></tr><tr><th>28</th><td>28</td><td>58</td><td>144</td><td>2800</td><td>M</td><td>0</td><td>51</td></tr><tr><th>29</th><td>29</td><td>64</td><td>142</td><td>2500</td><td>F</td><td>1</td><td>40</td></tr><tr><th>30</th><td>30</td><td>65</td><td>144</td><td>2350</td><td>F</td><td>0</td><td>42</td></tr></tbody></table>



このデータファイルを使って，肺活量を血圧と体重で予測する。

指定方法は単回帰分析の場合と同じである。


```julia
using GLM
results = lm(@formula(肺活量 ~ 血圧 + 体重), df);
println(results)
```

    StatsModels.DataFrameRegressionModel{LinearModel{LmResp{Vector{Float64}}, DensePredChol{Float64, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}}}, Matrix{Float64}}
    
    Formula: 肺活量 ~ 1 + 血圧 + 体重
    
    Coefficients:
    ─────────────────────────────────────────────────────
                  Estimate  Std.Error   t value  Pr(>|t|)
    ─────────────────────────────────────────────────────
    (Intercept)  3533.07    830.819     4.25251    0.0002
    血圧            -14.2164    5.57863  -2.54837    0.0168
    体重             30.7853    5.09687   6.04005    <1e-05
    ─────────────────────────────────────────────────────


`println()` で表示される総括表は `coeftable()` で表示されるものと同じである。


```julia
coeftable(results)
```




    ─────────────────────────────────────────────────────
                  Estimate  Std.Error   t value  Pr(>|t|)
    ─────────────────────────────────────────────────────
    (Intercept)  3533.07    830.819     4.25251    0.0002
    血圧            -14.2164    5.57863  -2.54837    0.0168
    体重             30.7853    5.09687   6.04005    <1e-05
    ─────────────────────────────────────────────────────



係数だけを取り出すのは `coef()` である。


```julia
coef(results)
```




    3-element Vector{Float64}:
     3533.065265426281
      -14.216447838209147
       30.785336447492494



変数名を取り出すのは `coefnames()` である。


```julia
coefnames(results)
```




    3-element Vector{String}:
     "(Intercept)"
     "血圧"
     "体重"



偏回帰係数（定数項）標準誤差は `stderror()` である。


```julia
stderror(results)
```




    3-element Vector{Float64}:
     830.8185624272833
       5.578634802056308
       5.096867359589727



R の `summary()` で出力されるもののうち，まず残差標準誤差 Residual standard error は以下のようにして求める。<br>
残差の自由度は `dof_residual(results)` である。


```julia
sqrt(sum(residuals(results) .^ 2) / dof_residual(results))
```




    368.26786392925806




```julia
dof_residual(results) # 残差の自由度
```




    27.0



重相関係数の二乗（決定係数），および自由度調整済みの重相関係数の二乗は以下のようにして求める。


```julia
r2(results) # 重相関係数の二乗（決定係数）
```




    0.6748009330930996




```julia
adjr2(results) # 重相関係数の二乗（決定係数）
```




    0.6507121133222181



AIC，AICC，BIC は `aic()`，`aicc()`，`bic()` で求める。`StatsBase` が必要である。


```julia
using StatsBase
aic(results)
```




    444.5041303923531




```julia
aicc(results)
```




    446.1041303923531




```julia
bic(results)
```




    450.1089199190017



そのほか，有効データ数は `nobs()`, 偏回帰係数の信頼区間は `confint()` , 自由度は `dof`，分散・共分散行列は `vcov()`，デビアンスは `deviance()`，対数尤度は `loglikelihood()`，予測値は `predict()`，残差は `rediduals()`，残差の自由度は `dof_residual()` で求める。


```julia
nobs(results) # 有効データ数
```




    30.0




```julia
confint(results) # 偏回帰係数（定数項）の信頼区間
```




    3×2 Matrix{Float64}:
     1828.37    5237.76
      -25.6629    -2.77003
       20.3274    41.2432




```julia
stderror(results) # 偏回帰係数（定数項）の標準誤差
```




    3-element Vector{Float64}:
     830.8185624272833
       5.578634802056308
       5.096867359589727



dof(results) # モデルで消費される自由度（説明変数の個数 + 1 + 1）


```julia
vcov(results) # 偏回帰係数の分散・共分散行列
```




    3×3 Matrix{Float64}:
         6.90259e5  -4341.97     -2496.46
     -4341.97          31.1212       7.97855
     -2496.46           7.97855     25.9781




```julia
deviance(results) # デビアンス
```




    3.6617729292815e6




```julia
nulldeviance(results) # null モデルのデビアンス
```




    1.126009666666667e7




```julia
loglikelihood(results) # 対数尤度
```




    -218.25206519617655




```julia
nullloglikelihood(results) # null モデルの対数尤度
```




    -235.10183175521513




```julia
predict(results) # 予測値
```




    30-element Vector{Float64}:
     4401.297582575182
     3714.4068112225223
     3686.1775219696315
     3326.1632476840186
     4460.5158146990925
     3281.0596601865527
     3065.5623050541053
     3747.7481948646164
     3202.9200986822843
     3335.5730107683153
     3854.2188488335396
     3648.1312567853884
     3503.8179440557506
        ⋮
     2828.6893765584628
     4531.496250678374
     3366.358347215808
     3685.973915546104
     3527.342351766493
     2928.1027082141627
     3162.4195895152034
     3955.984621260314
     3764.215280262136
     3055.948935546281
     2745.743130300282
     2778.8809075188483




```julia
residuals(results) # 残差
```




    30-element Vector{Float64}:
     -101.29758257518188
      785.5931887774777
      213.82247803036853
     -326.16324768401864
      339.4841853009075
      518.9403398134473
      734.4376949458947
      252.25180513538362
      397.07990131771567
       64.42698923168473
     -254.2188488335396
      251.86874321461164
     -403.81794405575056
        ⋮
     -128.68937655846275
       18.503749321625946
     -466.35834721580795
      264.02608445389615
     -427.3423517664928
       71.8972917858373
      237.58041048479663
     -555.9846212603138
     -294.2152802621358
     -255.94893554628106
     -245.7431303002818
     -428.8809075188483




```julia
scatter(df.肺活量, predict(results), smooth=true, xlabel="Lung capacity", ylabel="Predicted lung capacity", label="")
```




![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/0c959af1-13b1-f8a0-a8ec-9316772279df.png)