untitled

1

Rによる分析事例

本稿では，『計量政治分析入門』（増山幹高・山田真裕，東京大学出版会，2004年）の付録 として，同書における分析をRで実施する手順を説明します．RはR-projectが開発している統 計ソフトであり，http://www.R-project.org/を通じて自由にダウンロードできます（インターネッ ト上にはRの操作方法を解説しているサイトが多数存在します）．ここでは2005年12月20日に リリースされた2.2.1版を利用します1_{．プログラムをセットアップする際に日本語を利用できるよ} うインストールしておくと，メニューや操作画面が図R-1のように日本語表記になります2_． 図 R-1 R の操作画面

1 _{R Development Core Team (2005). R: A language and environment for statistical computing. R}

Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org. 2 _{フォントのデフォルトは Courier New ですので，最初は操作画面が文字化けします．日本語を表} 記するには，メニューの「編集」から「GUI プリファレンス」を選択し，Font を変更します（図 R-1 は MS Gothic に設定しています）．「反映」をクリックすると，文字化けが解消されます．また「保存」を クリックし，Rconsole を R-2.2.2 フォルダ内の etc フォルダにあるファイルと置換すれば，フォントの 設定を保存することができます．

ト・エディタなどで事前に準備しておけば，それらを操作画面にコピー・ペーストすることによっ て分析が可能となります．

【第 2 章】

ここでは，第 2 章の分析のうち，フリーダム・ハウスのデータに関して，2003 年のデータに限 って分析を再現します．R は多様な形式のデータ･ファイルに対応しますが，ここでは最もオー ソドックスなテキスト・ファイルを用いる手順を示しておきます．まず MS-Excel ファイルにおい て，文字列に含まれるスペースを _ などに置き換え，変数名を国名（Country），2003 年の政 治的権利（pr03），市民的自由（cl03）とします．また欠損値も NA としたうえで，テキスト形式の データ・ファイルとして保存します（たとえば，ch2.txt）． R を起動します（図 R-1）．単純化のため，データ・ファイルのディレクトリは C:/としましょう． 操作画面の>に続いて， ch2<-read.table ("C:/ch2.txt", header=T) attach(ch2) とすることでデータが読み込まれます．ch2 はデータ行列に便宜的につけた名前です．また header=T はデータの 1 行目に変数名があることを示しています．R でデータを表示するには， たとえば，list(ch2)とすれば，データ全体が出力されます．attach(ch2)としない場合，たとえば， 個別の変数は list(ch2$pr03)のように ch2 の pr03 であることを指定する必要がありますが，こ こでは attach(ch2)としていますので，list(pr03)とすることで pr03 のみが出力されます． R においては，mean(x)とすれば，x に関する平均値が求められます（メディアン：median，分 散：var，標準偏差：sd）．mean(ch2)とすると，ch2 に含まれる変数すべてについて平均値を求 めます（数値でない変数は計算できませんので NA となります）．以下は，2003 年の政治的権 利と市民的自由に関して，これらの記述統計量を示す方法の一例です．

Mean<-c(mean(pr03, na.rm=T), mean(cl03, na.rm=T)) Median<-c(median(pr03, na.rm=T), median(cl03, na.rm=T)) Variance<-c(var(pr03, na.rm=T), var(cl03, na.rm=T)) StDev<-c(sd(pr03, na.rm=T), sd(cl03, na.rm=T)) FD<-rbind(Mean, Median, Variance, StDev) colnames(FD) <- c("PR03","CL03") print(FD)

まず 1 行目では，mean(pr03, na.rm=T)によって 2003 年の政治的権利（pr03）の平均を求め ています．ここで na.rm=T は欠損値を除外することを意味します（事前に欠損値を除いている

