Homec4science
Diffusion hcv_modeling (HCV_with_sensitivity_reinfection) Browse

R3744/baseline_Cohort.R
a82971ce17f4HCV_with_sensitivity_re…

Last Change
Raw File

baseline_Cohort.R
View Options

######################################################################
### This code provides all the required baseline chracteristics
### and Specification of the distribution of baseline characteristics
### It includes 4 parameters:
### p_HIVMSM, p_IDU_start, p_IDU_stop, P_HIV
### Also it include the proportion of patients in each fibrosis stage (start)
### We calculate time to get HIV, which depends on the probability of getting HIV after HCV, and it can be a number between 3 month to maxTime, or it never happens.
### We calculate the time of starting and stopping IDU, which depend on the probability of starting and stoping IDU respectively.
### We also have two function:
# - sr.fun to calculate the spontaneous recovery time
# - pw.eval.ext to calculate the background mortality
#####################################################################
setwd("~/HCV_project_2017/hcv_git/Required_Data")
bl_number = 13 # number of baseline characteristics
statesNum = 52 # number of states
#proportion of patients in each fibrosis stages (F0, F1, F2, F3, F4) at simulation start (note: simulation starts at time of HCV infection)
start <- sample(c(1,2,3,4,5), cohortSize, c(0.8590, 0.1509, 0.0184, 0.00239, 0.00015), replace = T)
p_HIVMSM = 0.5
p_IDU_start = 0.1
p_IDU_stop = 1 # probability of stoping IDU
p_HIV = 0.01 # probability of getting HIV after HCV
################################Functions#########################################
#Baseline treatment time( Note: it can be used in scenarios where time between treatment and diagnosis is fixed.)
# It shows how long does it take to be treated after diagnosis and being eligible to be treated
##################################################
#Baseline spontaneous recovery time
###################################
sr.fun = function()
{
logisticfunc <- function(min, max, infl, slo) {min + ((max - min)/(1 + (teta / infl)^ slo))}
sr.t <- rep(0, cohortSize)
for (i in 1:cohortSize)
{
teta = runif( 1, 0, 1)
p = logisticfunc( min = 0, max = 1, infl = 0.25, slo = 2.23)
w <- runif(1, 0, 1)
ifelse( w <= p, sr.t[i] <- teta, sr.t[i] <- 999)
}
sr.t
}#End Of Spontaneous Recovery
#####################################################################################
#BACKGROUND MORTALITY
# will be used in the hazard_fun file to calculate the harzard function for mortality.
#To match age category with death rate (cuts is the value of the start of the intervall eg for eg 20-24 => cuts=20)
#########################################################
pw.eval.ext = function(Cuts, x, func.values, tol = 0.0001)
{
lengthCuts = length(Cuts)
unlist(lapply(x,
function(xp) {
if (xp < max(Cuts))
func.values[ which((Cuts[c(1:(lengthCuts - 1))] - rep(xp, lengthCuts - 1)) *
(Cuts[c(2:lengthCuts)] - tol - rep(xp, lengthCuts - 1)) <= 0 , arr.ind = TRUE)[1] ]
else func.values[lengthCuts]
}
))
}
###################################################################################
bl <- matrix(nrow = cohortSize, ncol = bl_number)
{
colnames(bl) <- c("Alcohol","Gender","HIV","HIV_time","MSM","IDU","IDU_time_start","IDU_time_stop", "Genotype" , "Origin", "Age", "BirthYear", "sr.time")
################################################################################
################################################################################################
# bl[,"Alcohol"] = sample(0:2, cohortSize, c(0.5381, 0.1623, 0.2994), replace=T)
# bl[,"Genotype"] = sample(0:3, cohortSize, c(0.667,0.19,0.016,0.127), replace=T)
# bl[,"Origin"] = sample(0:3, cohortSize, c(0.667,0.19,0.016,0.127), replace=T)
# bl[,"Gender"] = sample(0:3, cohortSize, c(0.667,0.19,0.016,0.127), replace=T)
# bl[,"HIV"] = rbinom(cohortSize, 0, 0.739)
# bl[,"IDU"] = rbinom(cohortSize, 0, 0.435)
# bl[,"MSM"] = rep(0, cohortSize)
# bl[,"Age"]<-rweibull(n=cohortSize, shape=4.23, scale=40.22) # baseLine age (SHCS data)
# #bl[,"MSM"] = rep(1, cohortSize)
########################################################
age = as.matrix(data$age)
year = as.matrix(data$BirthYear)
bl[,"Alcohol"] = as.matrix(data$alcohol)
bl[,"Genotype"] = as.matrix(data$genotype)
bl[,"Origin"] = as.matrix(data$origin)
bl[,"Gender"] = as.matrix(data$gender)
bl[,"HIV"] = as.matrix(as.numeric(data$HIV))
bl[,"IDU"] = as.matrix(as.numeric(data$IDU))
bl[,"MSM"] = rep(0, cohortSize)
bl[,"Age"] = rnorm(cohortSize, mean(age), sd(age))
#bl[,"BirthYear"] = rnorm(cohortSize, mean(year), sd(year))
bl[,"BirthYear"] = as.matrix(as.numeric(data$BirthYear))
#bl[,"MSM"] = rep(1, cohortSize)
########################################################
# setwd("~/HCV_project_2017/hcv_git/Required_Data/Swiss_data_prepared")
# AOB = read.csv("Data_Prepared_SCCS")
#####################################
bl[,"HIV_time"] <- (bl[,"HIV"])
#######################
for (i in 1:cohortSize)
{
bl[i,"HIV_time"] <- 999
w = runif(1, 0, 1)
if (bl[i,"HIV"] == 1)
bl[i,"HIV_time"] = 0 # time having hiv
if (bl[i,"HIV"] == 0)
{
ID = ifelse( w < p_HIV, w, -1)
bl[i,"HIV_time"] = ifelse( ID == -1, 999, runif(1, 0.25, maxTime)) # time having hiv
}
}
############################################################################################
bl[,"IDU_time_start"] <- bl[,"IDU"]
bl[,"IDU_time_stop"] <- bl[,"IDU"]
for (i in 1:cohortSize)
{
if (is.na(bl[i,"IDU"]))
{
bl[i,"IDU"] = 2
bl[,"IDU_time_start"] = 1000
bl[,"IDU_time_stop"] = 1000
}
w = runif(1, 0, 1)
if (bl[i,"IDU"] == 0)
{
IDU_ID = ifelse( w < p_IDU_start, w, 0)
bl[i,"IDU_time_start"] = ifelse( IDU_ID == 0, 999, maxTime * IDU_ID + runif(1, 0, maxTime) * rbinom(1, 1, 0.5)) # time starting idu
}
if (bl[i,"IDU"] == 1)
bl[i,"IDU_time_start"] = 0 # time having hiv
teta = runif(1, 0, 1)
bl[i,"IDU_time_stop"] = ifelse( bl[i,"IDU_time_start"] < maxTime && teta < p_IDU_stop, bl[i,"IDU_time_start"] + 0.2 * teta * maxTime , 999)
}
####################################################
# bl[,"Alcohol"] = rbinom(cohortSize, 1, 55.33) # Bl alcohol consumption (SHCS data)
# bl[,"Genotype"] = sample(0:3, cohortSize, c(59.5, 8.4, 26.4, 5.7), replace=T) # Genotype distibution: bl[,2]=(0,1,2,3)=>gen(1,4,2,3)
# bl[,"Gender"] <- rbinom(cohortSize, 1, 0.5)
# bl[,"HIV"] <- rbinom(cohortSize, 1, 0.13)
##########################################################################################
bl[,"sr.time"] <- sr.fun()
}
###################################################################################
#bm<-read.table("Mort_CH.txt", header=T) # Read file containing mortality rates by age groups in Switzerland
bm <- read.table("Mort_CH_swiss.txt", header = T) # Read file containing mortality rates by age groups in Switzerland
bm_M <- c(bm[,3],rep(1,70)) # to set background mortality of patients >75 years old as = background mortality of patients 75 years old
bm_F <- c(bm[,4],rep(1,70)) # to set background mortality of patients >75 years old as = background mortality of patients 75 years old
c4science · Help