比例ハザードモデル

#Cox比例ハザードモデル　中村剛著：医学統計学シリーズ

library(dplyr);library(dummies);library(survival)

data(kidney)

kidney=data.frame(kidney)

delta=kidney$status

kidney=kidney[order(kidney$time),]

disease=dummy(kidney$disease)

kidney=kidney%>%select(-c(disease,id))

kidney=data.frame(kidney,disease)


ite=1000

eta=0.0001

beta=rep(0.1,ncol(kidney)-2)

num=sort(unique(kidney$time))

log_lik=0;log_lik_pre=-100

for(l in 1:ite){

if(abs(log_lik-log_lik_pre)>0.01){

log_lik_pre=log_lik  

dl=rep(0,length(beta))

dl2=array(0,dim=c(length(beta),length(beta)))

log_lik=0;

for(i in num){

kidney_sub=kidney%>%filter(time<=i)%>%select(-time)

status=kidney_sub$status

kidney_sub=kidney_sub%>%select(-status)

dv1=apply(as.matrix(kidney_sub)*c(exp(as.matrix(kidney_sub)%*%beta)),2,sum)

dv2=sum(exp(as.matrix(kidney_sub)%*%beta))

dl=dl+apply(t(t(as.matrix(kidney_sub))-c(dv1/dv2))*c(status) ,2,sum)


dv1=array(0,dim=c(length(beta),length(beta)));dv2=rep(0,length(beta))

for(j in 1:nrow(kidney_sub)){

vec=as.numeric(kidney_sub[j,])

dv1=dv1+t(t(vec))%*%t(vec)*as.numeric(exp(t(vec)%*%t(t(c(beta)))))*status[j]

dv2=dv2+((vec)*exp(sum(vec*beta)))*status[j]

}

dl2=dl2+dv1/sum(exp(as.matrix(kidney_sub)%*%beta))-t(t(c(dv2)))%*%t(c(dv2))/(sum(exp(as.matrix(kidney_sub)%*%beta))^2)

log_lik=log_lik+sum(as.matrix(kidney_sub)%*%beta-log(sum(exp(as.matrix(kidney_sub)%*%beta))))



}

dl2=-dl2

beta=beta-eta*((svd(dl2)$u%*%diag(1/svd(dl2)$d)%*%t(svd(dl2)$v))
%*%dl)

}

print(log_lik)

}

#H:beta=0

var_beta=-((svd(dl2)$u%*%diag(1/svd(dl2)$d)%*%t(svd(dl2)$v)))

W_C=t(c(beta-rep(0,length(beta))))%*%solve(var_beta)%*%t(t(c(beta-rep(0,length(beta)))))    

qchisq(1-0.01,df=ncol(kidney_sub))

#個別に比較

(beta-rep(0,length(beta)))/sqrt(diag(var_beta))

qnorm(1-0.01)