person-year型の架空データ（その2）の結果を使ってお絵かき

w<-t(matrix(b[1:6],2,3))
rho<-b[7]
X<-cbind(1,seq(1, 20, 1),min(birthY))
U<-cbind(X%*%w,0)

ival<-log(exp(U[,1]/rho)+exp(U[,2]/rho))

#クラスター選択確率
CL1<-exp(rho*ival)/(exp(rho*ival)+exp(U[,3]))
CL2<-exp(U[,3])/(exp(rho*ival)+exp(U[,3]))#無職選択確率

#正規、非正規選択確率の計算
P1<-CL1*exp(U[,1]/rho)/(exp(U[,1]/rho)+exp(U[,2]/rho))
P2<-CL1*exp(U[,2]/rho)/(exp(U[,1]/rho)+exp(U[,2]/rho))

hazard<-cbind(P1,P2,CL2)
matplot(hazard,type="s")
surv<-matrix(NA,20,1)
for(k in 1:20){
  surv[k,]<-prod(CL2[1:k])
}

event<-cbind(surv*hazard[,1],surv*hazard[,2])
matplot(event,type="s")

surv<-cbind(matrix(NA,20,2),surv)
for(k in 1:20){
  surv[k,1]<-sum(event[1:k,1])
  surv[k,2]<-sum(event[1:k,2])
}

matplot(surv,type="s")