Rの基本的操作
西山
2016-10-3
「将来予測の技術」の第1回は統計分析ツールRの基本的操作にあてる。昨年度の基礎科目「統計分析の基本」はRコマンダーで利用したのでコマンドによる操作には（特に最初は）戸惑うかもしれない。Rコマンダーでも予測計算ができるツールがあることはあるのだが、最も有用な機能が備わっていない — 具体的には、タから最善の予測モデルを自動探索する機能がない。そのため、使いこなすための基礎的な勉強量がかえって増えてしまう。
「変化の規則性」をとらえるのにどんな統計的思考をしているのか？もちろん基礎的な勉強は重要だが、このために全ての時間を費やするのは余りにも非効率だろう。予測の計算では専用のパッケージ“forecast”を活用して、ルーティン的に作業を進めるので、Rの基本的操作で求められることは（実は）あまり多くはない。登場する予測コマンドも（最終的には）これ以上期待できないほど単純である。最終的には、すべての履修者がかなり高度の将来予測が出来るようになるといま断言しておきたい。
とはいえ、予測計算を行う前にデータに日付を付けたり、変化の規則性を図にするなど、時系列データにつきまとう前処理がある。こちらのほうがむしろ面倒臭いもので、しかもRコマンダーでは出来ない。所詮はルーティン作業であるので、コマンドを記憶するというより、作業の流れを頭に入れるつもりでフォローしてほしい。使われるコマンドは、まず覚えるのではなく、「進行シナリオ」を手元に置いて、作業の流れを覚えつつ、コマンドは記憶違いもあるので確かめながら、自然に覚えていく（覚えなくともよい）。このくらいのスタンスが「強迫感」もなく、ちょうど良いと思う（実際、西山もそうしてきました）。
特に必要な事柄は、データファイルの入力、Rによる分析結果をRDataファイルに保存すること、さらにRで値を保存する「変数」の形の一つであるデータフレームの特徴を知っておくことである。
第1回では授業ノートや参考書を使わず、資料「統計分析ソフトRの使用法」を用いる。この資料の1頁から13頁までは、Rの基本の中でも特に基本的なことを解説しているので、まだよく知らない人は必ず精読しておいて―できれば自分でも文章に沿って作業をしておいて―ほしい。授業では「データフレームの活用」（14頁）から始めたい。
ただ資料には一点だけ大変重要な勘所（というかコツ）に触れている。それは（作業用）ディレクトリーの変更である。Windows PCであれば、メニューの＜ファイル→ディレクトリーの変更…＞を選ぶ。Macであれば＜その他→作業ディレクトリの変更＞を選ぶ。
データファイルは、作業用ディレクトリーを指定せずとも読み込めるのであるが、各モジュールで使うファイルや結果の保存は（やはり）フォルダーを別々にしたほうが分かりやすいし、便利だろう。作業用ディレクトリーを一番初めに指定しておくことは良い習慣である。
データファイルの読み込みとデータフレーム
資料でも使っているデータファイル“gdpcons.txt”を読み込もう。データファイルは、普通はR以外の（例えば）エクセルなどで作成、保存することが多い。最近は、エクセルで保存した“xls”や“xlsx”ファイルをRからそのまま読めるようになっているが、これには専用のパッケージが必要であったり、また漢字が混在したりすると読めないことがある。データファイルは、エディターでも読めるテキストファイルで保存したほうがよい。またCSVファイルにしておけばエクセルで保存、修正ができるので最も便利である。
ファイル“gdpcons.txt”はダブルクリックすると（普通は）ノートパッド（Macではテキストエディット）が立ち上がり内容が表示される。第1行には保存されている個別データ名（Year、Gdp、Cons）があり、値は第2行以下に並んでいる—各行3個の値がブランクで区切られている。このファイルを読むには“read.table”を使う。
まずデータファイルの名前を確認するためディレクトリ内のファイル一覧を出しておく。
dir()## [1] "forecast-m1-2015-utf8.html" "forecast-m1-2015-utf8.Rmd"
## [3] "gdpcons-1994-2012.csv" "gdpcons.csv"
## [5] "gdpcons.txt"次いで、ファイルを読み込んで内容を“data0”という名前で保存する。
data0 <- read.table(file="gdpcons.txt",header=TRUE)
data0## YEAR GDP CONS
## 1 1980 312712.7 174382.7
## 2 1981 321490.5 177074.9
## 3 1982 331710.7 184799.3
## 4 1983 339823.8 189292.0
## 5 1984 353436.2 194237.4
## 6 1985 368184.1 201627.8
## 7 1986 379895.7 209050.0
## 8 1987 399442.3 217356.6
## 9 1988 424657.3 229129.5
## 10 1989 445468.8 238784.9
## 11 1990 469780.5 248840.1
## 12 1991 481660.7 256905.6
## 13 1992 483375.6 261560.2
## 14 1993 485498.4 266385.2
## 15 1994 490730.7 272342.2
## 16 1995 502794.3 277906.5
## 17 1996 520053.8 284766.8
## 18 1997 521315.1 281393.7
## 19 1998 518380.7 285094.0
## 20 1999 525695.8 289454.2
## 21 2000 530312.8 288981.1コマンドをキーボードから入力してリターンキーを打つと実行される。保存した変数名をうつと、その名前で保存されている内容が表示される。コマンドの中の“header=TRUE”は、ファイルの第1行が個別データ名であることを指定している。もし個別データ名のない、値だけのファイルを読み込むなら、“header=FALSE”とするのだが、個別データ名のないデータファイルは普通は作らないので、こういうことはまずないだろう。
別のやり方もあるので試されたい。
data0 <- read.table(file=file.choose(),header=TRUE)
data0上では“file.choose”を使っているのでファイル選択画面が出てくる。ここで読み込みたいファイルを選択すればよい。このfile.chooseを使うと、最初に作業用ディレクトリーを指定せずとも必要なデータは読み込める。便利であるが、作業用ディレクトリーを指定しないと、データファイルがあるフォルダーと分析結果を保存するフォルダ-が違うことになるので注意されたい。作業用ディレクトリーを指定しなければマイ・ドキュメント・フォルダーがデフォールトで作業用ディレクトリーになる。
同じ内容を保存したデータファイルとして“gdpcons.csv”も提供している。これはCSVファイルである。CSVというのは“Comma Separated Values”の略称であって、普通はエクセルなど表計算ソフトで編集、保存する。ファイルの中身は、エクセルのワークシートから色や文字の大きさなどの情報を除き、値だけを保存したものになっている。エクセルでなくともエディターでも表示できる。内容をみると、ブランクでなく、コンマで値が区切られている。CSVファイルを読むには“read.csv”を使う。
次を試されたい。
data0 <- read.csv(file="gdpcons.csv",header=TRUE)データの分析とグラフの作成
読みこんだデータは年次別のGDPと民間消費支出である。名前“data0”は、エクセルにたとえればシート名が“data0”ということであり、“YEAR”、GDP“、”CONS“（＝消費）などの個別データは”data0“と入力することでとり出せるわけではない。
たとえば、2列目の“Gdp”をとり出すには次のような方法がある。どれも結果は同じである。
data0$GDP## [1] 312712.7 321490.5 331710.7 339823.8 353436.2 368184.1 379895.7
## [8] 399442.3 424657.3 445468.8 469780.5 481660.7 483375.6 485498.4
## [15] 490730.7 502794.3 520053.8 521315.1 518380.7 525695.8 530312.8data0[,2]## [1] 312712.7 321490.5 331710.7 339823.8 353436.2 368184.1 379895.7
## [8] 399442.3 424657.3 445468.8 469780.5 481660.7 483375.6 485498.4
## [15] 490730.7 502794.3 520053.8 521315.1 518380.7 525695.8 530312.8attach(data0)
GDP## [1] 312712.7 321490.5 331710.7 339823.8 353436.2 368184.1 379895.7
## [8] 399442.3 424657.3 445468.8 469780.5 481660.7 483375.6 485498.4
## [15] 490730.7 502794.3 520053.8 521315.1 518380.7 525695.8 530312.8一番上では、データ全体の名前と個別データ名がドル記号“$”で区切られている。このドル記号を使った区切り方式は、分析結果の一部を取り出すときにも使うことがあり、Rを特徴づける「習慣」である。
2番目は、縦横に表形式で並んでいるデータから2列目をとりだす、という式である。詳しくは資料を呼んでほしい。
上の最後の2行で実行しているが、“attach(変数名)”とすると、アタッチされたデータに含まれている個別データ名がそのまま使えるようになる。
ここで横軸を年、縦軸をGDPとして折れ線グラフを描こう。
plot(YEAR,GDP,type="l",main="GDPの推移")
同じように横軸を年、縦軸を消費(CONS)をとった折れ線グラフを描け。
次に、横軸をGDP、縦軸をCONSとした散布図を描く。
plot(GDP,CONS,pch=16,col="red")
上から分かるように図を描くときは“plot”コマンドを使う。そして、最初に横軸（X軸）に割り当てるデータ名、次に縦軸（Y軸）に割り当てるデータ名を指定する。
あとはオプションで指定する必要はない。もし“pch”を指定すると、プロットされる記号を変えることができる。また“col”を指定すると色を変更できる。
【留意事項】Rで利用できる機能をフルに発揮するには、たとえば上のpchの確認など、ヘルプ機能を簡単に使いたい。それにはRStudioを利用すると遥かに便利になる。R本体をインストール済みであれば、フリーで追加インストールできるので活用されたい。授業でもRStudioで説明したほうが分かりやすい時期が来ると予想している。 入手先：https://www.rstudio.com/
GDPや消費の平均値を出しておこう。
mean(GDP); mean(CONS)## [1] 438401## [1] 239493.6セミコロンで区切れば、1行に複数のコマンドを書いてもよい。
データ全体の要約をしておくなら次のようにするのが定石である。
summary(data0)## YEAR GDP CONS
## Min. :1980 Min. :312713 Min. :174383
## 1st Qu.:1985 1st Qu.:368184 1st Qu.:201628
## Median :1990 Median :469780 Median :248840
## Mean :1990 Mean :438401 Mean :239494
## 3rd Qu.:1995 3rd Qu.:502794 3rd Qu.:277906
## Max. :2000 Max. :530313 Max. :289454一度attachで名前を登録した変数はdetach(変数名)で登録をはずすまで、その中の個別変数名をドル記号($)で区切らずに使うことができる。一度、登録から外すとGDPやCONSは名前としては使えない。
次のコマンドを実行して確認せよ。
detach(data0)
GDP
mean(data0$GDP)回帰分析と予測値の計算
データ全体か回帰直線\(CONS = a + b \times GDP\)を求めよう。求めた関係式に基づいて、GDP＝530.0の時の消費（CONS)を予測しよう。
Rでは単回帰分析と重回帰分析でやり方は変わらない。
kekka1 <- lm(CONS ~ GDP,data=data0)
summary(kekka1)##
## Call:
## lm(formula = CONS ~ GDP, data = data0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7363.8 -1883.4 542.3 2120.9 4981.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.035e+03 4.299e+03 1.404 0.177
## GDP 5.325e-01 9.666e-03 55.095 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3319 on 19 degrees of freedom
## Multiple R-squared: 0.9938, Adjusted R-squared: 0.9935
## F-statistic: 3035 on 1 and 19 DF, p-value: < 2.2e-16得られた結果は\(CONS = 6034.5028232+0.532524\times GDP\)である。
結果の要約をみると、決定係数は十分に高く、また有意である。というより、民間の消費がGDP（≠景気）に影響されることは先刻ご承知であろう(^^;;)
散布図に回帰直線を重ねて描いておこう。
plot(data0$GDP,data0$CONS)
abline(kekka1,lty=2,col="red")
まず元のデータから散布図を描き、回帰分析から得られた直線を点線(lty=2)、色は赤（col=“red”）で描き足している。
予測値の計算
回帰分析から得られた結果式の右辺にあるGDPに530兆円を代入すれば、その時の消費の予測値が出るわけだが、以下のようにすれば簡単に求められる。データの値は10億円単位になっていることに注意する。
newdata <- data.frame(GDP=530000)
predict(kekka1,new=newdata)## 1
## 288272.2もし、GDPが530.0、535.0という二つの値を代入して消費の予測値を出したいなら以下のように“newdata”を指定する。コマンド“c”を使っているが、これは複数の値を一つにまとめる最も基本的なコマンドである。資料を読まれたい。
newdata <- data.frame(GDP=c(530000,535000))
predict(kekka1,new=newdata)## 1 2
## 288272.2 290934.9ある区間全体について予測したいなら以下のようにする。
predict(kekka1,new=data.frame(GDP=seq(530000,535000,1000)))## 1 2 3 4 5 6
## 288272.2 288804.8 289337.3 289869.8 290402.3 290934.9上では、予測値を出したGDPの値として530兆円から535兆円まで1兆円ずつ刻んだ全ての値（6個）を指定している。故に予測値も順に6個計算されることになる。
コマンド“seq”はこのような数字列をつくるためのコマンドである。
【留意点】基本的なコマンドは知っておくにこしたことはないが、（案外多くの人が序盤に思うようなのだが）基本的なコマンドを十分理解してから、しかる後にRによる統計分析に進もうという姿勢は、理には適っているが『いつまでたっても使えない』という結果に至りがちである。時間のあるときに提供した資料を読んでおく程度が統計の勉強には最良だと感じる。何事によらず「勉強時間は十分にはない」ものである。
予測区間を確認したいときもある。確率95％で可能性が含まれる95%予測区間を上と同じ区間について出しておく。
predict(kekka1,new=data.frame(GDP=seq(530000,535000,1000)),interval="confidence")## fit lwr upr
## 1 288272.2 285878.2 290666.2
## 2 288804.8 286395.1 291214.5
## 3 289337.3 286911.8 291762.7
## 4 289869.8 287428.5 292311.1
## 5 290402.3 287945.2 292859.5
## 6 290934.9 288461.7 293408.0上のように“interval”を指定すれば予測区間を出してくれるわけである。もしinterval=“prediction”とすれば、その他要因による残差だけではなく、回帰直線自体に含まれる誤差も考慮した予測誤差を評価してくれる。
以下を実行せよ。
predict(kekka1,new=data.frame(GDP=seq(530000,535000,1000)),interval="prediction")新しい説明変数の作成と新たな結果、新たな予測値
1987年から1990年まではバブル景気だった。
元のデータには、どの年がバブル期にあったかという情報が含まれていないが、バブルで資産価格が高騰していたことにより民間の消費活動は他の期間と異なっていた可能性がある。そこで元のデータにその年がバブル期にあたっていたかどうかの情報を追加して回帰分析をやり直すことにしよう。
まず、その年がバブル期であったかどうかは“YEAR >= 1987”という不等式と“YEAR <= 1990”いう不等式が同時に満たされるかどうかで区別される。上の条件がTRUEならバブル期、FALSEであれば平常期であると設定しよう。
attach(data0)## The following objects are masked from data0 (pos = 3):
##
## CONS, GDP, YEARbubble <- (YEAR >= 1987) & (YEAR <= 1990)
bubble## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE
## [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEbubble <- factor(bubble,level=c(FALSE,TRUE),label=c("NO","YES"))
bubble## [1] NO NO NO NO NO NO NO YES YES YES YES NO NO NO NO NO NO
## [18] NO NO NO NO
## Levels: NO YESdata0$bubble <- bubble
data0## YEAR GDP CONS bubble
## 1 1980 312712.7 174382.7 NO
## 2 1981 321490.5 177074.9 NO
## 3 1982 331710.7 184799.3 NO
## 4 1983 339823.8 189292.0 NO
## 5 1984 353436.2 194237.4 NO
## 6 1985 368184.1 201627.8 NO
## 7 1986 379895.7 209050.0 NO
## 8 1987 399442.3 217356.6 YES
## 9 1988 424657.3 229129.5 YES
## 10 1989 445468.8 238784.9 YES
## 11 1990 469780.5 248840.1 YES
## 12 1991 481660.7 256905.6 NO
## 13 1992 483375.6 261560.2 NO
## 14 1993 485498.4 266385.2 NO
## 15 1994 490730.7 272342.2 NO
## 16 1995 502794.3 277906.5 NO
## 17 1996 520053.8 284766.8 NO
## 18 1997 521315.1 281393.7 NO
## 19 1998 518380.7 285094.0 NO
## 20 1999 525695.8 289454.2 NO
## 21 2000 530312.8 288981.1 NO上で作った変数（＝新たに作ったデータでもある）“bubble”の値は、本来はTRUE（＝真）、FALSE（＝偽）という論理値である。
値がこのまま表示されても意味はわかるはずだが、更に見やすくするためにTRUE、FALSEという値（levelと呼ぶ）に表示名（label）を対応させることができる。それが上のコマンドの中の“factor”である。
“YES”や“NO”という値は文字列だが、毎年を平常期かバブル期かで分類するための項目名になっている。こういうデータをRでは「因子（ファクター）」と呼んでいるのだが、こんなデータを説明変数に含めるときはなるべく上のようにして、分かりやすい項目名にしておくとよい。そのためのコマンドが“factor”である。
変数“bubble”は1987年から1990年までがYES、それ以外の年はすべてNOになっている。このように作った変数“bubble”をデータセット“data0”に追加するには、新しい名前をドル記号($)で区切って値を代入すればよい。
その後、回帰分析の説明変数に“bubble”を加えて回帰分析を行う。
kekka2 <- lm(CONS ~ GDP + bubble, data=data0)
summary(kekka2)##
## Call:
## lm(formula = CONS ~ GDP + bubble, data = data0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6556.2 -1189.8 787.1 1275.5 4055.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.259e+03 3.465e+03 2.095 0.0506 .
## GDP 5.319e-01 7.752e-03 68.620 <2e-16 ***
## bubbleYES -5.028e+03 1.479e+03 -3.399 0.0032 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2661 on 18 degrees of freedom
## Multiple R-squared: 0.9962, Adjusted R-squared: 0.9958
## F-statistic: 2367 on 2 and 18 DF, p-value: < 2.2e-16バブル期には財布の紐が緩んで消費にカネを使いがちである。そんな予想に反して、バブル期に該当した年には民間消費が下押しされるという結果になっている — 説明変数“bubble”がYESになると、民間の消費支出が抑えられるという係数推定値になっている。しかも誤差が小さく、有意である。
【考察】 予想しない結果が得られるときは、なぜそうなったかという理由を考える必要がある。PCや分析ソフトが計算で間違うことは（経験上）まずないと考えたほうがよい。今回の場合は、説明変数のGDPを本来は家計の可処分所得にするべきであったところ、GDPにしているのが原因である。GDPは消費によって増えるだけではなく、住宅投資や企業の設備投資によっても増える。バブル期には不動産投資や企業設備投資が非常に増えた。ということは、バブル期にあっては同じGDPであっても投資でもたらされていたのであり、平常期に当てはまっていた関係から予測されるより消費は低めに出ていた。そう解釈できる。
上の考察をさらに確認するためGDPに占める民間消費の割合をみておこう。
plot(data0$YEAR,data0$CONS/data0$GDP*100,xlab="年",ylab="消費対GDP")
abline(h=mean(data0$CONS/data0$GDP*100),lty=2)
図中に点線で引いた水平線は期間を通した消費割合の平均値（＝54.6819948を表す。図から分かるように、バブル景気が深まるにつれて、GDPに占める消費支出の割合は低下している。
『GDPの高さ（＝景気）の割には民間消費はそれほどでもなかったんだね・・・』というのがバブル期だった。
バブル期であるかどうかを説明変数に加えた上の結果に基づいて、前に行った予測を修正しよう。
次の作業を行え。 各行のコマンドでやっていることを自分でやりながら、よく考えること
newGdp <- seq(530000,535000,1000)
newBubble <- rep("NO",length(newGdp))
newData <- data.frame(GDP=newGdp,bubble=newBubble)
newData
predict(kekka2,new=newData,interval="confidence")前に行った予測計算の結果とどこが違っているか？
重回帰分析結果を利用した予測計算はよく分かったか？使えそうか？
結果の保存と分析の再開
これまでの分析によって得られた結果は、どんな名前で保存しているか。これを確かめるにはコマンド“ls”を用いる — lsはリスト(list)の意味である。コマンドの後ろは必ず丸括弧がつくことに注意せよ。
ls()## [1] "bubble" "data0" "kekka1" "kekka2" "newdata"これらの結果を保存しておくにはメニューから＜ファイル→作業スペースの保存…＞を選んでからファイル名を入力する。ファイル名には必ずファイル修飾子“RData”をつけておく — 例えば“TodayKekka.RData”という名前で保存すると、このファイルも作業用ディレクトリーに保存される。そしてRという文字が浮き出たアイコンになるので、Rが作成したファイルであることが一目で分かる。
RDataファイルをダブルクリックすると、直ちにRが起動され、作業用ディレクトリーはそのRDataファイルがあったフォルダーになる。
エクササイズ
（１）次のコマンドを実行せよ。どんな結果が得られるか？
dim(data0)
data1 <- data0[-1,]
data1
dim(data1)（２）データファイル“gdpcons-1994-2012.csv”を読み込み、上と同じ回帰分析を行え。但し、このデータファイルの期間は1994年から2012年までであるから、バブル景気の影響を調べるにはあたらない。