R と RStudio でのデータ分析 入門
太郎丸 博
2022-10-19 13:23:23
- 1 read.csv データの読み込み
- 2 Summerying Data データの概観
- 3 xtabs() 度数分布表/クロス表
- 4 Plot グラフを作る
- 5 Manipulating Variables 変数の操作
- 6 Regression Analysis 回帰分析
データの分析には、
- データの読み込み（あるいは呼び出し）、
- 分析できるようにデータを整形、
- 探索的分析、最終的な分析結果の確定、
- 図や表の整形、レポート／論文の執筆
といったプロセスがある。これらを Windows 上の R と RStudio で行う方法を概説しよう。
1 read.csv データの読み込み
データはふつうEXCELのような表計算ソフトなどで行い、それを R のような統計解析ソフトに読み込むことが多い。分析の第一歩はデータを統計解析ソフトに読み込むことであるが、けっこう失敗しやすいので注意。
1.1 setwd 作業ディレクトリの指定
データを読み込む前に、作業ディレクトリ（作業フォルダともいう）を指定する。作業ディレクトリとは、Rがデフォルトでファイルの読み書きをするフォルダであり、特にパスを指定しなければ、R は作業フォルダのファイルを読みに行くし、ファイルを出力する場合もこのフォルダに出力のファイルを保存する。分析課題ごとに新しいフォルダを作ってそこに関連するファイルをまとめて保存しておくとよい。
1.1.1 GUI での作業フォルダの指定
作業フォルダを指定するには、RStudio のメニューの [Session] [Set Working Directory] [Choose Directory…] を選ぶとフォルダを選択するダイアログボックスが出てくるので、ここで選べばよい。
1.1.2 フルパスで指定
しかし、毎回これを行うのはけっこう面倒なので、 setwd() という関数をスクリプトの最初に書いてフォルダを指定したほうが簡単な場合も多い。 setwd() は
setwd(“作業ディレクトリのフルパス”)
という書式で指定する。例えば、以下のように。
setwd("C:/Users/Hiroshi/Dropbox/15_IJ_Courses/R_Practice")スクリプトはすべて半角で書くこと。また R は大文字と小文字を区別しているので、SETWDと書いてもエラーになる。また Windows の場合ディレクトリの区切りは “\” ではなく “/” を使うので注意。
作業ディレクトリの指定を誤って、保存したファイルがどこに行ったのか、わからなくなる人がたまにいるが、そういう場合は、以下のように引数なしで getwd() を使えば作業ディレクトリの場所がわかる。
getwd()## [1] "C:/Users/taroh/Dropbox/22G_Course/Rmd"1.2 Sample Data サンプル・データの呼び出し
1.2.1 Preinstalled Sample Data いきなり使えるサンプル・データ
R には様々なサンプル・データがあらかじめ付属しており、これらを使って、分析の練習や様々な関数の挙動や性質を確認することができる。いくつかのパッケージのデータは最初からロードされており、データの名前をタイプするだけでいきなり使える。例えば、USArrests というサンプル・データをすべてそのまま表示するには、以下のようにすればよい。
USArrests## Murder Assault UrbanPop Rape
## Alabama 13.2 236 58 21.2
## Alaska 10.0 263 48 44.5
## Arizona 8.1 294 80 31.0
## Arkansas 8.8 190 50 19.5
## California 9.0 276 91 40.6
## Colorado 7.9 204 78 38.7
## Connecticut 3.3 110 77 11.1
## Delaware 5.9 238 72 15.8
## Florida 15.4 335 80 31.9
## Georgia 17.4 211 60 25.8
## Hawaii 5.3 46 83 20.2
## Idaho 2.6 120 54 14.2
## Illinois 10.4 249 83 24.0
## Indiana 7.2 113 65 21.0
## Iowa 2.2 56 57 11.3
## Kansas 6.0 115 66 18.0
## Kentucky 9.7 109 52 16.3
## Louisiana 15.4 249 66 22.2
## Maine 2.1 83 51 7.8
## Maryland 11.3 300 67 27.8
## Massachusetts 4.4 149 85 16.3
## Michigan 12.1 255 74 35.1
## Minnesota 2.7 72 66 14.9
## Mississippi 16.1 259 44 17.1
## Missouri 9.0 178 70 28.2
## Montana 6.0 109 53 16.4
## Nebraska 4.3 102 62 16.5
## Nevada 12.2 252 81 46.0
## New Hampshire 2.1 57 56 9.5
## New Jersey 7.4 159 89 18.8
## New Mexico 11.4 285 70 32.1
## New York 11.1 254 86 26.1
## North Carolina 13.0 337 45 16.1
## North Dakota 0.8 45 44 7.3
## Ohio 7.3 120 75 21.4
## Oklahoma 6.6 151 68 20.0
## Oregon 4.9 159 67 29.3
## Pennsylvania 6.3 106 72 14.9
## Rhode Island 3.4 174 87 8.3
## South Carolina 14.4 279 48 22.5
## South Dakota 3.8 86 45 12.8
## Tennessee 13.2 188 59 26.9
## Texas 12.7 201 80 25.5
## Utah 3.2 120 80 22.9
## Vermont 2.2 48 32 11.2
## Virginia 8.5 156 63 20.7
## Washington 4.0 145 73 26.2
## West Virginia 5.7 81 39 9.3
## Wisconsin 2.6 53 66 10.8
## Wyoming 6.8 161 60 15.6ただし、このようにデータをすべて表示させるのは、しばしば非常に鬱陶しいので、head()という関数を使うとよい。この関数は、
head(データフレーム名または行列名)
でデータの最初の数行だけを表示してくれるので、意外と便利である。例えば、以下のように。
head(USArrests)## Murder Assault UrbanPop Rape
## Alabama 13.2 236 58 21.2
## Alaska 10.0 263 48 44.5
## Arizona 8.1 294 80 31.0
## Arkansas 8.8 190 50 19.5
## California 9.0 276 91 40.6
## Colorado 7.9 204 78 38.71.3 Loading Sample Data サンプル・データの呼び出し
大半のパッケージのデータはいったんロード（呼び出し）しないと使えるようにならない。 例えば、AER というパッケージに含まれる CPS1985 というサンプル・データを使いたいとしよう。 いきなり以下のようにしてもうまくいかない。
library(AER)
CPS1985## Error in eval(expr, envir, enclos): オブジェクト 'CPS1985' がありませんちなみに library() はパッケージを呼び出すための関数で、
library(パッケージ名)
でパッケージを呼び出せる。パッケージは RStudio の Packages タブで使いたいパッケージにチェックをいれることでも呼び出せる。CPS1985 のようにあらかじめロードされていないサンプル・データをつかうためには、
data(データ名)
とすればよい。例えば以下のように。
data(CPS1985)Rがエラー・メッセージを出さず、無反応ならばうまく行っているはずである。この後ならば以下のようにするなどして、CPS1985が使えるようになる。
head(CPS1985)## wage education experience age ethnicity region gender occupation sector union married
## 1 5.10 8 21 35 hispanic other female worker manufacturing no yes
## 1100 4.95 9 42 57 cauc other female worker manufacturing no yes
## 2 6.67 12 1 19 cauc other male worker manufacturing no no
## 3 4.00 12 4 22 cauc other male worker other no no
## 4 7.50 12 17 35 cauc other male worker other no yes
## 5 13.07 13 9 28 cauc other male worker other yes noR があらかじめもっているサンプル・データを確認するためには、
data()とすればよい。サンプル・データのリストが表示される。ただし、ロードしていないパッケージのデータは表示されないので注意。
1.3.1 Data Frame データフレーム
複数の変数からなるデータはふつうの行列でもよいが、行列は数値と文字の要素を混在させることができない。数値なら数値だけ、文字なら文字だけしか扱えない。データフレームという形式は、様々なタイプのデータをひとまとめにして扱えるので、便利である。R のサンプル・データのほとんどもデータフレームである。データフレームを作るには、data.frame() 関数を使うのが基本である。
data.frame(変数名 = ベクトル, 変数名 = ベクトル, …)
例えば以下のように。
dat1 <- data.frame(
sex = rep(c("Female", "Male"), 5), # スクリプトはコンマのあとなど途中で改行できる
age = rep(18:22, each = 2),
score = c(46, 7, 88, 54, 66, 44, 59, 35, 22, 70)
)
dat1## sex age score
## 1 Female 18 46
## 2 Male 18 7
## 3 Female 19 88
## 4 Male 19 54
## 5 Female 20 66
## 6 Male 20 44
## 7 Female 21 59
## 8 Male 21 35
## 9 Female 22 22
## 10 Male 22 70ベクトルの長さが変数によって異なると、エラーが出るので注意。
データフレームの中の特定の変数は、行列と同じように何列目の変数なのかを [, ] を使って取り出すこともできるし、
データフレーム名$変数名
としても同じように取り出せる。例えば以下のように。
dat1[, 3]## [1] 46 7 88 54 66 44 59 35 22 70dat1$score## [1] 46 7 88 54 66 44 59 35 22 701.4 read.csv 外部データの読み込み
外部データの読み込みには、read.table() や read.csv() といった関数を使うが、これらが読み込めるデータの形式は、テキストファイルだけである。EXCELやSPSSなどのファイルを読み込める関数も gdata や foreign のようなパッケージにあるが、個人的には R でデータを読み込む場合は、CSV形式かタブ区切り形式のテキスト・ファイルにデータをいったん変換してから読み込むのが簡単だと感じている。一番最初の行に変数名を書いたCSV形式のデータ（１行目は変数名）を作業ディレクトリに保存した上で
read.csv(“ファイル名”)
とすればよい。例えば “sample.csv” というCSV形式のファイルを作業ディレクトリに保存し、
data1 <- read.csv("sample.csv")とすれば、data1 という名前で R の中でそのデータを分析できるようになる。もちろんデータは作業ディレクトリ以外の場所においておいて、以下のようにフルパスを指定して読み込んでもよい。
data1 <- read.csv("c:/myfolder/temp/sample.csv")また、以下の関数を引数無しでファイル名の代わりに書いておくと、ダイアログボックスが開いて、読み込むファイルを選ぶことができる。
choose.files()
data1 <- read.csv(choose.files()) # ファイルを選ぶダイアログボックスが開く2 Summerying Data データの概観
2.1 Help of Sample Data サンプル・データのヘルプ
サンプル・データの場合、ヘルプを呼び出すとサンプル・データについてある程度は解説してあるので、参考にするとよい。例えば、以下のように。
?CPS1985ちなみに “?” は
?関数名
で関数のヘルプを呼び出すことができる。例えば、以下のように。
?head2.2 head() データの中身を見る
データの中身を見るには、すでに述べたようにデータの名前をタイプするという方法がもっとも単純だが、データのサイズが大きいと表示しきれず、かえって不便である。これもすでに述べたように head() でも見られるが、変数の数が多いとやはり見にくい。R の場合、変数の数があまり多いデータを読み込むと使いにくいので、必要な変数だけを含んだデータを読み込んだほうがいいように思われる。
2.3 summery() 記述統計
簡単な記述統計は、
summary(データフレーム名)
で表示できる。例えば、以下のように。
summary(CPS1985)## wage education experience age ethnicity region
## Min. : 1.000 Min. : 2.00 Min. : 0.00 Min. :18.00 cauc :440 south:156
## 1st Qu.: 5.250 1st Qu.:12.00 1st Qu.: 8.00 1st Qu.:28.00 hispanic: 27 other:378
## Median : 7.780 Median :12.00 Median :15.00 Median :35.00 other : 67
## Mean : 9.024 Mean :13.02 Mean :17.82 Mean :36.83
## 3rd Qu.:11.250 3rd Qu.:15.00 3rd Qu.:26.00 3rd Qu.:44.00
## Max. :44.500 Max. :18.00 Max. :55.00 Max. :64.00
## gender occupation sector union married
## male :289 worker :156 manufacturing: 99 no :438 no :184
## female:245 technical :105 construction : 24 yes: 96 yes:350
## services : 83 other :411
## office : 97
## sales : 38
## management: 55数値型の変数と文字型の変数（あるいは因子、論理型の変数）では要約の仕方が異なり、文字型の場合は、各カテゴリの度数を示す度数分布表が示され、数値型の変数の場合は平均値などの統計が示される。各項目の意味は、下記の通り。
| 項目 | 意味 |
|---|---|
| Min. | 最小値 |
| 1st Qu. | 第1四分位点 |
| Median | 中央値 |
| Mean | 平均値 |
| 3rd Qu. | 第3四分位点 |
| Max. | 最大値 |
| NA’s | 欠損値の数 |
この他にも標準偏差や分散が知りたい場合は、
sd(データ名$変数名)
var(データ名$変数名)
で計算できる。例えば以下のように。
sd(CPS1985$wage)## [1] 5.139097var(CPS1985$wage)## [1] 26.41032もちろん sd や var の引数は、データフレーム中の変数だけでなく数値ベクトルであれば何でもよい。
x <- c(5, 4, 6, 7, 20)
sd(x)## [1] 6.5802742.3.1 NA 欠損値
データが得られない場合などは欠損値を用いることが多い。欠損値には、NA を用いる。例えば、5人分の身長と体重のデータを取ったが、2人目の身長と 5人目の人の体重がわからなかった場合、以下のように入力する。
身長 <- c(160, NA, 180, 155, 172)
体重 <- c(50, 55, 70, NA, 70)デフォルトでは欠損値があるとエラーが出て計算されない関数も多い。sd( ), var( ), mean( ), sum( ) といった関数の場合、欠損値を除外して計算する場合は、“, na.rm = TRUE” というオプションを付ける。
x <- c(1, 4, 6, 5, NA)
sd(x)## [1] NAsd(x, na.rm = T)## [1] 2.1602472.3.2 apply() を使って一度に標準偏差を計算する
複数の変数について一つ一つ標準偏差を計算するのは面倒なので、一度に計算してしまいたい場合もある。そのような場合はいくつかやり方があるが、ここでは apply() という関数を使ったやり方を紹介しよう。
apply(標準偏差を計算したいデータフレーム名, 2, sd)
で計算できる。例えば以下のように。
d <- CPS1985[,1:4] # [, 1:4]はデータフレームや行列の1:4行目だけを取り出す指定
head(d)## wage education experience age
## 1 5.10 8 21 35
## 1100 4.95 9 42 57
## 2 6.67 12 1 19
## 3 4.00 12 4 22
## 4 7.50 12 17 35
## 5 13.07 13 9 28apply(d, 2, sd)## wage education experience age
## 5.139097 2.615373 12.379710 11.726573やはり、欠損値がある場合は、以下のように “,na.rm = T”とオプションを付ければよい。
apply(d, 2, sd, na.rm = T)var() の場合はそのままデータフレームや行列を引数とすれば、分散・共分散行列を返すので、対角線の値がそれぞれの変数の分散である。
(var.cps <- var(d) ) # () でくくると代入するだけでなく中身も見られる。## wage education experience age
## wage 26.410316 5.133282 5.538773 10.664731
## education 5.133282 6.840174 -11.418801 -4.601001
## experience 5.538773 -11.418801 153.257222 141.972170
## age 10.664731 -4.601001 141.972170 137.512508diag(var.cps) # diag(行列名)で正方行列の主対角線の要素を取り出してくれる。## wage education experience age
## 26.410316 6.840174 153.257222 137.5125083 xtabs() 度数分布表/クロス表
度数分布表はグラフで見るのが基本であるが、数や比率をきちんと出したい場合もある。そのような場合、xtabs() 関数を使う。書式は以下の通り。
xtabs(~ 変数名, data = データフレーム名)
以下は使用例。
xtabs(~ ethnicity, data = CPS1985)## ethnicity
## cauc hispanic other
## 440 27 673.1 table()
xtabs() でできることはほとんど table() でもできるが、スクリプトの書き方が異なる。
table(ベクトル)
以下が使用例
table(CPS1985 $ ethnicity)##
## cauc hispanic other
## 440 27 673.2 prop.table() 相対度数
相対度数（個々のカテゴリの比率）は prop.table() で計算できる。
prop.table(度数分布表)
以下が使用例。
tab.ethnicity <- xtabs(~ ethnicity, data = CPS1985)
prop.table(tab.ethnicity)## ethnicity
## cauc hispanic other
## 0.8239700 0.0505618 0.12546823.3 contingency table 2変数のクロス表
2つ以上のカテゴリカルな変数を組み合わせてサンプルを分割する表をクロス表とか分割表と呼ぶ。クロス表も xtabs() や table() で作れる。
xtabs(~ 変数1 + 変数2, data = データフレーム名)
以下が使用例。
xtabs(~ gender + ethnicity, data = CPS1985)## ethnicity
## gender cauc hispanic other
## male 236 14 39
## female 204 13 28上の表は male かつ cauc の人は 236人、female かつ other の人は 28 人、といった意味になる。同じことは table() でもできる。これらの 236、 28 といった数はセル度数と呼ばれたりする。
table(変数1, 変数2)
table(CPS1985$gender, CPS1985$ethnicity)##
## cauc hispanic other
## male 236 14 39
## female 204 13 283.4 クロス表の比率
クロス表の各セル度数の比率を計算する場合も prop.table() を使うが、分母を何にするかによって、3通り（3つ以上の変数を使ったクロス表（後述）ならもっとたくさん）の計算の仕方がある。
- 一番簡単なのは、すべてのセル度数をセル度数の総和で割るやり方である。この場合、
prop.table(クロス表)
とすればよい。以下使用例。
tab1 <- xtabs(~ gender + ethnicity, data = CPS1985)
prop.table(tab1)## ethnicity
## gender cauc hispanic other
## male 0.44194757 0.02621723 0.07303371
## female 0.38202247 0.02434457 0.052434461行1列目の値は、 tab1 の1行1列目のセル度数 236 を総度数 534 で割った値である。ほかのセルも、すべて tab1 のセル度数を 534 で割った値になる。これに決まった名前はないが、総パーセントと呼んでおく。
- 次に、各行の和でセル度数を割るというセル度数の比率の計算もある。例えば以下のように。
rowSums(tab1)## male female
## 289 245tab1 [1, ] / rowSums(tab1) [1] # 1行目の 3つのセル度数は、1行目のセル度数の総和で割る## cauc hispanic other
## 0.81660900 0.04844291 0.13494810tab1 [2, ] / rowSums(tab1) [2] # 2行目の 3つのセル度数は、2行目の総和で割る## cauc hispanic other
## 0.83265306 0.05306122 0.114285711つ目のベクトルは、男性の総数に占める白人男性（male かつ cauc）の比率、ヒスパニック男性（male かつ hispanic) の比率、その他の男性の比率 (male かつ other) の比率をあらわし、2つ目のベクトルも男性と同じように女性の総数に占める各エスニシティの比率を示す。これを 行パーセント と呼んでおく。また、男性の総数やヒスパニックの総数のように、行または列の総和を 周辺度数 (marginal frequency) と呼ぶ。これらを周辺度数と呼ぶのは、以下のようにクロス表の周辺に記載するのが一般的だからかもしれない。
addmargins(tab1) # addmargins(クロス表) で周辺度数を書き加えたクロス表を返す## ethnicity
## gender cauc hispanic other Sum
## male 236 14 39 289
## female 204 13 28 245
## Sum 440 27 67 534行パーセントは prop.table() でも計算できる。
prop.table(クロス表, margin = 1) # margin = は省略可
以下は使用例。
prop.table(tab1, margin = 1)## ethnicity
## gender cauc hispanic other
## male 0.81660900 0.04844291 0.13494810
## female 0.83265306 0.05306122 0.11428571prop.table(tab1, 1)## ethnicity
## gender cauc hispanic other
## male 0.81660900 0.04844291 0.13494810
## female 0.83265306 0.05306122 0.11428571- さらに各列の総和で対応するセル度数を割るやり方もある。これを列パーセントと呼んでおく。
colSums(tab1)## cauc hispanic other
## 440 27 67tab1 [, 1] / colSums(tab1) [1]## male female
## 0.5363636 0.4636364tab1 [, 2] / colSums(tab1) [2]## male female
## 0.5185185 0.4814815tab1 [, 3] / colSums(tab1) [3]## male female
## 0.5820896 0.4179104これらは、各エスニシティ内部での男女の比率を示す。列パーセントも以下のようにして prop.table() で計算できる。
prop.table(クロス表, margin = 2) # margin = は省略可
以下は使用例。
prop.table(tab1, margin = 2)## ethnicity
## gender cauc hispanic other
## male 0.5363636 0.5185185 0.5820896
## female 0.4636364 0.4814815 0.41791043.5 3変数以上を組み合わせたクロス表
以下のように 3つ以上の変数を組み合わせてクロス表を作ることもできる。
xtabs(~ 変数1 + 変数2 + 変数3, data = データフレーム名)
以下使用例。
xtabs(~ gender + union + ethnicity, data = CPS1985)## , , ethnicity = cauc
##
## union
## gender no yes
## male 186 50
## female 181 23
##
## , , ethnicity = hispanic
##
## union
## gender no yes
## male 10 4
## female 12 1
##
## , , ethnicity = other
##
## union
## gender no yes
## male 25 14
## female 24 4上の表では白人かつ男性かつ組合に加入していない人が 186人、ヒスパニックかつ女性かつ組合に加入している人が 1人、という具合に読めばよい。さらに4つ目の変数を組み合わせたければ、
xtabs(~ 変数1 + 変数2 + 変数3 + 変数4, data = データフレーム名)
という具合にすればよいが、煩雑なのであまり使うことはないだろう。3変数を組み合わせたクロス表を三重クロス表、四変数を組み合わせたクロス表を四重クロス表というが、二重クロス表という言い方は聞いたことがない。
3.5.1 三重クロス表の比率
やはり prop.table() を使えば計算できるが、7通りの計算の仕方があり、ぜんぶ示すのは煩雑である。7通りあると言っても、三種類に分類できるので、計算例を示す。
tab2 <- xtabs(~ gender + union + ethnicity, data = CPS1985)
# 総パーセント：各セル度数をそれらの総和で割る
prop.table(tab2)## , , ethnicity = cauc
##
## union
## gender no yes
## male 0.348314607 0.093632959
## female 0.338951311 0.043071161
##
## , , ethnicity = hispanic
##
## union
## gender no yes
## male 0.018726592 0.007490637
## female 0.022471910 0.001872659
##
## , , ethnicity = other
##
## union
## gender no yes
## male 0.046816479 0.026217228
## female 0.044943820 0.007490637# 一行目のセル度数は、男性の人数の総和で割り、二行目のセル度数は女性の人数の総和で割る
round( # 数字が多すぎるので、四捨五入
prop.table(tab2, 1),
2
)## , , ethnicity = cauc
##
## union
## gender no yes
## male 0.64 0.17
## female 0.74 0.09
##
## , , ethnicity = hispanic
##
## union
## gender no yes
## male 0.03 0.01
## female 0.05 0.00
##
## , , ethnicity = other
##
## union
## gender no yes
## male 0.09 0.05
## female 0.10 0.02一行目（male) は、男性のセル度数を男性の総数でわったもので、それらをすべて足し合わせると \(0.64 + 0.17 + 0.03 + 0.01 + 0.09 + 0.05 \fallingdotseq 1\) になる。
# 性別とエスニシティを組み合わせた周辺度数で、対応するセル度数を割る。
round(
prop.table(tab2, margin = c(1, 3)),
2
)## , , ethnicity = cauc
##
## union
## gender no yes
## male 0.79 0.21
## female 0.89 0.11
##
## , , ethnicity = hispanic
##
## union
## gender no yes
## male 0.71 0.29
## female 0.92 0.08
##
## , , ethnicity = other
##
## union
## gender no yes
## male 0.64 0.36
## female 0.86 0.144 Plot グラフを作る
変数の分布を探索的にチェックする場合、グラフの作成が非常に有益である。
4.1 Histogram ヒストグラム
連続変数の場合はヒストグラムをまずは作るのが基本であろう。ヒストグラムは hist() 関数で作る。
hist(データフレーム名$変数名)
例えば以下のように。
hist(CPS1985$wage)
4.2 Bar Plot 棒グラフ
離散変数（文字型のデータ）の場合は度数分布表を単純にグラフにすればよい。そのためには、plot()関数で
plot(データフレーム名$変数名)
とすればよい。例えば以下のように。
plot(CPS1985$ ethnicity)
4.3 Scatter Plot 散布図
連続変数どうしの関係は散布図で見る。これにも plot() を使う。
plot(y軸に配する変数名 ~ x軸に配する変数名, data = データフレーム名)
例えば、以下のようにする。
plot(wage ~ experience, data = CPS1985)
すべての変数の組み合わせの散布図を作るのは、いちいちスクリプトを書いていてはたいへんなので、以下のようにする。
plot(データフレーム名)
例えば、以下のようにする。
plot(CPS1985[, 1:4])
上の図で1行2列目のグラフは縦軸が wage、 横軸が education の散布図である。一番右上のグラフは縦軸が wage で横軸が age である。一般にあるグラフと同じ列（上か下）にある文字がそのグラフの横軸の変数を示し、同じ行（右か左）にある文字がそのグラフの縦軸の変数を示す。
教育年数のようにいちおう数値型であっても実際には細かい数値をとりえない場合、散布図の点が完全に重なってしまい、どこに点が集中しているのかわからない。そのような場合、jitter() という関数を使うとほんの少し点をずらして散布図を見やすくすることができる。
jitter(ベクトル)
こうすると、元のベクトルに平均が 0 でちいさな分散の一様分布の乱数を足しあわせたベクトルを返してくれる。具体的には、以下のように使うとよい。
plot(wage ~ jitter(education), data = CPS1985)
こうするとどのあたりの事例が多いのかよくわかる。
4.4 Box Plot 箱ひげ図
グループ別に数値型の変数の分布を見たい場合は、箱ひげ図を使うのが基本である。箱ひげ図は、
plot(数値型の変数名 ~ グループ分けするためのカテゴリカル変数名, data = データフレーム名)
または
boxplot(数値型の変数名 ~ グループ分けするためのカテゴリカル変数名, data = データフレーム名)
とすれば作れる。
plot(wage ~ ethnicity, data = CPS1985)
4.5 Mosaic Plot モザイクプロット
離散変数どうしの関係は、モザイクプロットを上手に作れれば見やすい。これも散布図や箱ひげ図と同じように plot() で作れる。例えば以下のように。
plot(union ~ gender, data = CPS1985)
この例では、棒の幅が男女の比率に比例するように作られている。そして濃いグレーの部分の高さが男女それぞれの組合未加入率を示している。
spineplot() でも同じものが作れる。
5 Manipulating Variables 変数の操作
変数の操作はベクトルを使った計算なので、ベクトルを使った計算の応用出来さえすれば、新しい知識はあまり必要ない。
5.1 Simple Arithmetic 簡単な計算
CPS1985 の wage という変数は時給なので、これから１日８時間働いた場合の日給を計算したいなら、以下のようにすればよい。
CPS1985$wage.d <- CPS1985$wage * 8
head(CPS1985)## wage education experience age ethnicity region gender occupation sector union married
## 1 5.10 8 21 35 hispanic other female worker manufacturing no yes
## 1100 4.95 9 42 57 cauc other female worker manufacturing no yes
## 2 6.67 12 1 19 cauc other male worker manufacturing no no
## 3 4.00 12 4 22 cauc other male worker other no no
## 4 7.50 12 17 35 cauc other male worker other no yes
## 5 13.07 13 9 28 cauc other male worker other yes no
## wage.d
## 1 40.80
## 1100 39.60
## 2 53.36
## 3 32.00
## 4 60.00
## 5 104.56CPS1985の最後の列に、新しく wage.d という変数が付け加えられる。新しい変数を作ったときは、必ず意図したとおりに新しい変数が作られているかチェックすること。エラーが出ないからといって失敗していないという保証はない。
また、自然対数などの関数も使える。時給を対数変換した変数を新しく作りたいならば、以下のようにすれば良い。
CPS1985$wage.log <- log(CPS1985$wage)
head(CPS1985)## wage education experience age ethnicity region gender occupation sector union married
## 1 5.10 8 21 35 hispanic other female worker manufacturing no yes
## 1100 4.95 9 42 57 cauc other female worker manufacturing no yes
## 2 6.67 12 1 19 cauc other male worker manufacturing no no
## 3 4.00 12 4 22 cauc other male worker other no no
## 4 7.50 12 17 35 cauc other male worker other no yes
## 5 13.07 13 9 28 cauc other male worker other yes no
## wage.d wage.log
## 1 40.80 1.629241
## 1100 39.60 1.599388
## 2 53.36 1.897620
## 3 32.00 1.386294
## 4 60.00 2.014903
## 5 104.56 2.570320変数をZ得点（平均 0 標準偏差 1 の変数）に変換したいならば、scale() 関数を使えばよい。
scale(ベクトル)
例えば年齢をｚ以下のように。
CPS1985$age.z <- scale(CPS1985$age)
mean(CPS1985$age.z)## [1] -2.084775e-16sd(CPS1985$age.z)## [1] 15.2 Recode リコード
変数のいくつかのカテゴリや数値をまとめて新しい値やカテゴリを作りたい場合がある。やり方はいくつかあるが、SPSS の recode ライクな関数から紹介しよう。 SPSSの recode に近い働きをする関数をいくつかのパッケージが用意しているが、ここでは memisc パッケージの recode() 関数を紹介しよう（同じ名前の関数を複数のパッケージがもっている場合、最後に読み込んだパッケージの関数が用いられるので注意。recode() は car というパッケージにも含まれている）。これは以下のような書式で使う。
recode(ベクトル, 新しい値 <- 古い値（のベクトル）, 新しい値 <- 古い値（のベクトル）, … )
例えば、CPS1985 の occupation という変数は5つのカテゴリからなるが、これを upper, middle, lower の 3カテゴリにまとめなおしたいとしよう。その場合、以下のようにすればよい。
library(memisc)## 要求されたパッケージ lattice をロード中です## 要求されたパッケージ MASS をロード中です##
## 次のパッケージを付け加えます: 'memisc'## 以下のオブジェクトは 'package:car' からマスクされています:
##
## recode## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## contr.sum, contr.treatment, contrasts## 以下のオブジェクトは 'package:base' からマスクされています:
##
## as.arrayCPS1985$occ3 <- recode(CPS1985$occupation,
"upper" <- c("technical", "management"),
"middle" <- c("services", "office", "sales"),
"lower" <- "worker"
)
xtabs(~occupation + occ3, CPS1985) # もとの変数と新しい変数のクロス表## occ3
## occupation upper middle lower
## worker 0 0 156
## technical 105 0 0
## services 0 83 0
## office 0 97 0
## sales 0 38 0
## management 55 0 0xtabs() はクロス表を作る関数で、以下のような書式で書く。
xtabs(~ 行の変数名 + 列の変数名, データ名)
5.3 cut() 連続変数のリコード
連続変数をいくつかのカテゴリにまとめ直す場合、cut() という関数が便利である。書式は以下の通り。
cut(データ名$変数名, 区切り位置を示すベクトル)
ただし区切り位置は上限値と下限値も含まなければならない。
例えば、労働経験年数が 0～9年、10～19年、20～29年、30年以上という4つのカテゴリからなる変数を作りたい場合、以下のようにすればよい。
CPS1985$exp4 <- cut(CPS1985$experience, c(-1, 9, 19, 29, 60))
xtabs(~ exp4, CPS1985)## exp4
## (-1,9] (9,19] (19,29] (29,60]
## 156 183 91 104xtabs(~ experience + exp4, CPS1985) [1 : 25, ] # 長すぎるので最初の２５行だけ表示## exp4
## experience (-1,9] (9,19] (19,29] (29,60]
## 0 11 0 0 0
## 1 12 0 0 0
## 2 15 0 0 0
## 3 18 0 0 0
## 4 16 0 0 0
## 5 15 0 0 0
## 6 17 0 0 0
## 7 18 0 0 0
## 8 19 0 0 0
## 9 15 0 0 0
## 10 0 23 0 0
## 11 0 11 0 0
## 12 0 18 0 0
## 13 0 23 0 0
## 14 0 28 0 0
## 15 0 18 0 0
## 16 0 22 0 0
## 17 0 15 0 0
## 18 0 11 0 0
## 19 0 14 0 0
## 20 0 0 13 0
## 21 0 0 7 0
## 22 0 0 10 0
## 23 0 0 6 0
## 24 0 0 11 0この例の場合、experienceを \(x\) と表記すると、新しく作った exp4 という変数の4つのカテゴリは \(-1<x \leq 9\), \(9 < x \leq 19\), \(19 < x \leq 29\), \(29 < x \leq 60\) に対応している。注意すべきなのは、中間の区切り値 9, 19, 29 は値の小さい方のカテゴリに含まれるという点と、-1 と 60 のように上限と下限の値も区切り値として指定する必要があるということである。cut() で作ったカテゴリの名前を指定したい場合は、, labels = カテゴリの名前のベクトル という引数をつける。例えば以下のように。
CPS1985$exp4 <- cut(CPS1985$experience, c(-1, 9, 19, 29, 60),
labels = c("0～9年", "10～19年", "20～29年", "30年以上")
)
xtabs(~ exp4, CPS1985)## exp4
## 0～9年 10～19年 20～29年 30年以上
## 156 183 91 1045.4 If Condition 条件を指定して代入
上の労働経験年数を 4カテゴリに分類する例は、
データフレーム名$ 新しい変数名 [条件式] <- 新しい値
という書式でもできる。具体的には以下のようにすればよい。
x <- CPS1985$experience # 何回も参照するので短い名前をつけておく
CPS1985$ exp4 [x <= 9] <- "0～9年"
CPS1985$ exp4 [x >= 10 & x <= 19] <- "10～19年"
CPS1985$ exp4 [x >= 20 & x <= 29] <- "20～29年"
CPS1985$ exp4 [x >= 30] <- "30年以上"
xtabs( ~ exp4, CPS1985)## exp4
## 0～9年 10～19年 20～29年 30年以上
## 156 183 91 104この方法は煩雑に思える場合もあるが、柔軟で汎用性が高いのでぜひ習得してほしい。なおこのままではただの文字列で因子になっていないので、以下のように因子に変換しておいた方がいい。
CPS1985$exp4 <- factor(CPS1985$ exp4)5.5 Making Binary Variable from Logical Operation 論理値を使った二値変数の作成
論理演算子を使うと比較的簡単に二値変数が作れる。例えば教育年数が 13年以上のとき 1、13年未満のとき 0 をとるような二値変数を作りたい場合、以下のようにすればよい。
CPS1985$ edu13 <- CPS1985$ education >= 13
xtabs(~ edu13, CPS1985)## edu13
## FALSE TRUE
## 302 232xtabs(~ education + edu13, CPS1985)## edu13
## education FALSE TRUE
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 3 0
## 7 5 0
## 8 15 0
## 9 12 0
## 10 17 0
## 11 27 0
## 12 219 0
## 13 0 37
## 14 0 56
## 15 0 13
## 16 0 71
## 17 0 24
## 18 0 311, 0 ではなく TRUE と FALSE という2つの値をとるが、数値計算をすると TRUE と FALSE はそれぞれ 1 と 0 に置き換えて計算されるので、1, 0 の二値変数と同じである。
sum(CPS1985$edu13) # 1, 0 に変換されるので、TRUEの数と同じ値が返される## [1] 2325.6 subset() データの選択
特定の条件を満たすケースだけを取り出したデータを作りたい場合がある。その場合、データフレームでも行列と同じように処理できる。
データフレーム名 [条件式, ]
例えば CPS1985 から女性のケースだけを取り出したデータを作りたい場合、以下のようにすればよい。
dat2 <- CPS1985 [CPS1985$ gender == "female", ]
head(dat2)## wage education experience age ethnicity region gender occupation sector union married wage.d
## 1 5.10 8 21 35 hispanic other female worker manufacturing no yes 40.80
## 1100 4.95 9 42 57 cauc other female worker manufacturing no yes 39.60
## 20 5.71 12 16 34 cauc other female worker manufacturing yes yes 45.68
## 28 3.35 12 8 26 cauc other female worker manufacturing no yes 26.80
## 31 4.00 12 46 64 cauc south female worker other no no 32.00
## 33 5.00 17 1 24 cauc south female worker other no no 40.00
## wage.log age.z occ3 exp4 edu13
## 1 1.629241 -0.1563401 lower 20～29年 FALSE
## 1100 1.599388 1.7197409 lower 30年以上 FALSE
## 20 1.742219 -0.2416165 lower 10～19年 FALSE
## 28 1.208960 -0.9238278 lower 0～9年 FALSE
## 31 1.386294 2.3166758 lower 30年以上 FALSE
## 33 1.609438 -1.0943806 lower 0～9年 TRUEsummary(dat2$ gender)## male female
## 0 245また、年齢が30歳未満のケースだけを取り出したデータを作りたい場合、以下のようにすればよい。
dat3 <- CPS1985 [CPS1985$ age < 30, ]
head(dat3)## wage education experience age ethnicity region gender occupation sector union married wage.d
## 2 6.67 12 1 19 cauc other male worker manufacturing no no 53.36
## 3 4.00 12 4 22 cauc other male worker other no no 32.00
## 5 13.07 13 9 28 cauc other male worker other yes no 104.56
## 7 19.47 12 9 27 cauc other male worker other no no 155.76
## 9 8.75 12 9 27 cauc other male worker other no no 70.00
## 22 3.75 12 9 27 cauc other male worker other no no 30.00
## wage.log age.z occ3 exp4 edu13
## 2 1.897620 -1.5207626 lower 0～9年 FALSE
## 3 1.386294 -1.2649334 lower 0～9年 FALSE
## 5 2.570320 -0.7532749 lower 0～9年 TRUE
## 7 2.968875 -0.8385513 lower 0～9年 FALSE
## 9 2.169054 -0.8385513 lower 0～9年 FALSE
## 22 1.321756 -0.8385513 lower 0～9年 FALSEsummary(dat3$ age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 22.00 25.00 24.45 27.00 29.00R の場合、分析に使う変数だけからなるデータフレームを作ったほうが、しばしば分析は簡単になるので、必要な変数だけを取り出して新しいデータフレームを作ることがおすすめである。そのような場合、
データフレーム名[, 取り出す変数が何列目かを示すベクトル]
または、
データフレーム名[, 取り出す変数名のベクトル]
とすればよい。後者のほうがオススメ。なぜなら、データの操作を試行錯誤している段階では、変数の順番が入れ替わったり、途中に新しい変数が挿入されたりすることもあるから。そうすると、取り出す変数の列番号がずれるが、それを忘れて取り出す変数の列番号を修正しないと、思わぬ失敗につながる。私は何度かやらかした。
スクリプトの例。CPS1985 の 1, 4, 7 列目の wage, age, gender の 3 つの変数だけからなるデータフレームを作る場合、以下のようにすればよい。
dat4 <- CPS1985 [, c(1, 4, 7)]
head(dat4)## wage age gender
## 1 5.10 35 female
## 1100 4.95 57 female
## 2 6.67 19 male
## 3 4.00 22 male
## 4 7.50 35 male
## 5 13.07 28 maledat4 <- CPS1985 [, c("wage", "age", "gender")]
head(dat4)## wage age gender
## 1 5.10 35 female
## 1100 4.95 57 female
## 2 6.67 19 male
## 3 4.00 22 male
## 4 7.50 35 male
## 5 13.07 28 male変数名をいちいち ” ” でくくったり、使いたい変数が何列目か確認するのが意外と面倒なので、subset( ) という関数が一部のケースと変数を取り出したいとき意外と重宝する。
subset(データフレーム名, ケースを選択する条件, 取り出す変数名のベクトル)
例えば CPS1985 から30歳未満の男性のケースだけを取り出し、wage, experience, education の3変数だけを用いたい場合は、以下のようにすればよい。
dat5 <- subset(CPS1985, gender == "male" & age < 30, c(wage, experience, education))
head(dat5)## wage experience education
## 2 6.67 1 12
## 3 4.00 4 12
## 5 13.07 9 13
## 7 19.47 9 12
## 9 8.75 9 12
## 22 3.75 9 125.6.1 na.omit() 欠損値のあるケースをリストワイズで削除する
欠損値のあるケースの扱いは、実際のデータ分析では重要である。伝統的にはリストワイズ法（欠損値のあるケースは分析から除外）がよく用いられてきたが、リストワイズ法を用いる場合、あらかじめ欠損値を除外したデータフレームを作っておいたほうが便利な場合がしばしばある。その場合、na.omit() 関数を使う。以下のような書式でリストワイズで欠損値を含むケースを除外したデータフレームを返してくれる。
na.omit(欠損値を含んだデータフレーム)
例えば以下のようにする。
# 欠損値を含んだ架空のデータを作る
dat6 <- data.frame(sex = gl(2, 5, labels = c("female", "male")),
score = 1:10)
dat6$ sex [2] <- NA
dat6$ score [6:7] <- NA
dat6## sex score
## 1 female 1
## 2 <NA> 2
## 3 female 3
## 4 female 4
## 5 female 5
## 6 male NA
## 7 male NA
## 8 male 8
## 9 male 9
## 10 male 10# リストワイズで欠損値を含むケースを除外
dat7 <- na.omit(dat6)
dat7## sex score
## 1 female 1
## 3 female 3
## 4 female 4
## 5 female 5
## 8 male 8
## 9 male 9
## 10 male 106 Regression Analysis 回帰分析
回帰分析は、もっとも重要かつ基本的な分析法の一つであるが、lm( ) 関数で計算できる。
lm(従属変数名 ~ 独立変数名, data = データフレーム名)
例えば以下のように。
lm(wage ~ age, data = CPS1985)##
## Call:
## lm(formula = wage ~ age, data = CPS1985)
##
## Coefficients:
## (Intercept) age
## 6.16747 0.07755しかし、これでは係数が有意か、決定係数はどの程度かなど詳しいことがわからないので、さらに summary( ) で詳しい結果を引き出す。
summary(回帰分析の結果)
例えば以下のように。
lm1 <- lm(wage ~ age, data = CPS1985)
summary(lm1)##
## Call:
## lm(formula = wage ~ age, data = CPS1985)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.425 -3.503 -1.106 2.276 36.704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.16747 0.72280 8.533 < 2e-16 ***
## age 0.07755 0.01870 4.147 3.92e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.063 on 532 degrees of freedom
## Multiple R-squared: 0.03132, Adjusted R-squared: 0.0295
## F-statistic: 17.2 on 1 and 532 DF, p-value: 3.917e-056.1 Drawing a Regression Line 回帰直線の描画
散布図に回帰直線を引きたい場合は、まず散布図を作り、単回帰分析をし、その結果から abline( ) 関数で、
abline(回帰分析の結果)
とする。例えば以下のように。
plot(wage ~ age, data = CPS1985)
abline(lm1)
6.2 Multiple Regression 重回帰分析
重回帰分析の場合は、複数の独立変数を “+” でつないで指定する。例えば以下のように。
lm2 <- lm(wage ~ age + education + experience, data = CPS1985)
summary(lm2)##
## Call:
## lm(formula = wage ~ age + education + experience, data = CPS1985)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.351 -2.857 -0.599 1.994 36.336
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.76987 7.04271 -0.677 0.499
## age -0.02241 1.15475 -0.019 0.985
## education 0.94833 1.15524 0.821 0.412
## experience 0.12756 1.15571 0.110 0.912
##
## Residual standard error: 4.604 on 530 degrees of freedom
## Multiple R-squared: 0.202, Adjusted R-squared: 0.1975
## F-statistic: 44.73 on 3 and 530 DF, p-value: < 2.2e-166.2.1 Categorical Variable カテゴリカル変数
カテゴリカル変数もそのまま独立変数として指定すれば、自動的に最初のカテゴリを基準とするダミー変数としてモデルに投入される。例えば以下のように。
lm3 <- lm(wage ~ age + education + ethnicity, data = CPS1985)
summary(lm3)##
## Call:
## lm(formula = wage ~ age + education + ethnicity, data = CPS1985)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.302 -2.900 -0.655 2.087 36.212
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.26002 1.30800 -4.021 6.63e-05 ***
## age 0.10479 0.01722 6.085 2.24e-09 ***
## education 0.81050 0.07799 10.393 < 2e-16 ***
## ethnicityhispanic -0.40965 0.92267 -0.444 0.657
## ethnicityother -0.85043 0.60450 -1.407 0.160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.599 on 529 degrees of freedom
## Multiple R-squared: 0.2051, Adjusted R-squared: 0.1991
## F-statistic: 34.13 on 4 and 529 DF, p-value: < 2.2e-166.2.2 Squared Term 二乗項
x という変数の二乗項をモデルに加えたい場合、
I(x ^ 2)
という項を独立変数として指定すればよい。例えば以下のように。
lm4 <- lm(wage ~ age + I(age ^ 2), data = CPS1985)
summary(lm4)##
## Call:
## lm(formula = wage ~ age + I(age^2), data = CPS1985)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.396 -3.273 -1.016 1.935 38.157
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.672133 2.274396 -1.615 0.107
## age 0.618912 0.120294 5.145 3.77e-07 ***
## I(age^2) -0.006761 0.001485 -4.554 6.54e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.971 on 531 degrees of freedom
## Multiple R-squared: 0.06772, Adjusted R-squared: 0.06421
## F-statistic: 19.29 on 2 and 531 DF, p-value: 8.207e-09一般にモデルの式の中で “+, -, /, *” といった二項演算子を使う場合は I() でくくってやる必要がある。例えば、age をそのまま独立変数とするのではなく、平均値でセンタリング（この変数の値からこの変数の平均値を引くこと）して使う場合、以下のようにする。
lm4.1 <- lm(wage ~ I(age - mean(age)), data = CPS1985)
summary(lm4.1)##
## Call:
## lm(formula = wage ~ I(age - mean(age)), data = CPS1985)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.425 -3.503 -1.106 2.276 36.704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.02406 0.21909 41.190 < 2e-16 ***
## I(age - mean(age)) 0.07755 0.01870 4.147 3.92e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.063 on 532 degrees of freedom
## Multiple R-squared: 0.03132, Adjusted R-squared: 0.0295
## F-statistic: 17.2 on 1 and 532 DF, p-value: 3.917e-056.2.3 Log Transformation 対数変換
変数を自然対数に変換して使う場合は、I() は必要なく、 log() をモデルの式の中で指定してやればよい。例えば以下のように。
lm4.2 <- lm(log(wage) ~ occupation + sector, data = CPS1985)
summary(lm4.2)##
## Call:
## lm(formula = log(wage) ~ occupation + sector, data = CPS1985)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.35040 -0.34500 -0.01735 0.33188 1.46785
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.06602 0.05102 40.491 < 2e-16 ***
## occupationtechnical 0.42008 0.06547 6.417 3.11e-10 ***
## occupationservices -0.18889 0.07222 -2.615 0.00917 **
## occupationoffice -0.01626 0.06768 -0.240 0.81028
## occupationsales -0.06339 0.09077 -0.698 0.48524
## occupationmanagement 0.41048 0.07975 5.147 3.74e-07 ***
## sectorconstruction 0.04826 0.10940 0.441 0.65930
## sectorother -0.12610 0.06063 -2.080 0.03802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4793 on 526 degrees of freedom
## Multiple R-squared: 0.1861, Adjusted R-squared: 0.1753
## F-statistic: 17.18 on 7 and 526 DF, p-value: < 2.2e-166.2.4 Interaction Effects 交互作用項
複数の変数の交互作用もモデルに投入できる。x と y という変数の交互作用項をモデルに投入するには、独立変数を
x * y
とすれば主効果と交互作用効果がモデルに投入される。例えば以下のように。
lm5 <- lm(wage ~ gender * experience, data = CPS1985)さらに高次の交互作用を指定するには、“x * y * z” のようにすればよい。
6.3 複数の回帰分析の結果を表にまとめる
複数の回帰分析の結果をまとめて表にするには、memisc パッケージの mtable( ) 関数や texreg パッケージの screenreg() 関数など、いくつかの関数がいくつかのパッケージで提供されている。
mtable(回帰分析の結果1, 回帰分析の結果2, 回帰分析の結果3, …)
screenreg(list(回帰分析の結果1, 回帰分析の結果2, 回帰分析の結果3, …))
例えば以下のように
library(memisc)
mtable(lm1, lm2, lm3)##
## Calls:
## lm1: lm(formula = wage ~ age, data = CPS1985)
## lm2: lm(formula = wage ~ age + education + experience, data = CPS1985)
## lm3: lm(formula = wage ~ age + education + ethnicity, data = CPS1985)
##
## ==============================================================
## lm1 lm2 lm3
## --------------------------------------------------------------
## (Intercept) 6.167*** -4.770 -5.260***
## (0.723) (7.043) (1.308)
## age 0.078*** -0.022 0.105***
## (0.019) (1.155) (0.017)
## education 0.948 0.811***
## (1.155) (0.078)
## experience 0.128
## (1.156)
## ethnicity: hispanic/cauc -0.410
## (0.923)
## ethnicity: other/cauc -0.850
## (0.604)
## --------------------------------------------------------------
## R-squared 0.031 0.202 0.205
## N 534 534 534
## ==============================================================
## Significance: *** = p < 0.001; ** = p < 0.01; * = p < 0.05library(texreg)## Version: 1.38.6
## Date: 2022-04-06
## Author: Philip Leifeld (University of Essex)
##
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").screenreg(list(lm1, lm2, lm3))##
## ==================================================
## Model 1 Model 2 Model 3
## --------------------------------------------------
## (Intercept) 6.17 *** -4.77 -5.26 ***
## (0.72) (7.04) (1.31)
## age 0.08 *** -0.02 0.10 ***
## (0.02) (1.15) (0.02)
## education 0.95 0.81 ***
## (1.16) (0.08)
## experience 0.13
## (1.16)
## ethnicityhispanic -0.41
## (0.92)
## ethnicityother -0.85
## (0.60)
## --------------------------------------------------
## R^2 0.03 0.20 0.21
## Adj. R^2 0.03 0.20 0.20
## Num. obs. 534 534 534
## ==================================================
## *** p < 0.001; ** p < 0.01; * p < 0.056.4 Plotting Predictions モデルからの予測値のプロット
複雑なモデルを推定した場合、推定結果の意味は係数を見ているだけではよくわからない場合がある。その場合、モデルからの予測値をプロットするとわかりやすい。
pred(回帰分析の結果、予測値を出したい独立変数の値を入力したデータフレーム)
必ず「予測値を出したい独立変数の値を入力したデータフレーム」を作らなければならないので、自分がどんな予測値をプロットしたいのか、よく考える必要がある。例えば、先ほどの age の二乗項を投入した分析結果から予測される18から64歳までの賃金をプロットしたいとしよう。その場合、まずは以下のように計算し、プロットする。
data.for.pred <- data.frame(age = 18 : 64) # 年齢の最小値と最大値の範囲で指定
head(data.for.pred) # 作ったデータフレームの中身を確認## age
## 1 18
## 2 19
## 3 20
## 4 21
## 5 22
## 6 23pred1 <- predict(lm4, data.for.pred)
plot(18 : 64, pred1, type = "l") # type = "l" は点をプロットするのではなく線を描けという指定
gender と experience の交互作用を推定したモデル (lm5 と名づけた）の予測値をプロットしたいとしよう。この場合、男女それぞれの労働経験年数に対応する賃金の予測値をプロットしたいのだから、プロットしたい性別 (gender) と労働経験年数 (experience) の組み合わせをすべて持つデータフレームをあらかじめ作ってやる必要がある。以下のスクリプトがそれである。
data.for.pred2 <- data.frame(
experience = rep(0 : 55, 2), # 0～55 の数値を２回繰り返すベクトル
gender = gl(2, 56, labels = c("male", "female") # male を 56回、次に female を 56回繰り返す因子
)
)
head(data.for.pred2)## experience gender
## 1 0 male
## 2 1 male
## 3 2 male
## 4 3 male
## 5 4 male
## 6 5 malepred2 <- predict(lm5, data.for.pred2)
plot(
x = data.for.pred2$ experience,
y = pred2,
pch = 19,
xlab = "労働経験年数",
ylab = "予測される賃金"
)
text(25, c(8.2, 11), c("女性", "男性")) # text() は図に文字を書き加えるための関数
説明変数はたくさんあるが、特に興味のある変数によって予測値がどう変化するかだけプロットしたい場合もある。しかし、pred() 関数には、必ずモデルに投入したすべての独立変数を含むデータフレームを指定してやる必要があるので少し工夫が必要である。以下の例を見よ。
lm6 <- lm(wage ~ age + I(age^2) + education + occupation, data = CPS1985)
data.for.pred3 <- data.frame(
age = 18 : 64, # 年齢だけ 18～64歳まで変化させる
education = mean(CPS1985$ education), # 学歴は平均教育年数で固定
occupation = "worker" # 職種も worker に固定
)
head(data.for.pred3)## age education occupation
## 1 18 13.01873 worker
## 2 19 13.01873 worker
## 3 20 13.01873 worker
## 4 21 13.01873 worker
## 5 22 13.01873 worker
## 6 23 13.01873 workerpred3 <- predict(lm6, data.for.pred3)
plot(18 : 64, pred3, type = "l", xlab = "年齢", ylab = "モデルから予測される賃金")
上の例は、年齢 (age) の効果を図示しようとしているが、モデルからの予測値を計算するためには年齢だけではなく、教育年数と職業の値も指定してやる必要がある。そこで、教育年数には教育年数の平均値を、職業は最頻値の “worker” を一律に指定している。
6.5 ANOVA 分散分析
分散分析は lm() で引数をカテゴリカル変数にすればできるが、まれに伝統的な分散分析表を出力したい場合がある。その場合、aov() という関数を使う。書式は lm とだいたい同じである。以下は使用例
aov1 <- aov(wage ~ gender, data = CPS1985)
summary(aov1)## Df Sum Sq Mean Sq F value Pr(>F)
## gender 1 594 593.7 23.43 1.7e-06 ***
## Residuals 532 13483 25.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.6 Model Diagnostics モデルの診断
回帰分析の際には、モデルの診断は重要である。モデルの診断とは、回帰分析が満たすべき条件を、データやモデルが本当に満たしているのか検討することである。論文などにはモデルの診断の結果は示されないことが多いので、軽視されがちであるが、誤った分析を避ける上で非常に重要である。ふつうは以下の点をチェックする。
6.6.1 モデルが正しく特定 (specify) されているか
モデルの特定とは、独立変数と従属変数にどのような関係があるのか特定（仮定）することである。例えば、本当は曲線的な関係なのに線形の関係を仮定していないか、本当は交互作用効果があるのにそれを見落としていないか、除去したり追加したりすべき独立変数はないか、従属変数は対数変換などする必要はないか、といった点を検討する。
6.6.2 残差は正規分布しているか
残差が正規分布していないと、回帰分析の結果は最尤推定値にならない。最尤推定値でなくても Best Linear Unbiased Estimation (BLUE) であることには変わりないので、必ずしも深刻に受け止める必要はないが、従属変数の対数変換などでモデルの適合度を下げずに正規分布に近づけることができる場合もあるので、その場合は対数変換などを検討すると良い。
6.6.3 残差分散は一定か
回帰分析では、残差の分散は独立変数の値とは独立で、一定であることが仮定されている。この仮定が満たされないと標準誤差を過小に推定してしまう。このような場合は、ロバスト標準誤差が計算されることが多い。
6.6.4 残差は独立か
時系列データの場合や地域データの場合、残差が相関していることがある。また、横断調査のデータであっても、同じグループに属する対象者が近い残差を取ることがあり（マルチレベル・データ）、このような場合は、残差が独立とは言えず、やはり標準誤差が過小に推定されてしまうことが多い。時系列データや地域データ、マルチレベル・データの分析はここでは取り扱わない。診断法についてもこれらの分析法の解説をあたれ。
6.6.5 多重共線性はないか
共線性とは、重回帰分析において、ある独立変数がその他の独立変数によって予測できてしまうことを言う。例えば、横断データで出生年と年齢を独立変数としてモデルに投入する場合、両者には共線性が生じる。なぜなら、\(\mathrm{調査年} - \mathrm{年齢} = \mathrm{出生年}\) という関係がほぼ成り立つので、誕生日が調査時点の前か後かによって多少のズレはあるが、出生年がわかれば調査時点での年齢もほぼわかる。このような場合、調査年と年齢の傾きの標準誤差は非常に大きくなり、推定結果も 1, 2 の事例を取り除いただけで大きく変化してしまうほど不安定になる。多重共線性とは、共線性が 3変数以上の変数のあいだに生じている状態である。例えば、夫婦二人だけの世帯の調査で、夫収入、妻収入、世帯収入を独立変数とした回帰分析をしたいとしよう。\(\mathrm{夫収入} + \mathrm{妻収入} = \mathrm{世帯収入}\)