[[ノート/ノート]]~
[[ノート/R]]~
訪問者数 &counter();      最終更新 &lastmod();

*Rで(重)回帰分析をする (2017-01-11) [#y69279b6]
線形での依存性を見つけたい。Excelだと説明変数が16個までで具合が悪い。

手順を忘れるので、メモしておく。

**Excelで計算する [#if2bf644]
(あとで)

------------------------------
**Rで計算する [#ta9d989e]

***CSVファイルを読み込む [#z1e046b1]
まず、作業ディレクトリの変更が必要なら、[[作業ディレクトリの変更:http://cse.naro.affrc.go.jp/takezawa/r-tips/r/06.html]]
 > setwd("ディレクトリ指定")

CSVファイルの読み込みは、
 > x0 <- read.csv("全体検討.csv")
で読み込む。[[ファイルからデータを読み込む:http://cse.naro.affrc.go.jp/takezawa/r-tips/r/40.html]]、 [[Rメモ:read.csv/write.csvで読み込んだり書き出したりするときの注意点:http://statsbeginner.hatenablog.com/entry/2014/08/31/230758]]

この時の元データは、ExcelからCSV形式で書き出したもの。

デフォルトで1行目をタイトル行(列タイトルの指定)と見なすが、漢字もOK。
但し、漢字を含む場合は、Linux上のRを使うのであれば、CSVファイル上でUTF-8に
変換しておくのがよい。(でないと後で引っかかる)

行番号(行タイトル)は、今回は指定しないで上から順に打たせる。行タイトルを
含む欄を指定したい場合は、row.names="列名"かrow.names=列番号を指定する。[[Rメモ:read.csv/write.csvで読み込んだり書き出したりするときの注意点:http://statsbeginner.hatenablog.com/entry/2014/08/31/230758]]

列名の変更は、行列に対してもデータフレームに対しても、
 > colnames(dat) <- c("A","B","C","D")
のような形で良いらしい。[[行列やデータフレームの列名・変数名の変更:http://webbeginner.hatenablog.com/entry/2014/05/09/234219]]

***データフレームに変換する [#a674cb43]
 > x1 <- as.data.frame(x0)
で変換する。[[データ型とデータ構造:http://cse.naro.affrc.go.jp/takezawa/r-tips/r/25.html]]

念のため、x1の内容を確認しておくこと。

***標準化 [#i00fb8ee]
要は、「各列が変量となっているデータ行列 x の各変量の単位が異なる場合に,各変量を平均が 0 ,分散が 1 になるように変換すること.規準化または標準化という.」ということ。
-[[基本統計量の算出:http://cse.naro.affrc.go.jp/takezawa/r-tips/r/59.html]]
-[[R - データの標準化:http://d.hatena.ne.jp/ykchat/20110908/1315486610]]
-[[データの標準化 R:http://upaksta.hatenablog.jp/entry/2014/09/23/232050]]

で、やり方は、
 > x2 <- scale(x1)

***回帰分析 lm [#f2d6b6d5]
あとは、回帰分析の関数lmに食わせればよいだけ。
 > x3 <- lm(目的変数~説明変数, data=データフレーム)

重回帰分析(複数の説明変数の線形結合で説明したい)の場合は、説明変数の指定を
 > x3 <- lm(目的変数~説明変数1+説明変数2+... , data=データフレーム)
のように書けばよい。[[R による重回帰分析:http://stat.biopapyrus.net/nlm/mlm-r2.html]]

**どう動くか試してみる テキストにあるwomenのデータでやってみよう [#l06072d3]
まず、womenのデータ自体は、Rを起動すると変数womenとして既にあるので、そのまま使ってしまおう。
 > women
    height weight
 1      58    115
 2      59    117
 3      60    120
 4      61    123
 5      62    126
 6      63    129
 7      64    132
 8      65    135
 9      66    139
 10     67    142
 11     68    146
 12     69    150
 13     70    154
 14     71    159
 15     72    164
次に標準化する
 > w0 <- scale(women)
 > w0
           height     weight
  [1,] -1.5652476 -1.4022687
  [2,] -1.3416408 -1.2732255
  [3,] -1.1180340 -1.0796608
  [4,] -0.8944272 -0.8860962
  [5,] -0.6708204 -0.6925315
  [6,] -0.4472136 -0.4989668
  [7,] -0.2236068 -0.3054021
  [8,]  0.0000000 -0.1118374
  [9,]  0.2236068  0.1462489
 [10,]  0.4472136  0.3398136
 [11,]  0.6708204  0.5978998
 [12,]  0.8944272  0.8559861
 [13,]  1.1180340  1.1140723
 [14,]  1.3416408  1.4366802
 [15,]  1.5652476  1.7592880
 attr(,"scaled:center")
   height   weight
  65.0000 136.7333
 attr(,"scaled:scale")
    height    weight
  4.472136 15.498694

このままlmに食わせると、データフレームでない、というエラーが出るので、データフレームに変換する。
 > w1 <- scale(w0)
 > w1
        height     weight
 1  -1.5652476 -1.4022687
 2  -1.3416408 -1.2732255
 3  -1.1180340 -1.0796608
 4  -0.8944272 -0.8860962
 5  -0.6708204 -0.6925315
 6  -0.4472136 -0.4989668
 7  -0.2236068 -0.3054021
 8   0.0000000 -0.1118374
 9   0.2236068  0.1462489
 10  0.4472136  0.3398136
 11  0.6708204  0.5978998
 12  0.8944272  0.8559861
 13  1.1180340  1.1140723
 14  1.3416408  1.4366802
 15  1.5652476  1.7592880

これで、lmに食わせると、
 > w2 <- lm(weight~height, data=w1)
 > w2
 
 Call:
 lm(formula = weight ~ height, data = w1)
 
 Coefficients:
 (Intercept)       height
   8.597e-16    9.955e-01

更に、summaryを出力すると
 > summary(w2)
 
 Call:
 lm(formula = weight ~ height, data = w1)
 
 Residuals:
      Min       1Q   Median       3Q      Max
 -0.11184 -0.07312 -0.02473  0.04785  0.20109
 
 Coefficients:
              Estimate Std. Error t value Pr(>|t|)
 (Intercept) 8.597e-16  2.541e-02    0.00        1
 height      9.955e-01  2.630e-02   37.85 1.09e-14 ***
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
 Residual standard error: 0.0984 on 13 degrees of freedom
 Multiple R-squared:  0.991,     Adjusted R-squared:  0.9903
 F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

次に、漢字のカラム名が動作するかどうかを確認しておきたい。
 > colnames(w0) <- c("身長","体重")
でカラム名を書き換えられるのだが(うまく書き換えられる)、どうも裏に英語のカラム名
が潜んでいそうなので、全面的にデータを作り直してみる。まずこの書き変えたデータ(行列)をcsvファイルに書き出し、それを再度読み込んでみる。
 > write.csv(w1, "women-j.csv", quote=TRUE, row.names=TRUE)
 
 > w11 <- read.csv("women-j.csv",header=TRUE, row.names=1)
 > w11
          身長       体重
 1  -1.5652476 -1.4022687
 2  -1.3416408 -1.2732255
 3  -1.1180340 -1.0796608
 4  -0.8944272 -0.8860962
 5  -0.6708204 -0.6925315
 6  -0.4472136 -0.4989668
 7  -0.2236068 -0.3054021
 8   0.0000000 -0.1118374
 9   0.2236068  0.1462489
 10  0.4472136  0.3398136
 11  0.6708204  0.5978998
 12  0.8944272  0.8559861
 13  1.1180340  1.1140723
 14  1.3416408  1.4366802
 15  1.5652476  1.7592880

これでよさそうなので、データフレームに変換
 > w21 <- as.data.frame(w11)
 > w21
          身長       体重
 1  -1.5652476 -1.4022687
 2  -1.3416408 -1.2732255
 3  -1.1180340 -1.0796608
 4  -0.8944272 -0.8860962
 5  -0.6708204 -0.6925315
 6  -0.4472136 -0.4989668
 7  -0.2236068 -0.3054021
 8   0.0000000 -0.1118374
 9   0.2236068  0.1462489
 10  0.4472136  0.3398136
 11  0.6708204  0.5978998
 12  0.8944272  0.8559861
 13  1.1180340  1.1140723
 14  1.3416408  1.4366802
 15  1.5652476  1.7592880

これでlmに食わせる
 > w31 <- lm(体重~身長,data=w21)
 > w31
 
 Call:
 lm(formula = 体重 ~ 身長, data = w21)
 
 Coefficients:
 (Intercept)         身長
   1.663e-15    9.955e-01
 
 > summary(w31)
 
 Call:
 lm(formula = 体重 ~ 身長, data = w21)
 
 Residuals:
      Min       1Q   Median       3Q      Max
 -0.11184 -0.07312 -0.02473  0.04785  0.20109
 
 Coefficients:
              Estimate Std. Error t value Pr(>|t|)
 (Intercept) 1.663e-15  2.541e-02    0.00        1
 身長        9.955e-01  2.630e-02   37.85 1.09e-14 ***
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
 Residual standard error: 0.0984 on 13 degrees of freedom
 Multiple R-squared:  0.991,     Adjusted R-squared:  0.9903
 F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14
 
これで、lmの引数に、漢字のコラム名を直接書いてOKらしいことが分かる。

***やりたかった重回帰分析 [#zb4de94a]
いよいよ、目的とする「社会的強み」の各項目と成績(単位数)との重回帰分析をしてみる。

データの読み込みは、csvファイルからとし、UTF-8に変換して置くのを忘れないこと。
 > x0 <- read.csv("全体検討.csv")

あとで標準化(scale)するときに、学生番号は具合が悪い(平均が取れない?)ので、覗いてしまう。
 x0[,2:(ncol(x0))]

更に、データは部分的に穴があいている(欠落)と具合が悪いらしいので、na.omitで取り除く。その始末はRに任せてしまった。要確認。
 na.omit(x0[,2:(ncol(x0))])

その上で、標準化(scale)
 > x1 <- scale(na.omit(x0[,2:(ncol(x0))]))


データフレームに変換
 > x2 <- as.data.frame(x1)

これで、回帰分析lmを適用する。
 > x3 <- lm(単位数 ~ 社会的強み平均(偏差値)+ ... ,data=x2)

ここでエラー。
 エラー:  想定外の入力です  in "x3 <- lm(単位数 ~ 社会的強み平均・"

どうも、欄名 「社会的強み平均(偏差値)」 のカッコ (偏差値) がいけないらしい。
面倒なので(どうせこの平均との相関はあまりないことが分かっているので)説明変数から外してしまおう、ということで、やってみると今度はOK。
 > x3 <- lm(単位数 ~ 意欲_偏差値+ ... ,data=x2)

で、結果が得られた。この先は[[結果>YAMA/IR/Rで履修・調査データの回帰分析]]を参照

&ref(./mytest.png,999%);


トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS