パネル分析でcohort effect，age effect, period effectを同時に分析するには？ (Andreß et al. 2013: Chaper 4)
Hans-Jürgen Andreß, Katrin Golsch and Alexander W. Schmidt, 2013，"Applied Panel Analysis for Economic and Social Surveys", Springer, Heidelberg.
4.2.3 Analysis of Trends の記述をメモ的に書きます．
間違い等あれば，ご指摘いただけると幸いです．
パネル分析のトレンド効果には，3つのタイプがある
| 名称 | 概要 | 人が対象の場合の変数 |
| cohort effect | いつ成立したか(生まれたか) | 生まれた年 (the year of birth: 文中ではgenerationとも) |
| age effect | 開始時点からどれだけ経過したか | 年齢(the age) |
| period effect | 測定した時点そのもの | 調査を行った年 (the time point of the current panel wave ) |
3つの変数は同時に回帰モデルに投入することはできない(連続・カテゴリカル問わず)．
なぜなら， age = t - year_of_birth なので，完全な共線性が発生する．
解決方法としては，
- 特定の重要な期間に生まれた個人のさまざまなコホートを区別する(80~85年生まれ的なことだろうか？)，ageでも同様．
- 年齢を月単位で調査時に聞く⇒244か月(20歳4か月)という値をとらせる．ほかの2つが年単位で聞いてるなら， age = t - year_of_birth は成立しない
- コーホート，年齢，時点のどれかを代替となる変数で表す．(例: 当時のインフレーション率を変数として代用)
ただし上記の解決方法で，完全な共線性は回避できるが，変数間のちょっと相関が高すぎる
Klein and Pötschke(2004)では，7-12年のくくりで，6つのコーホートダミー(d)を作り，モデルは走らせたが，コーホートダミーとageの相関が高すぎた．
そこで，さらにt (period)を時点に対応するインフレ率(i_rate)で代用．((((加えて，maturationの効果を検討するために，ageの代わり，t-1984と(t-1984)二乗も使用．1984はfirst wave))
固定効果モデルでは，コーホートの効果(実際のモデルでは，gender とeducationも)が推定できないので，ML estimation of random effects modelsで推定．
それ以降は，回帰係数βをさらに，fixな部分とrandomな部分に分解して・・・っとしていましたが，ちょっと勉強不足だったので，いったんここまで．
カプランマイヤー曲線をggplot系の図で作りたい【R】【ggsurvplot】【 survminerパッケージ】
Rでカプランマイヤー曲線の図を作る機会があった。
せっかくなので、デフォルトのsurvivalパッケージだけではなく， survminerパッケージのggsurvplotというggplot系コードを使って鮮やかにしてみたので、備忘録的に
今回使うデータとは異なるデータですが、こんな感じのが作れたりします
使うデータと分析
データはMASSパッケージのgehanを使います。
このページで示されてるのと同じ分析をggsurvplotでやってみます
https://www1.doshisha.ac.jp/~mjin/R/Chap_36/36.html
通常のsurvival パッケージでのカプランマイヤー
#パッケージの読み込み
library(survival) #生存時間の関数を作るために呼びだし
library(survminer) # ggsurvplotのパッケージ
library(MASS) # gehanデータを使うために呼び出し
data(gehan) #データの呼び出し
#変数リスト
#pair: Pair index, time: Remission time (weeks), cens: 0=censored, treat: Treatment: control or 6-MP
ge.sf <- survival::survfit(Surv(time,cens)~treat,data=gehan) # モデル
plot(ge.sf,lty=1:2) #プロット
legend("topright" ,c("6-MP投与群","対照群"), lty=c(1,2))
ggsurvplotでカプランマイヤーやってみる
シンプルなカプランマイヤー曲線
ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,,xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray())
# sizeは曲線の太さ， ggthemeでは背景を様々に変えることができます。
ちょっとオプションを足すと・・・
#信頼区間やLog-rank testの結果を示す pval や conf.int をTRUEに
ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,
pval =T ,conf.int = T, xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray())
図の設定をいじりたい
#ggplot系ですので、どんどん後ろにスクリプトを足していくことで、いろいろな設定を変えていけます。
#たとえば さきほどのシンプルな図を temp に代入し，
temp <- ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,,xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray())
temp$plot + scale_x_continuous(breaks=seq(0,40,4))
#と付け足すと
x軸を4weekごとの目盛りを入れることができます。
ちょっとした注意点として，tempに代入後，temp$plot の続きに足していく要素を書かないとエラーが出てきてしまいます。
Risk Table?
あと社会学系では見たことないのですが、リスクテーブルなるものも表示できるようです。
ggsurvplot(ge.sf,dat=gehan, legend.title = "legend",legend=c(.8,.8), size=1,,xlab="xlab", ylab="ylab",censor=T, ggtheme = theme_gray(),risk.table=T)
おまけ
最初の図はこんなコードでつくりました
fit <- survfit(Surv(month,event)~age311_cat,data=dat)
temp <- ggsurvplot(fit, data = dat, size=2,legend.title = "世帯主年齢", legend.labs=c("20未満","20代","30代","40代","50代","60代","70代","80歳以上"),legend=c(.2,.85),conf.int = F, ylab="すまい再建率", xlab="時間", censor=F,pval = F, ggtheme = theme_gray(), fun="event")
temp <- temp$plot + theme(axis.text=element_text(size=30))
temp <- temp + theme(axis.title.x = element_text(size=40),axis.title.y = element_text(size=40))
temp <- temp + scale_x_continuous(breaks=c(0,12,24,36,48,60,72), labels=c("2011/3", "2012/3","2013/3","2014/3","2015/3","2016/3","2017/3"))
temp <- temp+ylim(0,1)
temp + theme(legend.title = element_text(size=38),legend.text = element_text(size=38))
# fun= "event" で y軸を逆転できます。 あとはこのデータは，2011年3月観察開始で，仮設住宅を退居するまでの時間を分析する図なので，x軸は12ヶ月後ごとに設定，そのラベルとして年を示しています。
論文メモ :Henderson et al., 2012, Measuring Economic Growth from Outer Space
個人メモ用ですので、論文の全ポイントというより、個人的に気になるポイントを挙げています
夜間光からGDPの変動の予測法について検討
1)発展途上国などでは，公式に記録された以外の取引などがおおい
2)また統計を取るためのインフラも整っていないこともある
3)都市単位でのGDPは公表されてないことがある
4)場合によってはもっと小さい範囲での総生産がほしいこともある
5) 1)~4)を解決する方法として、夜間光からGDPが予測できればいいよね
6)ただし以下の点は注意が必要
・衛生の夜間光の観測精度は，磨耗・技術的な発展によって異なる可能性がある
・文化によって，GDPに対する夜間光量が異なる可能性がある(あまり電気を使わない文化など)
・時代によってGDPに対する夜間光量が異なる可能性がある
↑のような点を調整するためには，fixed effectでyearの効果を調整しつつ，国間のGDP成長(変動)を検討するのがベスト
ここまでが前提条件
【R: mice: 傾向スコア】 miceの代入済みデータごとの回帰モデル予測値を統合したベクトルを得る
多重代入法による欠損値処理後に，傾向スコア分析をしたい．
ただし予測値を得るのは，代入済みデータセットごとに予測値があるので、なかなかうまくいかずに苦戦しました．
備忘録的にコードを残してきます
間違いがあれば指摘していただけばうれしいです。
なお変数は
z_dummy 割り当て変数 欠損無し
x1 共変量1(連続変数) 欠損あり
x2 共変量2(ダミー変数) 欠損あり
サンプルサイズ(N)を1000
とします
#miceの呼び出し
library("mice")
#多重代入を実行
imputed.data <- mice(data=data, m=100,method=c("",norm","logreg"),seed =1, maxit = 50)
#代入の概要や変数ごとのモデルを確認
summary(imputed)
imputed$method
#回帰モデルの実施
modelMI <- with(imputed.data, glm(z_dummy ~ x1+as.factor(x2), family = binomial(link = 'logit')))
pool(modelMI)
# p_hat に代入データとのモデル予測値を代入
p_hat <- with(imputed, glm(z_dummy ~ x1++as.factor(x2),family = binomial(link = 'logit'))$fitted.values)
#このp_hatは，modalMIと同じlarge miraという形式になってるが，pool()が使えない
#なので無理やり
p <- p_hat$analyses # 予測値だけ取り出す(large mira 形式から large listになる)
n <- 1000 #サンプルサイズをnに代入
DF <- structure(p, row.names = c(NA, -n), class = "data.frame") #listをデータフレームに変更
imputed_propensity_score <- rowMeans(DF[c(1:100)]) #代入済みデータごとの予測値を統合してベクトルに
imputed_propensity_scoreが統合後の傾向スコアになると思います．
論文サマリー: Bolin, 1986, Disaster Impact and Recovery: A Comparizon of Black and White Victims
Robert Bolin, 1986, Disaster Impact and Recovery: A Comparizon of Black and White Victims, "International Journal of Mass Energencies and Disaster"
復興の社会的な格差などの議論の際によく引用されている論文です。
自分の論文に引用するために読んだのでサマリーを作ってみました。
サマリー
Intro
アメリカにおけるBlack / White victims 間での復興の違いを分析
家族の復興はそれまで研究されてきたが、人種/エスニシティに着目した分析は行われてこなかった
Previous Research
災害時におけるマイノリティの研究では、
- 警報・避難行動
- マスメディアの影響
- 避難所での適応
- 住宅再建時の人種/エスニシティの居住地域上の統合/分離
また最近の研究ではソーシャルサポートの効果についても検討されてきている
Family Recovery
Family Recoveryは3つの要素で概念化される
- housing recoery 住宅再建
- financial recovery 経済的復興
- emotinal recovery 感情的な復興
本研究では 2と3に着目して分析
Method
判断分析で分析を行う
経済的復興と主観的ウェルビーイングについて5件法(0~4)で尋ねた
そして、4を復旧・復興済み(complete recovery), 3を復旧途上(intermediate recovery), 0~2を低復旧・復興(low recovery)とカテゴライズ
black, whiteのそれぞれのサブセットを使って、変数をステップワイズ投入することで、black, white victims の復旧・復興の予測因子の違いを検討
Analysis
2つの従属変数の3つの値(complete recovery, intermediate recovery, low recovery)を分ける決定因を分析
係数が0.5以上のものを考察
Discussion
- black victims において、経済的復興に対して、(公的な)金銭的・住宅援助が重要な予測要因となっていた。その一方で親戚間での援助は(経済復興に対して)負の効果を持ってた。
- 経済的復興については、black/white victims 間で予測因の違いはそれほど大きくなかった
- いくつかの違いが存在していおり
- 違いとしてはwhite victims においては、SBAローンは経済復興の重要な予測因となっているが、black victims においてはそうではない⇒black victims においてはローンが利用しにくい⇒black victims の復興の遅れを引き起こす
- black victims とは反対に、white victims では、primaly groupでの援助が経済的復興にプラスの影響を与えていた⇒支援する側の処分可能な所得の量の違いによるものである⇒社会的経済的な構造上の問題
- 主観的ウェルビーイングについては、black victimsにとっては、ソーシャルサポートと(被害による)心理的インパクトが大きな規定要因となっていた
- white victims にとっては、ソーシャルサポートが規定要因となっているのは少ない
- white victimsとっては、家族のケガ、知り合いのケガ、死亡は、心理的復興に大きい影響を与えていた
- 全体的に社会経済的変数は大きな影響は、直接的に影響はなかった
- その一方で、black victimsにとっては、小さい子どもの数が経済的復興に影響を与えてた。これらは経済的ニーズと関連していることから、社会経済的要因も災害からの復興を考える上で欠かせないだろう
感想
よく引用されている印象があったので、まとめてみました
1986年という、おそらく日本では災害に対するソーシャルなアプローチによる研究は非常に限られているであろう時期に、こうした研究(とくに計量的な分析を用いて)がされているのは、さすがって感じですね
非常におもしろいのは、現在の災害社会学の中心的概念であるVunerability（脆弱性）という言葉がでてきていない点でしょう。
この研究がある種 Vunerabilityなどの理論モデル構築の基礎となったということなんだろうと思います。
分析によってやりたいことも明確で、かつ分析結果と結論の間の距離が狭いので非常に読みやすいです
(ただフォントがやや読みにくい・・・)
疑問点
従属変数2つ(経済的復興と主観的ウェルビーイング)を恣意的にカテゴライズする意図が不明。よほど変な分布でもしていたのだろうか？連続変量でいいのでは？
Ｒの 4つのダミー変数化パッケージを使ったときの欠損値の扱いについて比較してみたよ 【4月14日改定】
makedummies関数がNAに対応してので，内容の改定を行いました
最近SPSS, Stata に加えてRを使い始めたんですが、makedummies関数でカテゴリカル変数をダミー変数に変換した際に、もともとNAだったケースにすべて0が入っていたことに腰を抜かしました。
そこで各パッケージのNAの処理を確認してみました。
ほぼ自分のためのメモですが、興味がある方がいらっしゃれば、どうぞ
R初心者なので、間違っている点やもっとシンプルなコードの書き方とかあれば教えてください。
今回検討するのは以下のパッケージに含まれるダミー変数化関数
(パッケージ名/関数名)
・makedummies / mekedummies
・fastDummies / dummy_cols
・dummies / dummy
・caret / dummyVars
他にも有名なダミー化変数があればご指摘ください(追記していきます)
caret以外はダミー変数作成のためのパッケージですね
今回は以下のような架空のデータセットを使います
9・10行目が欠損しています
Ｒのへの読み込み
data <- read.table("Book1.csv", header = T)
str(data)
V1が3水準のFactor型として読み込まれています
makedummies パッケージの makedummies
data.makedummies<- data #このパッケージ用のデータセット作成
install.packages("makedummies", dependencies = TRUE)
library("makedummies")
makedummies(data.makedummies, basal_level = T, col= data.makedummies$V1)
#basal_level=T で3水準すべてのダミーを作成を指定
detach("package:makedummies", unload=TRUE)
NAの8,9,行目はすべて0がはいっている
現在は，NAケースについては，各ダミーに欠損(NA)が入るように更新されました．
fastDummies パッケージの dummy_cols
data.fastDummies <- data #このパッケージ用のデータセット作成
install.packages("fastDummies", dependencies = TRUE)
library("fastDummies")
dummy_cols(data.fastDummies, remove_first_dummy = F, remove_most_frequent_dummy = F)
detach("package:data.fastDummies ", unload=TRUE)
ずれていて見にくいですが一番右の変数名の注目ですね
fasrDummyでは V1_NAが作成され、欠損値のダミーが作られるようです
dummies パッケージの dummy
data.dummies <- data #このパッケージ用のデータセット作成
install.packages("dummies", dependencies = TRUE)
library("dummies")
dummy.data.frame(data.dummies)
detach("dummies", unload=TRUE)
fastDummiesと同様に V1NAが作成されています
caret パッケージの dummyVars
install.packages("caret", dependencies = TRUE)
library("caret")
data.caret<- data #このパッケージ用のデータセット作成
temp <- dummyVars(~.,data=data.caret)
data.caret2 <- as.data.frame(predict(temp,data.caret)) #データフレームとして代入
NAだった8,9行目にはきちんとNAが入っていますね
SPSS, Stataユーザーだとこのダミー変換がもっとも慣れ親しんだものだと思うので、元SPSS, StataユーザーはcaretのdummyVarsを使うのが無難かもしれません
もちろん他の関数でも後からNAを代入すればもちろん問題ないです
私は今後はcaretを使うことにします。
まとめ
・mekedummies / mekedummies ⇒ 各ダミーにNAが入力されている【改定】
・fastDummies / dummy_cols ⇒ NAのダミー変数が作られる
・dummies / dummy ⇒ NAのダミー変数が作られる
・caret / dummyVars ⇒ 各ダミーにNAが入力されている
