トップ>R

Rについて

 Rは高機能な統計処理ソフトです(しかも無償!)。

 Rの詳細についてRjpWikiを、インストール方法についてはこちらをどうぞ。

 

 Rを使いやすくするためのインターフェイスには以下のものがあります。

 Rコマンダー

 RStudio

Rでの分析の前に

 Rで分析するにあたって、作業場(データファイル等の置き場所)を指定する必要があります。以下は32-bit版の画面で、「File」→「Change dir...」と進み、読み込みたいデータのあるフォルダを指定します。64-bit版では日本語対応していますので「ファイル」→「ディレクトリの変更...」と進んでフォルダを指定します。Mac版の場合は、「その他」→「作業ディレクトリの変更」でフォルダを指定します。

direc

 以下のようなコマンドを実行して作業場を指定することもできます。これは、デスクトップに読み込みたいデータのある「Rdata」というフォルダがある場合です。ユーザー名は各自のものになります。

> setwd("c:/Users/ユーザー名/Desktop/Rdata")

 Mac版の場合は以下のようになります。

> setwd("/Users/ユーザー名/Desktop/Rdata")

 

 分析にあたっては、コマンドを保存しながら進めるのが効率的です。32-bit版では「File」→「New script」でR-Editorが開くのでそこにコマンドを記述し、実行します。64-bit版では「ファイル」→「スクリプトを開く...」です。

scrip

 R-Editor上でコマンドを実行する場合は、実行したいコマンドのどこかにコマンドプロンプト(点滅している縦線)がある状態で、以下のボタンを押します。またショートカット「Ctrl+R」で実行できます。

run

Rのコマンド

注意

以下のコマンドでは

> attach(data1)

のようにデータフレームの指定を前提としているケースがあります。そのため、そのままコピペして実行すると失敗することがあります。attachを使用しないときは、変数名にデータフレーム名をつけて実行してください(var1→data1$var1)。あるいは使用するデータフレームを指定するのがよいでしょう(data=data1)。

ここに掲載のコマンドは授業の補足説明のためのメモ書きで、基本的に洗練されてないコマンドですがご容赦ください。

時間データを分に変換する

時間データで「時間.分」という形になっていることがあります。例えば2.54は2時間54分のことであり、2.54時間(約2時間30分)のことを意味しているのではありません。また、2.89というデータも存在しないので、そのまま使用するとおかしな散布図になります。そこで、分に換算して分析に用いることになります。以下のように関数を利用して換算します。

> floor(2.54)*60+(2.54%%1)*100
[1] 174

floor(x)はxを超えない整数を返す関数で、y%%z%%は、yをzで割り、その余りを得る演算子です。

データフレームの指定部分で記述統計量を計算2

データフレームの3~5列目と7列目の記述統計を計算したいときは以下のようにコマンドします。

> summary(data1[,c(3:5,7)])
> apply(data1[,c(3:5,7)],2,sd)

データフレームの3~5列目と7列目だけコンソールに表示させたい時は以下のとおりです。
> data1[,c(3:5,7)]

複数回答のデータを作り直す

複数回答(あてはまるものすべてに○系の設問)のデータで、選んだ選択肢の番号が1変数に1つずつ入っていることがあります。これを各選択肢についての1/0のダミー変数にしたい場合があります。例えば選択肢1を選んだダミー変数は以下のようにします。

> data1$c1[data1$q1_1 == 1 | data1$q1_2 ==1 | data1$q1_3 ==1 | data1$q1_4==1 ] <- 1

該当しない人はNAになっているので0にリコードしておきましょう。

> data1$c1 <- recode(data1$c1, '1=1;NA=0')

複数の変数から1つの新しい変数をつくる

例えば2つの変数の組み合わせで新しい変数をつくる場合、以下のようにコマンドします。

> data1$newvar[data1$var1==1 & data$var2==1] <- 1

> data1$newvar[data1$var1==1 & data$var2==2] <- 2

datファイルの読み込み

データファイルがxxx.datのように拡張子がdatになっている場合、以下のようにコマンドして読み込みます。

> data1 <- read.delim("xxx.dat")

なお、datファイルではデータがタブで区切られています。

データフレームの指定部分で記述統計量を計算

summary関数では標準偏差(sd)を計算してくれないので別途、計算する必要があります。ここではapply関数を利用して、複数の変数の標準偏差を一度に計算してみます。

> apply(data1[,2:5], 2, sd)

   flp80     tfr80    flp00    tfr00
6.6286828 0.1345359 4.0113553 0.1331380

apply関数の第1引数でデータフレームdata1の2列目から5列目を指定しています。第2引数で列単位で計算するよう指定しています。第3引数で計算する統計量をsd(標準偏差)に指定しています。

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)

Logitモデルの限界効果とオッズ比

