トップ 履歴 一覧 Farm ソース 検索 ヘルプ PDF RSS ログイン

Multiple Regression Analysis

*disclaimer
793841

R

重回帰分析


 サンプルデータとスクリプト

  • サンプルデータ:NICER 1.0より学習者のエッセイデータ287個
  • Criterionによるスコア、延べ語数、異なり語数、文数、TTR、ギローインデックス、平均単語長、平均文長、MATTR を各エッセイの特徴量として抽出
  • ファイル名とともに結果をテキストファイルに保存するスクリプト

myIndices4.R(342)

  1. NICER1.0を解凍したフォルダーの中のNICER_NNSに作業ディレクトリーを移動
  2. スクリプトを実行
  3. 結果を保存するファイルを作業ディレクトリーの外にファイル名をつけて「保存(作成)」例:jpn4.txt
  4. Rの中に読み込む
> jpn4 <- read.table(choose.files())
> class(jpn4)
[1] "data.frame"
> head(jpn4)
          V1 V2  V3  V4 V5        V6       V7        V8       V9      V10
1 JPN501.txt  4 319 135 30 0.4231975 7.558549 0.5921317 4.304075 10.63333
2 JPN502.txt  4 356 161 29 0.4522472 8.532983 0.6649157 4.233146 12.27586
3 JPN503.txt  3 201 121 13 0.6019900 8.534682 0.7170149 4.746269 15.46154
4 JPN504.txt  4 260 140 27 0.5384615 8.682431 0.6877692 4.761538  9.62963
5 JPN505.txt  4 420 175 25 0.4166667 8.539126 0.6341905 3.995238 16.80000
6 JPN506.txt  3 261 124 20 0.4750958 7.675407 0.6390038 4.072797 13.05000

 前処理

見出しをつける

> names(jpn4) <- c("file", "Score", "Token", "Type", "NoS", "TTR", "GI", "MATTR", "AWL", "ASL")
> head(jpn4)
        file Score Token Type NoS       TTR       GI     MATTR      AWL      ASL
1 JPN501.txt     4   319  135  30 0.4231975 7.558549 0.5921317 4.304075 10.63333
2 JPN502.txt     4   356  161  29 0.4522472 8.532983 0.6649157 4.233146 12.27586
3 JPN503.txt     3   201  121  13 0.6019900 8.534682 0.7170149 4.746269 15.46154
4 JPN504.txt     4   260  140  27 0.5384615 8.682431 0.6877692 4.761538  9.62963
5 JPN505.txt     4   420  175  25 0.4166667 8.539126 0.6341905 3.995238 16.80000
6 JPN506.txt     3   261  124  20 0.4750958 7.675407 0.6390038 4.072797 13.05000

データ数の確認

> nrow(jpn4)
[1] 287

> str(jpn4)
'data.frame':	287 obs. of  10 variables:
 $ file : Factor w/ 287 levels "JPN501.txt","JPN502.txt",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Score: int  4 4 3 4 4 3 4 3 4 3 ...
 $ Token: int  319 356 201 260 420 261 362 198 263 183 ...
 $ Type : int  135 161 121 140 175 124 151 98 104 99 ...
 $ NoS  : int  30 29 13 27 25 20 26 20 19 14 ...
 $ TTR  : num  0.423 0.452 0.602 0.538 0.417 ...
 $ GI   : num  7.56 8.53 8.53 8.68 8.54 ...
 $ MATTR: num  0.592 0.665 0.717 0.688 0.634 ...
 $ AWL  : num  4.3 4.23 4.75 4.76 4 ...
 $ ASL  : num  10.63 12.28 15.46 9.63 16.8 ...

欠損値の処理

  • 欠損値を調べる 欠損値があると計算できないことがあるので
> anyNA(jpn4)
[1] TRUE

> is.na(jpn4)
> sum(is.na(jpn4))
[1] 2
> jpn4.b <- na.omit(jpn4)
[1] 285

分析対象のオブジェクトattach()で「固定」しておくと便利

  • jpn4.bについて、各カラムを単位に操作する
