R
R
Embedded Files
R.pdfRstudio.pdf用語
データフレーム：RではCSVファイルなどのデータセットを読み込んで分析しますが、読み込んだデータセットをデータフレームと呼びます。このデータフレームは元のCSVファイルと同じ内容ですが、R専用の別ファイルとなっています。ですので、データフレームを加工しても元のCSVファイルは影響を受けません。
NA：Not Available。欠損値とも呼び、データがないことを意味しています。分析によって、NAのデータを持つ対象は自動的に分析に使われなかったり、エラーになります。
c( )：各コマンドで頻繁に出てくる関数です。基本的には、いくつかの数値や要素をベクトルやリストという一つの単位にまとめる機能を持ちます。combine values into a vector or listのcのようです。
注意
ここに掲載のコマンドは授業の補足のためのメモ書きであり、説明不足な部分が多いです。ご容赦ください。
以下では、出力を貼り付けているため > が表示されていますが入力しないでください。これはコマンド部分を示しているだけです。
欠損値のあるデータで行ごとに平均値を計算する
欠損値（NA）のある3つの変数からなるデータフレームを作成します。
> x1 <- c(1,NA,4,5,6)
> x2 <- c(2,NA,4,6,NA)
> x3 <- c(6,NA,3,7,8)
> data <- data.frame(x1,x2,x3)
> data
x1 x2 x3
1 1 2 6
2 NA NA NA
3 4 4 3
4 5 6 7
5 6 NA 8
行ごとに平均値を計算しようと以下を実行すると
> library(dplyr)
> data <- data %>%
+ mutate(mean1 = mean(c(x1,x2,x3),na.rm=TRUE))
> data
x1 x2 x3 mean1
1 1 2 6 4.727273
2 NA NA NA 4.727273
3 4 4 3 4.727273
4 5 6 7 4.727273
5 6 NA 8 4.727273
すべての数値から計算した平均値が1つだけ作成されます。na.rm=TRUEは欠損値を無視して計算するための命令です。
そこでrowwise()を使います。
> data <- data %>%
+ rowwise() %>%
+ mutate(mean2 = mean(c(x1,x2,x3),na.rm=TRUE))
> data
# A tibble: 5 × 5
# Rowwise:
x1 x2 x3 mean1 mean2
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1 2 6 4.73 3
2 NA NA NA 4.73 NaN
3 4 4 3 4.73 3.67
4 5 6 7 4.73 6
5 6 NA 8 4.73 7
データはtibble形式になっています。データフレームに戻したいときは以下のようにします。
> data <- data.frame(data)
> data
x1 x2 x3 mean1 mean2
1 1 2 6 4.727273 3.000000
2 NA NA NA 4.727273 NaN
3 4 4 3 4.727273 3.666667
4 5 6 7 4.727273 6.000000
5 6 NA 8 4.727273 7.000000
列（縦）方向のデータの結合
以下の列数が異なる2つのデータフレームがあったとします。
> data1 <- data.frame(id = 1:2, age = c(25,32), edu = c(1,5), sex=c(1,2))
> data1
id age edu sex
1 1 25 1 1
2 2 32 5 2
> data2 <- data.frame(id = 3:4, age = c(37,29), edu = c(3,4))
> data2
id age edu
1 3 37 3
2 4 29 4
rbindで列方向に足してみます。
> data3 <- rbind(data1,data2)
Error in rbind(deparse.level, ...) :
numbers of columns of arguments do not match
列数が異なるため、rbindではエラーになります。そこでパッケージdplyrを読み込み、bind_roｗsを使用します。
> library(dplyr)
> data3 <- bind_rows(data1,data2)
> data3
id age edu sex
1 1 25 1 1
2 2 32 5 2
3 3 37 3 NA
4 4 29 4 NA
行（横）方向のデータの結合
以下の2つのデータフレームがあったとします。
> data1 <- data.frame(id = 1:4, age = c(25,32,37,29), edu = c(1,5,3,4))
> data1
id age edu
1 1 25 1
2 2 32 5
3 3 37 3
4 4 29 4
> data2 <- data.frame(id = 1:3, sex = c(1,2,1))
> data2
id sex
1 1 1
2 2 2
3 3 1
変数idをキーにしてデータを結合します。
> data3 <- merge(data1,data2,by="id")
> data3
id age edu sex
1 1 25 1 1
2 2 32 5 2
3 3 37 3 1
data2にはid=4の人がいないので結合データから落ちています。id=4の人をデータに残したい場合は、以下のようにオプションをつけます。
> data4 <- merge(data1,data2,by="id",all=TRUE)
> data4
id age edu sex
1 1 25 1 1
2 2 32 5 2
3 3 37 3 1
4 4 29 4 NA
行数がそろっていればcbindも利用可能です。
パネルデータで個体数、個体ごとの観測数の確認
パネルデータ（data1）をRstudioのEnvironmentで見ると総観測数（個体数✕Wave数）がわかりますが、個体数は表示されません。Environmentに表示されているdata1の「◯◯obs.」の◯◯が総観測数です。そこで、以下のように個体（id）ごとにまとめたdata2をつくります。
> library(dplyr)
> data2 <- data1 %>%
+ group_by(id) %>%
+ summarize(count=n())
data2をEnvironmentで見ると個体数がわかりますし、データ自体を見ると個体ごとの観測数がわかります。
forによる繰り返し処理
iを1から10まで変化させながらvarという文字に、その数値をpasteしてprintしてみました。
> for(i in 1:10){
+ j <- paste("var",i,sep="")
+ print(j)
+ }
[1] "var1"
[1] "var2"
[1] "var3"
[1] "var4"
[1] "var5"
[1] "var6"
[1] "var7"
[1] "var8"
[1] "var9"
[1] "var10"
これを使ってtableなどの処理を自動で行いたいのですが、現状うまくいっていません。なお、sep="-"とするとvar-1のような文字ができます。
データフレームの特定の列を抜き出して新しいデータフレームをつくる
以下のようにv1からv10まで10個の変数のあるデータフレーム（data1）があるとします。
> data1
v1 v2 v3 v4 v5 v6 v7 v8 v9 v10
1 27 59 7 48 73 58 0 88 80 29
2 90 62 56 91 51 22 38 39 78 19
3 61 80 20 8 44 11 3 8 25 16
4 65 32 39 75 93 78 81 82 25 21
5 16 26 91 66 13 73 19 7 15 6
6 38 43 90 25 58 43 57 45 33 34
7 28 6 95 28 2 70 52 21 32 59
8 95 78 59 57 56 56 29 48 44 81
9 59 42 40 4 97 65 54 95 49 41
1〜3列、6列、9〜10列だけ抜き出してdata2をつくります。
> data2 <- data1[c(1:3,6,9:10)]
> data2
v1 v2 v3 v6 v9 v10
1 27 59 7 58 80 29
2 90 62 56 22 78 19
3 61 80 20 11 25 16
4 65 32 39 78 25 21
5 16 26 91 73 15 6
6 38 43 90 43 33 34
7 28 6 95 70 32 59
8 95 78 59 56 44 81
9 59 42 40 65 49 41
一行目に変数名のないデータを読み込む
一行目に変数名（id、waveなど）がないデータを読み込む場合には、以下のようにヘッダーをFにします。
> data1 <- read.csv("test.csv", header=F)
変数名は自動的に左からV1、V2のようにつけられます。
推定に使用した対象と変数だけで新たにデータフレームを作成する。
推定に使用した対象で記述統計を計算したい時がありますが、model.frameという関数を使うと、目的の対象だけ残したデータフレームができるようです。まだ使いこなしていないので、とりあえずメモ程度ですが、yとxに欠損値のあるdata1があったとします。
> data1
id y x w
1 1 21 3 1
2 2 34 5 0
3 3 15 1 1
4 4 NA 2 0
5 5 35 4 1
6 6 22 5 0
7 7 18 2 0
8 8 17 NA 0
9 9 31 6 1
10 10 28 4 1
そこでy~xの回帰を行い、使用された対象と変数だけで新たなデータフレームdata2を作成します。
> data2 <- model.frame(y~x, data=data1)
> data2
y x
1 21 3
2 34 5
3 15 1
5 35 4
6 22 5
7 18 2
9 31 6
10 28 4
idとwは推定に使用されなかったので、data2からは落ちています。このデータフレームを使って記述統計を計算すれば良いでしょう。
パネルデータでidごとに特定のwaveの値をすべてのwaveにコピーする。
パネルデータでは、wave1の調査でだけ性別など変わりにくい属性を質問することがあり、wave2以降でデータが欠損していることがあります。ここでは、wave１の値を他のすべてのwaveにコピーする方法を説明します。練習用に個体番号（id）、時点（wave）、変数（x1）からなるdata1を作成します。
> data1 <- data.frame(id = rep(c(1:3), each = 3), wave = c(1:3), x1 =c(1:9))
> data1
id wave x1
1 1 1 1
2 1 2 2
3 1 3 3
4 2 1 4
5 2 2 5
6 2 3 6
7 3 1 7
8 3 2 8
9 3 3 9
x1について、wave1時点の情報だけを残し、他のwaveの値を-9にした、x1aをdata1に追加します。
> data1$x1a <- ifelse(data1$wave==1, data1$x1, -9)
> data1
id wave x1 x1a
1 1 1 1 1
2 1 2 2 -9
3 1 3 3 -9
4 2 1 4 4
5 2 2 5 -9
6 2 3 6 -9
7 3 1 7 7
8 3 2 8 -9
9 3 3 9 -9
パイプ演算子（%>%）を使うために、パッケージdplyrを読み込みます。
> library(dplyr)
idごとにx1aの最大値（max）をとる新たな変数x1bを作成し、data1に追加します。
> data1 <- data1 %>% group_by(id) %>% mutate(x1b = max(x1a))
data1を見ると、x1bができているのが確認できます。
> data1
# A tibble: 9 × 5
# Groups: id [3]
id wave x1 x1a x1b
<int> <int> <int> <dbl> <dbl>
1 1 1 1 1 1
2 1 2 2 -9 1
3 1 3 3 -9 1
4 2 1 4 4 4
5 2 2 5 -9 4
6 2 3 6 -9 4
7 3 1 7 7 7
8 3 2 8 -9 7
9 3 3 9 -9 7
パネルデータでラグ付き変数を作成する。
練習用に個体番号（id）、時点（wave）、変数（x1）からなるdata1を作成します。
> data1 <- data.frame(id = rep(c(1:3), each = 3), wave = c(1:3), x1 =c(1:9))
> data1
id wave x1
1 1 1 1
2 1 2 2
3 1 3 3
4 2 1 4
5 2 2 5
6 2 3 6
7 3 1 7
8 3 2 8
9 3 3 9
パイプ演算子（%>%）を使うために、パッケージdplyrを読み込みます。
> library(dplyr)
2期前の値をとるx2を作成し、data1に追加します。
> data1 <- data1 %>%
+ group_by(id) %>%
+ mutate(x2 = lag(x1, n = 2, default = NA))
data1を見ると、x2ができているのが確認できます。
> data1
# A tibble: 9 × 4
# Groups: id [3]
id wave x1 x2
<int> <int> <int> <int>
1 1 1 1 NA
2 1 2 2 NA
3 1 3 3 1
4 2 1 4 NA
5 2 2 5 NA
6 2 3 6 4
7 3 1 7 NA
8 3 2 8 NA
9 3 3 9 7
場合分けで計算式を変えて新しい変数を作成する。
生年月（ybirth, mbirth）から年齢を計算するとき、生まれ月で年齢を変えたいときがあります。例えば、2020年5月に行われた調査データを使う時、4月生まれ以前と5月生まれ以降で年齢は異なるでしょう。ここでは、5月生まれは誕生日を迎えている人と迎えていない人が混在していますが、ここでは迎えていないとみなして計算してみます。パッケージdplyrとmutate関数を使います。
> library(dplyr)
> data1 <- data1%>%
+ mutate(age = ifelse(mbirth<=4, 2020-ybirth, 2020-ybirth-1))
> data1
id ybirth mbirth age
1 1 1980 1 40
2 2 1980 2 40
3 3 1980 3 40
4 4 1980 4 40
5 5 1980 5 39
6 6 1980 6 39
7 7 1980 7 39
8 8 1980 8 39
9 9 1980 9 39
10 10 1980 10 39
stargazerで複数の推定結果を比較する。
変数の組み合わせでいくつかのモデルを推定して結果を比較したいときがあります。その場合、stargazerを使用するのが便利です。例えば、以下のように2つの結果をreg1、reg2に格納します。
> reg1 <- lm(her~hinc+stu, data=data1)
> reg2 <- lm(her~hinc+stu+univ, data=data1)
パッケージstargazerを読み込んで、2つの結果を1つの表にまとめます。
> library(stargazer)
> stargazer(reg1,reg2,type="text")
=================================================================
Dependent variable:
---------------------------------------------
her
(1) (2)
-----------------------------------------------------------------
hinc 0.011 0.019
(0.012) (0.012)
stu 3.082*** 2.780***
(0.429) (0.411)
univ 8.685***
(3.003)
Constant 4.613 -0.511
(8.241) (7.830)
-----------------------------------------------------------------
Observations 47 47
R2 0.551 0.624
Adjusted R2 0.531 0.598
Residual Std. Error 4.516 (df = 44) 4.180 (df = 43)
F Statistic 27.024*** (df = 2; 44) 23.819*** (df = 3; 43)
=================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
なお、初期設定では有意水準は10％が下限となっていますが、以下のように指定することで有意水準を変更できます。これは下限を5％にする場合です。
> stargazer(reg1,reg2,type="text",star.cutoffs=c(0.05,0.01,0.001))
条件にあった対象のみで度数分布を確認する
あるデータフレーム（data1）で、男性（sex==1）だけで、配偶状態（mar）の度数分布を確認するには、以下のようにします。
> table(subset(data1,sex==1,c(mar)))
mar
1 2 4
337 156 6
女性（sex==2）だけで、最終学歴（educ）の度数分布を確認するには、以下のようにします。
> table(subset(data1,sex==2,c(edu)))
edu
1 2 3 4 5 6 9
2 118 109 118 139 13 2
いくつかの自由度（k）で見たカイ2乗分布のグラフ
自由度が大きくなると分布が右へ移るのがわかります。
> curve(dchisq(x,1),0,30,col="red",ylab="f(x)",family="serif")
> curve(dchisq(x,2),0,30,add=TRUE,col="blue")
> curve(dchisq(x,5),0,30,add=TRUE,col="green")
> curve(dchisq(x,10),0,30,add=TRUE)
> abline(h=0)
> text(2,0.6,labels="k=1",family="serif")
> text(2.6,0.29,labels="k=2",family="serif")
> text(4,0.18,labels="k=5",family="serif")
> text(10,0.12,labels="k=10",family="serif")
平均は異なるが等分散の正規分布のグラフ
平均0、分散1の正規分布（青）と平均1、分散1の正規分布（緑）のグラフを表示する。
> curve(dnorm(x,0,1),-6,6, ylim=c(0,0.5), ylab="", col="blue")
> curve(dnorm(x,1,1),-6,6, ylim=c(0,0.5), ylab="", col="green", add=TRUE)
> segments(qnorm(0.5,0,1), 0, qnorm(0.5,0,1), dnorm(qnorm(0.5,0,1),0,1))
> segments(qnorm(0.5,1,1), 0, qnorm(0.5,1,1), dnorm(qnorm(0.5,1,1),1,1))
> abline(h=0)
Macでグラフを作成する際、ラベルに漢字を表示させたいとき
Macでグラフに日本語を表示させたいときに□□□のようになってしまうことがあります。この場合、以下のようにpar()で使用フォントを設定すれば表示できるようになります。フォントはいくつか選択肢があります（HiraginoSans-W3, HiraMaruProN-W4, HirakakuProN-W3）。日本語以外の数字の部分も同じフォントになります。
> par(family= "HiraKakuProN-W3")
> plot(data$inc,data1$col,xlab="県民所得（千円）", ylab="大学進学率（％）")
時間データを分に変換する
時間データで「時間.分」という形になっていることがあります。例えば2.54は2時間54分のことであり、2.54時間（約2時間30分）のことを意味しているのではありません。また、2.60~2.99にはデータも存在しないので、そのまま使用するとおかしな散布図になります。そこで、分に換算して分析に用いることになります。以下のように関数を利用して換算します。
> floor(2.54)*60+(2.54%%1)*100
[1] 174
floor(x)はxを超えない整数を返す関数で、y%%zの%%は、ｙをｚで割り、その余りを得る演算子です。
データフレームの指定部分で記述統計量を計算
一部の変数だけで記述統計を計算した時があります。ここではapply関数を利用して、複数の変数の標準偏差を一度に計算してみます。
> data1
pref her hinc univ stu popd rcpi
1 北海道 43.3 525.9 0.69 12.32 240.5 100.3
2 青森 43.6 433.3 0.76 12.09 405.1 99.5
3 岩手 44.2 500.2 0.39 11.44 344.5 99.6
:
> apply(data1[,2:3],2,sd)
her hinc
6.593639 53.947523
apply関数の第1引数でデータフレームdata1の2列目から3列目を指定しています。[行,列]で場所をしています。第2引数で列単位で計算するよう指定しています。第3引数で計算する統計量をsd（標準偏差）に指定しています。
次に、データフレームの離れた複数箇所で記述統計を計算したい場合、例えば、2～3列目と5列目の記述統計を計算する場合は以下のようにします。
> summary(data1[,c(2:3,5)])
her hinc stu
Min. :39.20 Min. :395.8 Min. : 9.48
1st Qu.:45.40 1st Qu.:493.1 1st Qu.:12.18
Median :50.50 Median :533.2 Median :13.13
Mean :51.09 Mean :526.1 Mean :13.21
3rd Qu.:55.55 3rd Qu.:557.2 3rd Qu.:13.87
Max. :66.50 Max. :631.5 Max. :16.47
> apply(data1[,c(2:3,5)],2,sd)
her hinc stu
6.593639 53.947523 1.558138
変数名を指定してもできます。
> summary(data1[,c("her","hinc","stu")])
her hinc stu
Min. :39.20 Min. :395.8 Min. : 9.48
1st Qu.:45.40 1st Qu.:493.1 1st Qu.:12.18
Median :50.50 Median :533.2 Median :13.13
Mean :51.09 Mean :526.1 Mean :13.21
3rd Qu.:55.55 3rd Qu.:557.2 3rd Qu.:13.87
Max. :66.50 Max. :631.5 Max. :16.47
> apply(data1[,c("her","hinc","stu")],2,sd)
her hinc stu
6.593639 53.947523 1.558138
なお、データフレームの2～3列目と5列目だけコンソールに表示させたい時は以下のとおりです。
> data1[,c(2:3,5)]
her hinc stu
1 43.3 525.9 12.32
2 43.6 433.3 12.09
3 44.2 500.2 11.44
:
複数回答のデータからある選択肢該当ダミーを作成する
複数回答（あてはまるものに○系の設問）のデータで、選んだ選択肢の番号が1変数に1つずつ入っていることがあります。例えば、転職の理由として該当するもの3つに〇をするような場合です。以下は、5人が8つの選択肢から選んだものがq1_1~q1_3に入っています。1つだけ選んだ人もいれば3つ選んだ人もいます。
> data1
q1_1 q1_2 q1_3
1 1 2 3
2 3 5 .
3 2 . .
4 1 4 8
5 2 3 .
ここから、選択肢3を選んだというダミー変数（q1_c3）を作成したいとしましょう。選択肢3を選んだ人を1とする場合、以下のようにコマンドします。「|」は「または（or）」の演算子です。
> data1$q1_c3[data1$q1_1 == 3 | data1$q1_2 ==3 | data1$q1_3 ==3 ] <- 1
> data1
q1_1 q1_2 q1_3 q1_c3
1 1 2 3 1
2 3 5 . 1
3 2 . . NA
4 1 4 8 NA
5 2 3 . 1
該当しない人はNAになっているので0にリコードしておきましょう。パッケージcarを読み込んでから、リコードします。
> library(car)
> data1$q1_c3 <- recode(data1$q1_c3, '1=1;NA=0')
> data1
q1_1 q1_2 q1_3 q1_c3
1 1 2 3 1
2 3 5 . 1
3 2 . . 0
4 1 4 8 0
5 2 3 . 1
複数の変数から1つの新しい変数をつくる
2つの変数の回答の組み合わせで、新しい変数をつくりたい場合があります。以下のようなx1とx2の2つの変数があったとしましょう。
> data1
x1 x2
1 1 1
2 2 1
3 1 2
4 2 2
5 2 1
例えばx1が1かつx2が1の場合に1、x1が1かつx2が2の場合に2、x1が2かつx2が1の場合に3、x1が2かつx2が2の場合に4、となる新しい変数（x3）を作成するには以下のようにコマンドします。
> data1$x3[data1$x1==1 & data1$x2==1] <- 1
> data1$x3[data1$x1==1 & data1$x2==2] <- 2
> data1$x3[data1$x1==2 & data1$x2==1] <- 3
> data1$x3[data1$x1==2 & data1$x2==2] <- 4
> data1
x1 x2 x3
1 1 1 1
2 2 1 3
3 1 2 2
4 2 2 4
5 2 1 3
datファイルの読み込み
データファイルがxxx.datのように拡張子がdatになっている場合、以下のようにコマンドして読み込みます。
> data1 <- read.delim("xxx.dat")
なお、datファイルではデータがタブで区切られています。
2次関数による推定
独立変数に2次項を入れて推定したい場合、2次項を新たな変数として作成してモデルに入れます。例えば年齢（age）の2次項を使った推定をしたい場合は以下のようにします。
> data1$age2 <- data1$age^2
> lm(inc~age+age2, data=data1)
あるいは、I()関数を使うと2次項をデータフレームに追加することなく推定できます。
> lm(inc~age+I(age^2), data=data1)
Rの出力をExcelに貼り付け
Rの出力をそのままExcelに貼り付けると、図①のように変数名や係数、標準誤差などが1セルにまとめて入ってしまいます。これらの数値を1セルごとに入れたい場合、「データ」タブの「区切り位置」を押します。
すると、図②のような画面になりますので、「固定長―スペースに･･･」が選択されていることを確認して「次へ」を押します。
図③のように、どこで区切るか表示されますので、よければ「完了」を押してください。「次へ」を押すと表示形式を変更できますが、ここでは必要なさそうです。
図④のようにきれいにセルに入りました。
図①
図②
図③
図④
加工した数値変数が正しく作成されたか並べて確認
変数incを10倍してinc10を、incを100倍してinc100を作成した後、うまくできているか確認するために並べたいことがあります。以下のように実行します。新たな変数の作成方法として、＄でデータフレームをした場合とwith関数を使った2つの方法が示してありますが、どちらでもかまいません。with関数は使用するデータフレームをこの命令でだけ指定できます。attach()の場合はdetach()するまでデータフレームを指定し続けるのと対照的です。
> data1
x1 x2
1 1 1
2 2 1
3 1 2
4 2 2
5 2 1
> data1$inc10 <- data1$inc*10
> data1$inc100 <- with(data1,inc*100)
> data1[,c("inc","inc10","inc100")]
inc inc10 inc100
1 2476 24760 247600
2 2333 23330 233300
3 2309 23090 230900
:
因子の場合はクロス表でも確認できます。
サンプルサイズの確認
順序回帰モデルの推定などでは、使用したデータ数が表示されないことがあります。結果としては保存されているので表示することができます。例えば推定結果がOrdRegModel.1に格納されていたとします。この場合、以下のように入力し実行します。
> OrdRegModel.1$nobs
あるいは
> OrdRegModel.1$n
です。stagazerでも確認できます。
条件に合うデータだけ抜き出してデータフレームをつくる
まずはテストデータを作成します。
> w <- c(0,1,1,0,1,1,0,1)
> x <- c(6,8,10,5,7,8,4,3)
> y <- c(120,320,430,450,340,210,270,410)
> z <- c(1,2,3,2,3,3,2,1)
> data <- data.frame(w,x,y,z)
> data
w x y z
1 0 6 120 1
2 1 8 320 2
3 1 10 430 3
4 0 5 450 2
5 1 7 340 3
6 1 8 210 3
7 0 4 270 2
8 1 3 410 1
subset関数で条件を指定（z==3）し、抜き出す変数（x,y,z）を指定します。
> data1 <- subset(data,z==3,c(x,y,z))
> data1
x y z
3 10 430 3
5 7 340 3
6 8 210 3
zが3のデータだけ抜き出すことができました。ただし、attachして抜き出そうとしてもエラーになるようです。
> attach(data)
> data2 <- subset(z==3,c(x,y,z))
Error in subset.default(z == 3, c(x, y, z)) : 'subset' must be logical
> data2
Error: object 'data2' not found
データの四捨五入
小数部分のある変数incを、四捨五入で整数にして、新しい変数incrを作る場合は以下のようにコマンドします。
> data1$incr <- round(data1$inc,digits=0)
また、四捨五入してもとのincに上書きしてしまう場合は以下のようにコマンドします。
> data1$inc <- round(data1$inc,digits=0)
ベクトルをデータフレームに
以下のように2つのベクトルがあるとします。
> x <- c(1:10)
> y <- c(1,8,15,18,22,23,22,24,26,24)
これらを1つのデータフレームにまとめるには以下のようにコマンドします。
> data <- data.frame(var1=x,var2=y)
var1 var2
1 1 1
2 2 8
3 3 15
4 4 18
5 5 22
6 6 23
7 7 22
8 8 24
9 9 26
10 10 24
var1とvar2の変数名は任意です。
散布図のラベル表示
①都道府県データを使って通常の散布図を描きましたが、どのマークがどの都道府県なのかわかりません。
> plot(data1$univ,data1$her)
②そこで、データの並び番号で表示してみます。1行目のtype="n"はマークを表示しない命令です。2行目でデータの並び番号を表示させています。
> plot(data1$univ,data1$her,type="n")
> text(data1$univ,data1$her)
③番号ではどの都道府県かわからないので、都道府県名を表示したい場合は以下のようにします。1行目はデータフレームdata1の1列目にある都道府県名を抜き出し、labelという名前で保存しています。3行目ではMacのために日本語フォントの設定をしています。Windowsでは無視してください。4行目でlabelを使用して都道府県名を散布図内に表示しています。
> label <- data1[,1]
> plot(data1$univ,data1$her,type="n")
> par(family= "HiraKakuProN-W3")
> text(data1$univ,data1$her,label)
①通常の散布図
②都道府県番号を示した散布図
③都道府県名を示した散布図
変数の値にラベルをつける
テストデータをつくります。
> f1 <- c(0,1,1,0,1,1,0,1)
> q1 <- c(1,2,3,2,3,3,2,1)
> data <- data.frame(f1,q1)
> data
f1 q1
1 0 1
2 1 2
3 1 3
4 0 2
5 1 3
6 1 3
7 0 2
8 1 1
f1が性別、q1が何かに対する満足度だとして、以下のようにラベルをつけて、新変数として保存します。
> data$sex <- factor(data$f1, levels=c("0","1"), labels=c("女性","男性"))
> data$sat <- factor(data$q1, levels=c("1","2","3"), labels=c("満足","どちらでもない","不満足"))
以下のようにラベルのついた新変数ができました。
> data
f1 q1 sex sat
1 0 1 女性 満足
2 1 2 男性 どちらでもない
3 1 3 男性 不満足
4 0 2 女性 どちらでもない
5 1 3 男性 不満足
6 1 3 男性 不満足
7 0 2 女性 どちらでもない
8 1 1 男性 満足
この新変数を使うとクロス表もわかりやすいでしょう。
> table(data$sat)
満足 どちらでもない 不満足
2 3 3
> table(data$sex,data$sat)
満足 どちらでもない 不満足
女性 1 2 0
男性 1 1 3
対数をとる
10の自然対数をとるには
> log(10)
[1] 2.302585
底が2の場合
> log2(10)
[1] 3.321928
底が10の場合
> log10(10)
[1] 1
ある変数x1の自然対数をとってln_x1として保存する場合
> ln_x1 <- log(x1)
散布図
散布図でポイントを線でつなぐにはtypeで指定します。まずテストデータを生成します。
> x <- 0:10
> y <- x^2-5*x+5
> x;y
[1] 0 1 2 3 4 5 6 7 8 9 10
[1] 5 1 -1 -1 1 5 11 19 29 41 55
いくつかのタイプを示します。
> plot(x, y, type="b") #①線とポイント、重なりなし
> plot(x,y,type="l") #②線のみ
> plot(x,y,type="o") #③線とポイント、重なりあり
> plot(x,y,type="c") #④線のみ、ポイント部分空白
①線とポイント、重なりなし
②線のみ
③線とポイント、重なりあり
④線のみ、ポイント部分空白
データのソート
以下のようなデータフレームdata1があるとします。
> data1
pref her hinc
1 北海道 43.3 525.9
2 青森 43.6 433.3
3 岩手 44.2 500.2
:
45 宮崎 45.1 449.1
46 鹿児島 42.7 545.1
47 沖縄 39.2 427.6
これを変数herの昇順に並べ替えます。まずorder関数で順を作成しsortに格納します。
> sort <- order(data1$her)
> sort
[1] 47 35 46 41 1 31 2 3 5 42 6 45 7 43 15 44 32 39 20 4 30 8 24 33 37 36 9 16 38 10 22 40 17 25 21
[36] 12 18 19 11 23 29 34 27 28 14 26 13
sortには、上のように当初の都道府県の並び番号でherの順が示されます。この結果を使って並べ替えたデータフレームを新たに作成し、表示します。
> data1a <- data1[sort,]
> data1a
pref her hinc
47 沖縄 39.2 427.6
35 山口 42.7 576.4
46 鹿児島 42.7 545.1
:
14 神奈川 61.5 513.9
26 京都 66.4 495.3
13 東京 66.5 560.6
グラフを並べて表示する
デフォルトではグラフは1つだけ表示されます。並べて表示したい時はparのmfrowで指定します。横に２つ（１行で２列）並べる場合は以下のとおりです（左図）。
> par(mfrow=c(1,2))
> hist(data1$tfr80)
> hist(data1$tfr00)
縦に2つ（2行で1列）並べたいときは以下のようにコマンドします。ここでは、横軸の幅をxlimを使ってそろえてみました（右図）。
> par(mfrow=c(2,1))
> hist(data1$tfr80, xlim=range(1.0,2.5))
> hist(data1$tfr00, xlim=range(1.0,2.5))
横に並べた場合
縦に並べた場合
基本統計量
summaryでデータフレームを指定すると全変数の平均値などが表示されます。文字型変数（pref）については数値は出ません。
> summary(data1)
pref flp80 tfr80 flp00 tfr00
Length:47 Min. :39.10 Min. :1.440 Min. :50.30 Min. :1.070
Class :character 1st Qu.:49.15 1st Qu.:1.760 1st Qu.:57.50 1st Qu.:1.425
Mode :character Median :53.30 Median :1.820 Median :59.20 Median :1.470
Mean :53.37 Mean :1.829 Mean :59.79 Mean :1.473
3rd Qu.:58.40 3rd Qu.:1.885 3rd Qu.:62.75 3rd Qu.:1.560
Max. :67.00 Max. :2.380 Max. :68.20 Max. :1.820
csvファイルの書き出し
変数の加工などでデータフレームに修正を加えた場合、そのまま保存することがあります。一つの方法としてcsvファイルとして作業フォルダに書き出す命令があります。
以下は"data1"というデータフレームを"data2"というcsvファイルとして作業フォルダに書き出すコマンドです。
> write.csv(data1,file="data2.csv")
ただ、このままだと行番号を一つの変数として認識してしまうので
> write.csv(data1,file="data2.csv",row.names=FALSE)
とします。
ダミー変数2
東京か三重の場合に1、それ以外を0とするダミー変数をつくる場合、以下のようにします。ここでは都道府県名は文字型変数になっています。
> data1$tkymie1 <- data1$pref=="東京" | data1$pref=="三重"
データフレーム上、TRUE/FALSEになっていますが、1/0ダミーとして回帰分析に使用できます。1/0で表示したい場合はas.integerを使ってください。
> data1$tkymie2 <- as.integer(data1$pref=="東京" | data1$pref=="三重")
> data1[,c(1,8:9)]
pref tkymie1 tkymie2
:
12 千葉 FALSE 0
13 東京 TRUE 1
14 神奈川 FALSE 0
:
23 愛知 FALSE 0
24 三重 TRUE 1
25 滋賀 FALSE 0
データの一部を使って回帰
データの一部を使って分析したい時があります。例えばロング形式のパネルデータに1980年と2000年のデータがあったとして、1980年と2000年を別々に推定したい場合です。個票の分析では、男女別に推定したい場合もあります。以下のように、subsetで使用する年（year）を指定して回帰してみました。
> reg1 <- lm(tfr~flp, data=data1, subset=year==1980)
> reg2 <- lm(tfr~flp, data=data1, subset=year==2000)
> library(stargazer)
> stargazer(reg1,reg2,type="text")
==========================================================
Dependent variable:
----------------------------
tfr
(1) (2)
----------------------------------------------------------
flp 0.007** 0.018***
(0.003) (0.004)
Constant 1.441*** 0.383
(0.152) (0.248)
----------------------------------------------------------
Observations 47 47
R2 0.128 0.302
Adjusted R2 0.108 0.286
Residual Std. Error (df = 45) 0.127 0.112
F Statistic (df = 1; 45) 6.596** 19.432***
==========================================================
Note: *p<0.1; **p<0.05; ***p<0.01
時系列データ
データフレームdata1にある変数tfrを、1950年スタートの1年ごとの時系列データとして認識させます。
> tfr <- ts(data1$tfr,start=c(1950),frequency=1)
> tfr
Time Series:
Start = 1950
End = 2022
Frequency = 1
[1] 3.65 3.26 2.98 2.69 2.48 2.37 2.22 2.04 2.11 2.04 2.00 1.96 1.98 2.00 2.05 2.14 1.58 2.23 2.13 2.13 2.13
[22] 2.16 2.14 2.14 2.05 1.91 1.85 1.80 1.79 1.77 1.75 1.74 1.77 1.80 1.81 1.76 1.72 1.69 1.66 1.57 1.54 1.53
[43] 1.50 1.46 1.50 1.42 1.43 1.39 1.38 1.34 1.36 1.33 1.32 1.29 1.29 1.26 1.32 1.34 1.37 1.37 1.39 1.39 1.41
[64] 1.43 1.42 1.45 1.44 1.43 1.42 1.36 1.33 1.30 1.26
このデータを時系列プロットします（左図）。
> plot(tfr)
次に1階の階差をとってみます。diff関数を使います。
> tfr_d1 <- diff(tfr)
> tfr_d1
Time Series:
Start = 1951
End = 2022
Frequency = 1
[1] -0.39 -0.28 -0.29 -0.21 -0.11 -0.15 -0.18 0.07 -0.07 -0.04 -0.04 0.02 0.02 0.05 0.09 -0.56 0.65
[18] -0.10 0.00 0.00 0.03 -0.02 0.00 -0.09 -0.14 -0.06 -0.05 -0.01 -0.02 -0.02 -0.01 0.03 0.03 0.01
[35] -0.05 -0.04 -0.03 -0.03 -0.09 -0.03 -0.01 -0.03 -0.04 0.04 -0.08 0.01 -0.04 -0.01 -0.04 0.02 -0.03
[52] -0.01 -0.03 0.00 -0.03 0.06 0.02 0.03 0.00 0.02 0.00 0.02 0.02 -0.01 0.03 -0.01 -0.01 -0.01
[69] -0.06 -0.03 -0.03 -0.04
階差データをプロットします（右図）。
> plot(tfr_d1)
合計特殊出生率の推移
合計特殊出生率の1階階差の推移
交差項を使った回帰
交差項をデータとしてあらかじめ作成する必要はなく、*で結べば式に記述すれば単独項も交差項も両方回帰式に登場します。
> reg1 <- lm(her~hinc*univ, data=data1)
交差項のみを使いたい場合は:で結びます。
> reg2 <- lm(her~hinc:univ,data=data1)
> library(stargazer)
> stargazer(reg1,reg2,type="text")
===============================================================
Dependent variable:
-------------------------------------------
her
(1) (2)
---------------------------------------------------------------
hinc 0.023
(0.061)
univ 7.092
(53.072)
hinc:univ 0.013 0.026***
(0.099) (0.008)
Constant 30.958 43.146***
(33.051) (2.463)
---------------------------------------------------------------
Observations 47 47
R2 0.225 0.208
Adjusted R2 0.170 0.191
Residual Std. Error 6.006 (df = 43) 5.931 (df = 45)
F Statistic 4.150** (df = 3; 43) 11.849*** (df = 1; 45)
===============================================================
Note: *p<0.1; **p<0.05; ***p<0.01
観測値と予測値のプロット
観測値と単回帰による予測値を散布図にプロットします。
> reg1 <- lm(her~univ, data=data1)
予測値のマーカーを赤にしています（col=2）。
> plot(data1$univ,predict.lm(reg1),ylim=c(0,70),col=2,xlab="大学数",ylab="進学率（観測値と予測値）")
> par(new=T)
縦軸の範囲を指定しています（ylim=c(0,70)）。観測値と予測値の範囲が異なる場合、ずれてしまうので。
> plot(data1$univ,data1$her,ylim=c(0,70),pch=4,xlab="",ylab="")
recode
最終学歴（1：中学、2：高校、3：短大・高専、4：大学、5：大学院）を修学年数（9：中学、12：高校、14：短大・高専、16：大学、18：大学院）にリコードします。
recodeコマンドを利用するために、パッケージcarをDLしてインストールしておく必要があります。
> library(car)
> data1$yedu=recode(edu,'1=9;2=12;3=14;4=16;5=18')
散布図の任意の場所に任意の凡例を表示する
学歴（eduy：就学年数）別に年齢（age）と年収（ainc）の散布図を作成します。
> plot(data1$age,data1$ainc,type="n")
> text(data1$age,data1$ainc,data1$eduy)
続いて以下を実行して、グラフ上の任意の場所に左クリックで凡例を置きます。
> legend(locator(1),c("9:中学","12:高校","14:短大・高専","16:大学","18:大学院"))
一部の変数で別のデータフレームを作成
data1にいろいろな変数がある場合、以下のようにして抜き出すことができます。
> attach(data1)
> data2 <- data.frame(inc05,stu05,unv05)
残差プロット
時系列データの分析で、残差をプロットします。ols1に推定結果が格納されているとします。
> plot(year,residuals.lm(ols1))
また、以下のように予測値を計算して実現値から引いてもできます。
> data1$ptfr <- predict(ols1)
> data1$res <- tfr-ptfr
> attach(data1)
> plot(year,res)
> abline(h=0)
定数項なしの回帰
階差データの回帰などでは、定数項なしで回帰したいときがあります。その場合、式の部分に「0」と記述すればOKです。
> ols1 <- lm(tfr~0+uem)
また、「1」と記述すれば、通常の定数項ありの回帰になります。
> ols2 <- lm(tfr~1+uem)
t分布のグラフ
以下のコマンドで作成しました。もう少しシンプルにできるのかもしれません。「+」の部分は無視してつなげてOKです。
> curve(dt(x,45),-5,5,xlab="t")
> abline(h=0)
> segments(qt(0.025,45),0,qt(0.025,45),dt(qt(0.025,45),45));segments(qt(0.975,45),0,qt(0.975,45),
+ dt(qt(0.975,45),45))
> segments(qt(0.005,45),0,qt(0.005,45),dt(qt(0.005,45),45));segments(qt(0.995,45),0,qt(0.995,45),
+ dt(qt(0.995,45),45))
> text(0, 0.05, labels="95%")
> arrows(-0.6, 0.05,qt(0.025,45),0.05);arrows(0.6,0.05,qt(0.975,45),0.05)
> text(0,0.01,labels="99%")
> arrows(-0.6,0.01,qt(0.005,45),0.01);arrows(0.6,0.01,qt(0.995,45),0.01)
データの接合
data1とdata2を変数prefで一致するように接合し、新たなdata3としています。
> data3 <- merge(data1,data2,by="pref",all=TRUE)
相関係数の有意性検定
2つの変数x1とx2の相関係数の有意性検定を行います。
> cor.test(x1,x2)
偏相関係数
偏相関係数のコマンドが見当たらないので関数作成の練習として示しておきます。
3変数x、y、zがあったとき、zの影響をコントロールした後のxとyの偏相関係数は以下の関数で求められます。
> pcor <- function(x,y,z){
> rxy <- cor(x,y)
> rxz <- cor(x,z)
> ryz <- cor(y,z)
> rxyz <- (rxy-rxz*ryz)/((sqrt(1-(rxz)^2))*(sqrt(1-(ryz)^2)))
> rxyz
> }
と打ち込んで実行したらx、y、zに変数を指定します。都道府県データinc、enr、univがあり、univの影響をコントロールした場合のincとenrの偏相関係数について計算すると
> attach(data1)
> pcor(inc,enr,univ)
[1] 0.6897585
となります。この偏相関係数について有意性の検定をするときは、まずt値を計算します。
> pcor(inc,enr,univ)*sqrt(47-1-2)/sqrt(1-pcor(inc,enr,univ)^2)
[1] 6.319187
自由度44の両側10％点は
> qt(0.95,44)
[1] 1.68023
となるので、10％水準の境界値を超えていることがわかります。p値を求めると
> 2*pt(6.319187,44,lower.tail=FALSE)
[1] 1.144776e-07
となり、十分に低い値になっていることがわかります。
ダミー変数1
東京ダミーを作ります。
都道府県が文字列型の場合
> data6$tokyo <- data6$pref=="東 京 都"
都道府県が数値型の場合
> data6$tokyo <- data6$pref==13
相関行列
データフレームdata1にあるすべての変数間の相関係数を行列形式で出力する。
> cor(data1)
attach(データフレーム名)
> attach(data1)
としておくと、以下のようにデータフレームの指定を省略できます。
> ols1 <- lm(sch06~inc05+stu05+unv05)
ただし、新たな変数をデータフレームに追加した時は、再度attachするのを忘れないようにしましょう。
Page updated
Google Sites
Report abuse