Logitモデルの推定係数をそのまま解釈することは難しいです。そこで限界効果を計算する必要があります。ただし、非線形の回帰分析のため、限界効果を見たい独立変数およびその他の独立変数の水準(値)によって限界効果が異なります。

あらかじめmfxパッケージをインストールしておき、呼び出します。

> library(mfx)

平均的な限界効果(独立変数に実際の値を使い、各対象で限界効果を計算し推定サンプルで平均化したもの)を求める場合は

> logitmfx(Y~X1+X2+X3, data=data1, atmean=F)

平均値での限界効果(独立変数に平均値を代入して計算したもの)を求める場合は

> logitmfx(Y~X1+X2+X3, data=data1)

また、社会学等ではオッズ比が使われることが多いので参考までに。

> logitor(Y~X1+X2+X3, data=data1)

Rの出力をExcelに貼り付け

Rの出力をそのままExcelに貼り付けると、下の図のように変数名や係数、標準誤差などが1セルにまとめて入ってしまいます。

これらの数値を1セルごとに入れたい場合、「データ」タブの「区切り位置」を押します。

r_excel1

すると、以下のような画面になりますので、「スペースに・・・」が選択されていることを確認して「次へ」を押します。

r_excel2

どこで区切るか表示されますので、よければ「次へ」か「完了」を押してください。

次の画面では表示形式を操作できますが、ここでは必要なさそうです。

r_excel3

以下のようにきれいにセルに入りました。

r_excel4

加工した数値変数が正しく生成されたか並べて確認

例えば数値変数xを100倍してx2を作成し並べるには以下を実行します。

Dataset$x2 <- with(Dataset, x*100)
Dataset[,c("x","x2")]

因子の場合はクロス表でも確認できます。

サンプルサイズの確認

順序回帰モデルの推定などでは、使用したデータ数が表示されないことがあります。結果としては保存されているので表示することができます。例えば推定結果がOrdRegModel.1に格納されていたとします。この場合、以下のように入力し実行します。

> OrdRegModel.1$nobs

あるいは

> OrdRegModel.1$n

です。

条件に合うデータだけ抜き出してデータフレームをつくる

まずはテストデータを作成します。

> 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))
以下にエラー subset.default(z == 3, c(x, y, z)) :
'subset' は論理値でなければなりません
> data2
エラー: オブジェクト 'data2' がありません

データの四捨五入

小数部分のある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の変数名は任意です。

散布図のラベル表示

ここでは都道府県データを使って通常の散布図を描きました。ただし、どのマークがどの都道府県なのかわかりません。

> plot(ch3$inc05,ch3$sch05)

散布図1

そこで、データの並び番号で表示してみます。

1行目のtype="n"はマークを表示しない命令です。

2行目で並び番号を表示させています。

> plot(ch3$inc05,ch3$sch05,type="n")
> text(ch3$inc05,ch3$sch05)

散布図2

都道府県名を表示したい場合は以下のようにします。

1行目はデータフレームch3の1列目にある都道府県名を抜き出し、labelという名前で保存しています。

3行目でlabelを使用しています。

> label <- ch3[,1]
> plot(ch3$inc05,ch3$sch05,type="n")
> text(ch3$inc05,ch3$sch05,label)

散布図3

都道府県名は横に広いため、実際のマークの位置との対応を確認します。

上の散布図が表示された状態でさらに以下を実行します。

> text(ch3$inc05,ch3$sch05)
散布図4

見にくいですが、都道府県名の真ん中にマークがあることがわかります。愛知県であれば「知」がマークということです。

データ上、都道府県名は「愛 知 県」というようになっていますが「愛知」とすれば、もう少しすっきりと表せるでしょう。

対数

10の自然対数をとるには

> log(10)
[1] 2.302585

底が2の場合
> log2(10)
[1] 3.321928

底が10の場合
> log10(10)
[1] 1

ある変数varの自然対数をとってlnvarとして保存する場合

> lnvar <- log(var)

散布図

散布図でポイントを線でつなぐにはtypeで指定します。

まず適当なデータを生成します。

> x <- 0:10
> y <- x^2-5*x+5
> x
[1] 0 1 2 3 4 5 6 7 8 9 10
> y
[1] 5 1 -1 -1 1 5 11 19 29 41 55

次にplot命令でtypeを指定します。

> plot(x, y, type="b") #線とポイント、重なりなし

散布図b

そのほかにもいくつかバリエーションがあります。

> plot(x,y,type="l") #線のみ
> plot(x,y,type="o") #線とポイント、重なりあり

> plot(x,y,type="c") #線のみ、ポイント部分空白

データのソート

以下のようなデータフレーム(data2)があるとします。

> data2

       pref crm00 crm05 crm08

1  北 海 道 15.27 12.98 10.79

