multicoxのテスト

データの読み込み

library(KMsurv)

There were 18 warnings (use warnings() to see them)

data(bmt)

multicox functionを定義

multicox<-function(d,ls,digits=2,time,event, cutoff=0.1, topn=NA){
# eval
data<-eval(parse(text = d)) 
time<-eval(parse(text=paste0(d,"$",time)))
event<-eval(parse(text=paste0(d,"$",event)))
# function makep
makep<-function(x){
  pvalue<-ifelse(x>=0.1,as.character(round(x,digits=2)),
                 ifelse(x<0.1 & x>=0.01, as.character(round(x,digits=3)),
                        ifelse(x<0.01 & x>=0.001,as.character(round(x,digits=3)),
                               ifelse(x<0.001,"< 0.001",NA))))
  return(pvalue)}
# function univariate cox
f.cox<-function(x)
  {
  # library
  library(survival)
  library(dplyr)
  library(broom)
  library(stringr)
  library(purrr)
    # univariate cox
  factor<-eval(parse(text = paste0(d,"$",x)))
  u.cox<-coxph(formula = Surv(time, event)~ factor, data = data)
  u.cox<-u.cox %>% 
    tidy %>% 
    dplyr::mutate(HR=exp(estimate), HR_conf.low=exp(conf.low), HR_conf.high=exp(conf.high))%>%
    dplyr::select(term,starts_with("HR"),p.value)%>%
    dplyr::mutate(term=str_replace(term,"factor",x))%>%
    dplyr::mutate(P.value=makep(p.value))%>%
    dplyr::mutate_at (vars(starts_with("HR")),funs (round(.,digits = digits)))%>%
    dplyr::rename(Factor=term)%>%
    dplyr::select(Factor,HR,HR_conf.low,HR_conf.high,P.value,p.value)
  }
# purrr
res.u<-map_dfr(ls,f.cox) 
res.u<-cbind(res.u,ls)%>%
dplyr::mutate(ls=as.character(ls))
# making formula
res.ls<-res.u %>%
  dplyr::filter(p.value< cutoff) %>%
  dplyr::arrange(p.value)
#topn
ifelse(is.na(topn) == TRUE , res.ls<-res.ls,
       res.ls <- res.ls %>%
        dplyr::top_n(p.value, n = topn )
       )
  
ls.factor<-as.list(res.ls$ls) %>%
  paste0(.,collapse = " + ")
# multivariate lear regression
object<-paste0(
"m.cox<-coxph(formula = Surv(time, event)~",ls.factor , ",data = data)"
)
eval(parse(text=object))
# tidy and cleaning
res.multi<-m.cox %>%
  tidy%>% 
  dplyr::mutate(HR=exp(estimate), HR_conf.low=exp(conf.low), HR_conf.high=exp(conf.high)) %>%
  dplyr::select(term,starts_with("HR"),p.value)%>%
  dplyr::mutate(P.value=makep(p.value))%>%
  dplyr::mutate_at (vars(starts_with("HR")),funs (round(.,digits = digits)))%>%
  dplyr::rename(Factor=term)%>%
  dplyr::select(Factor,HR,HR_conf.low,HR_conf.high,P.value,p.value)
# cbind uni and multi
res.final<-res.u%>%
  dplyr::left_join(res.multi,by="Factor",suffix=c("_Univariate","_Multivariate"))%>%
  dplyr::select(-(p.value_Univariate),-(p.value_Multivariate),-(ls))
res.final
}

multicoxの実行

リストの作成

• factorsにはcoxの単変量を行う因子のリストを作成

multicox関数

• d データフレームの名前
• ls 単変量を行うリストの指定
• event イベントの有無
• time イベントまでの時間
• cutoff 単変量でpvalueがこの値より小さい場合は多変量に組み込む
• HR の丸めの単位 初期設定では小数点第2位
• topn　イベント数が少ないときや、多変量に投入する因子数を制限するときに使用　もし因子数を3つまでであれば topn=3などと指定 下の実行例ではcutoffを切る変数がもともと2個なので何も効果はありませんが

#coxにいれる因子のリスト
factors<-c("z1","z2","z3","z4","z5","z6","z7","z9","z10")
tab<-multicox(d = "bmt",ls = factors, event = "d1",time = "t1", cutoff = 0.1, digits = 2, topn = 3)
print(tab)