> attach(jpn4.b)
  • こうしておけば、いちいち jpn4.b$Score としなくても Score だけでよい。
  • 終わったら、detach()するのを忘れずに

 重回帰分析 lm()

lm(従属変数 ~ 独立変数1 + 独立変数2 + 独立変数3 + ... )

> lm(Score ~ Token + Type + NoS + TTR + GI + MATTR + AWL + ASL)
> lm(Score ~ Token + Type + NoS + TTR + GI + MATTR + AWL + ASL, data=jpn4.b)

  • データの中の一つを従属変数にして、残りはすべて独立変数として全部入れるときは、省略して書ける
lm(従属変数 ~ ., data=データ名)

Call:
lm(formula = Score ~ Token + Type + NoS + TTR + GI + MATTR + AWL + ASL)

Coefficients:
(Intercept)        Token         Type          NoS          TTR           GI        MATTR          AWL          ASL  
  -0.145765     0.004008    -0.021620     0.003506    -9.147468     1.004200    -1.306580     0.556470     0.021069  

  • 結果を詳しく見るには、一旦保存して、そのsummary()を見る。
> jukaiki <- lm(Score ~ Token + Type + NoS + TTR + GI + MATTR + AWL + ASL)
> summary(jukaiki)

Call:
lm(formula = Score ~ Token + Type + NoS + TTR + GI + MATTR + 
    AWL + ASL)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.9080 -0.2508  0.0035  0.2581  0.7614 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.145765   0.642548  -0.227 0.820705    
Token        0.004008   0.002144   1.870 0.062603 .  
Type        -0.021620   0.010421  -2.075 0.038953 *  
NoS          0.003506   0.014089   0.249 0.803665    
TTR         -9.147468   1.328107  -6.888 3.83e-11 ***
GI           1.004200   0.257697   3.897 0.000122 ***
MATTR       -1.306580   1.362782  -0.959 0.338519    
AWL          0.556470   0.064559   8.620 5.32e-16 ***
ASL          0.021069   0.025071   0.840 0.401424    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.334 on 276 degrees of freedom
Multiple R-squared:  0.8145,    Adjusted R-squared:  0.8091 
F-statistic: 151.5 on 8 and 276 DF,  p-value: < 2.2e-16


  • Estimateが各変数の係数
  • t値とその右のp値は、「それぞれの変数の係数が0(ゼロ)、つまり、その変数の影響はない」という帰無仮説が起きる確率
    • 「影響がない」確率が5%とか1%未満で低ければ、そういうことは起きえない、つまり「影響がないとは言えない」という話になる。
    • というわけで、有意(星がついている)であれば、その変数は影響力がある、と考えられえる

標準化

  • 変数の規模・単位が違うので、影響を等しく比べられるように「標準化(正規化)」する。

https://scrapbox.io/nishio/%E6%AD%A3%E8%A6%8F%E5%8C%96%E3%81%A8%E6%A8%99%E6%BA%96%E5%8C%96
関数は scale()
ただし、データは数値でないといけないので、まず、一番左にあるファイル名のカラムを除く。

> jpn4.c <- jpn4.b[,-1]
  • 様子を見てみる
> pairs(jpn4.c)
  • 標準化する
> jpn4.z <- scale(jpn4.c)
  • 型がmatrixになっているので、データフレームに変換
> jpn4.z <- as.data.frame(jpn4.z)
  • これをもとに、改めて、重回帰分析
> attach(jpn4.z)
> jukaiki.z <- lm(Score ~ Token + Type + NoS + TTR + GI + MATTR + AWL + ASL)
> summary(jukaiki.z) 

Call:
lm(formula = Score ~ Token + Type + NoS + TTR + GI + MATTR + 
    AWL + ASL)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.18748 -0.32797  0.00458  0.33760  0.99580 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.587e-17  2.588e-02   0.000 1.000000    
