Gerber and Green(2012) Field Experiment, Ch.4

勉強会のメモの続き．実験で共変量をどう使うかについての章である．

4.1. Using Covariates to Rescale Outcomes

共変量とはPotential Outcome(以下PO)を予測すると考えられる観察された変数であり，処置の割り当てに影響を受けないと仮定される．つまり共変量とは基本的に処置の前に決定しているpre-treatment variableである．以下の議論では共変量がこの仮定を満たすとする．割り当て確率が等しいとき共変量について以下が成り立つ．
{ \displaystyle E[X_i]=E[X_i|D_i=1]=E[X_i|D_i=0] }
{ \displaystyle X_i }がプレテストスコアで{ \displaystyle Y_i }が処置後のポストテストスコアだとするとDID推定量は，
{ \displaystyle E(\widehat{ATE})=E[Y_i-X_i|D_i=1]-E[Y_i-X_i|D_i=0] \\ =E[Y_i|D_i=1]-E[X_i|D_i=1]-E[Y_i|D_i=0]+E[X_i|D_i=0] \\ =E[Y_i(1)]-E[Y_i(0)] }
となる．ランダム割り当てがなされていればDIM(Differences-in-means)推定量{ \displaystyle E[Y_i(1)]-E[Y_i(0)] }とDID推定量はともに不偏推定量となるが，小標本の場合に両者は異なる推定値となる（証明略）．結論を述べると，DID推定量の標本分散がDIM推定量の標本分散より小さくなるのは，共変量とPOの関連が強い時である．フォーマルには，
{ \displaystyle \frac{Cov(Y_i(0), X_i)}{Var(X_i)}+\frac{Cov(Y_i(1), X)}{Var(X_i)} > 1 }
と書ける．つまり実際にフィールド実験を行う際には，a)サンプルサイズを大きくする，b)共変量の情報を集める，が予算制約下でのトレードオフとなっている．共変量を集めることによって生じうる深刻な問題としては，共変量を集めること自体が対象者の処置への反応の仕方を変えてしまい，それによってアウトカムも変化してしまうことであり，これでは除外制約等を満たさない．ここでは共変量を利用することでうまくいったケースについてMuralidharan and Sundararaman(JPE 2011)の結果を再現してみよう．アウトカムは学校の平均テストスコア，処置は教員へのボーナス，共変量は事前の平均テストスコアである．POがわかっているとして，Y_i(0)をX_iに回帰した時のR^2は0.87となり関連が強いことがわかる．その上で，10万回の実験シミュレーションを行ってみると以下のような結果となる．
f:id:analyticalsociology:20160325150723p:plain
DIDEの方がATEの標準誤差が小さいことがわかるだろう．DIMとDIDの平均ATEともに4.0で同じだが，標準誤差は4.77と1.53でありDIDの方が約68%(1-(1.53/4.77))小さい．ちなみにDIMでDIDと同程度に標準誤差を小さくするためにはサンプルサイズを40から392に大きくする必要がある(1/(1-0.68)^2)．このRスクリプトは以下．

# clear workspace
rm(list = ls())

# Load in library to read stata files
library(foreign)

# Load in Data from website
teacherout <- read.dta("http://hdl.handle.net/10079/wwpzgz8")

# Or load in Data as a .csv file from website
#teacherout <- read.csv(file="http://hdl.handle.net/10079/1ns1rzk",head=TRUE,sep=",")

# Attach Data
attach(teacherout)

## Reproduce Figure 4.1
# Create space to draw Figure 4.1
par(family="Gill Sans MT",font.main=1)
layout(matrix(c(1,2),2,1,byrow=TRUE))

# Graph the top part of Figure 4.1
hist(diffinmean,xlim=c(-10,20),freq=FALSE,ylim=c(0,.25),main="Sampling Distributions",xlab="Difference-in-Means")
lines(density(diffinmean))

# Graph the bottom part of Figure 4.1
hist(diffinchangemeans,xlim=c(-10,20),freq=FALSE,ylim=c(0,.25),main=NULL,xlab="Difference-in-Differences")
lines(density(diffinchangemeans))

# Detach data so can reload data with identical variable names in examples below.
detach(teacherout)

4.2. Adjusting Covariates Using Regression

POと回帰の関係を整理してみよう．
{ \displaystyle Y_i=Y_i(0)(1-d_i)+Y_i(1)d_i=Y_i(0)+(Y_i(1)-Y_i(0)d_i) \\ =\mu_{Y(0)}+[\mu_{Y(1)}-\mu_{Y(0)}]d_i+Y_i(0)-\mu_{Y(0)}+[(Y_i(1)-\mu_{Y(1)})-(Y_i(0)-\mu_{Y(0)})]d_i \\ =a+bd_i+u_i }
以上のようにYをdに回帰すればDIMによるATEと同じものを得る．共変量によってアウトカムを調整した場合も同様で，
{ \displaystyle Y_i-X_i=Y_i(0)(1-d_i)+Y_i(1)d_i-X_i\\ =a+bd_i+u_i-X_i }
となる．回帰式の右側に共変量を1つ以上含む場合にも，
{ \displaystyle Y_i=Y_i(0)(1-d_i)+Y_i(1)d_i \\ =a+bd_i+cX_i+(u_i-cX_i) }
cが1に近い時に上式は{ \displaystyle =a+bd_i+u_i-X_i }に近い値となる．
以上は単なる整理だが，実験ではアウトカムと共変量は相関であり因果ではない．ランダム割り当てを伴う実験で共変量を用いるのは，誤差項から共変量を引いてやることでアウトカムの説明されない部分を小さくし標準誤差を小さくするためである．但し問題としては，小標本の場合に回帰調整ではバイアスを持つ可能性があることに注意が必要である．

4.3. Covariate Imbalance and the Detection of Administrative Errors

さて，実験をやったことある人ならば経験済みだろうが，小標本である場合にはランダム割り当てでも処置群と対照群でバランスが取れない＝共変量と割り当てに相関が生じることがある．まずやるべきなのは2群のバランスをチェックすることだ．バランスしていない場合には，割り当てを共変量に回帰しF統計量を計算し，p-valueを求めるためにrandomization inferenceをする．偶然でないインバランスであればコントロールしてやれば良い．この点については以下のように著者らがそんな心配するなと回答している．

