-
Notifications
You must be signed in to change notification settings - Fork 0
/
inferKt_DiamondPrincess.Rev
executable file
·142 lines (115 loc) · 3.98 KB
/
inferKt_DiamondPrincess.Rev
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
############################################################################################
#
# RevBayes Ancestral Population Estimation - Empirical data - COVID-19 aboard the diamond princess cruise ship - Skyline occurrence birth-death process
#
# Model: Tree is drawn from a piecewise constant birth-death process with occurrences.
# We use the tree and parameter traces from an MCMC run to infer the instantaneous epidemiological
# incidence for the Diamond Princess SARS-CoV2 outbreak.
#
#
# authors: Antoine Zwaans.
#
############################################################################################
seed(12345)
#######################
# Reading in the Data #
#######################
# Choose your burnin percentage
burninPercentage <- 0.10
### Read in data
occurrence_ages <- readMatrix(file="data_diamond/data_occurrences_2.csv", delimiter=";")[1]
trace = readTrace("OBDP2_diam_timeline_temp/trace.out", burnin=0.10)
trees = readTreeTrace("OBDP2_diam_timeline_temp/trees.trees", burnin=burninPercentage)
# Trace[] vector with 9 values
# ============================
#
# [1]
# Trace (Iteration)
#
#
# [2]
# Trace (alpha_epi)
#
#
# [3]
# Trace (infectious[1])
#
#
# [4]
# Trace (infectious[5])
#
#
# [5]
# Trace (kappa)
#
#
# [6]
# Trace (obd_tree)
#
#
# [7]
# Trace (pi)
#
#
# [8]
# Trace (psiPlusOmega[1])
#
#
# [9]
# Trace (psiPlusOmega[4])
file_path = "output_diam_timeline_temp_infer/"
nb_trees <- trees.getNumberSamples()
burnin <- trees.getBurnin()
# Fixed OBDP parameters
rho <- 0.0
N <- 40
nb_time_points <- 500
interval_times <- v(9.0,12.0,16.0,22.0,23.0)
start_time <- 38.0
trace_start_time <- rep(start_time,nb_trees)
mu <- 1/20
verb <- TRUE
pm <- v(0,(71 / (257 + 71)),0,0,0,0)
rm <- 1.0
# variable OBDP parameters
trace_psiom_1 = trace[8]
trace_psiom_1.setBurnin(burnin)
trace_psiom_4 = trace[9]
trace_psiom_4.setBurnin(burnin)
trace_lambda_1 = trace[3]
trace_lambda_1.setBurnin(burnin)
trace_lambda_5 = trace[4]
trace_lambda_5.setBurnin(burnin)
# Define the times points at which density is computed according the the oldest start_time
max_age <- 38.0
time_points <- seq(0,38*(1+0.5/(nb_time_points-1)), max_age/(nb_time_points-1))
write("Start ages\n", trace_start_time, "\n", filename=file_path + "start_time_trace.txt", append=FALSE)
write("Kt\n", filename=file_path + "Kt_trace.txt", append=FALSE)
for (i in 1:(nb_trees-burnin)){
print(i)
# Tree: all trees are offset to account for no extant sampling
obd_tree <- trees.getTree(burnin+i)
obd_tree.offset(9.12784147309139)
# Variable OBDP parameters
lambda <- v(trace_lambda_1.getValues()[i],trace_lambda_1.getValues()[i],trace_lambda_1.getValues()[i],trace_lambda_1.getValues()[i],trace_lambda_5.getValues()[i],trace_lambda_5.getValues()[i])
psiom <- v(trace_psiom_1.getValues()[i],trace_psiom_1.getValues()[i],trace_psiom_1.getValues()[i],trace_psiom_4.getValues()[i],trace_psiom_4.getValues()[i],0.0)
psi <- psiom*pm
omega <- abs( psiom - psi )
Kt = fnInferAncestralPopSize( originAge=start_time,
lambda=lambda,
mu=mu,
psi=psi,
omega=omega,
rho=rho,
removalPr=rm,
maxHiddenLin=N,
occurrence_ages=occurrence_ages,
time_points=time_points,
timeTree=obd_tree,
timeline=interval_times,
verbose=verb
)
write(";", Kt, "\n", filename=file_path + "Kt_trace.txt", append=TRUE)
}
#you may want to quit RevBayes now
q()