Token        5.179e-01  2.770e-01   1.870 0.062603 .  
Type        -9.364e-01  4.514e-01  -2.075 0.038953 *  
NoS          3.299e-02  1.326e-01   0.249 0.803665    
TTR         -7.898e-01  1.147e-01  -6.888 3.83e-11 ***
GI           1.282e+00  3.290e-01   3.897 0.000122 ***
MATTR       -7.775e-02  8.110e-02  -0.959 0.338519    
AWL          2.423e-01  2.811e-02   8.620 5.32e-16 ***
ASL          7.988e-02  9.505e-02   0.840 0.401424    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4369 on 276 degrees of freedom
Multiple R-squared:  0.8145,    Adjusted R-squared:  0.8091 
F-statistic: 151.5 on 8 and 276 DF,  p-value: < 2.2e-16

VIF

  • ところが、独立変数間でお互いに相関が高いものがある。
    • 複数の相関が高いものを重ねて入れてしまうと、その影響が重複して評価されてしまう。
    • VIF(Variance Inflation factor)をチェック(10以上はダメ、ほとんど同じもの)
    • パッケージcarをインストールする
> library(car)
  • 関数 vif()の実行
> vif(jukaiki.z)
     Token       Type        NoS        TTR         GI      MATTR        AWL        ASL 
114.203516 303.150320  26.150429  19.565310 161.087442   9.787097   1.175870  13.443703 

10以上は外すべき。2以下が望ましい。

  • しかし、一つをはずすと、それに関連したものが変動して、値も変わるので、一概に10以上はすべてダメというわけではない。何を選ぶかは、判断による。で、やってみる。

  • VIFの計算は自分でやってもよい (VIFのページ参照)

VIF = 1/(1 - 相関係数^2)

変数選択 関数 step()

> jukaiki.z.result <- step(jukaiki.z)
> summary(jukaiki.z.result)

Call:
lm(formula = Score ~ Token + Type + TTR + GI + AWL + ASL)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.17036 -0.34503 -0.00159  0.34247  0.98419 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.508e-17  2.583e-02   0.000   1.0000    
Token        4.732e-01  2.183e-01   2.167   0.0311 *  
Type        -7.950e-01  4.082e-01  -1.948   0.0525 .  
TTR         -8.081e-01  1.057e-01  -7.648 3.36e-13 ***
GI           1.145e+00  2.709e-01   4.227 3.22e-05 ***
AWL          2.389e-01  2.768e-02   8.630 4.83e-16 ***
ASL          5.048e-02  2.918e-02   1.730   0.0848 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4361 on 278 degrees of freedom
Multiple R-squared:  0.8138,    Adjusted R-squared:  0.8098 
F-statistic: 202.6 on 6 and 278 DF,  p-value: < 2.2e-16
  • Token + Type + TTR + GI + AWL + ASL が変数の候補として残った。
    • ただし、TypeとASLは重要度が低い
  • ここで、また、VIFを見てみる。
> vif(jukaiki.z.result)
     Token       Type        TTR         GI        AWL        ASL 
 71.198905 248.805646  16.675124 109.594573   1.144067   1.271753 
  • 欠損値のためエラーが起きる場合は、欠損値を除いておく na.omit(データフレーム)

試行錯誤

  • TypeとASLを除いて、また、step()してみて、VIFを見てみる、、、という風に試行錯誤していく。
> jukaiki.z2 <- lm(Score ~ Token + TTR + GI + AWL)
> vif(jukaiki.z2)
    Token       TTR        GI       AWL 
12.765221 10.876050  8.233309  1.081001 
> jukaiki.z2.result <- step(jukaiki.z2)
> summary(jukaiki.z2.result)

Call:
lm(formula = Score ~ TTR + GI + AWL)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.24601 -0.34767  0.02472  0.32806  1.00389 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.446e-16  2.608e-02   0.000        1    
TTR         -7.848e-01  2.817e-02 -27.861   <2e-16 ***
GI           7.205e-01  2.752e-02  26.178   <2e-16 ***
AWL          2.539e-01  2.716e-02   9.346   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4403 on 281 degrees of freedom
Multiple R-squared:  0.8082,    Adjusted R-squared:  0.8061 
F-statistic: 394.6 on 3 and 281 DF,  p-value: < 2.2e-16
  • TTRとGIとAWLが残る。
