Homec4science
Diffusion hcv_modeling (master) Browse

R3744/HCV_model_exe3.R
5f1b74105180master

Last Change
Raw File

HCV_model_exe3.R
View Options

#############Executable file ######################################
#########################################
### Scenario: We need to choose the treatment scenario that we want to apply
### Simulation starts at time of HCV infection
### This file include all the required packages, and call all the required sources
### In this file the Cohort Size and the duration of simulation are defined.
#####################################################################
cat("\014")
cat('\n')
cat('***********************************************************\n')
cat('Mathematical Modeling of Hepatitis C Virus Progression Disease\n')
cat('\t \t 2017')
cat('\n')
cat('\n \t Version: 10')
cat('\n Requires files: \n')
cat('simulate_cohort.r: To do the simulation \n')
cat('baseline.r: Identify the baseline chracteristics \n')
cat('hazard_fun.r: Calls file specifying the hazard functions \n')
cat('params_cov.r: To choose parameters file \n')
cat('scenario_treatFx.r: To choose the treatment scenario \n')
cat('simulate_cohort.r: To do the simulation \n')
cat('TransitionTime.r: Shows all the possible transitions\n')
cat('***********************************************************\n')
#############################################################
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
rm(list = ls())
setwd("~/HCV_project_2017/hcv_git/Tables")
wd <- setwd("~/HCV_project_2017/hcv_git/Tables")
mainDir <- wd
today <- Sys.Date()
dd = format(today, format = "%B%d%Y")
subDir <- sprintf("outputDirectory%s",dd)
dir.create(file.path(mainDir, subDir), showWarnings = FALSE)
wd = setwd(file.path(mainDir, subDir))
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
######################### Required Packages #################
require(MASS)
require(msm)
require(mstate)
require(survival)
require(splines)
require(plyr)
require(gems)
#########################################################
start.time <- Sys.time()
########################################################
cohortSize <- 50 # specify how many patients to simulate
maxTime <- 100 # specify duration of simulation in years
setwd("~/HCV_project_2017/hcv_git/Required_Data/Swiss_data_prepared")
# data <- read.table("SCCS_cases_weigth.txt", header = T)
load("~/HCV_project_2017/hcv_git/Required_Data/Swiss_data_prepared/work_space_catagory_8172017.RData")
data <- categories
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
all_cases = matrix(NA,1,52 + 13)
#########################################################
outcome_averages = data
outcome_averages$mortality_background <- 0
outcome_averages$mortality_IDU <- 0
outcome_averages$mortality_HIV <- 0
outcome_averages$mortality_Liver <- 0
outcome_averages$HCC <- 0
outcome_averages$DC <- 0
outcome_averages$LT <- 0
########################################################
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
for (MSM_status in 0:1)
{
#for (case in 1:nrow(data))
for (case in 1:2)
{
#################################################
print(case)
#################################################
birth_ID = as.numeric(data$BirthYear[case])
data$yearOfbirth[case] <- switch(birth_ID,
sample(1937:1947, cohortSize, replace = T),
sample(1948:1957, cohortSize, replace = T),
sample( 1958:1967, cohortSize, replace = T),
sample(1968:1977, cohortSize, replace = T),
sample( 1978:1987, cohortSize, replace = T),
sample(1988:1997, cohortSize, replace = T),
sample(1998:2007, cohortSize, replace = T),
sample(2008:2016, cohortSize, replace = T)
)
#####################Source Codes ###########################
source("baseline1.r") # call file specifying baseline characteristics
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
cat('***********************************************************\n')
#ans <- readline(prompt="Please enter the treatment scenario you are interested to study: Treating patients in (and after) stage F0, F1, F2, F3, F4 or BL (baseline) \t")
ans <- "Fg"
source(sprintf("scenario_treat%s.r",ans))
source("hazard_fun.r") # call file specifying the hazard functions
source("params_cov.r") #choose parameters file
source("TransitionTime.r")
source("sim_cohort.r") # call the cohort simulator
cat('***********************************************************\n')
######################################################
setwd("~/HCV_project_2017/hcv_git/Tables")
wd <- setwd("~/HCV_project_2017/hcv_git/Tables")
subDir <- sprintf("outputDirectory%s",dd)
dir.create(file.path(mainDir, subDir), showWarnings = FALSE)
wd = setwd(file.path(mainDir, subDir))
mainDir <- wd
setwd(file.path(mainDir, subDir))
write.table(cohort@time.to.state, sprintf("treatFP%s_%s",case,dd))
#########################################################
dc_count = 0
hcc_count = 0
lt_count = 0
for (j in 1:nrow(Cohort))
{
if (!is.na(Cohort[j,31]) | !is.na(Cohort[j,32]) | !is.na(Cohort[j,33]) | !is.na(Cohort[j,34]) | !is.na(Cohort[j,35]) | !is.na(Cohort[j,36]))
dc_count = dc_count + 1
if (!is.na(Cohort[j,37]) | !is.na(Cohort[j,38]) | !is.na(Cohort[j,39]) | !is.na(Cohort[j,40]) | !is.na(Cohort[j,41]) | !is.na(Cohort[j,42]))
hcc_count = dc_count + 1
if (!is.na(Cohort[j,43]) | !is.na(Cohort[j,44]) | !is.na(Cohort[j,45]) | !is.na(Cohort[j,46]) | !is.na(Cohort[j,47]) | !is.na(Cohort[j,48]))
lt_count = lt_count + 1
}
########################################################
Patient = Cohort
Patient[is.na(Patient)] <- 0
outcome_averages$mortality_background[case] = sum(Patient[,49])/sum(!is.na(Cohort[,49]))
outcome_averages$mortality_IDU[case] = sum(Patient[,50])/sum(!is.na(Cohort[,50]))
outcome_averages$mortality_HIV[case] = sum(Patient[,51])/sum(!is.na(Cohort[,51]))
outcome_averages$mortality_Liver[case] = sum(Patient[,52])/sum(!is.na(Cohort[,52]))
outcome_averages$HCC[case] <- hcc_count
outcome_averages$DC[case] <- dc_count
outcome_averages$LT[case] <- lt_count
#######################################################################################
subDir <- sprintf("outputDirectory_with_Baseline%s",dd)
dir.create(file.path(mainDir, subDir), showWarnings = FALSE)
wd_bl = setwd(file.path(mainDir, subDir))
write.table(cbind(cohort@time.to.state, bl), sprintf("treatFP_baseline%s_%s",case, dd))
all_cases = rbind(as.matrix(all_cases),as.matrix(read.table(sprintf("treatFP_baseline%s_%s",case,dd))))
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
}
}
###################################
output_all = 0
for (i in 1:nrow(outcome_averages))
{
output_all = output_all + outcome_averages$mortality_background[i] * outcome_averages$weights$numdup[i] # mortality for Swiss
}
output_all_Background = output_all
print("Background mortality:")
print(output_all_Background)
output_all = 0
for (i in 1:nrow(outcome_averages))
{
output_all = output_all + (!is.na(outcome_averages$mortality_IDU[i])) * outcome_averages$weights$numdup[i] # mortality for Swiss
}
output_all_IDU = output_all
print("IDU mortality:")
print(output_all_IDU)
output_all = 0
#########################################
end.time <- Sys.time()
cat('***********************************************************\n')
print("Required time: ")
print(end.time - start.time)
cat('***********************************************************\n')
cat('You can find all the required data in the Table "all_cases"!! \n')
cat('all_cases are saved in Required_Data folder with suffix of date \n')
setwd("../Required_Data")
today <- Sys.Date()
dd = format(today, format = "%B%d%Y")
write.table(all_cases, "saa.csv")
write.table(all_cases, sprintf("All_cases%s",dd))
write.table(all_cases, sprintf("All_cases%s.csv",dd))
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
# cat('***********************************************************\n')
##################################################################################
hist(outcome_averages$mortality_background + outcome_averages$age + outcome_averages$yearOfbirth,
breaks = 10, col = "red", xlab = "Age of Death", xlim = c(0,100),main = "Histogram with Normal Curve")
hist(outcome_averages$mortality_background, xlim = c(0,100), ylim = c(0,100),main = "mortality_background")
outcome_averages$YearDeath <- 0
for (i in 1:nrow(outcome_averages))
{
outcome_averages$YearDeath[i] = outcome_averages$mortality_background[i] + outcome_averages$age[i] + outcome_averages$yearOfbirth[i]
}
hist(outcome_averages$YearDeath, main = "mortality_background")
par(new = TRUE)
d <- density(outcome_averages$YearDeath)
plot(d, col = "green")
c4science · Help