パーソン・イヤー型の架空データをRを使用し作成
架空データなので、
とりあえずざっくりと以下のような状況を仮定してみる
再就業の分析
・正規雇用からの離職を20-29歳の間に経験した女性が分析対象
・分析対象となるイベントは「再就業」
・就業状態は「正規」「非正規」「自営業」「無職」の四つ
・分析対象の定義より、リスク開始は「正規雇用からの離職」
・「無職」から何らかの就業形態に移行することでイベント発生とし、観察終了
・25年以上移行が発生していないケースは打ち切り
・共変量は出生年、リスク継続年数（線形関係を仮定）
→以上のような状況を、離散時間多項ロジットモデルで表現する
モデル定義
MNLmodel<-function(x,time,birthY){
b1_1<-x[1]#切片
b1_2<-x[2]#切片
b1_3<-x[3]#切片
b2_1<-x[4]#リスク継続年数
b2_2<-x[5]#リスク継続年数
b2_3<-x[6]#リスク継続年数
b3_1<-x[7]#出生年
b3_2<-x[8]#出生年
b3_3<-x[9]#出生年
#効用の計算
U1<-b1_1+b2_1*time[,1]+b3_1*birthY[,1]
U2<-b1_2+b2_2*time[,1]+b3_2*birthY[,1]
U3<-b1_3+b2_3*time[,1]+b3_3*birthY[,1]
U4<-0
#選択確率の分母を定義
d<-exp(U1)+exp(U2)+exp(U3)+exp(U4)
#選択確率の計算
P1<-exp(U1)/d
P2<-exp(U2)/d
P3<-exp(U3)/d
P4<-exp(U4)/d
Pr<-cbind(P1,P2,P3,P4)
return(Pr) #PrをMNLmodelの出力として返せ
}モデルからデータを発生させる
#①独立変数とパラメータを定義する
#乱数の結果を保存
set.seed(999)
#パラメタの設定
#正規
beta1_1<--3.8
beta1_2<--0.02
beta1_3<-0.008
#非正規
beta2_1<--4
beta2_2<-0.03
beta2_3<-0.02
#自営業
beta3_1<--4.5
beta3_2<-0.001
beta3_3<-0.002
betas<-c(beta1_1,beta2_1,beta3_1,beta1_2,beta2_2,beta3_2,beta1_3,beta2_3,beta3_3)
n<-1000 #サンプルサイズ
t<-25 #観測年数
py<-n*t
ID<-matrix(0,py,2) #個人ID,リスク継続年数
Y<-matrix(0,py,4) #就業状態
birthY<-matrix(0,py,1) #出生年
for(i in 1:n){
WKbirthY<-floor(rnorm(n=1,mean=65,sd=10))
for(j in 1:t){
r<-(i-1)*t+j #1からn*tまでのカウント変数として解釈
ID[r,1]<-i #IDのr行1列目にiを代入せよ
ID[r,2]<-j #IDのr行2列目にjを代入せよ
#出生年の定義
birthY[r,1]<-WKbirthY
}
}
time<-matrix(0,py,1)
time[,1]<-ID[,2]
#②MNLmodelを用いて選択確率を発生させる
PPr<-MNLmodel(betas,time,birthY)
#③発生させた選択確率をもとに各時点の就業状態を定義する
for(i in 1:py){
CSPPr<-cumsum(PPr[i,]) #i行目のケースの累積確率ベクトル
rn<-runif(1)
PPM<-which.max(CSPPr>=rn)
Y[i,PPM]<-1 #i行PPM列目のYに1を代入
}
data<-cbind(ID,birthY,Y)
#④再就業発生以降のデータは消去する
for(i in 2:py){
if(data[i,1]==data[i-1,1] && data[i-1,4]==1) data[i,4]<-9
if(data[i,1]==data[i-1,1] && data[i-1,4]==9) data[i,4]<-9
if(data[i,1]==data[i-1,1] && data[i-1,5]==1) data[i,5]<-9
if(data[i,1]==data[i-1,1] && data[i-1,5]==9) data[i,5]<-9
if(data[i,1]==data[i-1,1] && data[i-1,6]==1) data[i,6]<-9
if(data[i,1]==data[i-1,1] && data[i-1,6]==9) data[i,6]<-9
}
data<-as.data.frame(data)
data<-subset(data,data$V4!=9)
data<-subset(data,data$V5!=9)
data<-subset(data,data$V6!=9)
data<-as.matrix(data)
colnames(data)<-c("ID","time","birthY","seiki","hiseiki","ziei","mushoku")発生させたデータからパラメータ推定値を得る
#ベクトルとして認識しなおす
ID=data[,c("ID")]
time=data[,c("time")]
birthY=data[,c("birthY")]
Y=data[,c("seiki","hiseiki","ziei","mushoku")]
#MNLmodelの推定。
#尤度関数の定義
fr<-function(x,Y,time,birthY){
b1_1<-x[1]#切片
b1_2<-x[2]#切片
b1_3<-x[3]#切片
b2_1<-x[4]#リスク継続時間
b2_2<-x[5]#リスク継続時間
b2_3<-x[6]#リスク継続時間
b3_1<-x[7]#出生年
b3_2<-x[8]#出生年
b3_3<-x[9]#出生年
#効用の計算
U1<-b1_1+b2_1*time+b3_1*birthY
U2<-b1_2+b2_2*time+b3_2*birthY
U3<-b1_3+b2_3*time+b3_3*birthY
U4<-0
#選択確率の分母を定義
d<-exp(U1)+exp(U2)+exp(U3)+exp(U4)
#選択確率の計算
P1<-exp(U1)/d
P2<-exp(U2)/d
P3<-exp(U3)/d
P4<-exp(U4)/d
#これらより、尤度関数は
likelihood<-Y[,1]*P1+Y[,2]*P2+Y[,3]*P3+Y[,4]*P4
L<-prod(likelihood) #尤度関数
#よって、対数尤度関数は
loglik<-log(Y[,1]*P1+Y[,2]*P2+Y[,3]*P3+Y[,4]*P4)
LL<-sum(loglik) #対数尤度関数
return(LL)
}
b0<-c(0,0,0,0,0,0,0,0,0)#パラメータの初期値
#frを最大化する
res<-optim(b0,fr,gr=NULL,Y,time,birthY,method="BFGS",
hessian=T,control=list(fnscale=-1))
#パラメータ推定値
b<-res$par
#t値と情報量基準
tval<-b/sqrt(-diag(solve(res$hessian)))
AIC<--2*res$value+2*length(res$par) #-2*対数尤度+2*自由度
BIC<--2*res$value+log(nrow(Y))*length(res$par) #-2*対数尤度+log(サンプルサイズ)*自由度
