GLM パッケージによる定数項を含まない重回帰分析

# GLM パッケージによる定数項を含まない重回帰分析の結果に一部不備がある

以下のデータを用いて検証する。


```julia
data = [25.43  78.11  34.35   98.88  53.28
        18.56  76.34  35.76  100.48  56.91
        28.66  79.05  44.00  103.21  54.52
        13.52  72.24  35.26   97.97  50.13
        30.52  77.86  41.97   97.32  59.89
        30.16  77.88  50.14   98.35  53.43
        29.30  79.26  35.15   98.55  57.35
        31.28  79.38  49.57   97.74  55.13
        31.46  79.14  46.80  101.98  53.71
        29.21  77.63  37.07  102.35  36.10
        22.70  76.80  37.95   88.83  47.18
        34.47  78.88  53.90  100.21  59.90
        37.38  78.90  54.16  104.02  67.05
        13.38  73.97  42.15   92.01  42.44
        16.98  75.59  30.25   92.59  32.98];
using DataFrames
df = DataFrame(data, [:x1, :x2, :x3, :x4, :y]);
first(df, 2)
```




<table class="data-frame"><thead><tr><th></th><th>x1</th><th>x2</th><th>x3</th><th>x4</th><th>y</th></tr><tr><th></th><th>Float64</th><th>Float64</th><th>Float64</th><th>Float64</th><th>Float64</th></tr></thead><tbody><p>2 rows × 5 columns</p><tr><th>1</th><td>25.43</td><td>78.11</td><td>34.35</td><td>98.88</td><td>53.28</td></tr><tr><th>2</th><td>18.56</td><td>76.34</td><td>35.76</td><td>100.48</td><td>56.91</td></tr></tbody></table>



## Julia の GLM.lm( ) では，定数項を含まない場合の決定係数，調整済み決定係数が普通の定義（たとえば R での結果）とは異なる結果となる。

### Julia での分析結果


```julia
using GLM
result2 = lm(@formula(y ~ 0 + x1 + x2 + x3 + x4), df);
println(result2)
```

    StatsModels.DataFrameRegressionModel{LinearModel{LmResp{Vector{Float64}}, DensePredChol{Float64, Cholesky{Float64, Matrix{Float64}}}}, Matrix{Float64}}
    
    Formula: y ~ x1 + x2 + x3 + x4
    
    Coefficients:
    ─────────────────────────────────────────────
         Estimate  Std.Error    t value  Pr(>|t|)
    ─────────────────────────────────────────────
    x1   0.338163   0.377074   0.896807    0.3890
    x2  -0.216267   0.687945  -0.314368    0.7591
    x3   0.412138   0.367005   1.12298     0.2854
    x4   0.433544   0.54367    0.797439    0.4421
    ─────────────────────────────────────────────



```julia
println("R2      = $(r2(result2))")
println("R2Adj.  = $(adjr2(result2))")
```

    R2      = 0.4549168173881398
    R2Adj.  = 0.3062577675849053


### R での実行結果


```julia
using RCall
R"""
result2 = lm(y ~ 0 + x1 + x2 + x3 + x4, $df)
summary2 = summary(result2)
print(summary2, digits=6)
""";

R"""
options(digits=16)
cat(sprintf("R2     = %.15g\n", summary2$r.squared))
cat(sprintf("R2Adj. = %.15g\n", summary2$adj.r.squared))
""";
```

    
    Call:
    lm(formula = y ~ 0 + x1 + x2 + x3 + x4, data = `#JL`$df)
    
    Residuals:
           Min         1Q     Median         3Q        Max 
    -16.640026  -3.271851  -0.356926   4.361203   8.843050 
    
    Coefficients:
        Estimate Std. Error  t value Pr(>|t|)
    x1  0.338163   0.377074  0.89681  0.38903
    x2 -0.216267   0.687945 -0.31437  0.75912
    x3  0.412138   0.367005  1.12298  0.28536
    x4  0.433544   0.543670  0.79744  0.44207
    
    Residual standard error: 7.57951 on 11 degrees of freedom
    Multiple R-squared:  0.984853,	Adjusted R-squared:  0.979344 
    F-statistic:   178.8 on 4 and 11 DF,  p-value: 6.29747e-10
    
    R2     = 0.984852626318541
    R2Adj. = 0.979344490434375


R での結果において，決定係数は 0.9848526263185413，自由度調整済み決定係数は 0.9793444904343745 である。<br>
Julia の GLM.ml() ではそれぞれ，0.4549168173881398，0.3062577675849053 と大きく異なる。

### Julia での R2, R2adj の計算法

結果が違う原因は，Julia の GLM.ml() での決定係数と自由度調整済み決定係数は定数項を含む場合も含まない場合も同じ計算法（平均値を中心とする平方和に基づいている）をとっていることである。<br>[^1]

[^1]: 定数項を使う場合は $(n - p)$ の項は $(n - p - 1)$ になる。


$$
R^2 = 1 - \frac{S_e}{S_t}, \ \ \ \ 
R^{2}_{adj} = 1 - \frac{S_e / (n - p)}{S_t / (n - 1)}
$$


ここで，$n$ はサンプルサイズ，$p$ は独立変数の個数，$S_t$ は従属変数 $y$ の平方和，$S_e$ は残差平方和である。
$$
S_t = \displaystyle \sum_{i=1}^n \left (y_i-\bar{y}\right )^2, \ \ \ \ 
S_e = \displaystyle \sum_{i=1}^n \left (y_i-\hat{y} \right )^2
$$

計算をしてみると確かに GLM.lm( ) と同じ結果になる。


```julia
n, p = size(df)
p -= 1
println("n = $n, p = $p")
St = var(df.y) * (n - 1)
Se = sum(residuals(result2) .^ 2)
R² = 1 - Se / St
R²adj = 1 - (Se / (n - p)) / (St / (n - 1))
println("R²     = $R²")
println("R²adj  = $R²adj")
```

    n = 15, p = 4
    R²     = 0.4549168173881398
    R²adj  = 0.3062577675849052


### R などでの R2, R2adj の計算法

原点を中心とした平方和に基づく。

$$
R^2 = \frac{S_{predict}}{S_{observe}}, \ \ \ \ 
R^{2}_{adj} = 1 - \left (1 - R^2 \right )\frac{n}{n - p}
$$


$$
S_{observe} = \sum_{i=1}^n y_i^2, \ \ \ \ 
S_{predict} = \sum_{i=1}^n \hat{y}_i^2
$$


```julia
Spredict = sum(predict(result2) .^ 2)
Sobserve = sum(df.y .^ 2)
R2 = Spredict / Sobserve
R2adj = 1 - (1 - R2) * (n / (n - p))
println("R2    = $R2")
println("R2adj = $R2adj")
```

    R2    = 0.9848526263185415
    R2adj = 0.9793444904343748