Created
February 18, 2020 11:00
-
-
Save kstawiski/1ff37522ae8ce1f604bc97782cd63951 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
suppressMessages(library(rmdformats)) | |
suppressMessages(library(dplyr)) | |
suppressMessages(library(plyr)) | |
suppressMessages(library(tidyverse)) | |
suppressMessages(library(reshape2)) | |
suppressMessages(library(xlsx)) | |
suppressMessages(library(psych)) | |
suppressMessages(library(knitr)) | |
suppressMessages(library(gplots)) | |
suppressMessages(library(car)) | |
suppressMessages(library(kableExtra)) | |
suppressMessages(library(RVAideMemoire)) | |
#suppressMessages(library(rms)) | |
suppressMessages(library(survival)) | |
suppressMessages(library(cmprsk)) | |
# install.packages("biostatUZH", repos="http://R-Forge.R-project.org") | |
suppressMessages(library(biostatUZH)) | |
suppressMessages(library(cr17)) | |
suppressMessages(library(ggpubr)) | |
suppressMessages(library(crrstep)) | |
suppressMessages(library(mice)) | |
suppressMessages(library(PairedData)) | |
suppressMessages(library(Hmisc)) | |
suppressMessages(library(corrplot)) | |
suppressMessages(library(profileR)) | |
suppressMessages(library(bestNormalize)) | |
suppressMessages(library(mgcv)) | |
suppressMessages(library(lm.beta)) | |
suppressMessages(library(olsrr)) | |
suppressMessages(library(inflection)) | |
suppressMessages(library(data.table)) | |
#suppressMessages(library(usdm)) | |
suppressMessages(library(survminer)) | |
# devtools::install_github("zabore/ezfun") | |
suppressMessages(library(ezfun)) | |
library(tibble) | |
library(janitor) | |
# Funkcje generyczne | |
ks.wykreskorelacji = function(var1, var2, labvar1, labvar2, title, yx = T, metoda = 'pearson', gdzie_legenda = "topleft") { | |
var1 = as.numeric(as.character(var1)) | |
var2 = as.numeric(as.character(var2)) | |
plot(var1, var2, pch = 19, cex=0.5, xlab=labvar1, ylab=labvar2, main=title) | |
if (yx == T) { abline(0,1, col='gray') } | |
abline(fit <- lm(var2 ~ var1), col='black', lty = 2) | |
temp = cor.test(var1, var2, method = metoda) | |
if (metoda=='pearson') { | |
if (temp$p.value < 0.0001) { | |
legend(gdzie_legenda, bty="n", legend=paste("r =",round(cor(var1,var2,use="complete.obs"),2),", p < 0.0001\nadj. R² =", format(summary(fit)$adj.r.squared, digits=4))) | |
} else { | |
legend(gdzie_legenda, bty="n", legend=paste("r =",round(cor(var1,var2,use="complete.obs"),2),", p =",round(temp$p.value, 4),"\nadj. R² =", format(summary(fit)$adj.r.squared, digits=4))) } } | |
if (metoda=="spearman") | |
{ if (temp$p.value < 0.0001) { | |
legend(gdzie_legenda, bty="n", legend=paste("rho =",round(cor(var1,var2,use="complete.obs"),2),", p < 0.0001")) | |
} else { | |
legend(gdzie_legenda, bty="n", legend=paste("rho =",round(cor(var1,var2,use="complete.obs"),2),", p =",round(temp$p.value, 4))) } } | |
} | |
ks.correct01 = function(x, eng = T) { | |
if (eng == T) { | |
y = as.factor(ifelse(x==1, "1 Yes", "0 No")) } | |
else { y = as.factor(ifelse(x==1, "1 Tak", "0 Nie")) } | |
y | |
} | |
ks.correctvar = function(data, nvar, func) { | |
newdata = do.call(func, list(data[, as.numeric(nvar)])) | |
return(newdata) | |
} | |
ks.excelnumericdate = function(x) | |
{ | |
excel_numeric_to_date(as.numeric(as.character(x)), date_system = "modern") | |
} | |
ks.plot_with_save = function(templot, tempname, height = 8.27, width = 11.69) | |
{ | |
#tiff(paste0("Plots/",tempname,".tiff"), res = 300, width = width, height = height) | |
#templot | |
#dev.off() | |
#ggsave(paste0("Plots/",tempname,".tiff"), plot = templot, width = width, height = height, dpi = 300) | |
templot | |
} | |
ks.competingrisk.univariate.nominal = function(daneOS, variable, variable_name, survival_name = "OS", time_name = "Months", our_cencode = "Alive", our_failcode = "Related", adjust_to = NA) { | |
tempdaneOS = dplyr::select(daneOS, paste0(survival_name,c("","time"))) | |
colnames(tempdaneOS) = c("Event","Time") # czyli OS -> OS i OStime | |
#tempdaneOS$Time = as.numeric(as.character(tempdaneOS$Time)) | |
tempdaneOS = cbind(tempdaneOS, data.frame(selected = as.factor(variable), adjust_to = adjust_to)) | |
tempdaneOS$selected[tempdaneOS$selected == ""] = NA | |
tempdaneOS$selected = as.factor(as.character(tempdaneOS$selected)) | |
tempdaneOS = tempdaneOS[complete.cases(tempdaneOS),] | |
resCumIncByDis <- cuminc(ftime = tempdaneOS$Time, # failure time variable | |
fstatus = tempdaneOS$Event, # variable with distinct codes for different causes of failure | |
group = tempdaneOS$selected, # estimates will calculated within groups | |
## strata = , # Tests will be stratified on this variable. | |
rho = 0, # Power of the weight function used in the tests. | |
cencode = our_cencode, | |
## subset = , | |
## na.action = na.omit | |
) | |
plot1 = ggcompetingrisks(resCumIncByDis, ggtheme = theme_bw(), ylab = paste0(survival_name," probability"), xlab = time_name, main = paste0(survival_name, " vs. ", variable_name), legend.title = "", multiple_panels = F) | |
resCumIncByDis | |
plot2 = ggcompetingrisks(resCumIncByDis, ggtheme = theme_bw(), ylab = paste0(survival_name," probability"), xlab = time_name, main = paste0(survival_name, " vs. ", variable_name), legend.title = "") | |
if(!is.na(adjust_to)) { | |
covs1 <- model.matrix(~ as.factor(selected) + as.numeric(adjust_to), data = tempdaneOS)[, -1] | |
} else { | |
covs1 <- model.matrix(~ as.factor(selected), data = tempdaneOS)[, -1] } | |
shr_fit <- | |
crr( | |
ftime = tempdaneOS$Time, | |
fstatus = tempdaneOS$Event, | |
cov1 = covs1, | |
failcode = our_failcode, | |
cencode = our_cencode | |
) | |
shr_fit | |
tempdaneOS2 = as.data.frame(tempdaneOS) | |
tempsurv = survfit(Surv(Time,Event == "Related") ~ as.factor(selected), data=tempdaneOS2, type = "kaplan-meier", error = "greenwood", conf.type = "log-log") | |
names(tempsurv$strata) = levels(tempdaneOS$selected) | |
plot_kaplan = ggsurvplot( | |
fit = tempsurv,data=tempdaneOS2, | |
xlab = time_name, | |
ylab = paste0(survival_name," probability"), | |
risk.table = "abs_pct", # absolute number and percentage at risk. | |
risk.table.y.text.col = T,# colour risk table text annotations. | |
risk.table.y.text = T,# show bars instead of names in text annotations | |
# in legend of risk table. | |
ncensor.plot = F, # plot the number of censored subjects at time t | |
surv.median.line = "hv", # add the median survival pointer. | |
) | |
plot_kaplan | |
median_survival = surv_median(tempsurv) | |
res = cbind("Parametr" = c(paste0(variable_name, " (", levels(tempdaneOS$selected)[2:length(levels(tempdaneOS$selected))], " vs. ", levels(tempdaneOS$selected)[1], ")"),"Adjustment") | |
,ezfun::mvcrrres(shr_fit), "Gray's test (p-value)" = resCumIncByDis$Tests[1, 2]) | |
list(result_table = res[-nrow(res),], model = shr_fit, cuminc = resCumIncByDis, plot1 = plot1, plot2 = plot2, plot_kaplan = plot_kaplan, median_surv = median_survival, used_data = tempdaneOS) | |
} | |
ks.competingrisk.univariate.continous = function(daneOS, variable, variable_name, survival_name = "OS", time_name = "Months", our_cencode = "Alive", our_failcode = "Related", adjust_to = daneOS$CzasDoRT) { | |
tempdaneOS = dplyr::select(daneOS, paste0(survival_name,c("","time"))) | |
colnames(tempdaneOS) = c("Event","Time") # czyli OS -> OS i OStime | |
tempdaneOS$Time = as.numeric(as.character(tempdaneOS$Time)) | |
tempdaneOS = cbind(tempdaneOS, data.frame(selected = as.character(variable), adjust_to2 = adjust_to)) | |
tempdaneOS$selected[tempdaneOS$selected == ""] = NA | |
tempdaneOS$selected = as.numeric(tempdaneOS$selected) | |
tempdaneOS = tempdaneOS[complete.cases(tempdaneOS),] | |
if(!is.na(adjust_to)) { | |
covs1 <- model.matrix(~ as.numeric(selected) + as.numeric(adjust_to2), data = tempdaneOS)[, -1] | |
} else { | |
covs1 <- model.matrix(~ as.numeric(selected), data = tempdaneOS)[, -1] } | |
shr_fit <- | |
crr( | |
ftime = tempdaneOS$Time, | |
fstatus = tempdaneOS$Event, | |
cov1 = covs1, | |
cencode = our_cencode, | |
failcode = our_failcode | |
) | |
shr_fit | |
res = cbind("Parametr" = variable_name,ezfun::mvcrrres(shr_fit), "Gray's test (p-value)" = "NA (continous variable)") | |
list(result_table = res[-nrow(res),], model = shr_fit, used_data = tempdaneOS) | |
} | |
ks.competingrisk.univariate = function(daneOS, nominal_variables = NA, continous_variables = NA, nominal_variables_names = nominal_variables, continous_variables_names = continous_variables, ...) | |
{ | |
univariate = list() | |
if(!is.na(nominal_variables)) | |
{ | |
for(i in 1:length(nominal_variables)) | |
{ | |
tryCatch({ | |
univariate[[nominal_variables[i]]] = ks.competingrisk.univariate.nominal(daneOS, as.factor(unlist(dplyr::select(daneOS, nominal_variables[i])[,1])), nominal_variables_names[i], ...)}) | |
} | |
} | |
if(!is.na(continous_variables)) | |
{ | |
for(i in 1:length(continous_variables)) | |
{ | |
tryCatch({univariate[[continous_variables[i]]] = ks.competingrisk.univariate.continous(daneOS, as.numeric(unlist(dplyr::select(daneOS, continous_variables[i])[,1])), continous_variables_names[i], ...)}) | |
} | |
} | |
return(univariate) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment