## Warning: package 'knitr' was built under R version 3.0.3
plm パッケージの Cigar というパネルデータを使って、米国の各州におけるタバコの消費量を規定する要因を分析したいとする。最終的には観察されない異質性を統制したモデルを推定したいが、その前にやるべき探索的分析をせよ。
特に人口規模の及ぼす影響について検証するとし、コントロール変数としてタバコの価格、週の一人当たり可処分所得、周辺の州のタバコの最低価格を用いる。
library(plm)
data(Cigar)
head(Cigar)
## state year price pop pop16 cpi ndi sales pimin
## 1 1 63 28.6 3383 2236.5 30.6 1558.305 93.9 26.1
## 2 1 64 29.8 3431 2276.7 31.0 1684.073 95.4 27.5
## 3 1 65 29.8 3486 2327.5 31.5 1809.842 98.5 28.9
## 4 1 66 31.5 3524 2369.7 32.4 1915.160 96.4 29.5
## 5 1 67 31.6 3533 2393.7 33.4 2023.546 95.5 29.6
## 6 1 68 35.6 3522 2405.2 34.8 2202.486 88.4 32.0
tab.state <- xtabs(~state, Cigar)
tab.state
## state
## 1 3 4 5 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 35 36 37
## 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30
## 39 40 41 42 43 44 45 46 47 48 49 50 51
## 30 30 30 30 30 30 30 30 30 30 30 30 30
tab.year <- xtabs(~year, Cigar)
tab.year
## year
## 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
## 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46 46
このデータは米国の46の州に関する1963年から1992年までの29年分のパネルデータでバランスが取れている。用いる変数の記述統計は以下の通りである。
d0 <- Cigar[, -5:-6] # いらない変数を削除したデータを作る
summary(d0[, -1:-2]) # state と year の記述統計はいらない
## price pop ndi sales pimin
## Min. : 23.40 Min. : 319 Min. : 1323 Min. : 53.4 Min. : 23.40
## 1st Qu.: 34.77 1st Qu.: 1053 1st Qu.: 3328 1st Qu.:107.9 1st Qu.: 31.98
## Median : 52.30 Median : 3174 Median : 6281 Median :121.2 Median : 46.40
## Mean : 68.70 Mean : 4537 Mean : 7525 Mean :124.0 Mean : 62.90
## 3rd Qu.: 98.10 3rd Qu.: 5280 3rd Qu.:11024 3rd Qu.:133.2 3rd Qu.: 90.50
## Max. :201.90 Max. :30703 Max. :23074 Max. :297.9 Max. :178.50
# tapply(d0$sales, d0$year,summary) # やたら長いので出力は省略
タバコの消費量の級内相関係数は、以下のようになる。あまり大きな値ではないが、無視できるというほどでもない。
aov1 <- aov(sales ~ state, d0)
library(multilevel)
ICC1(aov1)
## [1] 0.01980365
次にタバコの消費量の変化をグラフにしたのが、以下の図1 である。70年ごろ一度落ち込み、その後70年代後半に上昇し、その後は下降る傾向が続いている。
plot(sales ~ year, data=d0)
id.state <- unique(d0$state) # 州IDのリストを作る
for(i in 1:length(id.state)){ # 州別に線を引く
d1 <- subset(d0, state==id.state[i]) # 各州のデータフレームを作成
lines(d1$year, d1$sales) # 州別のデータフレームで線を描画
}
次に各州の人口規模とタバコの一人当たり消費量の散布図を示したのが、図2 である。人口規模が大きいとタバコの一人当たり消費量は少なくなるようにも見えるが、明らかに一部の例外的に人口規模の大きい州(5と17)とタバコ消費量の多い州がそのような印象を強めている点に注意が必要だろう。
plot(sales ~ pop, data=d0)
for(i in 1:length(tab.state)){
d1 <- subset(d0, state==id.state[i])
lines(d1$pop, d1$sales)
}
sort(tapply(d0$pop, d0$state, mean)) # どの州が外れ値かチェック
## 51 46 8 35 9 42 29 27 13
## 414.8433 489.1133 589.6800 650.4133 689.9800 692.2533 733.8267 761.7433 867.2900
## 30 40 20 32 45 28 49 4 17
## 874.4333 946.9000 1091.5667 1244.3700 1352.0033 1546.7067 1853.2833 2172.1400 2362.7433
## 25 3 16 37 41 7 18 1 48
## 2448.9667 2499.1700 2848.6900 2889.5333 2998.1133 3080.8600 3482.0233 3767.5933 3894.0500
## 24 19 21 43 50 26 47 11 15
## 3991.1967 4020.2000 4149.6500 4379.5733 4586.4533 4851.0100 5220.7100 5305.5800 5329.9000
## 22 31 23 10 36 14 39 44 33
## 5704.9200 7330.1033 8986.7567 9083.4900 10691.2733 11259.3400 11832.7567 13651.5733 17941.1567
## 5
## 23149.2733
sort(tapply(d0$sales, d0$state, mean))
## 45 32 48 42 35 13 1 28 25 24
## 67.87333 92.06000 95.78333 103.27667 104.35333 104.81667 107.47667 107.64667 108.13333 108.34000
## 50 27 44 16 17 3 49 5 39 4
## 108.58000 111.33000 111.51000 111.72333 112.66333 113.24667 114.75667 114.85000 115.01000 115.33667
## 7 43 11 33 41 22 37 31 19 21
## 118.06667 118.96000 119.38667 119.48333 119.51333 119.82333 120.48000 120.59000 124.02333 124.42333
## 10 14 36 23 26 20 47 40 51 46
## 126.04667 126.24333 126.39667 129.85333 132.48000 133.56000 134.34333 135.23000 136.49667 140.69000
## 15 8 9 29 18 30
## 142.12333 150.62667 164.22000 173.41333 179.97667 236.52333
次に州内平均からの偏差をプロットしたのが図3である。図3をみると、やはり州平均を統制すると、ほとんど相関がなくなってしまうのがわかる。
p.d0 <- pdata.frame(d0) # pdata.frame 形式にする。これやると偏差の計算が簡単
plot(Within(p.d0$sales) ~ Within(p.d0$pop)) # Within は p.d0$変数名という形で引数指定
for(i in 1:length(tab.state)){
d1 <- subset(p.d0, state==id.state[i])
lines(Within(d1$pop), Within(d1$sales))
}
plm パッケージの Crime というパネルデータを使って、ノースカロライナ州の各郡における犯罪の発生率を規定する要因を分析したいとする。最終的には観察されない異質性を統制したモデルを推定したいが、その前にやるべき探索的分析をせよ。