Sometimes researchers worry that observed imbalance is symptomatic of a broader problem of imbalance among other unobserved factors that affect outcomes. So long as imbalance is solely due to random chance(as opposed to administrative error) and so long as we control for the covariate that is imbalanced, there is no reason to expect imbalance on other covariates or on unmeasured causes of the outcome variable.

4.4. Blocked Randomization and Covariate Adjustment

偶然のインバランスが生じた時に共変量をどう用いるかによってATEの推定値は異なってくることを確認したが，事前にブロックランダマイゼーションという方法もある．実験家の格言としては，

Block what you can, and randomize what you cannot

というものがあるそうだ．準備期間が十分にあるのであれば，複数の共変量についてブロックすることができるが，そうでない場合には先行研究を参考にPOを強く予測する変数を選択すべきである．Block randomizationの場合，Yをdとxに回帰するかYをdに回帰するかに関わらず同じATEを得る．共変量調整を行ったrandomizationの場合，偶然にインバランスして共変量と割り当てが相関する可能性があり，不偏推定量を得てもATEの推定値の標準誤差が大きくなってしまうことがある．これをcollinearity penaltyと呼ぶ．但し大標本では2つの標本分布は近いものとなる．改めてComplete randomizationとblock randomizationの違いをみてみよう．
f:id:analyticalsociology:20160326121651p:plain
POを強く予測する共変量でブロックしたATEが最もばらつきが小さいことがわかる．

4.5. Analysis of Block Randomized Experiments with Treatment Probabilities That Vary by Block

これまではブロックごとの割り当て確率が等しいと仮定してきたが，ブロックごとに割り当て確率が異なる場合は全てをプールしてATEを推定するとバイアスが生じる．この場合，各ブロックごとにATEを推定してそこに当該ブロックにおける割り当て確率をかけてやり，各ブロックから得るこの値を足し合わせてやれば良い．以上を簡単にまとめると，Block randomizationの手順は概ね以下となろう．

  1. どの共変量がPOを強く予測するか確認する→POを共変量に回帰してR^2を確認
  2. 共変量についてブロックランダマイゼーション→RのblockToolsパッケージが便利→処置群と対照群のサンプルサイズを決定→各ブロックで2群のバランスがとれているかチェック
  3. ブロックがATEのばらつきを小さくするのにどの程度うまくいったかの確認→上の図でみたようにcomplete randomizationとblock randomizationの標本分布を比較
  4. block randomizationとblockに用いた共変量で調整した回帰のATE標本分布を比較
  5. ばらつきの小さいものを選ぶ！

Gerber and Green(2012) Field Experiment, Ch.3

メモの続き．第3章はSampling Distributions, Statistical Inference, and Hypothesis Testingである．

3.1. Sampling Distribution

The term sampling distribution refers to the collection of estimates that could have been generated by every possible random assignment.

Rosenbaum(1984)ではrandomized distributionと呼んでいるが，標本分布のことを意味している．Ch.2の例では7つの村のうち2つをtreatmentとして扱ったが，2つのtreatmentの選び方は21通りあり，それぞれの組み合わせに応じてATEは変化する（つまり実験の結果はバイアスはないが必ずしも正確とは言えない）．しかし，ランダム割り当て，除外制約，SUTVAを満たしていれば，推定されたATEの平均は真のATEに一致する．

3.2. The Standard Error as a Measure of Uncertainty

こうしたsampling variabilityの指標となるのが標準誤差である．標準誤差はばらつきの指標なので小さいとありがたい．ではどうしたら標準誤差を小さくできるのか？換言すれば，どうしたらより正確なATEを得るような実験をデザインできるのだろうか．Nをサンプルサイズ，mをtreatment unitのサイズとすると，推定されたATEの標準誤差は，
{ \displaystyle SE(\widehat{ATE})=\sqrt{\frac{1}{N-1}\Bigl\{\frac{mVar(Y_i(0))}{N-m}+\frac{(N-m)Var(Y_i(1))}{m}+2Cov(Y_i(0), Y_i(1))\Bigr\}} }
となることから，標準誤差を小さくするには以下４点がわかる．

  1. Nを大きくすると標準誤差は小さくなる
  2. { \displaystyle Y_i(1), Y_i(0)}の分散を小さくすると標準誤差は小さくなる
  3. { \displaystyle Cov(Y_i(1), Y_i(0))}を小さくすると標準誤差は小さくなる
  4. { \displaystyle m \approx 2/N}の時に標準誤差は最も小さくなる

3.3. Estimating Sampling Variability

真の標準誤差はわからないので以下の式で推定しなければならない．
{ \displaystyle \widehat{SE}=\sqrt{\frac{\widehat{Var}(Y_i(0))}{N-m}+\frac{\widehat{Var}(Y_i(1))}{m}} }

3.4. Hypothesis Testing

代表的な仮説検定は以下である．
Sharp Null Hypothesis of No Effect
全てのユニットについて因果効果がゼロ：{ \displaystyle Y_i(1)=Y_i(0)}
Null Hypothesis of No Average Effect
ATEがゼロ：\mu_{Y(1)}=\mu_{Y(0)}
こうした仮説検定をrandomization inferenceと呼び，以下で説明される．

Randomization Inference: The Sampling distribution of the test statistic under the null hypothesis is computed by simulating all possible random assignments. When the number of random assignments is too large to simulate, the sampling distribution may be approximated by a large sample of possible assignments. p-values are calculated by comparing the observed test statistic to the distribution of test statistics under the null hypothesis.

t値による検定もあるが，推定値の分布が歪んでいるときにも応用できるのはp-valueであり本書では後者を使っている．

3.5. Confidence Intervals

ATEがどのくらいの範囲にあるのかを知りたい時に区間推定を参照することがある．この場合には，推定されたATEを用いてPotential Outcomeに代入して信頼区間を推定する方法がある．Clingingsmith et al.(QJE 2009)ではPOの表を完成させ100000回のランダム割り当てをするシミュレーションをして信頼区間を求めている．

3.6. Sampling Distributions for Experiments that Use Block or Cluster Random Assignment

Block Random Assignmentとは対象者をサブグループに分けた上でランダム割り当てすることである．ランダム割り当てをしても処置群と対照群でバランスがとれていないことが起こる．例えば20人の対象者がいて，男女が半々とし，ランダム割り当てをすると，処置群に男性ばっかりもしくは女性ばっかりという可能性もなくはない．この場合，あらかじめ男女でサブグループに分けてランダム割り当てをすれば処置群と対照群で男女の比が同じになる．さらにBlock Random AssignmentはPOの分散が小さくなるので標準誤差が小さくなるという利点もあり，著者らは全てをプールした状態でのランダム割り当て＝Complete Random Assignmentと比べてBlock Random Assignmentのほうが望ましいとしている．ちなみに全体のATEはブロックjのサイズを考慮してやればよいので以下となる．
{ \displaystyle ATE=\sum_{j=1}^J\frac{N_j}{N}ATE_j }
Olsen(JPE 2007)の例ではBlock Random AssignmentとComplete Random AssignmentのATEの標準誤差の違いがはっきりと現れており分りやすい．以下のようにブロックしたATEの推定値のばらつきのほうが明らかに小さい．
f:id:analyticalsociology:20160324183512p:plain
ちなみにRスクリプトを以下．

rm(list=ls(all=TRUE))
library(ri)
set.seed(1234567)

Y0 <- c(0,1,2,4,4,6,6,9,14,15,16,16,17,18)
Y1 <- c(0,0,1,2,0,0,2,3,12,9,8,15,5,17)

Z <- c(1,1,0,0,0,0,0,0,0,0,0,0,1,1)

# generate all permutations of Z under _complete_ random assignment
# note that default is to do every possible permutation if less than 10,000 permutations

compperms <- genperms(Z)
numperms <- ncol(compperms)

# create empty vector
compmeans <- rep(NA,numperms)

# loop to create average treatment effect estimates for each randomization
for (i in 1:numperms) compmeans[i] <- mean(Y1[compperms[,i]==1]) - mean(Y0[compperms[,i]==0])

# randomize within blocks
block <- c(1,1,1,1,1,1,1,1,2,2,2,2,2,2)

# generate all permutations of Z under block random assignment

blockperms <- genperms(Z,block)
numperms <- ncol(blockperms)

# create empty vector
blockmeans <- rep(NA,numperms)

# loop to create average treatment effect estimates for each randomization
for (i in 1:numperms) blockmeans[i] <- weighted.mean(Y1[blockperms[,i]==1],c(8/2,8/2,6/2,6/2)) - weighted.mean(Y0[blockperms[,i]==0],c(8/6,8/6,8/6,8/6,8/6,8/6,6/4,6/4,6/4,6/4))

save(compmeans,blockmeans,file="figure3.1.Rdata")

# Draw histograms

par(mfrow=c(2,1))
hist(compmeans,main="Sampling Distribution under Complete Randomization",xlim=c(-15,10),xlab="ATE Estimates",freq=FALSE,ylim=c(0,.3))
hist(blockmeans,main="Sampling Distribution under Blocked Randomization",xlim=c(-15,10),xlab="ATE Estimates",freq=FALSE,ylim=c(0,.3))

# calculate the proportion of esitmates that are above zero

length(compmeans[compmeans > 0])
length(compmeans[compmeans > 0])/length(compmeans)

length(blockmeans[blockmeans > 0])
length(blockmeans[blockmeans > 0])/length(blockmeans)

ブロックしたい共変量が明確な場合は以上の手順で実行可能だが，実際には共変量が多い場合などは何でブロックするかは悩ましいことがある．その場合にはどうしたら良いのだろうか．この点について，ルービンの講義を受けた時に彼は以下の会話をスライドで紹介していた．

Rubin: What if, in a randomized experiment, the chosen randomized allocation exhibited substantial imbalance on a prognostically important baseline covariate?
Cochran: Why didn't you block on that variable?
Rubin: Well, there were many baseline covariates, and the
correct blocking wasn't obvious; and I was lazy at that time.
Cochran: This is a question that I once asked Fisher, and his reply was unequivocal:
Fisher (recreated via Cochran): Of course, if the experiment had not been started, I would rerandomize.

実際に多くの場合にはフィールド実験の前に共変量がわかっていることはあまりないと思うので，実験前にバランスチェックするのが困難だと思うが，実験環境によってはルービンやフィッシャーがいうようにランダム割り当てをもう一回すれば良いのだろう．

Gerber and Green(2012) Field Experiment, Ch.2

先日からフィールド実験の勉強会に参加しているので，毎回の内容について簡単にメモをとることにした．テキストは以下．

Field Experiments: Design, Analysis, and Interpretation

Field Experiments: Design, Analysis, and Interpretation

2章は基本的なノテーションの確認なので，一般的な因果推論のテキストの冒頭と類似している．

2.1. Potential Outcome

　本章で取り上げられる論文がChattopadhyay and Duflo(ECTA 2004)であり，この論文は，村長の性別によって予算配分に違いがあるかを検討している．女性が村長になった場合には，男性に比べてより多くを水を綺麗にするために設備投資等に費やすと想定されるそうだ． 「処置=女性村長選出」とすればPotential Outcome(以下PO)は以下で定義される．
{ \displaystyle Y_i(1)}: 村iが女性村長の時のアウトカム
{ \displaystyle Y_i(0)}: 村iが男性村長の時のアウトカム

2.2 Average Treatment Effects

　各村における平均処置効果は，
{ \displaystyle \tau_i=Y_i(1)-Y_i(0)}
となり，各村の観察されるアウトカムは，
{ \displaystyle Y_i=d_iY_i(1)+(1-d_i)Y_i(0)}
と書ける．そして全体の平均処置効果(平均因果効果=Average Treatment Effects)は，
{ \displaystyle ATE\equiv\frac{1}{N}\sum_{i=1}^N\tau_i}