2  青 森 県 11.34 10.28  7.91

3  岩 手 県  9.79  8.15  6.74

4  宮 城 県 19.96 14.13 12.21

5  秋 田 県 10.22  7.51  5.54

:   :   :   :   :

43 熊 本 県 14.45 12.05  9.57

44 大 分 県 12.45 10.78  8.20

45 宮 崎 県 13.59 10.43  9.77

46 鹿児島県 10.91  8.52  7.95

47 沖 縄 県 16.30 14.15 11.16

 

これを変数crm00の昇順に並べ替えます。

まずorder関数でcrm00の順を作成し、その順を表示してみます。

> sort<-order(data2$crm00)

> sort
すると以下のように当初のサンプル番号で順が示されます。

[1] 42 6 3 32 5 15 16 17 31 46 2 44 41 36 7 37 18 45 24 19 10 43 20 35
[25] 1 39 22 47 33 38 9 28 29 8 34 21 25 4 14 26 30 11 13 23 12 27 40

この結果を使って並べ替えたデータフレームを新たに作成し、表示します。

> data2a<-data2[sort,]

> data2a

       pref crm00 crm05 crm08

42 長 崎 県  8.90  8.68  7.17

6  山 形 県  9.78  8.51  6.67

3  岩 手 県  9.79  8.15  6.74

32 島 根 県  9.89 10.22  8.00

5  秋 田 県 10.22  7.51  5.54

:   :   :   :   :

13 東 京 都 24.15 20.19 16.52

23 愛 知 県 25.01 27.42 19.54

12 千 葉 県 25.74 21.64 16.47

27 大 阪 府 28.66 28.30 22.92

40 福 岡 県 30.63 21.15 17.88

ヒストグラムの同時表示

デフォルトではヒストグラムは1つだけ出力されます。同時に表示したい時は以下のようにコマンドします。

> par(mfrow=c(1,2))
> hist(data2$crm00)
> hist(data2$crm05)

 

ヒストグラム横

 

縦に並べたいときは以下のようにコマンドします。ここでは、横軸をそろえてみました。

> par(mfrow=c(2,1))
> hist(data2$crm00,xlim=range(0,40))
> hist(data2$crm05,xlim=range(0,40))

 

ヒストグラム縦

 

基本統計量

データフレームを指定すると平均などが出ます。

標準偏差はでません。

> summary(data1)

csvファイルの書き出し

変数の加工などでデータフレームに修正を加えた場合、保存する必要があります。一つの方法としてcsvファイルとして作業フォルダに書き出す命令があります。

以下は"data1"というデータフレームを"data2"というcsvファイルとして作業フォルダに書き出すコマンドです。

> write.csv(data1,file="data2.csv")

ただ、このままだと行番号を一つの変数として認識してしまうので

> write.csv(data1,file="data2.csv",row.names=FALSE)

とします。

ダミー変数3

東京or三重ダミーの作り方。

都道府県が文字型変数の場合です。

> data6$tkymie <- data6$pref=="東 京 都" | data6$pref=="三 重 県"

カイ2乗分布(自由度4)

> curve(dchisq(x,df=4),xlim=c(0,15))
> abline(h=0)

chi

サブセットでの回帰

講義用データ"panel"を使用します。

> data <- read.csv("panel.csv")
> attach(data)

2000年のデータを1、それ以外を0とするyearダミー(d00)を作成します。

> data$d00 <- as.integer(year==2000)

> attach(data)

全サンプル、1980年サンプル、2000年サンプルの順に推定してみます(結果表示は省略)。

全サンプル

> reg <- lm(sch~rinc)

1980年(d00=0)のサンプル

> reg3 <- lm(sch~rinc, subset=d00==0)

2000年(d00=1)のサンプル

> reg4 <- lm(sch~rinc, subset=d00==1)

そのあと、以下のようにコマンドし、散布図の好きなところで左クリックします。

> plot(rinc,sch,type="n")
> text(rinc,sch,d00)
> abline(reg,lty=1)
> abline(reg3,lty=2)
> abline(reg4,lty=3)
> legend(locator(1),c("all","1980","2000"),lty=c(1,2,3))

サブセット

プールド回帰(全サンプル)と単年回帰(1980年or2000年)の傾きの違いが視認できます。

ダミー変数2

講義用データ"panel"を使用します。

> data <- read.csv("panel.csv")
> attach(data)

2000年のデータを1、それ以外を0とするダミー変数の作成方法は以下のようにしますが、こうするとTRUE、FALSEという表示になります。推定では1,0として扱われますが。

> data$d00 <- year==2000

データ上でも1,0表示にしたければ以下のようにコマンドします。

> data$d00 <- as.integer(year==2000)

時系列データ

データフレーム(data1)にある変数(tfr)を1975年スタート、1年ごとの時系列データとして認識させます。