> jukaiki.z3 <- lm(Score ~ TTR + GI + AWL)
> vif(jukaiki.z3)
     TTR       GI      AWL 
1.162202 1.109615 1.080821 
  • これで、VIFも問題ない。
  • R2も0.8とかなり高い。(8割のデータがこのモデルで説明できる)

結論

エッセイのスコアはTTRとGIとAWLで、8割が説明できる。

興味深い点

  • TTRもGIも語彙の多様性を表す指標といわれているのに、、、
  • 相関を見てみる
> cor.test(TTR, GI)

        Pearson's product-moment correlation

data:  TTR and GI
t = 5.3729, df = 283, p-value = 1.62e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1949506 0.4060786
sample estimates:
      cor 
0.3042462 
  • 相関は弱い。

もう一度やってみる。

  • 語彙の多様性指標は、MATTRにしてみる。
  • Tokenはいろいろと相関が出るので外してみる。
> head(jpn4.b)
        file Score Token Type NoS       TTR       GI     MATTR      AWL      ASL
1 JPN501.txt     4   319  135  30 0.4231975 7.558549 0.5921317 4.304075 10.63333
2 JPN502.txt     4   356  161  29 0.4522472 8.532983 0.6649157 4.233146 12.27586
3 JPN503.txt     3   201  121  13 0.6019900 8.534682 0.7170149 4.746269 15.46154
4 JPN504.txt     4   260  140  27 0.5384615 8.682431 0.6877692 4.761538  9.62963
5 JPN505.txt     4   420  175  25 0.4166667 8.539126 0.6341905 3.995238 16.80000
6 JPN506.txt     3   261  124  20 0.4750958 7.675407 0.6390038 4.072797 13.05000
> cor(jpn4.b[3:10])
           Token         Type         NoS        TTR        GI      MATTR          AWL        ASL
Token  1.0000000  0.889026009  0.75740319 -0.6440682 0.4822118 0.12764105 -0.108067282  0.4091337
Type   0.8890260  1.000000000  0.67377228 -0.2703160 0.8237958 0.51028459  0.002673065  0.3775746
NoS    0.7574032  0.673772281  1.00000000 -0.5133734 0.3599018 0.02273514 -0.209605424 -0.2531863
TTR   -0.6440682 -0.270316024 -0.51337343  1.0000000 0.3042462 0.59836349  0.261516052 -0.2370583
GI     0.4822118  0.823795754  0.35990182  0.3042462 1.0000000 0.86155201  0.155691920  0.2309099
MATTR  0.1276411  0.510284594  0.02273514  0.5983635 0.8615520 1.00000000  0.286936810  0.1777219
AWL   -0.1080673  0.002673065 -0.20960542  0.2615161 0.1556919 0.28693681  1.000000000  0.1645096
ASL    0.4091337  0.377574643 -0.25318632 -0.2370583 0.2309099 0.17772193  0.164509608  1.0000000
> jpn4.b.cor <- cor(jpn4.b[3:10])
> jpn4.b.vif <- 1/(1-jpn4.b.cor^2)
> jpn4.b.vif
         Token     Type      NoS      TTR       GI    MATTR      AWL      ASL
Token      Inf 4.770247 2.345544 1.708887 1.302979 1.016562 1.011817 1.201043
Type  4.770247      Inf 1.831398 1.078831 3.111770 1.352065 1.000007 1.166266
NoS   2.345544 1.831398      Inf 1.357870 1.148804 1.000517 1.045953 1.068494
TTR   1.708887 1.078831 1.357870      Inf 1.102008 1.557727 1.073411 1.059543
GI    1.302979 3.111770 1.148804 1.102008      Inf 3.880058 1.024842 1.056322
MATTR 1.016562 1.352065 1.000517 1.557727 3.880058      Inf 1.089720 1.032615
AWL   1.011817 1.000007 1.045953 1.073411 1.024842 1.089720      Inf 1.027816
ASL   1.201043 1.166266 1.068494 1.059543 1.056322 1.032615 1.027816      Inf
> round(jpn4.b.vif, digits=2)
      Token Type  NoS  TTR   GI MATTR  AWL  ASL