2.3. Random Sampling and Expectations

　POの期待値はATEに等しい．
{ \displaystyle E[Y_i(1)-Y_i(0)]=E[Y_i(1)]-E[Y_i(0)]=\frac{1}{N}\sum_{i=1}^N[Y_i(1)-Y_i(0)] }

2.4. Random Assignment and Unbiased Inference

　ランダムに割り当てされた処置{ \displaystyle D_i}は，実際に観察されるアウトカム{ \displaystyle Y_i}に影響を与えるが，POや共変量とは独立．
{ \displaystyle Y_i(0), Y_i(1), X \perp D_i}
不偏推定量とは{ \displaystyle E(\hat\theta)=\theta}であるが，ランダム割り当てがなされている時のATEは不偏推定量となる．

2.5. The Mechanics of Random Assignment

　単純なランダム割り当てはコイントスやサイコロを振って行われるが，サンプルサイズが小さい場合には実際には特定の値が出やすくなったりなかなかランダムにならない．著者らが完全ランダム割り当てとしてあげている方法はふたつ．

  1. 全対象者=N人の中から1人を選ぶ→N-1人の中から1人を選ぶ→処置群に入る数=m回繰り返す
  2. N人に番号をランダムに割り当て，順番にm番目の人を選ぶ

2.6. The Threat of Selection Bias When Random Assignment is not Used

　セレクションバイアスとはPOと処置が独立でないことなので，
{ \displaystyle E[Y_i(1)|D_i=1]=E[Y_i(0)|D_i=0]=E[Y_i(1)-Y_i(0)|D_i=1]+E[Y_i(1)|D_i=1]-E[Y_i(0)|D_i=0] }
の中で{ \displaystyle E[Y_i(1)|D_i=1]-E[Y_i(0)|D_i=0] }
がセレクションバイアスになる．

2.7. Two Core Assumption about Potential Outcomes

　「POは本人が処置を受けるか受けないかだけに依存する」という仮定がある．他人の処置に影響を受けないというのがSUTVAで，本人の処置「だけ」に依存するのが除外制約である．IVの文脈と同じである．

Takaku(SSM 2016) 乳幼児医療費助成は子どもの健康を改善したか

　子どもが医療機関にかかる際に自己負担額が軽減される，乳幼児医療費助成制度が子どもの健康を改善したのか否かを検証した論文．乳幼児期の健康は，その後の健康のみならず様々な面で影響があることが知られており重要なテーマである．

Takaku, R. (2016). "Effects of reduced cost-sharing on children’s health: Evidence from Japan." Social Science & Medicine, 151: 46-55.

　日本の医療保険制度では，就学前児童の外来診療は2割負担となっているが，ほとんどの自治体が独自に2割負担をさらに軽減するための乳幼児医療費助成制度を拡充してきた． 著者が着目するのはこの点である．すなわち，乳幼児医療費助成の在り方が自治体によって異なるのであれば，それを自然実験とみなして識別戦略を練ることが可能になるということだ．社会学ではこうした分析デザインはあまりないが，経済学や近年のポリサイではよくみられる方法である．この論文の面白いところは，自治体によって異なる乳幼児医療費助成の拡充について，独自の調査で明らかにし，その結果を分析に用いているところである．本論文で扱われているのは1995~2010年であるが，この期間に乳幼児医療費助成が自治体ごとにどう変化したのかが分かる公式な資料はないそうだ．そこで全国市区町村に対して郵送調査を行い，乳幼児医療費助成の拡充過程をまず明らかにする．回収率は55%だったそうだが，15歳以下人口でウェイトをかけたところ回収率は75%まで上昇している．郵送調査によって得られた市区町村の医療費助成対象年齢から，都道府県ごとに子どもの年齢及び調査年別にみた「医療費助成の対象となる確率」を下式で算出している．
{ \displaystyle Prob_{npt} = \sum_{m=1}^{N} W_{mt}*Elig_{nmt} }
　t時点における都道府県pに居住するn歳のiさんが乳幼児医療費助成の対象となる確率{ \displaystyle Prob_{npt} }は「都道府県pにおける市区町村mに居住するn歳のiさんが乳幼児医療費助成の対象年齢であるか否か」に15歳以下人口のウェイトをかけた{ \displaystyle \sum_{m=1}^{N} W_{mt}*Elig_{nmt} }に等しいという設定である．そして，「国民生活基礎調査」のおける子どもの健康に関する指標と「医療費助成の対象となる確率」をマッチさせ，乳幼児医療費助成を拡充した都道府県で子どもの健康指標が改善したかを検証している．したがって，最終的な分析モデルは，アウトカム=子どもの健康指標を{ \displaystyle H_{it} }，統制変数を{ \displaystyle X_{it} }，誤差項を{ \displaystyle \upsilon_{it} }とすると，
{ \displaystyle H_{it}=G(Prob_{npt}, X_{it}, \upsilon_{it}) }
となる．知りたいのは{ \displaystyle Prob_{npt} }の因果効果であるため，幾つかのSpecification Testを行い，さらに測定誤差の問題に対処するために回収率の低かった都道府県を除いた分析を行っているが，結果に大きな違いはみられなかった．アウトカムである子どもの健康指標については，有訴(熱，だるさ，咳，頭痛，ぜいぜい，⻭痛，鼻づまり，便秘，下痢，胃痛，発疹，切り傷)，入院しているか否か，健康上の問題による日常生活の困難の有無，主観的健康観，病院外来の有無を設定している．また就学児と未就学児に分けて分析を行っている．ちなみにこ健康指標については本人が回答しているのではなく親が代わりに回答している．

　分析結果は以下の通りである．

  • 有訴確率について，「医療費助成の対象となる確率」は未就学児においてのみ有意に負の効果(=有訴確率を下げる)があるが，就学児については有意な効果はない．有訴の内容について細かくみると，未就学児フルサンプルでは熱，せき，歯痛，便秘が全体の負の効果を牽引している．就学児については健康上の問題による日常生活の困難の有無，主観的健康観も検討したが効果はなかった．
  • 入院確率について，「医療費助成の対象となる確率」は未就学児・就学児ともに有意な効果はない．
  • 病院外来にかかっている確率について，「医療費助成の対象となる確率」は未就学児においてのみ有意に正の効果があるが，回収率の低かった都道府県を除くと有意でなくなる．症状について細かくみると，主に重篤でないと思われる「せき」が有意に正の効果をもっており，先の入院確率において有意な効果がないという分析結果と整合的だとされている．