3 場合，これを指定しなくても構いません）．同様に，mean(cl03, na.rm=T)は市民的自由（cl03） の平均を求めています．c(x, y)は x と y を列結合するコマンドですので，Mean は二つの平均 の 1×2 行列となります．2∼4 行目で Median，Variance，StDev をそれぞれ中位，分散，標準 偏差の 1×2 行列とし，5 行目の rbind によって，それらを 4×2 行列の FD に行結合させてい ます．6 行目の colnames は列名を指定するコマンドであり(行名は rownames を用います)，こ こでは各列を PR03 と CL03 としています．最後に print によって以下が出力されます3_． PR03 CL03 Mean 3.369792 3.302083 Median 3.000000 3.000000 Variance 4.642643 3.321880 StDev 2.154679 1.822603 R における度数分布は table によります．たとえば， T<-rbind(table(pr03), table(cl03)) rownames(T) <- c("PR03", "CL03") T とすれば，以下のようにそれぞれの度数分布を表示します．（またヒストグラムを作成するに は hist を用います）． 1 2 3 4 5 6 7 PR03 59 27 23 17 16 33 17 CL03 41 37 30 26 32 17 9 3 _{単に FD としても同じ出力を得ます．また summary(x)は x の最小，第 1 四分位，中位，第 3 四分} 位，最大，欠損数を求めます．母分散，母標準偏差を求めるには以下を実行します． sum((pr03[!is.na(pr03)]-mean(pr03, na.rm=T))^2)/NROW(pr03[!is.na(pr03)]) sqrt(sum((pr03[!is.na(pr03)]-mean(pr03, na.rm=T))^2)/NROW(pr03[!is.na(pr03)])) ここで[!is.na(x)]は x が欠損値である場合には適用しないことを意味します．NROW は変数を 1 列 の行列とみなして行数を求めます（NCOL：列数）．これに対して，行列一般の行数（列数）には nrow（ncol）を用います．sum(x)は x の合計，sqrt(x)は x の平方根を求めます．

MS-Excel で Jes2-1b.csv を開き，必要な変数に限って変数名をアルファベット表記に直しま す．たとえば，miyazawa を「感情温度：宮澤喜一（93 衆院事前）」とし，party を「支持政党（93 衆院事前）」としましょう．欠損値は NA とし，これをテキスト形式のデータ・ファイルとして保存 します． データを読み込んだうえで，以下の手順によって，自民支持グループと非自民支持グルー プを区別する変数を作成し，miyazawa が 100 を超えるデータを欠損値にします． ldp<-match(party, c(1), nomatch=0) miyazawa[miyazawa>100]<-NA R における等分散の検定は var.test を用います． var.test(miyazawa[ldp == 1], miyazawa[ldp == 0]) を入力すると，以下が出力されます4_．

F test to compare two variances

data: miyazawa[ldp == 1] and miyazawa[ldp == 0]

F = 1.0918, num df = 932, denom df = 1219, p-value = 0.1522

alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval:

0.968138 1.232403 sample estimates: ratio of variances 1.09177 Variance ratio test

R における平均値の t 検定は t.test を用います． t.test(miyazawa[ldp == 1], miyazawa[ldp == 0])

を入力すると，等分散を仮定しない検定が行なわれます．これにオプション指定を追加して， t.test(miyazawa[ldp == 1], miyazawa[ldp == 0], var.equal=T)

とすれば，等分散を仮定する検定となります．それぞれの出力結果は以下のとおりです． Welch Two Sample t-test

data: miyazawa[ldp == 1] and miyazawa[ldp == 0] t = 13.9879, df = 1959.782, p-value < 2.2e-16

alternative hypothesis: true difference in means is not equal to 0

4 _{等分散検定のデフォルトは両側検定ですが，以下のようにオプション指定をすれば，片側検定}

を行います．

5 95 percent confidence interval:

10.17193 13.48937 sample estimates: mean of x mean of y 42.94212 31.11148 Two Sample t-test

data: miyazawa[ldp == 1] and miyazawa[ldp == 0] t = 14.0699, df = 2151, p-value < 2.2e-16

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:

