.Rhistory

# anyway that shouldn't take me too long, so I can start with that tonight. #DONEZOES. Now I've specified the detectability delays ACCORDING to the
# next on the list is to check the old figures with some mixing included. (imperfect correlation). ALSO, check that zero correlation is what they say it is....
#  Figure 1
################
# : produce a family of sensitivity curves and their average
#goto_1
#Parameters
n=10
# test_details
mean_delay_t1=55
mean_delay_t2=45       #we have a second test listed so that we can easily swap between them when generating the sensitivity curves
sd_size_t1 = 5
sd_size_t2 = 1
detail = 10 # 1/Step size
timeaxis = seq(0,100,1/detail)
scale_t1= 3
shape_t1 = 2
scale_t2= 3
shape_t2 = 2
#visuals
# col_negative <- rgb(27/255,158/255,119/255)
# col_positive <- rgb(217/255,95/255,2/255)
# col_mean <- rgb(231/255,41/255,138/255)
# col_truth <- rgb(117/255,112/255,179/255)
# col_likelihood <- col_truth
# col_dotted <- rgb(3/7,3/7,3/7) #color for dotted lines
#col_negative <- "green" #rgb(27/255,158/255,119/255)
#col_positive <- "red" #rgb(217/255,95/255,2/255)
col_negative <- rgb(27/255,158/255,119/255)
col_positive <- rgb(217/255,95/255,2/255)
col_mean <- rgb(231/255,41/255,138/255)
#col_truth <- "purple"  #rgb(117/255,112/255,179/255)
col_truth <- rgb(117/255,112/255,179/255)
col_likelihood <- col_truth
col_dotted <- rgb(3/7,3/7,3/7) #color for dotted lines
#y-axis limits
ylim_chosen <- c(-.002,1.007)
#Generate Data
sensitivity_family_1 <-  family_sensitivity_weibul(n=n, scale=scale_t1,shape=shape_t1,mean_delay = mean_delay_t1 , sd_size=sd_size_t1,times=timeaxis)
sensitivity_family_background <- family_sensitivity_weibul(n=n+150,scale=scale_t1,shape=shape_t1,mean_delay=mean_delay_t1,sd_size=sd_size_t1,times=timeaxis)
sensitivity_average <- generate_mean_of_family(sensitivity_family_background)
# Plot
# pdf(file='figure_1.pdf',width=7.3,height=4.7)
plot(timeaxis,sensitivity_family_1[,1],type='l',xaxt='n',xaxs='i',yaxs='i',xlim=c(mean_delay_t1-4*sd_size_t1,mean_delay_t1+shift_to_half_likelihood_weibul(shape=shape_t1,scale=scale_t1)+2.35*sd_size_t1),ylim=c(-.005,1.005),xlab="",ylab="",col='green',yaxt='n',bty='L' )
# plot(timeaxis,sensitivity_average,col=col_truth,lwd=5,type='l',xaxt='n',xaxs='i',yaxs='i',xlim=c(mean_delay_t1-4*sd_size_t1,mean_delay_t1+shift_to_half_likelihood_weibul(shape=shape_t1,scale=scale_t1)+2.35*sd_size_t1),ylim=c(-.005,1.005),xlab="",ylab="",yaxt='n',bty='L')
# todo:   title
# label sizes
# colors
# line widths
# scale
#  you're right (alex) about the size of the lines rendering correctly in the pdf
title(xlab="Time since infection", line=1.5, cex.lab=1.2)
title(ylab='Probability of Infection', line=2, cex.lab=1.2)
## goto_fix
# axis ticks, remove box, bring lower axis up to zero or thereabout
yaxis_pos <- c(0,1)
yaxis_names <- c('0','1')
zero_pos <- c(0)
zero_name <- c(expression('0'['']))
axis(side=2, at=yaxis_pos, labels= yaxis_names,tck=-0.037, padj=.237)
#axis(side=1, at=xaxis_pos, labels= xaxis_names,padj=-.35,hadj=-.137)
axis(side=1, at=zero_pos, labels=zero_name,padj=-0.45,hadj=0.37)
for (i in seq(1:n)){
lines(timeaxis,sensitivity_family_1[,i],col=col_negative,lwd=1.5)
}
lines(timeaxis,sensitivity_average,col=col_truth,lwd=5)
# Average delay arrows
Arrows(x0 = mean_delay_t1-4*sd_size_t1, y0 = 0.501 ,x1= timeaxis[which.min(abs(sensitivity_average-0.5))], y1 = 0.5,code=3 ,arr.adj=1)
text(x=(mean_delay_t1-4*sd_size_t1+timeaxis[which.min(abs(sensitivity_average-0.5))])/2, y=c(0.53), pos=4, labels=expression(italic("d")))
#Standard deviation arrows
Arrows(x1=timeaxis[which.min(abs(sensitivity_average-0.5))],y0=0.5, x0= timeaxis[which.min(abs(sensitivity_average-0.5))]+sd_size_t1,y1=0.5,code=3,arr.adj=1)
text(x=  timeaxis[which.min(abs(sensitivity_average-0.5))]+sd_size_t1/2-.5,y=0.53,pos=4,labels=expression(sigma))
# Arrows(c(0,1.7),c(1.3,-1.8),c(0.8,1.1),c(1.2,-1), lwd=2
setwd("C:\\Users\\jeremyb\\Documents\\A_Master\\Code")
}else if(Sys.info()['login']=='jeremy') {
setwd("C:/Users/JumpCo Vostro3700/infection-dating-tool/manuscripts/figures")
}else{
setwd(".") #what does this do?
}
source("part_Onefunctions.R")
file.create("C:\\Users\\jeremyb\\Documents\\A_Master\\Code\\runlog.txt")
runlog <- file("runlog.txt")
writeLines(capture.output(sessionInfo()),runlog)
close(runlog)
?writeLines
?date
##                    SCRIPT
# Now we are going towards risk estimation
# So we define the quick transition from non-infectious to infectious, and the probability of detection (the 'sensitivity' of the test), both with Weibul distributions
n = 5
detail <- 10
timeaxis <- seq(0,90,1/detail)
shape_infect <- 1
scale_infect <- 5
mean_delay_infect <- 7
population_sd_infect <- 10
shape_detect <- 5
scale_detect <- 5
mean_delay_detect <- mean_delay_infect + 14
population_sd_detect <- population_sd_infect
donation_time <- 77
# now we generate and plot the curves
infectious_delays <- generate_positions_cumulative_normal(n=n,mean_center_position = mean_delay_infect, sd_size = population_sd_infect)
detectable_delays <- infectious_delays + 14
infectious_curves <- family_positive_likelihood_weibul_pos(times = timeaxis, shape = shape_infect, scale = scale_infect, positions = infectious_delays, n = n, test_time = donation_time)
undetected_curves <- family_negative_likelihood_weibul_pos(times = timeaxis, shape = shape_detect, scale = scale_detect, positions = detectable_delays, n = n, test_time = donation_time)
# I want a function to take a set of infectious curves and undetected curves and do the three-layered analysis we're proposing
# or it could take the parameters required for all three calculations, and do them.
# then we can take sets of parameters gleaned from literature and stick them together...
#but now how do I calculate the background ones....... if I decide to shift the likelihood by some amount? I'll probably have to define a function or a block of
# code which shifts the infectious delays to get corresponding detectable delays, then just run that block for a much bigger infectiousness set.. no problem
# need to think about whether there is any difference between defining that everyone waits two weeks from infectious to detectable,but the total delay is distributed, and defining that the infectiousness
# infectiousness aand detectability as distributed, using the same distribution and different means (delta is two weeks). No difference in the non-mixed case, but say we want to
# allow the time between infectiousness and detectability to vary. In that case it is more transparent to define a bunch of infectiousness times and delay-between-infectiousness-and-detectability-times,
# then calculate detectability times accordingly. That way we can change the map from 1 - 1 to 1 - many.... Defining them conditionally seeeems to make more sense but I'm not sure.
# anyway that shouldn't take me too long, so I can start with that tonight. #DONEZOES. Now I've specified the detectability delays ACCORDING to the
# next on the list is to check the old figures with some mixing included. (imperfect correlation). ALSO, check that zero correlation is what they say it is....
#  Figure 1
################
# : produce a family of sensitivity curves and their average
#goto_1
#Parameters
n=10
# test_details
mean_delay_t1=55
mean_delay_t2=45       #we have a second test listed so that we can easily swap between them when generating the sensitivity curves
sd_size_t1 = 5
sd_size_t2 = 1
detail = 10 # 1/Step size
timeaxis = seq(0,100,1/detail)
scale_t1= 3
shape_t1 = 2
scale_t2= 3
shape_t2 = 2
#visuals
# col_negative <- rgb(27/255,158/255,119/255)
# col_positive <- rgb(217/255,95/255,2/255)
# col_mean <- rgb(231/255,41/255,138/255)
# col_truth <- rgb(117/255,112/255,179/255)
# col_likelihood <- col_truth
# col_dotted <- rgb(3/7,3/7,3/7) #color for dotted lines
#col_negative <- "green" #rgb(27/255,158/255,119/255)
#col_positive <- "red" #rgb(217/255,95/255,2/255)
col_negative <- rgb(27/255,158/255,119/255)
col_positive <- rgb(217/255,95/255,2/255)
col_mean <- rgb(231/255,41/255,138/255)
#col_truth <- "purple"  #rgb(117/255,112/255,179/255)
col_truth <- rgb(117/255,112/255,179/255)
col_likelihood <- col_truth
col_dotted <- rgb(3/7,3/7,3/7) #color for dotted lines
#y-axis limits
ylim_chosen <- c(-.002,1.007)
#Generate Data
sensitivity_family_1 <-  family_sensitivity_weibul(n=n, scale=scale_t1,shape=shape_t1,mean_delay = mean_delay_t1 , sd_size=sd_size_t1,times=timeaxis)
sensitivity_family_background <- family_sensitivity_weibul(n=n+150,scale=scale_t1,shape=shape_t1,mean_delay=mean_delay_t1,sd_size=sd_size_t1,times=timeaxis)
sensitivity_average <- generate_mean_of_family(sensitivity_family_background)
# Plot
# pdf(file='figure_1.pdf',width=7.3,height=4.7)
plot(timeaxis,sensitivity_family_1[,1],type='l',xaxt='n',xaxs='i',yaxs='i',xlim=c(mean_delay_t1-4*sd_size_t1,mean_delay_t1+shift_to_half_likelihood_weibul(shape=shape_t1,scale=scale_t1)+2.35*sd_size_t1),ylim=c(-.005,1.005),xlab="",ylab="",col='green',yaxt='n',bty='L' )
# plot(timeaxis,sensitivity_average,col=col_truth,lwd=5,type='l',xaxt='n',xaxs='i',yaxs='i',xlim=c(mean_delay_t1-4*sd_size_t1,mean_delay_t1+shift_to_half_likelihood_weibul(shape=shape_t1,scale=scale_t1)+2.35*sd_size_t1),ylim=c(-.005,1.005),xlab="",ylab="",yaxt='n',bty='L')
# todo:   title
# label sizes
# colors
# line widths
# scale
#  you're right (alex) about the size of the lines rendering correctly in the pdf
title(xlab="Time since infection", line=1.5, cex.lab=1.2)
title(ylab='Probability of Infection', line=2, cex.lab=1.2)
## goto_fix
# axis ticks, remove box, bring lower axis up to zero or thereabout
yaxis_pos <- c(0,1)
yaxis_names <- c('0','1')
zero_pos <- c(0)
zero_name <- c(expression('0'['']))
axis(side=2, at=yaxis_pos, labels= yaxis_names,tck=-0.037, padj=.237)
#axis(side=1, at=xaxis_pos, labels= xaxis_names,padj=-.35,hadj=-.137)
axis(side=1, at=zero_pos, labels=zero_name,padj=-0.45,hadj=0.37)
for (i in seq(1:n)){
lines(timeaxis,sensitivity_family_1[,i],col=col_negative,lwd=1.5)
}
lines(timeaxis,sensitivity_average,col=col_truth,lwd=5)
# Average delay arrows
Arrows(x0 = mean_delay_t1-4*sd_size_t1, y0 = 0.501 ,x1= timeaxis[which.min(abs(sensitivity_average-0.5))], y1 = 0.5,code=3 ,arr.adj=1)
text(x=(mean_delay_t1-4*sd_size_t1+timeaxis[which.min(abs(sensitivity_average-0.5))])/2, y=c(0.53), pos=4, labels=expression(italic("d")))
#Standard deviation arrows
Arrows(x1=timeaxis[which.min(abs(sensitivity_average-0.5))],y0=0.5, x0= timeaxis[which.min(abs(sensitivity_average-0.5))]+sd_size_t1,y1=0.5,code=3,arr.adj=1)
text(x=  timeaxis[which.min(abs(sensitivity_average-0.5))]+sd_size_t1/2-.5,y=0.53,pos=4,labels=expression(sigma))
# Arrows(c(0,1.7),c(1.3,-1.8),c(0.8,1.1),c(1.2,-1), lwd=2
n=10
detail=10
timeaxis=seq(0,100,1/detail)
#                       TEST 1 (negative)
#     Describe individual (person) test sensitivity form with population mean-delay and standard deviation of delay
# delay is the variable we distribute across the population - it could in principle be anything else of course
#so each individual has the same SHAPE of sensitivity, but different delays
## Visuals
lwd_means <- 4
lwd_ind <- 1.37
lwd_likelihood <- lwd_means - 3
#High scale causes slower swap
#high shape causes quicker and steeper swap AND more symetrical swap
scale_t1= 2
shape_t1 = 1.73
mean_delay_t1 = 12
sd_size_t1 = 10
#   Time of negative test (relative to arbitrary t=0)
test_time_1 = 48
##                      TEST 2 (positive)
scale_t2 = scale_t1*4.7
shape_t2 = shape_t1
mean_delay_t2 = mean_delay_t1*2+10
sd_size_t2 = 1.3*sd_size_t1
#   Time of positive test
test_time_2 = timeaxis[length(timeaxis)]-5
## Generate Data
# Data = the individual likelihood curves for the first (negative) and second (positive) test
#Test 1
plotdata_negative = family_negative_likelihood_weibul(n=n, scale=scale_t1, shape=shape_t1, mean_delay=mean_delay_t1, sd_size= sd_size_t1, times = timeaxis, test_time = test_time_1)
#for generating mean curve
plotdata_negative_background <- family_negative_likelihood_weibul(n=n+150, scale=scale_t1, shape=shape_t1, mean_delay=mean_delay_t1, sd_size= sd_size_t1, times = timeaxis, test_time = test_time_1)
#Test 2
plotdata_positive = family_positive_likelihood_weibul(n=n, scale=scale_t2, shape=shape_t2, mean_delay=mean_delay_t2, sd_size= sd_size_t2, times = timeaxis, test_time = test_time_2)
#for generating mean curve
plotdata_positive_background <- family_positive_likelihood_weibul(n=n+150, scale=scale_t2, shape=shape_t2, mean_delay=mean_delay_t2, sd_size= sd_size_t2, times = timeaxis, test_time = test_time_2)
##        COMMUNICATE
# pdf(file = "simple_case.pdf", width = 7.3, height = 4.7)
plot(timeaxis,plotdata_negative[,1],type='l',xlim=c(timeaxis[1],timeaxis[length(timeaxis)]),ylim=c(-.002,1.002),xaxt='n',yaxt='n',xaxs='i',yaxs='i',bty='l',xlab='',ylab='',col='green') #clarify label in comment
title(xlab="Hypothetical time of infection", line=1.5, cex.lab=1.4)
title(ylab=expression('Likelihood'), line=2, cex.lab=1.4)
yaxis_pos <- c(0,0.5,1)
yaxis_names <- c('0',"",'1')
xaxis_pos <- c(test_time_2)
xaxis_names <- c(expression('T'['Donation']))
zero_pos <- c(0)
zero_name <- c(expression('0'['']))
#shift axes
axis(side=2, at=yaxis_pos, labels= yaxis_names,tck=-0.037, padj=.437)
axis(side=1, at=xaxis_pos, labels= xaxis_names,padj=-.35,hadj=0)#-.137)
#axis(side=1, at=zero_pos, labels=zero_name,padj=-0.45,hadj=0.37)
# goto_do
#points(plotdata[,1],plotdata[,3])
for (i in seq(1:n)){
lines(timeaxis,plotdata_negative[,i],col=col_negative,lwd=lwd_ind)
}
#points(plotdata[,1],plotdata[,3])
for (i in seq(1:n)){
lines(timeaxis,plotdata_positive[,i],col=col_positive,lwd=lwd_ind)
}
positive_mean_background <- rowMeans(plotdata_positive_background)
negative_mean_background <- rowMeans(plotdata_negative_background)
lines(timeaxis,negative_mean_background, lwd=lwd_means, col=col_negative)
lines(timeaxis,positive_mean_background, lwd=lwd_means, col=col_positive)
naive_likelihood <- positive_mean_background*negative_mean_background
lines(timeaxis,naive_likelihood,col='grey42',lwd=4,lty=5)
true_likelihood <- likelihood_by_DDI(plotdata_positive_background,plotdata_negative_background,timeaxis)
lines(timeaxis,true_likelihood,col=col_truth,lwd=4)
# segments(x0=test_time_1,y0=0,x1=test_time_1,y1=1,lty=4)
segments(x0=test_time_2,y0=0,x1=test_time_2,y1=1,lty=4)
# mean delay arrows (shape package)
#
# Arrows(x0 = timeaxis[which.min(abs(negative_mean_background-0.5))],y0=0.5,x1=test_time_1,y1=0.5,code=3 ,arr.adj=1,arr.width = .1/2)
# text(x=(timeaxis[which.min(abs(negative_mean_background-0.5))]+test_time_1)/2-2.4,y=0.53,pos=4,labels=expression(italic('d')['1']))
#
# Arrows(x0 = timeaxis[which.min(abs(positive_mean_background-0.5))],y0=0.5,x1=test_time_2,y1=0.5,code=3 ,arr.adj=1,arr.width = .1/2)
# text(x=(timeaxis[which.min(abs(positive_mean_background-0.5))]+test_time_2)/2-2.5,y=0.53,pos=4,labels=expression(italic('d')['2']))
# Standard deviation arrows
Arrows(x0 = timeaxis[which.min(abs(negative_mean_background-0.5))]-sd_size_t1,y0=0.5,x1= timeaxis[which.min(abs(negative_mean_background-0.5))], y1=0.5, code=3,arr.adj=1,arr.length = .1,arr.width = .1/2)
text(x=timeaxis[which.min(abs(negative_mean_background-0.5))]-sd_size_t1/2-2.4,y=0.53,pos=4,labels=expression(sigma['1']))
Arrows(x0 = timeaxis[which.min(abs(positive_mean_background-0.5))]-sd_size_t2,y0=0.5,x1= timeaxis[which.min(abs(positive_mean_background-0.5))], y1=0.5, code=3,arr.adj=1,arr.length = .1,arr.width = .1/2)
text(x=timeaxis[which.min(abs(positive_mean_background-0.5))]-sd_size_t2/2-2.5,y=0.53,pos=4,labels=expression(sigma['2']))
# Delta arrow
if(Sys.info()['login']=='jeremyb'){
setwd("C:\\Users\\jeremyb\\Documents\\A_Master\\Code")
}else if(Sys.info()['login']=='jeremy') {
setwd("C:/Users/JumpCo Vostro3700/infection-dating-tool/manuscripts/figures")
}else{
setwd(".") #what does this do?
}
library("DescTools")
library("shape")
#Tool for exploring the effect of incorporating more realistic nuances of test `sensitivies' on prospective
# Initiated March 2018
# Jeremy Bingham
#Definitions
# All test functions [unless specified otherwise] are sensitivity functions ie probability of testing positive as a function of time since exposure
# The next section of code defines functions that are used to:
# choose a set of values for one parameter (usually position), with spacing defined using an inverse cumulative normal distribution.
# define and generate individual curves - a sensitivity curve for each type of test, and a description of how to include the population variability to generate a set of curves
# Generate families of n curves for sensitivity
# Generate families of n curves for likelihood of a certain test result. These functions take a 'test' object which specifies the parameter values
#   (including the set of values for the varying parameter(s) that represent the population-level variability) and the test result ('negative' or 'positive')
#   and a set of times. They return a matrix, each column of which contains the likelihood values for one of the `individual` curves; each curve specifies
#   the likelihood of the test result given for each hypothetical date of first detectable infection.
# Note: generate_family functions don't return the range, granularity etc of thetime axis - this should be defined somewhere more central in the script.
#Tool for generating plots that illustrate the value of the knowing the details of individual HIV detectability
#Initiated February 2018
#Code by Jeremy Bingham with Eduard Grebe and Alex Welte
#Definitions
# All test functions [unless specified otherwise] are sensitivity functions ie probability of testing positive as a function of time since exposure
# The next section of code defines functions that are used to:
# choose a set of values for one parameter (usually position), with spacing defined using an inverse cumulative normal distribution.
# define and generate individual curves - a sensitivity curve for each type of test, and a description of how to include the population variability to generate a set of curves
# Generate families of n curves for sensitivity
# Generate families of n curves for likelihood of a certain test result. These functions take a 'test' object which specifies the parameter values
#   (including the set of values for the varying parameter(s) that represent the population-level variability) and the test result ('negative' or 'positive')
#   and a set of times. They return a matrix, each column of which contains the likelihood values for one of the `individual` curves; each curve specifies
#   the likelihood of the test result given for each hypothetical date of first detectable infection.
# Note: generate_family functions don't return the range, granularity etc of thetime axis - this should be defined somewhere more central in the script.
generate_positions_for_individual_curves = function(n,scale,shape,mean_center_position,sd_size){
shift_to_half_likelihood = scale*(-1*log(1/2))^(1/shape)
myseq <- seq(1/(2*n),1-1/(2*n),1/n)
positions <- qnorm(myseq,mean=mean_center_position-shift_to_half_likelihood,sd=sd_size)
return(positions)
}
generate_positions_cumulative_normal = function(n,mean_center_position,sd_size){
myseq <- seq(1/(2*n),1-1/(2*n),1/n)
positions<-qnorm( myseq,mean=mean_center_position,sd=sd_size )
return(positions)
}
#a sensitivity is a probability of correctly detecting an infection as a function of time after infection (or fddi doesn't matter here)
#a sensitivity is a probability of a positive result as a function of time since infection
#in our graphs (of likelihood functions) the time since infection is the difference between the time we observe and the test-time. This increases as we move left, which is
# why we need to reverse the ``direction of time" when we convert from a sensitivity function to a likelihood. We should call the sensitivity function, reversed (comment this clearly so people aren't confused or angry),
#the probability of testing negative is just 1 -probability of testing positive. no indeterminate results here
# Hello and Welcome
# These comments will guide you gently and clearly throught this code
# if you have any questions don't hesistate to email us at jeremyb@sun.ac.za
#questions I'm wondering about: legends instead of long y-axis labels?
# SECTION ONE: THE SENSITIVITY CURVES
# We begin by defining any and all sensitivity functions we will use. make sure to include an adjustable position parameter
# Each function you define here will have a function of its own in each of the following sections
shift_to_half_likelihood_weibul <- function(scale,shape){
return(scale*(-1*log(1/2))^(1/shape))
}
sensitivity_weibul <- function(x,scale,shape,position){
return(ifelse(x-position<0,0,1-exp(-((x-position)/scale)^shape))) #Note that the zero in the ifelse implies that pre-infection test results are never positive.
#adjust to incorporate imperfect specificity JEREMY CHECK THIS
}
biomarker_weibul <- function(x,scale,shape,position,height){
return(ifelse(x-position<0,0,height*(1-exp(-((x-position)/scale)^shape))))
}
#SECTION TWO: THE LIKELIHOOD CURVES
# we now write functions to represent each sensitivity function as a likelihood of testing positive/negative, given an infection/detectability time.
# the functional-form name at the end of each function name (eg "weibul") refers to the shape of the SENSITIVITY function for a particular test
# position is determined by the delay and the time of the test (0 for sensitivity) otherwise test_time
individual_sensitivity_weibul <- function(times,scale,shape,delay){
return(sensitivity_weibul(x=times,scale=scale,shape=shape,position=delay-shift_to_half_likelihood_weibul(shape=shape,scale=scale)))
}
individual_biomarker_weibul <- function(times,scale,shape,delay,height){ # so delays are delays to half-likelihood
return(biomarker_weibul(x=times,scale=scale,shape=shape,position=delay-shift_to_half_likelihood_weibul(shape=shape,scale=scale),height=height))
}
individual_negative_likelihood_weibul  <- function(times, scale, shape, delay, test_time){ #times is a vector of all the times we consider
return(1-sensitivity_weibul(x=test_time-times,scale=scale,shape=shape,position=delay-shift_to_half_likelihood_weibul(scale=scale,shape=shape)))
}
individual_positive_likelihood_weibul <- function(times, scale, shape, delay, test_time){
return (sensitivity_weibul(x=test_time-times,scale=scale,shape=shape,position=delay-shift_to_half_likelihood_weibul(scale=scale,shape=shape)))
}
# We now have code  generate the likelihood curves for our chosen tests
#SECTION THREE: THE FAMILIES
#note that the factor which varies within the family can be any of the parameters or even a combination of parameters
# currently the family-generating functions have variability driven by one input
family_sensitivity_weibul <- function(times, scale, shape, mean_delay, sd_size, n){ #could also list all parameters (with default values) and use if else to select particular sensitivity shape
list_of_positions <- generate_positions_cumulative_normal(n=n, mean_center_position=mean_delay, sd_size=sd_size)
set_of_sensitivity_curves <- matrix(nrow=length(times),ncol=n)
for(i in seq(1,n)){
set_of_sensitivity_curves[,i] <- individual_sensitivity_weibul(times,scale,shape,list_of_positions[i])
}
return(set_of_sensitivity_curves)
}
family_negative_likelihood_weibul <- function(times, scale, shape, mean_delay, sd_size, n, test_time){ #could also list all parameters (with default values) and use if else to select particular sensitivity shape
list_of_positions <- generate_positions_cumulative_normal(n=n, mean_center_position=mean_delay, sd_size=sd_size)
set_of_negative_curves <- matrix(nrow=length(times),ncol=n)
for(i in seq(1,n)){
set_of_negative_curves[,i] <- individual_negative_likelihood_weibul(times=times,scale=scale,shape=shape,delay = list_of_positions[i],test_time=test_time)
}
return(set_of_negative_curves)
}
family_positive_likelihood_weibul <- function(times, scale, shape, mean_delay, sd_size, n, test_time){
list_of_positions = generate_positions_cumulative_normal(n=n, mean_center_position=mean_delay, sd_size=sd_size)
set_of_positive_curves <- matrix(nrow=length(times),ncol=n)
for(i in seq(1,n)){
set_of_positive_curves[,i] <- individual_positive_likelihood_weibul(times=times,scale=scale,shape=shape,delay = list_of_positions[i],test_time=test_time)
}
return(set_of_positive_curves)
}
family_biomarker_weibul<- function(times,set_of_scales,set_of_shapes,set_of_delays, set_of_heights){ #totally generic function that just
# creates a set of heights. The order of variable combinations, if any, must be determined via the
if(length(set_of_scales)!= length(set_of_shapes)| length(set_of_shapes) != length(set_of_delays) | length(set_of_delays)!= length(set_of_heights)){stop("Some sizes are not lining up!")}
n = length(set_of_scales)
set_of_positive_curves <- matrix(nrow = length(times),ncol=n)
for(i in seq(1,n)){
set_of_positive_curves[,i] <- individual_biomarker_weibul(times=times,scale=set_of_scales[i],shape=set_of_shapes[i],delay=set_of_delays[i],height = set_of_heights[i])
}
return (set_of_positive_curves)
}
n=10
detail=10
timeaxis=seq(0,100,1/detail)
#                       TEST 1 (negative)
#     Describe individual (person) test sensitivity form with population mean-delay and standard deviation of delay
# delay is the variable we distribute across the population - it could in principle be anything else of course
#so each individual has the same SHAPE of sensitivity, but different delays
## Visuals
lwd_means <- 4
lwd_ind <- 1.37
lwd_likelihood <- lwd_means - 3
#High scale causes slower swap
#high shape causes quicker and steeper swap AND more symetrical swap
scale_t1= 4
shape_t1 = 1.73
mean_delay_t1 = 12
sd_size_t1 = 5
#   Time of negative test (relative to arbitrary t=0)
test_time_1 = 28
##                      TEST 2 (positive)
scale_t2 = scale_t1
shape_t2 = shape_t1
mean_delay_t2 = mean_delay_t1
sd_size_t2 = sd_size_t1
#   Time of positive test
test_time_2 = timeaxis[length(timeaxis)]-10
## Generate Data
# Data = the individual likelihood curves for the first (negative) and second (positive) test
#Test 1
plotdata_negative = family_negative_likelihood_weibul(n=n, scale=scale_t1, shape=shape_t1, mean_delay=mean_delay_t1, sd_size= sd_size_t1, times = timeaxis, test_time = test_time_1)
#for generating mean curve
plotdata_negative_background <- family_negative_likelihood_weibul(n=n+50, scale=scale_t1, shape=shape_t1, mean_delay=mean_delay_t1, sd_size= sd_size_t1, times = timeaxis, test_time = test_time_1)
#Test 2
plotdata_positive = family_positive_likelihood_weibul(n=n, scale=scale_t2, shape=shape_t2, mean_delay=mean_delay_t2, sd_size= sd_size_t2, times = timeaxis, test_time = test_time_2)
#for generating mean curve
plotdata_positive_background <- family_positive_likelihood_weibul(n=n+50, scale=scale_t2, shape=shape_t2, mean_delay=mean_delay_t2, sd_size= sd_size_t2, times = timeaxis, test_time = test_time_2)
##        COMMUNICATE
# pdf(file = "figure_2.pdf", width = 7.3, height = 3.7)
plot(timeaxis,plotdata_negative[,1],type='l',xlim=c(timeaxis[1],timeaxis[length(timeaxis)]),ylim=c(-.002,1.002),xaxt='n',yaxt='n',xaxs='i',yaxs='i',bty='l',xlab='',ylab='',col='green') #clarify label in comment
title(xlab="Time of infection", line=1.5, cex.lab=1.2)
title(ylab=expression('Likelihood of test result'), line=2, cex.lab=1.05)
yaxis_pos <- c(0,0.5,1)
yaxis_names <- c('0',"",'1')
xaxis_pos <- c(test_time_1,test_time_2)
xaxis_names <- c(expression('t'['1']*'(-)'),expression('t'['2']*'(+)'))
zero_pos <- c(0)
zero_name <- c(expression('0'['']))
#shift axes
axis(side=2, at=yaxis_pos, labels= yaxis_names,tck=-0.037, padj=.437)
axis(side=1, at=xaxis_pos, labels= xaxis_names,padj=-.35,hadj=-.137)
#axis(side=1, at=zero_pos, labels=zero_name,padj=-0.45,hadj=0.37)
# goto_do
#points(plotdata[,1],plotdata[,3])
for (i in seq(1:n)){
lines(timeaxis,plotdata_negative[,i],col=col_negative,lwd=lwd_ind)
}
#points(plotdata[,1],plotdata[,3])
for (i in seq(1:n)){
lines(timeaxis,plotdata_positive[,i],col=col_positive,lwd=lwd_ind)
}
positive_mean_background <- rowMeans(plotdata_positive_background)
negative_mean_background <- rowMeans(plotdata_negative_background)
lines(timeaxis,negative_mean_background, lwd=lwd_means, col=col_negative)
lines(timeaxis,positive_mean_background, lwd=lwd_means, col=col_positive)
segments(x0=test_time_1,y0=0,x1=test_time_1,y1=1,lty=4)
segments(x0=test_time_2,y0=0,x1=test_time_2,y1=1,lty=4)
# mean delay arrows (shape package)
Arrows(x0 = timeaxis[which.min(abs(negative_mean_background-0.5))],y0=0.5,x1=test_time_1,y1=0.5,code=3 ,arr.adj=1,arr.width = .1/2)
text(x=(timeaxis[which.min(abs(negative_mean_background-0.5))]+test_time_1)/2-2.4,y=0.53,pos=4,labels=expression(italic('d')['1']))
Arrows(x0 = timeaxis[which.min(abs(positive_mean_background-0.5))],y0=0.5,x1=test_time_2,y1=0.5,code=3 ,arr.adj=1,arr.width = .1/2)
text(x=(timeaxis[which.min(abs(positive_mean_background-0.5))]+test_time_2)/2-2.5,y=0.53,pos=4,labels=expression(italic('d')['2']))
# Standard deviation arrows
Arrows(x0 = timeaxis[which.min(abs(negative_mean_background-0.5))]-sd_size_t1,y0=0.5,x1= timeaxis[which.min(abs(negative_mean_background-0.5))], y1=0.5, code=3,arr.adj=1,arr.length = .1,arr.width = .1/2)
text(x=timeaxis[which.min(abs(negative_mean_background-0.5))]-sd_size_t1/2-2.4,y=0.53,pos=4,labels=expression(sigma['1']))
Arrows(x0 = timeaxis[which.min(abs(positive_mean_background-0.5))]-sd_size_t2,y0=0.5,x1= timeaxis[which.min(abs(positive_mean_background-0.5))], y1=0.5, code=3,arr.adj=1,arr.length = .1,arr.width = .1/2)
text(x=timeaxis[which.min(abs(positive_mean_background-0.5))]-sd_size_t2/2-2.5,y=0.53,pos=4,labels=expression(sigma['2']))
# Delta arrow
Arrows(x0 = test_time_1,y0=0.75,x1=test_time_2,y1=0.75, code=3 ,arr.adj=1,arr.width = .1/2)
text(x = test_time_1 + (test_time_2 - test_time_1)/2 - 5, y=0.8,pos=4,labels=expression(delta))
individual_biomarker_weibul <- function(times,scale,shape,delay,height){ # so delays are delays to half-likelihood
return(biomarker_weibul(x=times,scale=scale,shape=shape,position=delay-shift_to_half_likelihood_weibul(shape=shape,scale=scale),height=height))
}
family_biomarker_weibul<- function(times,set_of_scales,set_of_shapes,set_of_delays, set_of_heights){ #totally generic function that just
# creates a set of heights. The order of variable combinations, if any, must be determined via the
if(length(set_of_scales)!= length(set_of_shapes)| length(set_of_shapes) != length(set_of_delays) | length(set_of_delays)!= length(set_of_heights)){stop("Some sizes are not lining up!")}
n = length(set_of_scales)
set_of_positive_curves <- matrix(nrow = length(times),ncol=n)
for(i in seq(1,n)){
set_of_positive_curves[,i] <- individual_biomarker_weibul(times=times,scale=set_of_scales[i],shape=set_of_shapes[i],delay=set_of_delays[i],height = set_of_heights[i])
}
return (set_of_positive_curves)
}
timeaxis <- seq(0,40,1/2)
plotdata_biomarker <- family_biomarker_weibul(timeaxis, set_of_scales = c(1,2,3,4,5,6,7) , set_of_shapes = c(2,3,4,5,6,7,8) , set_of_delays = generate_positions_cumulative_normal(7,10,5) , set_of_heights = seq(0.4,1,.1))
plot(timeaxis,plotdata_biomarker[,1],type = 'n',ylim=c(0,1))
for (i in 1:ncol(plotdata_biomarker)){
print(i)
lines(timeaxis,plotdata_biomarker[,i])
}
lines(timeaxis,proportions_recent(timeaxis,plotdata_biomarker,.4),lty=1,lwd=4,col=2)
n = 5
detail <- 10
timeaxis <- seq(0,90,1/detail)
shape_infect <- 1
scale_infect <- 5
mean_delay_infect <- 7
population_sd_infect <- 10
shape_detect <- 5
scale_detect <- 5
mean_delay_detect <- mean_delay_infect + 14
population_sd_detect <- population_sd_infect
donation_time <- 77
# now we generate and plot the curves
infectious_delays <- generate_positions_cumulative_normal(n=n,mean_center_position = mean_delay_infect, sd_size = population_sd_infect)
detectable_delays <- infectious_delays + 14
infectious_curves <- family_positive_likelihood_weibul_pos(times = timeaxis, shape = shape_infect, scale = scale_infect, positions = infectious_delays, n = n, test_time = donation_time)
undetected_curves <- family_negative_likelihood_weibul_pos(times = timeaxis, shape = shape_detect, scale = scale_detect, positions = detectable_delays, n = n, test_time = donation_time)