R と RStudio でのデータ分析 入門
太郎丸 博
Thursday, April 09, 2015
データの分析には、
- データの読み込み（あるいは呼び出し）、
- 分析できるようにデータを整形、
- 探索的分析、最終的な分析結果の確定、
- 図や表の整形、レポート／論文の執筆
といったプロセスがある。これらを Windows 上の R と RStudio で行う方法を概説しよう。
1 read.csv データの読み込み
データはふつうEXCELのような表計算ソフトなどで行い、それを R のような統計解析ソフトに読み込むことが多い。分析の第一歩はデータを統計解析ソフトに読み込むことであるが、けっこう失敗しやすいので注意。
1.1 setwd 作業ディレクトリの指定
データを読み込む前に、作業ディレクトリ（作業フォルダともいう）を指定する。作業ディレクトリとは、Rがデフォルトでファイルの読み書きをするフォルダであり、特にパスを指定しなければ、R は作業フォルダのファイルを読みに行くし、ファイルを出力する場合もこのフォルダに出力のファイルを保存する。分析課題ごとに新しいフォルダを作ってそこに関連するファイルをまとめて保存しておくとよい。
作業フォルダを指定するには、RStudio のメニューの [Session] [Set Working Directory] [Choose Directory…] を選ぶとフォルダを選択するダイアログボックスが出てくるので、ここで選べばよい。
しかし、毎回これを行うのはけっこう面倒なので、 setwd() という関数をスクリプトの最初に書いてフォルダを指定したほうが簡単な場合も多い。 setwd() は
setwd(“作業ディレクトリのフルパス”)
という書式で指定する。例えば、以下のように。
setwd("C:/Users/Hiroshi/Dropbox/15_IJ_Courses/R_Practice")スクリプトはすべて半角で書くこと。また R は大文字と小文字を区別しているので、SETWDと書いてもエラーになる。また Windows の場合ディレクトリの区切りは “\” ではなく “/” を使うので注意。
作業ディレクトリの指定を誤って、保存したファイルがどこに行ったのか、わからなくなる人がたまにいるが、そういう場合は、以下のように引数なしで getwd() を使えば作業ディレクトリの場所がわかる。
getwd()## [1] "C:/Users/Hiroshi/Dropbox/15_IJ_Courses/R_Practice"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)## Loading required package: car## Warning: package 'car' was built under R version 3.1.3## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: sandwich
## Loading required package: survival
## Loading required package: splinesCPS1985## 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 = gl(2, 5, labels=c("Female", "Male")),
age = rep(18:22, 2), # スクリプトはコンマのあとなど途中で改行できる
score =c(46, 7, 88, 54, 66, 44, 59, 35, 22, 70)
)
dat1## sex age score
## 1 Female 18 46
## 2 Female 19 7
## 3 Female 20 88
## 4 Female 21 54
## 5 Female 22 66
## 6 Male 18 44
## 7 Male 19 59
## 8 Male 20 35
## 9 Male 21 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")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人分の身長と体重のデータを取ったが、二人目の身長と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 674 Plot グラフを作る
4.1 Histogram ヒストグラム
変数の分布を探索的にチェックする場合、グラフの作成が非常に有益である。連続変数の場合はヒストグラムをまずは作るのが基本であろう。ヒストグラムは hist()関数で作る。
hist(データフレーム名$変数名)
例えば以下のように。
hist(CPS1985$wage)
4.1.1 Bar Plot 棒グラフ
離散変数（文字型のデータ）の場合は度数分布表を単純にグラフにすればよい。そのためには、plot()関数で
plot(データフレーム名$変数名)
とすればよい。例えば以下のように。
plot(CPS1985$ethnicity)
上の１行２列目の図は、縦軸に wage 横軸に education をとっている。１行３列目の図は縦軸に wage 横軸に experience をとっている。つまり、その図の横にある文字が縦軸の変数の名前であり、その図の上か下にある文字が、横軸の変数名である。
4.1.2 Scaptter 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.1.3 Box Plot 箱ひげ図
グループ別に数値型の変数の分布を見たい場合は、箱ひげ図を使うのが基本である。箱ひげ図は、
plot(数値型の変数名 ~ グループ分けするためのカテゴリカル変数名, data=データフレーム名)
または
boxplot(数値型の変数名 ~ グループ分けするためのカテゴリカル変数名, data=データフレーム名)
とすれば作れる。
plot(wage ~ ethnicity, data=CPS1985)
4.1.4 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 というパッケージにも含まれており、car はデフォルトでロードされている）。これは以下のような書式で使う。
recode(ベクトル, 新しい値 <- 古い値（のベクトル）, 新しい値 <- 古い値（のベクトル）, … )
例えば、CPS1985 の occupation という変数は5つのカテゴリからなるが、これを upper, middle, lower の 3カテゴリにまとめなおしたいとしよう。その場合、以下のようにすればよい。
library(memisc)## Warning: package 'memisc' was built under R version 3.1.3## Loading required package: lattice
## Loading required package: MASS
##
## Attaching package: 'memisc'
##
## The following object is masked from 'package:car':
##
## recode
##
## The following objects are masked from 'package:stats':
##
## contr.sum, contr.treatment, contrasts
##
## The following object is masked from '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(~ 行の変数名 + 列の変数名, データ名)
連続変数をいくつかのカテゴリにまとめ直す場合、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)## 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
## 25 0 0 10 0
## 26 0 0 10 0
## 27 0 0 8 0
## 28 0 0 8 0
## 29 0 0 8 0
## 30 0 0 0 6
## 31 0 0 0 4
## 32 0 0 0 8
## 33 0 0 0 15
## 34 0 0 0 5
## 35 0 0 0 3
## 36 0 0 0 4
## 37 0 0 0 5
## 38 0 0 0 9
## 39 0 0 0 5
## 40 0 0 0 3
## 41 0 0 0 4
## 42 0 0 0 7
## 43 0 0 0 7
## 44 0 0 0 5
## 45 0 0 0 6
## 46 0 0 0 2
## 47 0 0 0 2
## 48 0 0 0 1
## 49 0 0 0 1
## 54 0 0 0 1
## 55 0 0 0 1この例の場合、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.3 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.3.1 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.4 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.4.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.2.5 mtable 複数の回帰分析の結果を表にまとめる
複数の回帰分析の結果をまとめて表にするには、memisc パッケージの mtable( ) 関数を使う。
mtable(回帰分析の結果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
## adj. R-squared 0.029 0.198 0.199
## sigma 5.063 4.604 4.599
## F 17.199 44.727 34.129
## p 0.000 0.000 0.000
## Log-likelihood -1622.810 -1571.049 -1570.009
## Deviance 13635.855 11232.849 11189.190
## AIC 3251.620 3152.098 3152.019
## BIC 3264.461 3173.500 3177.701
## N 534 534 534
## =======================================================6.3 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),
gender = gl(2, 56, labels = c("male", "female")
)
)
pred2 <- predict(lm5, data.for.pred2)
plot(data.for.pred2$experience, 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,
education=mean(CPS1985$education),
occupation="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” を一律に指定している。
7 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 ' ' 17.1 Model Diagnostics モデルの診断
回帰分析の際には、モデルの診断は重要である。モデルの診断とは、回帰分析が満たすべき条件を、データやモデルが本当に満たしているのか検討することである。論文などにはモデルの診断の結果は示されないことが多いので、軽視されがちであるが、誤った分析を避ける上で非常に重要である。