　以上から，乳幼児医療費助成の拡充が子どもの健康改善に与える効果は限定的だったと結論づけている．分析結果をうけて，著者は幾つかの論点を挙げている．まず，今回の分析では乳幼児医療費助成が子どもの健康改善に与える効果は限定的だったが，医療の目的はそれだけではないので，効果が限定的だからという理由のみで制度は否定されないということ．また，今回の分析では乳幼児医療費助成の短効果に着目したが，長期的にはどうなるかわからないということ．例えばCurrie et al.(2008)では短期的な効果はないが長期的には効果があることが示されているそうだ．さらに，今回の分析とは直接関連しないが，日本では乳幼児医療費助成は拡充される一方で，妊娠中の女性に対する医療費助成を実施している自治体は少ないそうである．Currie and Gruber(1996)では妊娠中の女性に対する医療費助成が生まれてくる子どもの健康状態を改善したことを報告しているため，今後の乳幼児医療費助成の在り方については様々な選択肢があるだろう．実際に，子どもの医療については厚労省「子どもの医療制度の在り方等に関する検討会」でも議論されている．この結果を知ってショックを受ける関係者は多いかもしれないが，分析デザインが練られており，またテーマも重要で政策的インプリケーションも豊富なので，各方面で貢献のある分析だと思った．

Polavieja(ASR 2015) 伝統主義は女性労働参加を抑制するか：Synthetic IVによる推定

価値観が行動を規定することは社会科学でしばしば指摘されているが，本稿では伝統主義(traditionalism)が女性の就業/非就業に影響を与えているのかを分析している．

Polavieja, J. G. 2015. "Capturing Culture: A New Method to Estimate Exogenous Cultural Effects Using Migrant Populations." American Sociological Review 80(1): 166–191.

