準備：Ross2012のデータを読み込む。

ワーキングディレクトリは各自で設定すること。

load("Replication data for The Oil Curse - Ross 2012.RData")

浜中コードの確認

コメントアウト(“#”記号以下のコメント)は松尾による。

#ライブラリの読み込み
library(dplyr)
#グループの作成
IMFa<- filter(x, year==2003)#データの中から2003年だけ抽出
IMFaa <- filter(IMFa, gdpcap2000_sup<5000)#低所得国を抽出
IMFab <- filter(IMFa, gdpcap2000_sup>5000)#高所得を抽出
#1行目（低所得国のt検定)
t.test(IMFaa$IMFrev[IMFaa$oil_gas_value100==1],
       IMFaa$IMFrev[IMFaa$oil_gas_value100==0])
#2行目(高所得国のt検定)
t.test(IMFab$IMFrev[IMFab$oil_gas_value100==1],
       IMFab$IMFrev[IMFab$oil_gas_value100==0])
#3行目(すべてのt検定)
t.test(IMFa$IMFrev[IMFa$oil_gas_value100==1],
       IMFa$IMFrev[IMFa$oil_gas_value100==0])

t検定の比較対象が、表とは逆の順番(表では「非産油国と産油国」、浜中コードでは「産油国と非産油国」)になっている点に注意。t検定のコード(“t.test”)にはよりシンプルな書式がある。これらの点については以下に浜中コードを各行ごとに確認してゆく中で記述する。

1行目(低所得国）の確認

    Welch Two Sample t-test

data:  IMFaa$IMFrev[IMFaa$oil_gas_value100 == 1] and IMFaa$IMFrev[IMFaa$oil_gas_value100 == 0]
t = 2.9636, df = 35.106, p-value = 0.005431
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  2.060337 11.018722
sample estimates:
mean of x mean of y 
 27.74868  21.20915 

この結果はロスの表2.1の1行目と同じになるが、表に現れない細かい部分でロスとは異なる。この分析ではt値が2.9636、dfが35.106、p値が0.005431だが、ロスがSTATAで行った分析ではt値が-2.6513、dfが107、p値が0.0092である。この違いは、ロスが等分散を前提とした検定を行っているのに対して、浜中コードが等分散を前提としないWelchのt検定を行っている為である。また、浜中コードでは比較対象の配列順序が表とは逆なので、平均値の表基準も逆になっている点に注意。

浜中コードを等分散を前提としたコードに改め、比較対象の配列順序を表に揃えると、結果はロスと同じになる。

t.test(IMFaa$IMFrev[IMFaa$oil_gas_value100==0],
       IMFaa$IMFrev[IMFaa$oil_gas_value100==1],var.equal = T)

    Two Sample t-test

data:  IMFaa$IMFrev[IMFaa$oil_gas_value100 == 0] and IMFaa$IMFrev[IMFaa$oil_gas_value100 == 1]
t = -2.6513, df = 107, p-value = 0.009236
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -11.429108  -1.649951
sample estimates:
mean of x mean of y 
 21.20915  27.74868 

なお、上記のコードに含まれるt検定のコードは、以下のようによりシンプルに表記できる。結果は上記の改定コードと同じになるので、ここでは表記しない。

t.test(IMFrev~oil_gas_value100,data = IMFaa,var.equal = T)

2行目(高所得国のt検定)の確認

    Welch Two Sample t-test

data:  IMFab$IMFrev[IMFab$oil_gas_value100 == 1] and IMFab$IMFrev[IMFab$oil_gas_value100 == 0]
t = 0.9833, df = 11.595, p-value = 0.3455
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5.941541 15.646605
sample estimates:
mean of x mean of y 
 40.54099  35.68846 

これはロスとかなり違う結果が出る。p値が0.34なので、統計的には有意ではないが、表では５％水準で有意となっている。ロスの分析ではp値は0.0319であり、浜中コードとはかなり違う結果が出る。ただし、『石油の呪い』の議論は「高所得国では政府の規模に対する石油の効果が弱い」というものなので、浜中コードの結果がロスの議論を否定する、というほとではない。

浜中コードとロスの分析結果の違いは、等分散を前提としているかどうかとは別の問題である。二つ目のt検定でも、ロスは等分散を前提とした検定を行なっている。これに合わせて浜中コードを変更しても、分析結果がロスのものと一致することはない(等分散に関する修正と同時に、“t.test”のコードの書式も上記1行目の分析と同様に改めてある)。

t.test(IMFrev~oil_gas_value100,data = IMFab,var.equal = T)

    Two Sample t-test

data:  IMFrev by oil_gas_value100
t = -1.131, df = 26, p-value = 0.2684
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -13.671909   3.966845
sample estimates:
mean in group 0 mean in group 1 
       35.68846        40.54099 

ロスの分析では、平均値は32.819（非産油国）と44.576（産油国）だが、浜中コード(および上記の改定コード)では35.688（非産油国）、40.541（産油国）とかなり違う。平均値が違うということは、平均値を求めるサンプルが異なることを意味している。ロスの分析では高所得国のサンプル数は32になっているが、Rで確認するとサンプル数は46あり、この中でNAが18で、NAを除くと28になる。

table(is.na(IMFab$IMFrev))

FALSE  TRUE 
   28    18 

問題はSTATAが行うt検定におけるNAの処理だと思われる。STATA用のデータファイルとR用のデータファイルの内容が異なるかもしれないので、念の為STATA用のファイルをRで読み込んでR用のファイルとの違いを確認したが、高所得国の産油国・非産油国の部分については、違いは見られなかった。

3行目(すべてのt検定)のt検定

    Welch Two Sample t-test

data:  IMFa$IMFrev[IMFa$oil_gas_value100 == 1] and IMFa$IMFrev[IMFa$oil_gas_value100 == 0]
t = 3.4009, df = 41.302, p-value = 0.001502
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  3.919061 15.372301
sample estimates:
mean of x mean of y 
 33.17683  23.53115 

1行目のt検定とほぼ同じ結果。表にすると違いが現れないが、細かいところではロスの分析と違う結果が出る。浜中コードではp値が0.002なので1%水準で有意となり、表と変わらない（t値は表にない）。一方のロスではp値が0.0002となっている。この違いは、一つ目のt検定と同様に、ロスが等分散を前提としたt検定を実施しているためである。

t検定の比較対象の配列順を表に従って修正し、“t.test”のコードをシンプルなものに変更しつつ、浜中コードを等分散を前提としたものに変更すると、ロスの検定と結果が同じになる。

t.test(IMFrev~oil_gas_value100,data=IMFa,var.equal = T)

    Two Sample t-test

data:  IMFrev by oil_gas_value100
t = -3.8634, df = 139, p-value = 0.0001709
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -14.582126  -4.709236
sample estimates:
mean in group 0 mean in group 1 
       23.53115        33.17683 

これでロスの分析結果と同じになった。

念の為の等分散検定

ロスは３つの検定を全て等分散を前提に行なっている。これは妥当なのか、念の為確認する

var.test(IMFrev~oil_gas_value100,data = IMFaa)#低所得国

    F test to compare two variances

data:  IMFrev by oil_gas_value100
F = 1.4358, num df = 87, denom df = 20, p-value = 0.3612
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.6579836 2.6847364
sample estimates:
ratio of variances 
          1.435791 

var.test(IMFrev~oil_gas_value100,data=IMFab,)#高所得国

    F test to compare two variances

data:  IMFrev by oil_gas_value100
F = 0.45677, num df = 18, denom df = 8, p-value = 0.1597
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.1132362 1.3727109
sample estimates:
ratio of variances 
         0.4567677 

var.test(IMFrev~oil_gas_value100,data = IMFa)#全て

    F test to compare two variances

data:  IMFrev by oil_gas_value100
F = 0.63329, num df = 109, denom df = 30, p-value = 0.09288
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.3379181 1.0781788
sample estimates:
ratio of variances 
         0.6332903 

どれも５％水準で等分散の仮定を棄却できないので、等分散を前提とすることには問題はない。等分散検定を行わずにt検定を行うという方法も存在するので(例えば青木)、t検定という観点では浜中コードも間違いではない。ここでは、ロスの分析を正確に再現することが目的なので、ロスが行なっているのと同様に、等分散を前提としてt検定を行う。