10.18169 13.47960 sample estimates: mean of x mean of y 42.94212 31.11148

ここでは投票経験と性別の関係に限って R による分析手順を紹介します．まず SPSS で jeds2000r3.sav を読み込み，付属の欠損値処理プログラムを実行します（この手順に関しては 『計量政治分析入門』第 4 章の該当箇所を参照して下さい）．まず分析に必要な変数に限って タブ区切りのテキスト・ファイルとして保存します（ファイル形式として Tab-delimited(*.dat)を選 択します）．あるいは，必要な変数を MS-Excel ファイルにコピー・ペーストし，1 行目に変数名 を追加したうえで，それをテキスト形式で保存します． データを読み込うんだうえで，table(a26a, f48)とすると，以下の出力を得ます． f48 a26a 1 2 1 682 794 2 18 29 3 37 40 7 7 2 8 5 4 このように投票経験（a26a）に欠損値である 7 と 8 が数値として含まれていますので， T<-table(vote<-a26a[a26a<=3], sex<-f48[a26a<=3]) を実行し，欠損値を除く新たな変数 vote と sex によるクロス表 T を作成します．ここで T とす ると，上記出力の a26a が 3 までのクロス表を得ます（出力省略）． またピアソンのカイ 2 乗検定は chisq.test(vote, sex) chisq.test(T) のいずれかを実行します．以下は前者を実行した場合の出力です． Pearson's Chi-squared test

data: vote and sex

X-squared = 1.2754, df = 2, p-value = 0.5285 ピアソンのカイ 2 乗検定を求めるだけでなく，クラマーの V や尤度比カイ 2 乗検定を求める には，ｖｃｄパッケージを追加し，assocstats を利用します5_{．パッケージのダウンロードは R の} パッケージ・メニューから「パッケージのインストール」によって行います（一度ダウンロードす れば，この手続きを繰り返す必要はありません）．ｖｃｄパッケージをダウンロードしたうえで， library(vcd)

5 _{David Meyer, Achim Zeileis and Kurt Hornik (2005). vcd: Visualizing Categorical Data. R}

7 summary(assocstats(T))

を実行すると，以下が出力されます． Number of cases in table: 1600 Number of factors: 2

Test for independence of all factors:

Chisq = 1.2754, df = 2, p-value = 0.5285 X^2 df P(> X^2) Likelihood Ratio 1.2895 2 0.52481 Pearson 1.2754 2 0.52851 Phi-Coefficient : 0.028 Contingency Coeff.: 0.028 Cramer's V : 0.028

第 3 章と同様，MS-Excel で Jes2-1b.csv を開き，必要な変数に限って変数名をアルファベ ット表記に直します．たとえば，miyazawa を「感情温度：宮澤喜一（93 衆院事前）」とし，jimin を 「感情温度：自民党（93 衆院事前）」としましょう．欠損値は NA とし，これをテキスト形式のデー タ・ファイルとして保存します． データを読み込うんだうえで，以下の手順によって，感情温度が 100 を超えるデータを欠損 値にします． miyazawa[miyazawa>100]<-NA jimin[jimin>100]<-NA R において散布図を作成するには plot を用います．具体的には， plot(jimin, miyazawa) とすると，図 R-2 を得ます（ファイル・メニューから画像を保存したり，クリップボードにコピーす ることができます）． 図 R-2 R による散布図 0 20 40 60 80 100 0 2 04 0 6 08 0 1 0 0 jimin m iya za w a

9 また相関分析は cor.test を用います．

cor.test(jimin, miyazawa) とすると，以下が出力されます．

Pearson's product-moment correlation data: jimin and miyazawa

t = 32.5999, df = 2102, p-value < 2.2e-16

alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval:

0.5503868 0.6071883 sample estimates: cor

