Homec4science
Diffusion hcv_modeling (master) Browse

R3744/HCV_model_exe_reinfection.R
a82971ce17f4master

Last Change
Raw File

HCV_model_exe_reinfection.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('\n')
cat('***********************************************************\n')
cat('\n Re-infection \n')
cat('***********************************************************\n')
########################################################################
require(MASS)
require(msm)
require(mstate)
require(survival)
require(splines)
require(plyr)
require(gems)
#######################################################################
start.time <- Sys.time()
maxTime <- 100 # specify duration of simulation in years
#diag_ind <- readline(prompt = "Please enter the screening scenario: 1.Baseline,2.Birth, 3.IDU, 4.ex_IDU")
diag_ind = 1
#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 <- "All"
FP_scenario <- 1
#########################################################################
# setwd("~/HCV_project_2017/hcv_git/Tables/outputDirectorySeptember042017/")
# all_cases = matrix(NA,1,52 + 13)
# # for (case in 1:nrow(data))
# for (case in 1:10)
# {
# temp = read.table(sprintf("treatFP_baseline%s_September042017", case))
# all_cases = rbind(as.matrix(all_cases),as.matrix(temp))
# }
#setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
########################################################################
all_cases_reinfection = matrix(NA,1,52 + 16)
###################################################################################################################
bl_reinf <- data.frame(matrix(vector(), nrow(all_cases), 17, dimnames =
list(c(),c("Alcohol","Gender","HIV","HIV_time","MSM","IDU","IDU_time_start",
"IDU_time_stop", "Genotype" , "Origin", "Age", "BirthYear","initial_index_fs",
"patient_ID", "death_reason", "death_reason", "IDU_former"))), stringsAsFactors = F)
###################################################################################################################
# 1. run the origin code first
#setwd("../Required_Data")
#data <- read.table("all_casesAugust162017", header = T) # check folder hcv_modeling, file HCV_model_exe2.R to make sure
data = all_cases
data[1,] = 0
##################################################################################################################
# 2. pick some individuals randomly; choose their baseline characteristics, and the liver status after they had SVR
cohortSize = 0
for (case in 1:nrow(data)) ## we need to add up all the baselines to the individuals
{
## Status
status = 0
status_IDU = 0
status_svr = 0
status_HIV_MSM = 0
status_svr = 0
patient_id = case
##################################################
p = 0.5
stage = max(which(!is.na(data[case,1:48])))
if (is.element(stage,c(26,27,28,29,30,36,42,48)))
status_svr = 1 # experienced SVR
###################################################
## we need to know the first time that the person experienced SVR
#################################################
ss = data[case, stage]
svr_ind = c(26:30, 36, 42, 48)
if (status_svr == 1)
{
ind = which(complete.cases(data[case,svr_ind]))
ind = min(ind)
stage = svr_ind[ind]
ss = data[case, stage]
}
##########################################################################################################################
if (data[case, "IDU"] == 1 && ss <= data[case,"IDU_time_stop"])
status_IDU = 1 # current IDU
if (data[case, "MSM"] == 1 && ss >= data[case,"HIV_time"])
status_HIV_MSM = 1 # current IDU
########################
W = runif( 1, 0, 1)
status = ifelse(W < p && status_svr == 1 && (status_IDU == 1), 1, 0)
status = ifelse(W < p && status_svr == 1 && (status_HIV_MSM == 1), 1, status)
if (status == 0 && status_svr == 1)
{
W = runif( 1, 0, 1)
status = ifelse(W < 0.01, 1, 0)
}
######################################################################################################################
if (status == 1)
{
cohortSize = cohortSize + 1
blocks_size = 5
blocks_nb = 6
ref0 = cbind(rep(1:blocks_size, blocks_nb), 1:(ns - 22))
ref1 = cbind(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6), c(31,37,43,32,38,44,33,39,45,34,40,46,35,41,47,36,42,48))
ref2 = rbind(ref0,ref1)
#################################################################################################
index_fs.fun = function(state, ref = ref2)
{
ref[which(ref[,2] == state, arr.ind = TRUE),1]
}
index_fs = index_fs.fun(state = stage, ref = ref2)
########################################################
Death_ind = max(which(!is.na(data[case,49:52])), 0)
Death = data[case, 48 + Death_ind]
stage_Death = max(which(!is.na(data[case,1:52])))
if (Death_ind == 0)
Death = stage_Death
#############################################################################################################
# 3. We assume that x months later they will be reinfected again.
w = runif( 1, 1, 5)
bl_reinf$reinf_time = w # time of re-infecting
#############################################################################################################
### Picking the baseline characteristics
bl_reinf$BirthYear[cohortSize] <- data[case, "BirthYear"]
bl_reinf$Alcohol[cohortSize] <- data[case, "Alcohol"]
bl_reinf$Age[cohortSize] <- data[case, "Age"] + ss + w
bl_reinf$HIV[cohortSize] <- data[case, "HIV"]
if ( status_HIV_MSM == 1)
bl_reinf$HIV_time[cohortSize] <- 0
else
bl_reinf$HIV_time[cohortSize] <- ifelse(data[case, "HIV_time"] <= ss, 0, data[case, "HIV_time"] - ss)
bl_reinf$Origin[cohortSize] <- data[case, "Origin"]
bl_reinf$Genotype[cohortSize] <- data[case, "Genotype"]
bl_reinf$Gender[cohortSize] <- data[case, "Gender"]
bl_reinf$MSM[cohortSize] <- data[case, "MSM"]
bl_reinf$IDU[cohortSize] <- data[case, "IDU"]
bl_reinf$IDU_time_start[cohortSize] <- data[case, "IDU_time_start"] - ss # if it is negative it means that they started before this simulation
bl_reinf$IDU_former[cohortSize] <- ifelse(data[case, "IDU_time_start"] < 999 && data[case, "IDU_time_stop"] <= ss, 1, 0)
bl_reinf$IDU_time_stop[cohortSize] <- data[case, "IDU_time_stop"] - ss # if it is negative it means that they stoped before this simulation
bl_reinf$initial_index_fs[cohortSize] <- index_fs
bl_reinf$death_time[cohortSize] <- Death
bl_reinf$death_reason[cohortSize] <- Death_ind + 48
bl_reinf$patient_ID[cohortSize] <- patient_id
}
}
# 4. run the code for them for the second time, with the modified baseline characteristics and new parameters
wd = getwd()
#####################Source Codes ###########################
source("baseline_reinfection.R") # call file specifying baseline characteristics
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
#######################################################
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))
mainDir <- wd
setwd(file.path(mainDir, subDir))
all_cases_reinfection = cbind(cohort@time.to.state, bl)
write.table(all_cases_reinfection, sprintf("treatFP_baseline_reinfection_%s", dd))
##########################################################
c4science · Help