[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[ノート/ノート]]~
[[ノート/R]]~
訪問者数 &counter(); 最終更新 &lastmod();
*Ｒで（重）回帰分析をする （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形式で書き出したもの。
デフォルトで１行目をタイトル行（列タイトルの指定）と見なすが、漢字も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(目的変数~説明変数１+説明変数２+... , 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の引数に、漢字のコラム名を直接書いてＯＫらしいことが分かる。
***やりたかった重回帰分析 [#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(単位数 ~ 社会的強み平均・"
どうも、欄名 「社会的強み平均（偏差値）」 のカッコ （偏差値） がいけないらしい。
面倒なので（どうせこの平均との相関はあまりないことが分かっているので）説明変数から外してしまおう、ということで、やってみると今度はＯＫ。
> x3 <- lm(単位数 ~ 意欲_偏差値+ ... ,data=x2)
で、結果が得られた。この先は[[結果>YAMA/IR/Rで履修・調査データの回帰分析]]を参照
終了行:
[[ノート/ノート]]~
[[ノート/R]]~
訪問者数 &counter(); 最終更新 &lastmod();
*Ｒで（重）回帰分析をする （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形式で書き出したもの。
デフォルトで１行目をタイトル行（列タイトルの指定）と見なすが、漢字も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(目的変数~説明変数１+説明変数２+... , 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の引数に、漢字のコラム名を直接書いてＯＫらしいことが分かる。
***やりたかった重回帰分析 [#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(単位数 ~ 社会的強み平均・"
どうも、欄名 「社会的強み平均（偏差値）」 のカッコ （偏差値） がいけないらしい。
面倒なので（どうせこの平均との相関はあまりないことが分かっているので）説明変数から外してしまおう、ということで、やってみると今度はＯＫ。
> x3 <- lm(単位数 ~ 意欲_偏差値+ ... ,data=x2)
で、結果が得られた。この先は[[結果>YAMA/IR/Rで履修・調査データの回帰分析]]を参照
ページ名:
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