推定したいのは伝統主義T(0~5の値をとる)が女性の就労に与える効果である．ここでIVを使う．まず分析対象を移民に絞り，彼らの伝統主義を欠損値にする(実際には観察されている)．次に，その移民経験者の出身国の人たちの伝統主義を共変量に回帰させてimputationして得られるsynthetic trait(T')を得る．このT'をTの操作変数にするということだ．操作変数の仮定としては，独立性，関連性，除外制約があるが，著者はこれを満たしていると述べている．だが，明らかに独立性，除外制約は満たしていないだろうというのが読んでみての感想である．

独立性について，著者のディフェンスは「実際に移民が現在住んでいないところの人たち，つまり社会環境が違う人たちから合成されたT'なので現在の移民が住んでいる環境からは独立である」というものである．おかしいのは，第1に，IVが独立でなければならないのは著者がいう社会環境だけでなく，共変量を含めたあらゆる欠落変数に対してだ．したがって著者はそうした変数との関連がないのかをチェックすべきだが，その点については何も言及されていない．第2に，住んでいる国が違うからといって社会環境が全く異なるというわけではない．確かに国が違えば文化は異なるが，当然共通する部分だってある．共通する部分があれば著者のいう"By construction, synthetic traits satisfy the exogeneity condition"というのは成り立たない．著者曰く，このSynthetic IVは疫学で使われるらしいのだが，まずい気がする．

除外制約について，著者はかなり苦し紛れに書いているのが分かるが，結果として除外制約が成立しているとするのは無理がある．著者はまず以下のように述べる．

Unlike the exogeneity condition, the exclusion restriction is not automatically ensured by the use of migrating populations. Suppose, for instance, that we used donors' age and schooling as the sole predictors of traditional values in the imputation equation. In this case, we would impute synthetic values of traditionalism for migrants by imposing on them the same age-schooling-traditionalism covariance matrix as observed for donors. The problem in this example is that synthetic values would be computed as a linear product of two variables that do not satisfy exclusion.

このままでは除外制約を満たさないので，以下の対処をしている．

To ensure the exclusion restriction is met, the imputation regression should include at least one regressor that is (arguably) orthogonal to the error term in the structural equation of interest (i.e., one that has no effect on Y other than through T). I propose to use regressors that measure cultural transmission as a means to build exclusion into the imputation model.

ここで言われている誤差項に直交する説明変数として，移民の出身国における(おおまかな)一世代前の女性労働参加率である．除外制約はIVであるT'がTを通してのみアウトカムに効果があるということであるのに，imputationの回帰式に特定の説明変数を加えたからOKという話とは別である．さらに言えば，一世代前の女性労働参加率がほんとうに直交しているかどうかも怪しい．

ざっと読んだだけでも問題点と思われる箇所が目に付いたが，ASRに通っているし，実はSynthetic IVというのはかなり進んでいたりするということもあるのだろうか(ざっとググった感じはそうでもなさそう)．このあたりは詳しい人に教えて欲しい．アバディーらのSynthetic Controlはかなり使える手法だと思うが，著者が提唱しているSISTER=Synthetic IVというのは可能なのだろうか．分析の結論としては，伝統主義思想は女性労働参加をかなり抑制するようで，教育効果の約2倍ほどあるらしい．

Gelman(2006) マルチレベルモデリングで出来ること出来ないこと

マルチレベルモデリングの有用性と限界が明確に述べてある論文．著者はData Analysis Using Regression and Multilevel/Hierarchical Modelsで有名なAndrew Gelman．

Gelman, A. 2006. “Multilevel (Hierarchical) Modeling: What It Can and Cannot Do” TECHNOMETRICS 48(3): 432-435.

要旨は以下．

Multilevel (hierarchical) modeling is a generalization of linear and generalized linear modeling in which regression coefficients are themselves given a model, whose parameters are also estimated from data. We illustrate the strengths and limitations of multilevel modeling through an example of the prediction of home radon levels in U.S. counties. The multilevel model is highly effective for predictions at both levels of the model, but could easily be misinterpreted for causal inference.

結論から述べれば，マルチレベルモデリングは予測には役立つが因果推論として解釈するには困難があるということだ．ここでは肺がんを引き起こすといわれるラドンを例にとりマルチレベルモデリングを定式している．ラドンは地域によって偏りがあるらしく，本分析の目的はラドン分布の偏りを調べることであり，そうすれば各家庭はなんらかの対処ができるかもしれない．居住地(レベル1)は群(レベル2)にネストされていると考えよう．レベル1の説明変数は居住地(各家庭のこと)の地下で測定がされたか否か(ラドンは地下にあるため地下で測定していたらより高い値がでやすいと予想される)．レベル2の説明変数は群で観測される土壌ウランの値．定式を以下とする．
{ \displaystyle y_{ij}\sim N(\alpha_j+\beta x_{ij}, \sigma_y^2) }
{ \displaystyle \alpha_j\sim N(\gamma_0+\gamma_1 u_{j}, \sigma_\alpha^2) }
上式をミネソタ州データ(919house, 85counties)を用いた階層ベイズ(Hierarchical Bayes Methods)で分析する．

Data Reduction
ここでは上式のマルチレベルモデリングを，すべてプールした{ \displaystyle y=\alpha+\beta x}とプールなしの{ \displaystyle y=\alpha_j+\beta x}と比較している．Fig1.からマルチレベルモデリングが他の2モデルに比べてData Reductionの面で優れていることが一目瞭然である．

Prediction
Gelmanはマルチレベルモデリングが最も役に立つのはこの「予測」であると述べている．どこかの群で新しい家が建てられた時に計測されるラドンの値を予測するのに適しているのだ．変量効果のおかげである．例えば，

We can use cross-validation to formally demonstrate the benefits of multilevel modeling. We perform two cross-validation tests: first removing single data points and checking the pre- diction from the model fit to the rest of the data, then removing single counties and performing the same procedure. For each cross-validation step, we compare complete-pooling, nopooling, and multilevel estimates. Other cross-validation tests for this example were performed by Price et al. (1996).

のようにして色々と試行錯誤できる．

Causal Inference
マルチレベルモデリングは因果推論として使えるのだろうか．つまり，地下測定ダミーの係数を因果効果として解釈できるのだろうか．かなり難しいというのがGelmanの結論である．問題はレベル1の変数とレベル2の誤差項の相関である．上式における尤度と事前密度をかけた場合には，事後密度ではレベル2の誤差項とレベル1の変数は独立であることを仮定している．この独立性はまず成立しないだろう．経済学でマルチレベルモデリングがあまり使われないのは，こういうところにあるのだろう．

In other settings, especially in social science, individual averages used as group-level predictors are often interpreted as “contextual effects.” For example, the presence of more basements in a county would somehow have a radon-lowering effect. This makes no sense here, but it serves as a warning that, with identical data of a social nature (e.g., consider substituting “income” for “radon level” and “ethnic minority” for “basement” in our study), it would be easy to leap to a misleading conclusion and find contextual effects where none necessarily exist.

要するにこれは生態学的誤謬の話であり，何の用心もなしにレベル2の係数を因果効果のように解釈することは危険である(ちなみにレベル2変数にIVを使う方法もあることはある)．以上の点は，マルチレベルモデリングをよく理解している人にとっては「何をいまさら」という点であるが，分析目的に応じてマルチレベルモデリングを採用するか否かを決める必要がある(すべての分析手法について言えることですが)．Gelmanはブログでもこういう分かりやすい記事を書いてくれるから重宝する．

Charles et al. (SS 2015) ケアワークにおける世代間再生産

親がケアワーカーだと子もケアワーカーになりやすいのかを検証した論文．

Charles, M., Ellis, C., & England, P. 2015. "Is There a Caring Class? Intergenerational Transmission of Care Work." Sociological Science 2: 527–547.

親の社会経済状況が子の社会経済状況に大きな影響を与えることは広く知られている．社会学では伝統的に階級や階層の世代間継承が分析されてきたが，特定の職業の世代間継承については蓄積が少ない(Weeden and Grusky 2005; Jonsson et al. 2009)．そこで，著者らはケアワークという職業に着目し，親から子への世代間再生産が生じているのかを検証している．データはGSSで，1977年から2010年のなかの23waveを使用している．GSSを使用するのは，回答者の親の職業や従業上の地位に関する質問項目を含むからである．ケアワークの定義に関しては，先行研究で指摘されているものに則っている(看護や介護系はもちろん教師，保育士，ベビーシッター，チャイルドマインダー，在宅医療従事者等も含む)．現在ケアワークに従事しているかをアウトカムに，着目する独立変数には親(父と母を分けて)がケアワーカーであったかを設定している．基本的な属性や職業威信等が統制されている．ロジスティック回帰分析を行った結果は，同性の親がケアワークであった場合には子もケアワークに従事する確率が高くなっている．但し，女性に限っては異性である父親がケアワークに従事しても自らがケアワークに従事しやすいという結果を得ており，これは男性の場合に異性である母親ケアワークの効果がないことと異なる．世代間継承のメカニズムは様々あるが，著者らが挙げているのは利他傾向のあるケアワークに就く親は子に生活面でも利他心を伝えることがあり，子に継承された利他心は，子が職業選択の際に「利他傾向が求められるだろうケアワーク」を選択しているという図式だ．それに対して，そんな図式は存在せず，ただ単に子が親と全く同じ職業を選択している可能性もある．そこで著者らは親子で全くおなじ職業(親と子が両方とも幼稚園の職員，親と子が両方とも看護師など)の対象者を除いて再分析している．すると先ほど得られた親の効果はすべてなくなっており，このことから，著者らは親から子へのケアワーク継承は，ケアワークという職種全体で生じているというよりも，ケアワークにおける特定の職業を継承するという形で生じているだろうと述べる．なお，補助的な分析によって親がケアワークだと本人が(質問紙上で)利他的になりやすいという傾向はほぼない*1ことも確認している．分析の粗さはあるものの，この枠組みは日本でも応用可能なので非常に興味深いし，やってみる価値があると思う．著者らがあげている世代間継承のメカニズムについて一点付記すれば，今日のケアワーク選択理由としては積極的/消極的に二極化してる*2と思われるので，そのあたりをうまいマクロ変数を設定して分析することも可能だろう．今回の分析から著者らがいうような，

Care-working parents may transmit values, networks, and human capital resources that are specific to detailed occupations, rather than generic to care work.

こういうことが本当に生じているかは怪しい．

*1:父ケアワークの場合のみ息子が利他傾向

*2:本当はケアワークに就きたくないがその他に選択肢がないのでケアワークに就いたという人もいまや少なくない．

King and Roberts(PA 2015) 頑健標準誤差を無自覚に使ってはいけません

計量社会学のレクチャーで漸近理論の話が出てくることはあまりないが，これはまずいと思っている．漸近的性質については計量屋のフロンティアでいまもたくさん成果が出ているし，それらをフォローするのは難しいが，基本的なことは共有されるべきである．さて，推定値を正当化する漸近的近似が適切でない場合，頑健な標準誤差の推定値が過小になったりする．やみくもに頑健標準誤差やクラスター標準誤差を使えば良いというわけではないということは，モストリーハームレスの8章でも扱わているが，本論文もその類のものである．

King, R. and Robert, M. 2015. "How Robust Standard Errors Expose Methodological Problems They Do Not Fix, and What to Do About It" Political Analysis 23: 159-179.

要旨は以下である．

“Robust standard errors” are used in a vast array of scholarship to correct standard errors for model misspecification. However, when misspecification is bad enough to make classical and robust standard errors diverge, assuming that it is nevertheless not so bad as to bias everything else requires considerable optimism. And even if the optimism is warranted, settling for a misspecified model, with or without robust standard errors, will still bias estimators of all but a few quantities of interest. The resulting cavernous gap between theory and practice suggests that considerable gains in applied statistics may be possible. We seek to help researchers realize these gains via a more productive way to understand and use robust standard errors; a new general and easier-to-use “generalized information matrix test” statistic that can formally assess misspecification (based on differences between robust and classical variance estimates); and practical illustra- tions via simulations and real examples from published research. How robust standard errors are used needs to change, but instead of jettisoning this popular tool we show how to use it to provide effective clues about model misspecification, likely biases, and a guide to considerably more reliable, and defensible, inferences. Accompanying this article is software that implements the methods we describe.

主眼は，頑健標準誤差においてmisspecificationがおきていないかを検定する"generalized information matrix test"の提案である．頑健標準誤差はポリサイでも多用されているようで，著者らによれば，

Among all articles between 2009 and 2012 that used some type of regression analysis published in the American Political Science Review, 66% reported robust standard errors. In International Organization, the figure is 73%, and in American Journal of Political Science, it is 45%.

ほどである．社会学でもかなり無自覚に使っているユーザーは多いと思われるが，通常の標準誤差と頑健標準誤差に大きな違いがある場合には注意が必要となる．かつては不均一分散の検定に使われていたブロッシュ＝ペーガン検定やホワイト検定では，こうした問題をスルーしてしまう(モデルミススペシフィケーションなのか単なる不均一分散があるだけなのかの区別ができないということ)．キングらに言わせればこういうことである．

However, they are often used in applications as a default setting, without justification (sometimes even as an effort to inoculate oneself from criticism), and without regard to the serious consequences their use implies about the likely misspecification in the rest of one’s model. Moreover, a model for which robust and classical standard error estimates differ is direct confirmation of misspecification that extends beyond what the procedure corrects, which means that some estimates drawn from it will be biased—often in a way that can be fixed but not merely by using robust standard errors.

まず頑健標準誤差の成功例を確認しよう．均一分散で線形正規分布を仮定した回帰モデルで最尤推定を行った場合，データ生成過程で実際には不均一分散が生じていたとしても，頑健標準誤差を用いれば，有効性はないが一致性と不偏性をもつ推定量となる．こうした頑健標準誤差の利点に対して，キングらは問題点は2つあると指摘する．

First, even if the functional form, independence, and other specification assumptions of this regression are correct, only certain quantities of interest can be consistently estimated. For example, if the dependent variable is the Democratic proportion of the two-party vote, we can consistently estimate a regression coefficient, but not the probability that the Democrat wins, the variation in vote outcome, risk ratios, vote predictions with confidence intervals, or other quantities. In general, computing quantities of interest from a model, such as by simulation, requires not only valid point estimates and a variance matrix, but also the veracity of the model’s complete stochastic component (King, Tomz, and Wittenberg 2000; Imai, King, and Lau 2008).
Second, if robust and classical standard errors diverge—which means the author acknowledges that one part of his or her model is wrong—then why should readers believe that all the other parts of the model that have not been examined are correctly specified? We normally prefer theories that come with measures of many validated observable implications; when one is shown to be inconsistent with the evidence, the validity of the whole theory is normally given more scrutiny, if not rejected (King, Keohane, and Verba 1994). Statistical modeling works the same way: each of the standard diagnostic tests evaluates an observable implication of the statistical model. The more these observable implications are evaluated, the better, since each one makes the theory vulnerable to being proven wrong. This is how science progresses. According to the contrary philosophy of science implied by the most common use of robust standard errors, if it looks like a duck and smells like a duck, it is just possible that it could be a beautiful blue-crested falcon.

キングらは頑健標準誤差を使うべきではないと言っているのではない．モデル推定において通常の標準誤差と頑健標準誤差に大きな違いがある場合には，ミススペシフィケーションの可能性が高いため，手持ちデータに関する基礎的な診断をすべきだとしている．その診断に応じて適切な推定手法を選ぶべきだと．彼らが結論が述べているように，要はこういうことである．

Robust standard errors should be treated not as a way to avoid reviewer criticism or as a magical cure-all. They are neither. They should instead be used for their fundamental contribution—as an excellent model diagnostic procedure. We strongly echo what the best data analysts have been saying for decades: use all the standard diagnostic tests; be sure that your model actually fits the data; seek out as many observable implications as you can observe from your model. And use all these diagnostic evaluation procedures to respecify your model. If you have succeeded in choosing a better model, your robust and classical standard errors should now approximately coincide.

GIM(generalized information matrix)テストでモデルを改善させた証は，大きな乖離のあった通常の標準誤差と頑健標準誤差がだいたい一致するということである．

O'Brien et al.(SMR 2008) APC混合効果モデル

簡単なメモをしておく．アウトカムに対する年齢(age)，時代(period)，世代(cohort)の効果を識別するのは単純でない．なぜなら「世代＝時代ー年齢」が成立しているからだ．これらを工夫して何とか識別しようとするのがAPCモデルであるが，これは社会学での呼び方であり(おそらく)，経済学や政治学ではなんと呼ばれているか知らない．例えば，経済学で有名なDeaton and Paxson(JPE 1994)ではAPCモデルという語は使われていない．それはさておき，本論文はAPCの世代特性におけるランダム効果推定のメリットを提示するものである．

O’Brien, R. M., Hudson, K., & Stockard, J. (2008). "A Mixed Model Estimation of Age, Period, and Cohort Effects." Sociological Methods & Research, 36(3), 402–428.

APCを同時投入すると完全な多重共線性が生じるため，著者らによればこれまで以下の解決策がとられてきた．

1. APCの関連を線形から非線形にするために少なくとも1つの変数を適当なものに変換する
2. APCの2つ以上の変数についてアウトカムに対する効果が全く同じであると仮定する
3. 特定のコーホートの効果が時間の変化に応じて変化すると仮定する

私は第1の方法のものをよく見る気がするのだが，著者らによれば，上記で最も用いられるのは第2の方法らしい．さて，本稿では次式のように年齢，時代を固定効果，世代特性をランダム効果として推定している(Age Period Cohort Mixed Model: APCMM)．
{ \displaystyle Y_{ij}=\mu+\alpha_i+\pi_j+\chi_k+\upsilon_k+\epsilon_{ij} }
{ \displaystyle \mu }は切片，{ \displaystyle \alpha_i, \pi_j }は年齢と時代の固定効果(ダミー変数)，{ \displaystyle \chi_k}は世代特性，{ \displaystyle \upsilon_k }はランダム誤差項である．APCMMの利点として著者らは以下のように述べている．

They can (1) assess the amount of variance in the dependent variable that is associated with cohorts while controlling for the age and period dummy variables, (2) model the dependencies that result from the age-period-specific rates for a single cohort being observed multiple times, and (3) assess how much of the variance in observations that is associated with cohorts is explained by differences in the characteristics of cohorts.

世代の固定効果はコントロールされていないので，年齢と時代の効果については留意が必要となるのは言うまでもない．APCのうち1つ以上の効果が0である場合にはOLSで推定しれやれば良いが，APCすべての効果を無視できない場合には本論文のように混合効果モデルで推定するのもひとつの手である．

Western and Bloome(SM 2009) Variance Function Regressionによる不平等の分析

最近Bruce Westernがよく使っているVariance Function Regression(VFR)を本人が解説した論文．国内学会の報告等でもVFRは最近みかける．

Western, B., and Bloome, D. 2009. "VARIANCE FUNCTION REGRESSIONS FOR STUDYING INEQUALITY." Sociological Methodology 39(1): 293-326.

通常の回帰分析は女性より男性のほうが高い賃金を予測するとか，大卒者が非大卒より高い階層帰属意識であるとか，グループ間の比較(between-group)を主眼としている．それに対してVFRはグループ内(within-group)の分散に着目し，女性や男性の賃金の不平等のばらつきに何が影響を与えているのかに関心がある．簡潔に述べれば，回帰分析で推定された残差分散を分析するのがVFRである．手順としては，第1に被説明変数を説明変数に回帰して{ \displaystyle \beta_{ki} }を推定し残差を保存，第2に(対数)残差二乗をさらに説明変数に回帰し{ \displaystyle \lambda_{ki} }を得る．平均値と分散は逆の動きをすることが多いので{ \displaystyle \beta_{ki} }と{ \displaystyle \lambda_{ki} }はだいたいが逆の符号になるが，{ \displaystyle \beta_{ki} }と{ \displaystyle \lambda_{ki} }が同じ符号になっていればストーリとしては面白くなるだろう．推定方法について，本稿で取り上げられているのは，(1)GLM，(2)ML，(3) REML，(4)Bayesである．GLM以外は先の手順を繰り返す．著者らはモンテカルロシミュレーションを行い，バイアスが小さいのは(3) REMLと(4)Bayesでありこれらを用いることを推奨している．VFRを用いれば不平等変化の要因分解も可能である．不平等の変化は，a)between-group不平等，b)within-group不平等，c)グループ(集団)比率の変化に分解できる．この計算は簡単にできるし，少しのプログラムでベイズでもできる．要因分解についてはLemieux(AER 2006)に従っているのでそちらを参照するほうが分かりやすいし詳しい．VFRは記述的な分析なので，モデルを慎重にたてる必要がある．とりわけステップ2を意識して欠落変数バイアスにご注意あれ．ちなみにstataでのVFRの一例は以下．

reg Y ‘X’;
predict R, r;
gen R2=Rˆ2;
glm R2 ‘Z’, family(gamma) link(log); predict S2, mu;
gen LOGLIK=−.5∗(ln(S2)+(R2/S2)); egen LL0 = sum(LOGLIK);
di LL0;
∗Updating beta and lambda coefficients; gen DLL=1;
while DLL > .00001 {;
drop R;
quietly: reg Y ‘X’ [aw=1/S2]; drop S2;
predict R, r;
replace R2=Rˆ2;
est store BETA;
quietly:
glm R2 ‘Z’, family(gamma) link(log); predict S2, mu;
est store LAMBDA;
replace LOGLIK=−.5∗(ln(S2)+(R2/S2)); egen LLN = sum(LOGLIK);
di LLN;
replace DLL=LLN-LL0; replace LL0=LLN; drop LLN;
};
est table BETA LAMBDA, b se;