パネルデータの探索的分析
太郎丸 博
Thursday, May 21, 2015
- 1 Why Exploration Is Necessary 分布のチェックと探索的分析の重要性
- 2 Checking Descrptive Statistics 記述統計のチェック
- 3 Plot グラフで分布をチェック
1 Why Exploration Is Necessary 分布のチェックと探索的分析の重要性
- 最終目的は内生性や観察されない異質性をコントロールした正確な因果関係の推定
- それゆえ、線形加法モデル（固定効果モデルやランダム効果モデル）の結果が重要
- しかし、用いる変数の分布がどうなっているのか、知っておくことは重要
- モデルの特定の失敗を避けられるかも
- 結果の解釈の間違いを避けられるかも
1.1 An Example of Wrong Model Specification 探索的分析を怠ったせいで生じるモデル特定の失敗例
## Warning: package 'knitr' was built under R version 3.1.3Y という変数をXという連続変数とZというダミー変数で予測した。
summary(data.frame(Y, X, Z))## Y X Z
## Min. :-10.51 Min. : 1.0 京:50
## 1st Qu.:325.78 1st Qu.:13.0 阪:50
## Median :432.15 Median :25.5
## Mean :387.71 Mean :25.5
## 3rd Qu.:490.14 3rd Qu.:38.0
## Max. :629.14 Max. :50.0l1 <- lm(Y ~ X + Z)
summary(l1)##
## Call:
## lm(formula = Y ~ X + Z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -233.58 -86.42 17.37 90.02 197.43
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 225.5809 25.5180 8.840 4.27e-14 ***
## X 7.1080 0.7813 9.098 1.19e-14 ***
## Z阪 -38.2489 22.5496 -1.696 0.0931 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 112.7 on 97 degrees of freedom
## Multiple R-squared: 0.4689, Adjusted R-squared: 0.458
## F-statistic: 42.82 on 2 and 97 DF, p-value: 4.683e-14上の結果を見ると、Xは有意な正の効果を持ち、Zは有意ではない。Xが大きくなるのに比例してYも大きな値をとる、と考えていいだろうか。しかし、XとYの散布図を見てみると、放物線を描くような関係になっており、単純な線形関係ではないことがわかる。
plot(Y ~ X, pch=paste(Z))
legend("topleft", pch=levels(Z), legend=c(": Z = 京", ": Z = 阪"))
図1 架空の変数 Y と X の散布図
そこで、Xの二乗項をモデルに投入すると、以下のようにモデルの適合度は上昇し、Zの効果も有意になった。
l2 <- update(l1, ~. + I(X^2))
summary(l2)##
## Call:
## lm(formula = Y ~ X + Z + I(X^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -126.065 -36.666 1.953 40.428 140.238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.5972 18.6273 -0.139 0.88940
## X 33.4362 1.6047 20.836 < 2e-16 ***
## Z阪 -38.2489 11.3571 -3.368 0.00109 **
## I(X^2) -0.5162 0.0305 -16.923 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56.79 on 96 degrees of freedom
## Multiple R-squared: 0.8667, Adjusted R-squared: 0.8625
## F-statistic: 208 on 3 and 96 DF, p-value: < 2.2e-16散布図を見ておけば、最初のモデルを採択するような失敗は犯さないだろう。
1.2 An Example of Wrong Model Interpretation 探索的分析を怠ったせいで生じる結果の解釈の失敗例
Yを幸福感、Xを収入（単位は百万円）、Zを居住地とする。回帰分析すると以下のような結果が得られた。
summary(data.frame(Y, X, Z))## Y X Z
## Min. :2.00 Min. : 0.002 京:5000
## 1st Qu.:4.00 1st Qu.: 2.502 阪:5000
## Median :5.00 Median : 5.001
## Mean :5.16 Mean : 5.001
## 3rd Qu.:6.00 3rd Qu.: 7.500
## Max. :9.00 Max. :10.000lm3 <- lm(Y ~ X +Z)
summary(lm3)##
## Call:
## lm(formula = Y ~ X + Z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2872 -1.0278 -0.0697 0.7536 3.7333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.231518 0.023211 225.386 <2e-16 ***
## X 0.006215 0.003595 1.729 0.0839 .
## Z阪 -0.205200 0.020758 -9.885 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.038 on 9997 degrees of freedom
## Multiple R-squared: 0.009973, Adjusted R-squared: 0.009775
## F-statistic: 50.35 on 2 and 9997 DF, p-value: < 2.2e-16X は正の有意な効果を持つので、収入が高まるほど幸福感が上がる、という解釈は間違ってはいない。しかし、実はこの収入の効果は非常に小さいというべきである。散布図を見てみよう。
plot(jitter(Y) ~ X, pch=20, cex=0.3)
図2 収入 (X) と幸福度 (Y) の散布図
明確な右肩上がりの傾向など見られないことがわかるだろう。Xの最小値は、0.002 （0.2万円）、最大値は10（1000万円）、傾きが0.006なので、最低の収入の人と最高の収入の人の平均的な幸福度の差は、0.062にすぎない。幸福感が2から9の間で分布していることを考えると、Xの効果は小さいというべきだろう。
このような「解釈」は回帰分析の結果だけを見ていても出てこない。分析に用いた変数の分布、特にバラつきの大きさについて把握しておくことが重要である。
このような探索的分析はパネルデータの分析においても同じように重要である。
2 Checking Descrptive Statistics 記述統計のチェック
探索的分析は、記述統計の検討と、グラフなどを使った複数の変数の関連の検討に分けられるだろう。グラフが視覚的に分布を認識しようとするのに対して、記述統計のチェックはデータの基本的な特徴を数値で確認しておく作業である。
2.1 Checking N and T 個体数や時点数の確認
まず、パネルデータの基本構造は個体 X 時点であり、それぞれの数やデータがバランスしているか（すべての個体に関してすべての時点のデータが得られている場合バランスしているという）、を確認する。
library(plm)# パッケージのロード
data(Crime) # サンプルデータの呼び出し
head(Crime) # データの最初の数行を見てみる## county year crmrte prbarr prbconv prbpris avgsen polpc density taxpc region smsa pctmin
## 1 1 81 0.0398849 0.289696 0.402062 0.472222 5.61 0.0017868 2.307159 25.69763 central no 20.2187
## 2 1 82 0.0383449 0.338111 0.433005 0.506993 5.59 0.0017666 2.330254 24.87425 central no 20.2187
## 3 1 83 0.0303048 0.330449 0.525703 0.479705 5.80 0.0018358 2.341801 26.45144 central no 20.2187
## 4 1 84 0.0347259 0.362525 0.604706 0.520104 6.89 0.0018859 2.346420 26.84235 central no 20.2187
## 5 1 85 0.0365730 0.325395 0.578723 0.497059 6.55 0.0019244 2.364896 28.14034 central no 20.2187
## 6 1 86 0.0347524 0.326062 0.512324 0.439863 6.90 0.0018952 2.385681 29.74098 central no 20.2187
## wcon wtuc wtrd wfir wser wmfg wfed wsta wloc mix pctymle
## 1 206.4803 333.6209 182.3330 272.4492 215.7335 229.12 409.37 236.24 231.47 0.0999179 0.0876968
## 2 212.7542 369.2964 189.5414 300.8788 231.5767 240.33 419.70 253.88 236.79 0.1030491 0.0863767
## 3 219.7802 1394.8030 196.6395 309.9696 240.1568 269.70 438.85 250.36 248.58 0.0806787 0.0850909
## 4 223.4238 398.8604 200.5629 350.0863 252.4477 281.74 459.17 261.93 264.38 0.0785035 0.0838333
## 5 243.7562 358.7830 206.8827 383.0707 261.0861 298.88 490.43 281.44 288.58 0.0932486 0.0823065
## 6 257.9139 369.5465 218.5165 409.8842 269.6129 322.65 478.67 286.91 306.70 0.0973228 0.0800806?Crime # サンプルデータのヘルプを見て、どの変数が個体IDと時点なのかチェック
tab.county <- xtabs(~county, data=Crime) # 個体IDの度数分布表を作る
tab.county## county
## 1 3 5 7 9 11 13 15 17 19 21 23 25 27 33 35 37 39 41 45 47 49 51 53 55 57 59
## 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## 61 63 65 67 69 71 77 79 81 83 85 87 89 91 93 97 99 101 105 107 109 111 113 115 117 119 123
## 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## 125 127 129 131 133 135 137 139 141 143 145 147 149 151 153 155 157 159 161 163 165 167 169 171 173 175 179
## 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## 181 183 185 187 189 191 193 195 197
## 7 7 7 7 7 7 7 7 7tab.year <- xtabs(~year, data=Crime) # 時点の度数分布表を作る
tab.year## year
## 81 82 83 84 85 86 87
## 90 90 90 90 90 90 90このデータの場合、90の郡に関して7時点分のデータがあることがわかる。個体数や時点数が多い場合、上のようなやり方では表示しきれないので、時点や郡IDの度数分布表を作ったあと下記のようにすればいい。
length(tab.county) # 郡の数 = 表の長さ## [1] 90tab.county[1:10] # 度数分布表の最初の10個の数値を示す## county
## 1 3 5 7 9 11 13 15 17 19
## 7 7 7 7 7 7 7 7 7 72.2 Descrbing Panel Data パネルデータの記述統計
分析に用いる変数の平均値、最小値、最大値、標準偏差といった基本的な統計量を概観しておくことは、多変量解析のような複雑なモデルの推定の前に必ずやっておくべき基本である。パネルデータの場合、ロング形式で（つまりプール・データで）全体の統計量を計算するだけでなく、時点数が少なければ時点別に記述統計を見ておきたい。
2.2.1 Running R: R での計算
- まず、用いる変数のみからなるデータ・フレームを作る。
d0 <- subset(Crime, select=c(county, year, crmrte, smsa, pctmin, wcon, pctymle))
head(d0)## county year crmrte smsa pctmin wcon pctymle
## 1 1 81 0.0398849 no 20.2187 206.4803 0.0876968
## 2 1 82 0.0383449 no 20.2187 212.7542 0.0863767
## 3 1 83 0.0303048 no 20.2187 219.7802 0.0850909
## 4 1 84 0.0347259 no 20.2187 223.4238 0.0838333
## 5 1 85 0.0365730 no 20.2187 243.7562 0.0823065
## 6 1 86 0.0347524 no 20.2187 257.9139 0.0800806- 全体で計算した記述統計を計算
summary(d0)## county year crmrte smsa pctmin wcon
## Min. : 1.0 Min. :81 Min. :0.001812 no :574 Min. : 1.284 Min. : 65.62
## 1st Qu.: 51.0 1st Qu.:82 1st Qu.:0.018352 yes: 56 1st Qu.:10.005 1st Qu.: 201.66
## Median :103.0 Median :84 Median :0.028441 Median :24.852 Median : 236.46
## Mean :100.6 Mean :84 Mean :0.031588 Mean :25.713 Mean : 245.67
## 3rd Qu.:151.0 3rd Qu.:86 3rd Qu.:0.038406 3rd Qu.:38.223 3rd Qu.: 269.69
## Max. :197.0 Max. :87 Max. :0.163835 Max. :64.348 Max. :2324.60
## pctymle
## Min. :0.06216
## 1st Qu.:0.07859
## Median :0.08316
## Mean :0.08897
## 3rd Qu.:0.08919
## Max. :0.27436apply(d0[,c(-1, -2, -4)], 2, sd, na.rm=T) # 1,2列目と4列目の smsa は標準偏差を計算しても意味がないので除外## crmrte pctmin wcon pctymle
## 0.01812087 16.90353614 121.98372022 0.02434931- 時点数が少なければ時点別に記述統計を計算。被説明変数に関しては特にやっておきたい。
# tapply(d0$crmrte, d0$year, summary) # 長いので出力は省略
tapply(d0$crmrte, d0$year, sd, na.rm=T)## 81 82 83 84 85 86 87
## 0.01700389 0.01713244 0.01713540 0.01676258 0.01696727 0.02243483 0.01888699# 以下同様に他の変数についても計算すればよい2.3 Intraclass Correlation 級内相関
固定効果モデルを推定する場合、個人内の変数の分散も個人間の変数の分散も、両方ある程度あるのが大前提である。まずは、級内分散と級間分散（あるいは個人が単位のデータであれば個人内分散と個人間分散）について解説しよう。
2.3.1 Intraclass Variance 級内分散
個人 \(i \; (i=1, \; \ldots,\; N)\) の時点 \(t (t=1, \; \ldots, \; T)\) における変数 \(X\) の値を \(X_{it}\) とすると、個人 \(i\) の個人内平均 \(\bar X_i\) を
\(\bar X_i =\frac{\sum_{t=1}^T X_{it} }{T}\)
と定義する。
例えば、下の表のようなデータが得られた場合、
# 表1 架空のパネルデータ（4人 X 3時点)
data.frame(i, t, X, Y)## i t X Y
## 1 1 1 -1.8 -4.2
## 2 1 2 1.7 1.2
## 3 1 3 2.2 -2.8
## 4 2 1 1.6 -3.0
## 5 2 2 8.4 2.7
## 6 2 3 4.4 -2.1
## 7 3 1 3.5 -2.7
## 8 3 2 7.7 2.7
## 9 3 3 1.0 -4.5
## 10 4 1 3.6 -1.7
## 11 4 2 5.9 -0.5
## 12 4 3 8.1 -0.2\(i=1\) の人の \(X\) の個人内平均は、
\(\bar X_1 = \frac{-1.8 + 1.7 + 2.2}{3} = 0.7\)
である。個人 \(i\) の個人内分散 \(s^2_i\) は、
\(s^2_i = \frac{\sum_{t=1}^T (X_{it} - \bar X_i)^2}{T-1}\)
で定義する。例えば、上の表1の一人目の人の個人内分散は、
\(s_1^2 = \frac{(-1.8 - 0.7)^2 + (1.7 - 0.7)^2 + (2.2 - 0.7 )^2}{3 - 1} = 4.8\)
である。上の個々人の個人内分散をすべての個人に関して計算し、その平均値を計算したものを単に個人内分散\(s^2_\mathrm{within}\)と呼ぶ。表1の \(X\)の場合、4人の個人の個人内分散は、
ss.i <- tapply(X, i, var)
eq03 <- round(ss.i, 1)
eq03## 1 2 3 4
## 4.8 11.7 11.5 5.1である。それゆえ個人内分散は、
\(s^2_\mathrm{within}= \frac{4.8+11.7+11.5+5.1}{4}= 8.24\)
mean(ss.i)## [1] 8.239167である。
2.3.2 Between-Class Variance 個人間分散
次に個人間分散は個人内平均の分散である。すなわち、\(s^2_{\mathrm{between}}\)を個人間分散とすると、
\(s^2_\mathrm{between} = \frac{\sum_{i=1}^N (\bar X_i - \bar X)^2}{N-1}\)
で定義される。ただし、\(\bar X\) は全体の平均である。
\(\bar X =\frac{\sum X_{it}}{T \times N}\)
ただし、\(\sum\) はすべての個人、時点に関して足し合わせるという意味である。
例えば、表1の例でいえば、\(X\)の全体平均は、
\(\bar X =\frac{-1.8 + 1.7 + 2.2, \; \ldots, \; + 8.1}{3 \times 4} =3.9\)
である。また、4人の個人内平均は、
m.in <- tapply(X, i, mean)
m.in## 1 2 3 4
## 0.700000 4.800000 4.066667 5.866667であるから、\(X\)の個人間分散は
\(s^2_\mathrm{between}=\frac{0.7 - 3.9)^2 + (4.8 - 3.9)^2 + (4.1 - 3.9)^2 + (5.9 - 3.9 )^2}{4-1}\)
var(m.in)## [1] 4.979537である。
2.3.3 Definition of Intraclass Correlation 級内相関
級間分散と級内分散の和で級間分散を割った値を級内相関 (intraclass correlation) と呼ぶ。すなわち、
級内相関 \(=\frac{s^2_\mathrm{between}}{s^2_\mathrm{between} + s^2_\mathrm{within}}\)
である。表1の \(X\) の例では、
var(m.in) / (var(m.in) + mean(ss.i))## [1] 0.3767039もしも個人内の分散がなければ（つまり級内相関が 1 ならば）、一時点分のデータでふつうに回帰分析したほうがいいかもしれない。逆に個人間の分散がまったくなければ（つまり級内相関が 0 ならば）プール・データでふつうに回帰分析したほうがいいかもしれない。つまり、固定効果モデルやランダム効果モデルのような手法に意味があるのは、級内相関が 0 からも 1 からもある程度離れているような場合である。
上の定義が最近のパネルデータ分析の教科書やマルチレベル・モデルでは一般的であるが、級内相関係数にはほかにもいくつか定義がある。R でも multilevelパッケージや psych パッケージ、 ICCパッケージなど、級内相関係数を計算する関数がいくつかあるのだが、私の知る限りどれも上の定義とは異なる計算法であった。どれも級間分散が級内分散に比べてどれくらい大きいかを示す数値ではあるが、やはり最近の一般的な定義にしたがったほうがいいと思われる。そこで、以下のような関数を作ってみたので活用してほしい。
icc1 <- function(X, id){ # ここから icc1 という関数の定義
aov1 <- aov(X ~ -1 + factor(id))
v.btwn <- var(coef(aov1))
v.wthn <- sum(resid(aov1)^2)/aov1$df.residual
return(v.btwn/(v.wthn + v.btwn))
} # ここまでコピペしてまず実行
icc1(X, i) # 変数名 ~ -1 + 個体 id で計算できる。欠損値のないデータで## [1] 0.3767039上の icc1() という関数の定義をいちいちコピペしてから実行する必要があるので、面倒だが、いちいち自分で計算するよりはマシだろう。
3 Plot グラフで分布をチェック
分布のチェックはぜひグラフでやりたい。私が学生の頃は度数分布表やクロス表のチェックが重要だと習ったが、度数分布表やクロス表を見て分布のイメージをつかむのは、難しいことも多い。特に連続変数をあつかう場合、カテゴリのわけ方でずいぶん印象が変わる場合もあるので、可能な限り連続変数は連続変数のまま扱いたい。しかし、度数分布表やクロス表ではそういう訳にはいかない。コンピュータの性能が向上し、簡単にきれいなグラフが作れる今日、グラフを使わない手はない。
あつかう変数が本物の連続変数の場合、散布図や箱ひげ図は非常に有効であるし、カテゴリカル変数もモザイク・プロットを使ったほうがずっと見やすいので、強くおすすめしたい。
3.1 Scatter Plot with Multiple Lines 個体別の折れ線を加えた散布図
散布図は2つの連続変数の関連を見るためにとても役立つ。しかし、プール・データで散布図を作る場合、固体内の変化が見えない。そこで、同一個体の値を時間順に線で結んでやれば固体内の変化も、2変数間の関連も見えて、有益である。例えば以下のように。
plot(Y ~ t)
# 個体別に折れ線を散布図に付け加えるための関数
m.lines <- function(x, y, group){
id <- unique(group)
dat <- na.omit(data.frame(x, y))
for(i in 1:length(id)){ #
lines(dat[group==id[i], 1], dat[group==id[i], 2]) #
}
}
m.lines(t, Y, i)
図3 Y の時系列変化
上のように、時点と注目している変数の散布図を作り、個体別に折れ線を加えたものを 多重時系列プロット (multiple time series plot) と Free (2004) は呼んでいる。このデータの場合、特に明確な時点と Y の相関は見られない。このようなデータは時点による擬似相関を気にする必要がないので、固定効果モデルやランダム効果モデルを推定する上では好都合である。
以下は X と Y の折れ線付き散布図である。
plot(Y ~ X)
for(j in 1:4){ # 上の図3と同じようにスクリプトを書いてもいいがループを回したほうが簡、
lines(X[i==j], Y[i==j]) # j に 1～4を代入して4本の折れ線を個体別に描く
}
図4 X と Y の個体別の折れ線付き散布図
全体レベルでの相関はゆるやかだが、折れ線の傾きを見ると、明らかに右肩上がりなので、X と Y の間に何らかの因果関係が期待できそうである。
3.2 Scatter Plot of Within Deviation 固体内偏差の散布図
時間によって変化しない観察されない異質性 (以下では Time-Constant Unobserved Heterogeneity: TCUH と略称) を統制するのが、固定効果モデルやランダム効果モデルの目的であるから、散布図でも、TCUH を統制した上での関連を見たい。TCUH は固体内平均\(\bar X_i\) によって表されるので、これを \(X_{it}\) の個々の値から引いてやればよい。同じ計算を \(Y_{it}\) に関しても行い、\(X_{it} - \bar X_i\) と \(Y_{it} - \bar Y_i\)の散布図を作れば、2変数の TCUH を統制した上での関連を図示できる。 \(X_{it} - \bar X_i, \; Y_{it} - \bar Y_i\) を固体内偏差と呼んでおく。
固体内偏差の計算も何通りかやり方があろうが、plm パッケージの Within() 関数が簡単だと思う。Within() を使うためには、まずデータをデータフレームにまとめて、データフレームを pdata.frame という形式（R 的には「形式」ではなく class というべきだあろう）に変換する。pdata.frame は通常のデータフレームに、どの変数とどの変数が個体 ID と時点の変数かを示す情報（これも R 的には「情報」ではなく属性 (attribute) というべき）を付け加えたもので、Within 関数はその情報を参照して、固体内偏差を計算する。
library(plm)
d0 <- data.frame(i, t, X, Y)
d1 <- pdata.frame(d0) # 特に指定しなければ1列目の変数が個体ID、2列目が時点の変数とみなされる
tapply(d1$X, d1$i, mean) # 個人内平均の確認## 1 2 3 4
## 0.700000 4.800000 4.066667 5.866667cbind(d1$X, # cbind は複数のベクトルを横に、くっつけて行列を返す
Within(d1$X) # Within() が個人内偏差を返していることを確認せよ
) ## [,1] [,2]
## 1-1 -1.8 -2.50000000
## 1-2 1.7 1.00000000
## 1-3 2.2 1.50000000
## 2-1 1.6 -3.20000000
## 2-2 8.4 3.60000000
## 2-3 4.4 -0.40000000
## 3-1 3.5 -0.56666667
## 3-2 7.7 3.63333333
## 3-3 1.0 -3.06666667
## 4-1 3.6 -2.26666667
## 4-2 5.9 0.03333333
## 4-3 8.1 2.23333333plot(Within(d1$X), Within(d1$Y))
図5 固体内偏差の散布図
上の固体内偏差のデータを使って、OLS (Ordinary Least Square) で回帰直線の傾きを推定したのが、固定効果モデルである。重要なのは、この固体内偏差のプロットで、線形の関係が見られるかどうかである。曲線的な関係や閾値の存在が示唆されるならば、それに合わせたモデルを考える必要がある。
3.3 Plotting Discrete or Large N Data 個体数が多い場合や離散変数を扱う場合
散布図は連続変数を扱う場合は有効な探索的分析法であるが、個体数が多いと線や点が重なりあいすぎて分布の特徴が認識できないことが多い。同様の問題は、離散型の従属変数を連続変数とみなして分析する場合にも起きることがある。そのような場合は、時点と従属変数のモザイクプロット、時点数が少なければ、移動表を作って移動表のモザイクプロットを作ればよい。
以下の架空のデータを使って、モザイクプロットを使った探索的分析をやってみよう。
head(d2)## id time sex X Y
## 1 1 1 male 2 3
## 2 1 2 male 1 3
## 3 1 3 male 1 3
## 4 2 1 female 1 2
## 5 2 2 female 2 3
## 6 2 3 female 2 3tab.id <- xtabs(~id, d2)
tab.id[1:20]## id
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3length( tab.id )## [1] 1000xtabs(~time, d2)## time
## 1 2 3
## 1000 1000 1000summary(d2[, 3:4])## sex X
## male :1500 Min. :1.000
## female:1500 1st Qu.:1.000
## Median :2.000
## Mean :2.294
## 3rd Qu.:3.000
## Max. :4.000icc1(d2$X, d2$id)## [1] 0.4707172icc1(d2$Y, d2$id)## [1] 0.368257plot(factor(Y) ~ factor(time), data=d2) # factor で因子に変換してやると離散変数とみなされる
ちなみに散布図にすると、図7のようになる。
plot(jitter(Y) ~ jitter(X), data=d2)
図7 jitter() を使った離散変数の散布図
Y の分布に多少の変化はあるものの大きな変化ではないことがわかる。
移動表とは、ワイド形式のデータの t 時点の Y と t + 1 時点の Y のクロス表のことである。そのため、まずはロング形式のデータをワイド形式に変換する方法を紹介しよう。
3.3.1 Wide and Long Formats ワイド形式とロング形式のデータの変換
データをワイド形式とロング形式の間で変換するのには、 reshape() 関数を用いる。ロングからワイドに変換する場合、以下のように指定する。
reshape(データフレーム名, v.names = 時変変数名のベクトル, idvar=個体id変数名, timevar=時点変数名, direction=“wide”)
以下は使用例。
d2w <- reshape(d2, v.names = c("X", "Y"), idvar="id",
timevar="time", direction="wide")
head(d2w)## id sex X.1 Y.1 X.2 Y.2 X.3 Y.3
## 1 1 male 2 3 1 3 1 3
## 4 2 female 1 2 2 3 2 3
## 7 3 male 3 3 3 4 1 3
## 10 4 female 4 4 2 5 4 5
## 13 5 male 2 3 3 4 2 3
## 16 6 female 2 3 2 4 2 2nrow(d2w)## [1] 1000summary(d2w)## id sex X.1 Y.1 X.2 Y.2 X.3
## Min. : 1.0 male :500 Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.00
## 1st Qu.: 250.8 female:500 1st Qu.:2.000 1st Qu.:3.000 1st Qu.:2.000 1st Qu.:3.000 1st Qu.:1.00
## Median : 500.5 Median :2.000 Median :3.000 Median :2.000 Median :3.000 Median :2.00
## Mean : 500.5 Mean :2.299 Mean :3.351 Mean :2.292 Mean :3.411 Mean :2.29
## 3rd Qu.: 750.2 3rd Qu.:3.000 3rd Qu.:4.000 3rd Qu.:3.000 3rd Qu.:4.000 3rd Qu.:3.00
## Max. :1000.0 Max. :4.000 Max. :5.000 Max. :4.000 Max. :5.000 Max. :4.00
## Y.3
## Min. :1.000
## 1st Qu.:3.000
## Median :3.000
## Mean :3.409
## 3rd Qu.:4.000
## Max. :5.000逆にワイド形式のデータをロング形式に変換する場合は以下のような書式になる。
reshape(データフレーム名, varying=時変変数名のリストまたは行列, direction=“long”)
以下の例を見よ。
temp # 6人 X 4時点のワイド形式のサンプルデータ## id sex x1 x2 x3 x4 y1 y2 y3 y4
## 1 1 m 13 19 25 31 37 43 49 55
## 2 2 f 14 20 26 32 38 44 50 56
## 3 3 m 15 21 27 33 39 45 51 57
## 4 4 f 16 22 28 34 40 46 52 58
## 5 5 m 17 23 29 35 41 47 53 59
## 6 6 f 18 24 30 36 42 48 54 60names.v <- t(matrix(names(temp)[3:10], 4)) # t() は行列を転置する。
names.v # 時変変数の対応関係を示す行列。1行がロング形式ではひとつの変数になる## [,1] [,2] [,3] [,4]
## [1,] "x1" "x2" "x3" "x4"
## [2,] "y1" "y2" "y3" "y4"temp.long <- reshape(temp, varying=names.v, direction="long")
head(temp.long, 15) # 時点変数でまずソートされたデータが得られる## id sex time x1 y1
## 1.1 1 m 1 13 37
## 2.1 2 f 1 14 38
## 3.1 3 m 1 15 39
## 4.1 4 f 1 16 40
## 5.1 5 m 1 17 41
## 6.1 6 f 1 18 42
## 1.2 1 m 2 19 43
## 2.2 2 f 2 20 44
## 3.2 3 m 2 21 45
## 4.2 4 f 2 22 46
## 5.2 5 m 2 23 47
## 6.2 6 f 2 24 48
## 1.3 1 m 3 25 49
## 2.3 2 f 3 26 50
## 3.3 3 m 3 27 51temp.long2 <- temp.long[order(temp.long$id), ] # id で先にソートしたければこうする
head(temp.long2, 15)## id sex time x1 y1
## 1.1 1 m 1 13 37
## 1.2 1 m 2 19 43
## 1.3 1 m 3 25 49
## 1.4 1 m 4 31 55
## 2.1 2 f 1 14 38
## 2.2 2 f 2 20 44
## 2.3 2 f 3 26 50
## 2.4 2 f 4 32 56
## 3.1 3 m 1 15 39
## 3.2 3 m 2 21 45
## 3.3 3 m 3 27 51
## 3.4 3 m 4 33 57
## 4.1 4 f 1 16 40
## 4.2 4 f 2 22 46
## 4.3 4 f 3 28 523.3.2 Ploting Mobility Tables 移動表とそのプロット
それでは、ワイド形式のデータから移動表とモザイクプロットを作ってみよう。
tab.y12 <- xtabs(~ Y.1 + Y.2, data=d2w)
tab.y23 <- xtabs(~ Y.2 + Y.3, data=d2w)
tab.y12## Y.2
## Y.1 1 2 3 4 5
## 1 2 14 12 16 5
## 2 15 37 72 47 14
## 3 16 53 99 87 41
## 4 9 37 91 97 72
## 5 4 21 31 62 46tab.y23## Y.3
## Y.2 1 2 3 4 5
## 1 3 8 13 18 4
## 2 12 34 54 44 18
## 3 13 64 98 94 36
## 4 7 41 94 110 57
## 5 5 18 44 64 47par(mfrow=c(1, 2))
plot(factor(Y.2) ~ factor(Y.1), data=d2w)
plot(factor(Y.3) ~ factor(Y.2), data=d2w)
図8 Yという変数の移動表のモザイクプロット
緩やかな相関があるのがわかる。ちなみに以下が同じ移動表を散布図で作ったもの。モザイクプロットのほうがわかりやすいと私は思うのだが、いかがだろうか。
plot(jitter(Y.2) ~ jitter(Y.1), data=d2w) 
図9 変数 Y の移動表の散布図
3.3.3 Bubble Plot バブルプロット
移動表の各カテゴリの数を点の大きさで表し、線の太さで移動した人数を表すという方法も考えられる（あまり見たことはない）が、スクリプトを書くのが面倒なわりにあまりわかりやすいとはいえない。
# データを用意
tab.Yt <- xtabs(~ Y + time, d2)
tab.Yt## time
## Y 1 2 3
## 1 49 46 40
## 2 185 162 165
## 3 296 305 303
## 4 306 309 330
## 5 164 178 162tab.Yt.long <- as.data.frame.table(tab.Yt) # クロス表のかたちをロング形式みたいに変える
tab.Yt.long## Y time Freq
## 1 1 1 49
## 2 2 1 185
## 3 3 1 296
## 4 4 1 306
## 5 5 1 164
## 6 1 2 46
## 7 2 2 162
## 8 3 2 305
## 9 4 2 309
## 10 5 2 178
## 11 1 3 40
## 12 2 3 165
## 13 3 3 303
## 14 4 3 330
## 15 5 3 162# as,numeric() で因子を数値に変換してプロット
plot(as.numeric(Y) ~ as.numeric(time), data=tab.Yt.long, cex= Freq / 50)
tab.y12.long <- as.data.frame.table(tab.y12)
tab.y23.long <- as.data.frame.table(tab.y23)
for(i in 1: nrow(tab.y12.long) ){
lines(1:2, tab.y12.long[i, 1:2], lwd = tab.y12.long$Freq[i] /10)
lines(2:3, tab.y23.long[i, 1:2], lwd = tab.y23.long$Freq[i] /10)
}
3.4 Plotting Discrete and Continuous Variables 離散変数と連続変数の関連
これについてはもう書かないが、これまでの知識を適宜応用すればできるはずである。変数の分布などデータの特徴によってどんな分析法が適切かは異なる。データの分布をよく見ながらわかりやすい図を作って欲しい。