まず pref（都道府県）をローマ字表記の文字列とし，変数名をそれぞれ kofuzei（交付税）， fusoku（財源不足），kokko（国庫支出金），yoto（与党得票）とします．R における線形モデル の推定は lm によって行ないます．

ols <- lm(kofuzei fusoku+kokko+yoto)

とすると，推定結果が ols に保存され，summary(ols)とすれば，以下が出力されます． Call:

lm(formula = kofuzei fusoku + kokko + yoto) Residuals:

Min 1Q Median 3Q Max -321.574 -22.414 3.529 21.693 277.412 Coefficients:

Estimate Std. Error t value Pr(>¦t¦) (Intercept) 61.763536 32.443631 1.904 0.0637 . fusoku 1.023416 0.008415 121.623 <2e-16 *** kokko 0.039934 0.019185 2.082 0.0434 * yoto 0.595574 0.798382 0.746 0.4597 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 91.83 on 43 degrees of freedom

Multiple R-Squared: 0.9988, Adjusted R-squared: 0.9987 F-statistic: 1.188e+04 on 3 and 43 DF, p-value: < 2.2e-16 また説明変数の VIF を相関行列から求めるには，

diag(solve(cor(cbind(fusoku, kokko, yoto))))

を入力すると，以下が出力されます6_．

fusoku kokko yoto 2.310137 8.251854 5.911245 6 _{第 4 章で説明したパッケージ追加の方法により，ここでは car パッケージをダウンロードしたうえ} で，以下を入力しても同様の出力を得ます． library(car) vif(ols)

（John Fox. I am grateful to Douglas Bates, David Firth, Michael Friendly, Gregor Gorjanc, Georges Monette, Henric Nilsson, Brian Ripley, Sanford Weisberg, and Achim Zeleis for various suggestions and contributions. (2005). car: Companion to Applied Regression. R package version 1.0-18. http://www.r-project.org, http://socserv.socsci.mcmaster.ca/jfox/）

11

【第 7 章】

第 6 章で用いた変数に pop（人口），area（面積）を加えます．国庫支出金の OLS 推定を ols<-lm(kokko pop+area+yoto) とすると，これに続いて plot (ols) と入力すれば，一連の残差に関する散布図が作成されます． たとえば，図 R-3 は最初に作成される予測値に対する残差の散布図です． 図 R-3 R による残差の散布図 2000 4000 6000 8000 10000 -1 0 0 0 0 1000 2000 Fitted values R e s idual s

lm(kokko ~ pop + area + yoto) Residuals vs Fitted

47

40

27

R における WLS は，人口の多い地域のデータを割り引く場合， wls<-lm(kokko pop+area+yoto, weight=1/pop)

とウェイトをオプション指定として追加します． Call:

Min 1Q Median 3Q Max -56.99 -26.79 -12.64 14.95 145.20 Coefficients:

Estimate Std. Error t value Pr(>¦t¦) (Intercept) 815.383 174.088 4.684 2.83e-05 *** pop 7.919 1.509 5.249 4.47e-06 *** area 6.643 1.228 5.409 2.62e-06 *** yoto -7.457 8.481 -0.879 0.384 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 42.97 on 43 degrees of freedom

Multiple R-Squared: 0.8139, Adjusted R-squared: 0.8009 F-statistic: 62.68 on 3 and 43 DF, p-value: 9.643e-16

また不均一分散に一致させた標準誤差（White の Robust 標準誤差）を用いるには，先に定 義した推定モデル ols を前提として，以下を実行します（ここでもサンプルの有限性を考慮する 方法として

e

ˆ

_i2

n

(

n

−

k

)

を用いる手法を紹介しておきます）． X<-model.matrix(ols) s<-summary(ols) Xc<-solve(crossprod(X)) r<-sqrt(diag(Xc%*%(crossprod(X*resid(ols)))%*%Xc*nrow(X)/(nrow(X)-ncol(X)))) t<-abs(coef(ols)/r) p<-2*(1-pt(t, nrow(X)-ncol(X))) results<-cbind(coef(ols), r, t, p) dimnames(results)<-dimnames(s$coefficients) results を実行すると，以下の出力を得ます7_．