> tfr <- ts(data1$tfr,start=c(1975),frequency=1)

データを表示すると

> tfr
Time Series:
Start = 1975
End = 2007
Frequency = 1
[1] 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
[15] 1.57 1.54 1.53 1.50 1.46 1.50 1.42 1.43 1.39 1.38 1.34 1.36 1.33 1.32
[29] 1.29 1.29 1.26 1.32 1.34

のようになります。また、以下のコマンドで時系列プロットができるようになります。

> plot(tfr)

TFR

さらに1階の階差をとるコマンドはdiff( )です。

> tfr_d1 <- diff(tfr)

以下のようにデータが作成されます。

> tfr_d1
Time Series:
Start = 1976
End = 2007
Frequency = 1
[1] -0.06 -0.05 -0.01 -0.02 -0.02 -0.01 0.03 0.03 0.01 -0.05 -0.04 -0.03
[13] -0.03 -0.09 -0.03 -0.01 -0.03 -0.04 0.04 -0.08 0.01 -0.04 -0.01 -0.04
[25] 0.02 -0.03 -0.01 -0.03 0.00 -0.03 0.06 0.02

交差項を使った回帰

交差項をデータとしてあらかじめ作成する必要はなく、式に記述すればOK。

reg1 <- lm(Y~x1+x2+x1*x2)

実現値と予測値のプロット

実現値yと回帰(ols6)による予測値y_hatを散布図にプロットする。

> plot(Inc05r,predict.lm(ols6),ylim=c(0,65),xlim=c(0,50),col=2,xlab="x",ylab="y and y_hat")
> par(new=T)
> plot(Inc05r,sch06,ylim=c(0,65),xlim=c(0,50),pch=4,xlab="",ylab="")

実現値と予測値

系列相関のテスト

パッケージlmtestをDLしてインストールしておく。ここでは推定結果がols1に格納されているものとします。

> dwtest(ols1)

不均一分散のテスト

同じくlmtestがインストールされていれば実行可能。

> bptest(ols1)

recode

最終学歴(1:中学、2:高校、3:短大・高専、4:大学・大学院)を修学年数(9:中学、12:高校、14:短大・高専、16:大学・大学院)にリコードする。
recodeコマンドを利用するためにパッケージcarをDLしてインストールしておく。

> library(car)
> data1$yedu=recode(edu,'1=9;2=12;3=14;4=16')

散布図のラベル表示

散布図

学歴(edu)別に年齢(age)と年収(ainc)の対応をプロットする。

> plot(age,ainc,type="n")
> text(age,ainc,edu)

そのあと

> legend(locator(1),c("1:中学","2:高校","3:短大・高専","4:大学・大学院"))

を実行してグラフの任意の場所に右クリックで凡例を表示する。

一部の変数で別のデータフレームを作成

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)

read.csv("外部ファイル名.csv")

> read.table("外部ファイル名.csv",header=TRUE,sep=",")

で外部データを読み込めますが、

> read.csv()

ではheader以降の部分を省略することができます。

edit(データ名)

データエディタを開きます。データを直接、修正することができます。

分散

不偏分散を求めるコマンドはvar(var)です。標本分散を求めるには以下のどちらかでどうぞ。

> var(var)*(length(var)-1)/length(var)

> sum((var-mean(var))^2)/(length(var))

t分布のグラフ

t分布

以下のコマンドで作成しました。もっとシンプルなものもあると思うのですがとりあえず。

> 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つの変数var1var2の相関係数の有意性検定を行います。

> cor.test(var1,var2)

偏相関係数

偏相関係数のコマンドが見当たらないので関数作成の練習も込めて以下。

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に変数を指定して

> pcor(x,y,z)

を実行すると出ます。

都道府県データx、y、zがあり、xとyの偏相関係数について有意性の検定をするときは

> pcor(x,y,z)*sqrt(47-1-2)/sqrt(1-pcor(x,y,z)^2)

でt値を計算します。t値が1.28だったとしましょう。自由度44の両側10%点は

> qt(0.95,44)

より1.68023となるので、10%水準をクリアーできていないことがわかります。p値を求めると

> 2*pt(1.28,44,lower.tail=FALSE)

より0.207となります。

ダミー変数1

東京ダミーの作り方。

都道府県が文字型変数の場合

> data6$tokyo <- data6$pref=="東 京 都"

都道府県が数値の場合

> data6$tokyo <- data6$pref==13

相関行列

データフレームdata1にある変数の相関行列を出力する。

> cor(data1)

attach( )

> attach(data1)

としておくと、以下のようにデータフレームの指定を省略できる。

> ols1 <- lm(sch06~inc05+stu05+unv05)

ただし、新たな変数をデータフレームに追加した時は、再度attachするのを忘れないように。