Token   Inf 4.77 2.35 1.71 1.30  1.02 1.01 1.20
Type   4.77  Inf 1.83 1.08 3.11  1.35 1.00 1.17
NoS    2.35 1.83  Inf 1.36 1.15  1.00 1.05 1.07
TTR    1.71 1.08 1.36  Inf 1.10  1.56 1.07 1.06
GI     1.30 3.11 1.15 1.10  Inf  3.88 1.02 1.06
MATTR  1.02 1.35 1.00 1.56 3.88   Inf 1.09 1.03
AWL    1.01 1.00 1.05 1.07 1.02  1.09  Inf 1.03
ASL    1.20 1.17 1.07 1.06 1.06  1.03 1.03  Inf
> jukaiki.ml2 <- lm(Score ~ Type + NoS + MATTR + AWL + ASL, data=jpn4.b)
> cor(jpn4.b[c(4,5,8,9,10)])
             Type         NoS      MATTR          AWL        ASL
Type  1.000000000  0.67377228 0.51028459  0.002673065  0.3775746
NoS   0.673772281  1.00000000 0.02273514 -0.209605424 -0.2531863
MATTR 0.510284594  0.02273514 1.00000000  0.286936810  0.1777219
AWL   0.002673065 -0.20960542 0.28693681  1.000000000  0.1645096
ASL   0.377574643 -0.25318632 0.17772193  0.164509608  1.0000000
> jpn4.b.cor2 <- cor(jpn4.b[c(4,5,8,9,10)])
> jpn4.b.vif2 <- 1/(1-jpn4.b.cor2^2)
> jpn4.b.vif2
          Type      NoS    MATTR      AWL      ASL
Type       Inf 1.831398 1.352065 1.000007 1.166266
NoS   1.831398      Inf 1.000517 1.045953 1.068494
MATTR 1.352065 1.000517      Inf 1.089720 1.032615
AWL   1.000007 1.045953 1.089720      Inf 1.027816
ASL   1.166266 1.068494 1.032615 1.027816      Inf
> plot(jukaiki.ml2)
 次の図を見るためには <Return> キーを押して下さい: 
 次の図を見るためには <Return> キーを押して下さい: 
 次の図を見るためには <Return> キーを押して下さい: 
 次の図を見るためには <Return> キーを押して下さい: 
> step(jukaiki.ml2)
Start:  AIC=-576.28
Score ~ Type + NoS + MATTR + AWL + ASL

        Df Sum of Sq    RSS     AIC
<none>               36.174 -576.28
- MATTR  1    1.1418 37.316 -569.43
- NoS    1    2.6349 38.809 -558.24
- Type   1    2.6557 38.830 -558.09
- ASL    1    4.0863 40.260 -547.78
- AWL    1    7.6256 43.800 -523.77

Call:
lm(formula = Score ~ Type + NoS + MATTR + AWL + ASL, data = jpn4.b)

Coefficients:
(Intercept)         Type          NoS        MATTR          AWL          ASL  
   -0.64428      0.01119      0.04275     -2.65094      0.53249      0.09177  

 データは正規分布していないといけないか

plot(jukaiki.ml2)


 重回帰分析 lmres()

https://cran.r-project.org/web/packages/pequod/index.html

  • ランダム効果は入れられない

install.packages("pequod")
library(pequod)

model.lmres <- lmres(目的 ~ 説明1+ 説明2 + 説明3 + 説明2:説明3, データ)

交互作用の分析

  • 単純主効果
model.int <- simpleSlope(model.lmres, pred="説明2", mod1="説明3")
PlotSlope(model.int)


データの中心化

  • 中心化をしたい項目について、オプションで指定する
, centered = c("項目A", "項目B")

Reference