Estimate Std. Error t value Pr(>¦t¦) (Intercept) 817.721950 197.3887800 4.1426972 0.0001576243 pop 7.754188 1.8237001 4.2518988 0.0001120459 area 6.728679 0.3641458 18.4779803 0.0000000000 yoto -6.849713 10.5763026 0.6476473 0.5206562177 7 _{同様の標準誤差は car パッケージの hccm で type=c("hc1")と設定することによっても得られま} す（注 6 参照）．また Design パッケージの robcov はサンプルの有限性を考慮しない標準誤差を求 めます（注 10 参照）．

13

【第 8 章】

変数名を yr（年），doro（道路整備費），kokyo（その他の公共事業費），giseki（衆議院にお ける与党議席率）とします．まず線形モデルとして ols<-lm(doro kokyo+giseki) を定義します．R では lmtest パッケージの dwtest によってダービン・ワトソン検定ができます8_． 具体的には，lmtest パッケージをダウンロードしたうえで， library(lmtest) dwtest(ols) を実行すると以下の出力を得ます． Durbin-Watson test data: ols DW = 0.5123, p-value = 4.109e-09

alternative hypothesis: true autocorrelation is greater than 0

以下では，GLS 反復推定（Prais-Winsten 法）を pw.ar1 というプログラムにまとめています （SPSS による分析事例と同様にダービン・ワトソン比からρを求める方法です）．手順として は，pw.ar1 を実行したうえで，それを推定モデルに適用することになります． pw.ar1<-function(gls) { X<-model.matrix(gls) y<-model.response(model.frame(gls)) e<-resid(gls) n<-length(e) dw<-sum((e[1:(n-1)]-e[2:n])^2)/sum(e^2) r<-1-dw/2 #1 original<-append(r, dw) repeat { yg<-c(y[1]*(1-r^2)^0.5, y[2:n]-r*y[1:(n-1)]) Xg<-rbind(X[1,]*(1-r^2)^0.5, X[2:n,]-r*X[1:(n-1),]) gls<-lm(yg Xg-1) e<-y-X%*%coef(gls) dw<-sum((e[1:(n-1)]-e[2:n])^2)/sum(e^2) rg<-1-dw/2 #2

if (delta<1e-3) break } #3 eg<-resid(gls) dwg<-sum((eg[1:(n-1)]-eg[2:n])^2)/sum(eg^2) rd<-rbind(original, final=append(r, dwg)) colnames(rd)<-c("RHO", "DW") print(rd) summary(gls) } pw.ar1(lm(doro kokyo+giseki)) 出力は以下のとおりです（事前に定義した推定モデル，たとえば ols があれば，pw.ar1(ols) とすることもできます）． RHO DW original 0.7438477 0.5123046 final 0.8707564 1.6711165 Call: lm(formula = yg Xg - 1) Residuals:

Min 1Q Median 3Q Max -2248.0 -327.2 117.7 447.0 1160.0 Coefficients:

Estimate Std. Error t value Pr(>¦t¦) Xg(Intercept) 3049.78187 3224.69135 0.946 0.351 Xgkokyo 0.35837 0.01407 25.473 <2e-16 *** Xggiseki -33.12387 51.63346 -0.642 0.526 ---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 756.5 on 32 degrees of freedom

Multiple R-Squared: 0.9686, Adjusted R-squared: 0.9657 F-statistic: 329.1 on 3 and 32 DF, p-value: < 2.2e-16

またe_t =ρe_t−1+ε_tの推定からρを求める方法は，上記のダービン・ワトソン比からρを求 める手順の#1 と#2 を以下に代えます． r<-coef(lm(c(e[2:n]) c(e[1:(n-1)])-1)) #1 rg<-coef(lm(c(e[2:n]) c(e[1:(n-1)])-1)) #2 ここではρの変化が 0.001 を下回るまで反復推定するよう設定しています．反復終了基準 は#3 行の 1e-3 によって変更できます． News 2(3), 7-10. URL http://CRAN.R-project.org/doc/Rnews/

