|
@@ -0,0 +1,67 @@
|
|
|
+#' COX hazards
|
|
|
+#'
|
|
|
+#' @param x data frame
|
|
|
+#' @param vars a vector of variables used in Cox's survival analysis
|
|
|
+#'
|
|
|
+#' @return model as returned by coxph
|
|
|
+#'
|
|
|
+#' @export
|
|
|
+
|
|
|
+cox.hazard<-function(x,vars){
|
|
|
+ vlist=base::paste(vars,collapse='+')
|
|
|
+ f<-stats::as.formula(base::paste('Surv(followup,censored)',vlist,sep='~'))
|
|
|
+ s=survival::coxph(f,data=x)
|
|
|
+ s
|
|
|
+}
|
|
|
+
|
|
|
+#' COX for univariate variables
|
|
|
+#'
|
|
|
+#' @param x data frame containing followup and censored columns
|
|
|
+#' @param var the variable used to explain survival in format mode:name where mode drives unitChange selection, SD for m, 1 for r or logic and X for uX
|
|
|
+#'
|
|
|
+#' @return data frame with a single row for tested variable, with fields variableName,unitChange,HR,l95,u95,pWald,pLikelihood
|
|
|
+#'
|
|
|
+#' @export
|
|
|
+
|
|
|
+cox.univariate<-function(x,var){
|
|
|
+ data=base::c(var)
|
|
|
+ qv=base::strsplit(data,split=':')[[1]]
|
|
|
+ mode=qv[1]
|
|
|
+ v=qv[2]
|
|
|
+ #print(sprintf('Calculating for %s',var))
|
|
|
+ s=cox.hazard(x,v)
|
|
|
+ s1<-base::summary(s)
|
|
|
+
|
|
|
+ #calculate survival probability at mean or prescribed value, which is given as false for binary and as predifined value with categorical data
|
|
|
+ if (mode=='m'){
|
|
|
+ #use sd as a unit change
|
|
|
+ x.dx<-stats::sd(x[,v],na.rm=TRUE)
|
|
|
+ }
|
|
|
+ if (base::substr(mode,1,1)=='u'){
|
|
|
+ #use specified unit
|
|
|
+ x.dx=base::as.numeric(base::substring(mode,2))
|
|
|
+ }
|
|
|
+
|
|
|
+ if (mode=='logic'){
|
|
|
+ #unit is 1
|
|
|
+ x.dx=1
|
|
|
+
|
|
|
+ }
|
|
|
+ if (substr(mode,1,1)=='r'){
|
|
|
+ #the same as for logic
|
|
|
+ x.dx=1
|
|
|
+ }
|
|
|
+ #adjust output data
|
|
|
+ sd=base::sqrt(s$var)*x.dx
|
|
|
+ c=s1$coef[1]*x.dx
|
|
|
+ l95=base::round(base::exp(c-1.96*sd),2)
|
|
|
+ u95=base::round(base::exp(c+1.96*sd),2)
|
|
|
+ HR=base::round(base::exp(c),2)
|
|
|
+ pWald=base::round(s1$coefficients[1,5],2)
|
|
|
+ pLikelihood=base::round(s1$logtest['pvalue'],3)
|
|
|
+ unitChange=base::round(x.dx,2)
|
|
|
+#return as data frame with a single row
|
|
|
+ base::data.frame(varName=v,unitChange=unitChange,HR=HR,low95=l95,up95=u95,pWald=pWald,pLikelihood=pLikelihood)
|
|
|
+
|
|
|
+}
|
|
|
+
|