person-year型の架空データ(その1)をベースに、
仮定するデータ発生メカニズムを変更
就業選択に階層構造を仮定する
①就業するか否か
②就業するとしたら、正規雇用か、非正規雇用か
＊単純化のために、ここでは就業形態を正規、非正規の2つにしている
その他の条件はperson-year型の架空データ(その1)と同じ
→以上のような状況を、離散時間入れ子型ロジットモデルで表現する
モデル定義
NestedLogitModel<-function(x,time,birthY){
b1_1<-x[1]#切片
b1_2<-x[2]#切片
b2_1<-x[3]#リスク継続時間
b2_2<-x[4]#リスク継続時間
b3_1<-x[5]#出生コーホート
b3_2<-x[6]#出生コーホート
rho<-x[7]
#効用の計算
U1<-b1_1+b2_1*time[,1]+b3_1*birthY[,1]
U2<-b1_2+b2_2*time[,1]+b3_2*birthY[,1]
U3<-0
#クラスター1のログサム変数
#「正規」「非正規」の間に相関を仮定
ival<-log(exp(U1/rho)+exp(U2/rho))
#クラスター選択確率
CL1<-exp(rho*ival)/(exp(rho*ival)+exp(U3))
CL2<-exp(U3)/(exp(rho*ival)+exp(U3))#無職選択確率
#正規、非正規選択確率の計算
P1<-CL1*exp(U1/rho)/(exp(U1/rho)+exp(U2/rho))
P2<-CL1*exp(U2/rho)/(exp(U1/rho)+exp(U2/rho))
Pr<-cbind(P1,P2,CL2)
return(Pr)
}モデルからデータを発生させる
#①独立変数とパラメータを定義する
#乱数のseedを設定
set.seed(999)
#パラメタの設定
#正規
beta1_1<--3.5
beta1_2<--0.03
beta1_3<-0.022
#非正規
beta2_1<--4
beta2_2<-0.035
beta2_3<-0.02
#rho
rho<-0.2
betas<-c(beta1_1,beta2_1,beta1_2,beta2_2,beta1_3,beta2_3,rho)
n<-1000 #サンプルサイズ
t<-25 #観測年数
py<-n*t
ID<-matrix(0,py,2) #個人ID,リスク継続年数
Y<-matrix(0,py,3) #就業状態
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]
#②NestedLogitModelを用いて選択確率を発生させる
PPr<-NestedLogitModel(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
}
data<-as.data.frame(data)
data<-subset(data,data$V4!=9)
data<-subset(data,data$V5!=9)
data<-as.matrix(data)
colnames(data)<-c("ID","time","birthY","seiki","hiseiki","mushoku")発生させたデータからパラメータ推定値を得る
#ベクトルとして認識しなおす
ID=data[,c("ID")]
time=data[,c("time")]
birthY=data[,c("birthY")]
Y=data[,c("seiki","hiseiki","mushoku")]
#NestedLogitModelの推定
#尤度関数の定義
frnes<-function(x,Y,time,birthY){
b1_1<-x[1]#切片
b1_2<-x[2]#切片
b2_1<-x[3]#リスク継続時間
b2_2<-x[4]#リスク継続時間
b3_1<-x[5]#出生年
b3_2<-x[6]#出生年
rho<-x[7]
#効用の計算
U1<-b1_1+b2_1*time+b3_1*birthY
U2<-b1_2+b2_2*time+b3_2*birthY
U3<-0
#クラスター1のログサム変数
#「正規」「非正規」の間に相関を仮定
ival<-log(exp(U1/rho)+exp(U2/rho))
#クラスター選択確率
CL1<-exp(rho*ival)/(exp(rho*ival)+exp(U3))
CL2<-exp(U3)/(exp(rho*ival)+exp(U3))#無職選択確率
#正規、非正規選択確率の計算
P1<-CL1*exp(U1/rho)/(exp(U1/rho)+exp(U2/rho))
P2<-CL1*exp(U2/rho)/(exp(U1/rho)+exp(U2/rho))
#これらより、尤度関数は
likelihood<-Y[,1]*P1+Y[,2]*P2+Y[,3]*CL2
L<-prod(likelihood) #尤度関数
#対数尤度関数は
loglik<-Y[,1]*log(P1)+Y[,2]*log(P2)+Y[,3]*log(CL2)
LL<-sum(loglik) #対数尤度関数
return(LL)
}
b0<-c(0,0,0,0,0,0,0.15)#パラメータの初期値
#frを最大化する
res<-optim(b0,frnes,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(サンプルサイズ)*自由度
以下は実行結果。
> b
[1] -3.27779339 -3.68524722 -0.02670424 0.03344167 0.01846105 0.01539592 0.19054824
> tval
[1] -13.290818 -12.071085 -2.088800 2.899009 4.970577 4.074132 2.731473
> AIC
[1] 6814.842
> BIC
[1] 6864.681多項ロジットモデルで推定すると
#MNLmodelの推定。
#尤度関数の定義
fr<-function(x,Y,time,birthY){
b1_1<-x[1]#切片
b1_2<-x[2]#切片
b2_1<-x[3]#リスク継続時間
b2_2<-x[4]#リスク継続時間
b3_1<-x[5]#出生年
b3_2<-x[6]#出生年
#効用の計算
U1<-b1_1+b2_1*time+b3_1*birthY
U2<-b1_2+b2_2*time+b3_2*birthY
U3<-0
#選択確率の分母を定義
d<-exp(U1)+exp(U2)+exp(U3)
#選択確率の計算
P1<-exp(U1)/d
P2<-exp(U2)/d
P3<-exp(U3)/d
#これらより、尤度関数は
likelihood<-Y[,1]*P1+Y[,2]*P2+Y[,3]*P3
L<-prod(likelihood) #尤度関数
#よって、対数尤度関数は
loglik<-log(Y[,1]*P1+Y[,2]*P2+Y[,3]*P3)
LL<-sum(loglik) #対数尤度関数
return(LL)
}
b0<-c(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(サンプルサイズ)*自由度
結果
> b
[1] -3.382198060 -4.969953179 -0.116013794 0.120316407 0.023009270 0.008775567
> tval
[1] -11.164212 -12.533858 -11.328430 15.260224 5.128217 1.491521
> AIC
[1] 6855.353
> BIC
[1] 6898.073→入れ子型ロジットモデルのほうがあてはまりがよい