15

【第 9 章】

変数名を shiji（自民支持），ondo（感情温度），chiiki（地域改善）とします．R において最尤法 によるロジスティック回帰を行なうにはいくつかの方法がありますが9_{，ここでは Design パッケ} ージの lrm を用います10_{．まず推定モデルの構造を簡略化しておきます．} X<-model.matrix(ols<-lm(shiji ondo+chiiki)) y<-model.response(model.frame(ols)) Design パッケージをダウンロードしたうえで， library(Design) print(logit<-lrm(ols)) を実行すると，以下が出力されます． Logistic Regression Model

lrm(formula = ols)

Frequencies of Responses 0 1

32 18

Obs Max Deriv Model L.R. d.f. P C Dxy 50 7e-05 26.06 2 0 0.891 0.781 Gamma Tau-a R2 Brier

0.812 0.367 0.557 0.117 Coef S.E. Wald Z P Intercept -7.8085 2.23003 -3.50 0.0005 ondo 0.1151 0.03638 3.16 0.0016 chiiki 1.8550 0.86575 2.14 0.0321 推定係数の有意性に関しては Z 統計量による正規分布確率を報告します．また出力中央 にカイ 2 乗尤度比検定も示されていますが（ちなみに，その下段の R2 はNagelkerke R2_です）， 以下を続けて実行すると，定数項のみによる対数尤度（Lc），最大化対数尤度（Lk），カイ 2 乗 尤度比検定（PrX），擬似決定係数（PR2），50％を基準とする正判別率（PCP）を整理して出力 することができます． n<-length(y) r<-sum(y) 9 _{非線形モデル nlm や一般化線形モデル glm によってロジスティック回帰を行うこともできます．} 10 _{Frank E Harrell Jr (2005). Design: Design Package. R package version 2.0-12.}

PrX<-1-pchisq(cbind(logit$stats)[3], cbind(logit$stats)[4]) PR2<--(cbind(logit$stats)[3]/2)/Lc

yh<-c(round(1/(1+exp(-X%*%coefficients(logit))))) PCP<-sum(y==yh)/n

results<-cbind(Lc, Lk, PrX, PR2, PCP) rownames(results) <-rep("", times=1) results Lc Lk PrX PR2 PCP -32.67091 -19.63894 2.189208e-06 0.3988861 0.8 さらに，自民に対する感情温度が中立である有権者に関して，地域改善志向の有無による 自民支持確率は， X1<-cbind(1, 50, chiiki<-c(0,1)) P_shiji<-c(1/(1+exp(-X1%*%coefficients(logit)))) P1<-cbind(chiiki, P_shiji) rownames(P1)<-rep("", times=2) P1 とすることによって求められ，以下の出力を得ます． chiiki P_shiji 0 0.1136111 1 0.4503220 あるいは， X2<-cbind(1, ondo<-c(0:100), 1) p2<-c(1/(1+exp(-X2%*%coefficients(logit)))) X3<-cbind(1, ondo<-c(0:100), 0) p3<-c(1/(1+exp(-X3%*%coefficients(logit))))

plot(ondo, p2, type='l', lty=1,ylab='Pr(LDP shiji=1)', xlab='LDP ondo') points(ondo, p3, type='l', lty=2)

を実行すると，図 R-4 のように自民に対する感情温度に応じた自民支持の予測確率を地域 改善志向の有る有権者（実線）と無い有権者（破線）で図示することができます．

17 図 R-4 R による自民支持の確率変化 0 20 40 60 80 100 0. 0 0 .2 0 .4 0 .6 0. 8 1 .0 LDP ondo P r(L D P s h iji = 1 )