R Markdown
This vignette illustrates the analysis of access period datasets with the EpiPvr package. First for semi-persistently transmitted viruses and then for persistently transmitted viruses we use simulated AP datasets to demonstrate virus parameter estimation. This process leads to posterior distributions for viral parameters.
In the EpiPvr package we also include a function for calculating the probability of a local epidemic for a single introduced inoculum. This function requires virus parameters and local parameters. We show how to take samples from the posterior distributions for the virus parameters (as outlined above) and to use them to calculate epidemic probability in local conditions.
In practice, results from these stochastic functions will vary
between runs because of random sampling. We fix a seed
(set.seed(123)) at the start of this vignette to ensure
reproducibility for CRAN checks.
virus_type='SPT'
lsEst_in=40; # upper bound on per insect survival in the experiment in daysLoading the example dataset for an SPT virus. See package documentation entry for the ap_data_sim_SPT dataset for explanation of the dataset columns which correspond to AP experiment structure.
data("ap_data_sim_SPT", package = "EpiPvr")
print(ap_data_sim_SPT)#> $d_AAP
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> T_vec 2 3 3.5 4 4.5 5 6 8
#> R_vec 30 30 30.0 30 30.0 30 30 30
#> I_vec 16 19 20.0 15 18.0 15 16 16
#>
#> $d_IAP
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> T_vec 0.08333333 0.1666667 0.25 0.3333333 0.4166667 0.5 0.6666667
#> R_vec 30.00000000 30.0000000 30.00 30.0000000 30.0000000 30.0 30.0000000
#> I_vec 4.00000000 6.0000000 9.00 13.0000000 13.0000000 8.0 14.0000000
#> [,8] [,9]
#> T_vec 0.8333333 1
#> R_vec 30.0000000 30
#> I_vec 13.0000000 22
#>
#> $d_durations
#> [,1] [,2]
#> [1,] -1 6
#> [2,] 4 -1
#>
#> $d_vectorspp
#> [1] 20
#>
#> $d_virusType
#> [1] "SPT"
#>
#> attr(,"alpha")
#> [1] 0.1
#> attr(,"beta")
#> [1] 1
#> attr(,"mu")
#> [1] 1We now run the estimate_virus_parameters_SPT function for viral parameter estimation from the SPT virus dataset. Note: For runtime reasons, CRAN runs a precomputed fit. To rerun the full analysis, set Sys.setenv(NOT_CRAN = “true”) before rebuilding the vignette.
if (identical(Sys.getenv("NOT_CRAN"), "true")) {
EVSPT_sim=estimate_virus_parameters_SPT(ap_data_sim_SPT,lsEst_in,D_numPtsPdin=1,
mcmcOptions=c(4500,6000),numChainsIn=4,mc.parallel=0)
} else {
# Load precomputed object for CRAN
EVSPT_sim <- readRDS(system.file("extdata", "EVSPT_sim.rds", package = "EpiPvr"))
}#>
#> SAMPLING FOR MODEL 'APmodel_SPT_virus' NOW (CHAIN 1).
#> Chain 1: Rejecting initial value:
#> Chain 1: Error evaluating the log probability at the initial value.
#> Chain 1: Exception: binomial_lpmf: Probability parameter is -9.29352e+12, but must be in the interval [0, 1] (in 'string', line 124, column 2 to column 73)
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000146 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.46 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 6000 [ 0%] (Warmup)
#> Chain 1: Iteration: 600 / 6000 [ 10%] (Warmup)
#> Chain 1: Iteration: 1200 / 6000 [ 20%] (Warmup)
#> Chain 1: Iteration: 1800 / 6000 [ 30%] (Warmup)
#> Chain 1: Iteration: 2400 / 6000 [ 40%] (Warmup)
#> Chain 1: Iteration: 3000 / 6000 [ 50%] (Warmup)
#> Chain 1: Iteration: 3600 / 6000 [ 60%] (Warmup)
#> Chain 1: Iteration: 4200 / 6000 [ 70%] (Warmup)
#> Chain 1: Iteration: 4501 / 6000 [ 75%] (Sampling)
#> Chain 1: Iteration: 5100 / 6000 [ 85%] (Sampling)
#> Chain 1: Iteration: 5700 / 6000 [ 95%] (Sampling)
#> Chain 1: Iteration: 6000 / 6000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 1.408 seconds (Warm-up)
#> Chain 1: 0.514 seconds (Sampling)
#> Chain 1: 1.922 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'APmodel_SPT_virus' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 3.7e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.37 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 6000 [ 0%] (Warmup)
#> Chain 2: Iteration: 600 / 6000 [ 10%] (Warmup)
#> Chain 2: Iteration: 1200 / 6000 [ 20%] (Warmup)
#> Chain 2: Iteration: 1800 / 6000 [ 30%] (Warmup)
#> Chain 2: Iteration: 2400 / 6000 [ 40%] (Warmup)
#> Chain 2: Iteration: 3000 / 6000 [ 50%] (Warmup)
#> Chain 2: Iteration: 3600 / 6000 [ 60%] (Warmup)
#> Chain 2: Iteration: 4200 / 6000 [ 70%] (Warmup)
#> Chain 2: Iteration: 4501 / 6000 [ 75%] (Sampling)
#> Chain 2: Iteration: 5100 / 6000 [ 85%] (Sampling)
#> Chain 2: Iteration: 5700 / 6000 [ 95%] (Sampling)
#> Chain 2: Iteration: 6000 / 6000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 1.544 seconds (Warm-up)
#> Chain 2: 0.463 seconds (Sampling)
#> Chain 2: 2.007 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'APmodel_SPT_virus' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 7.1e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.71 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 6000 [ 0%] (Warmup)
#> Chain 3: Iteration: 600 / 6000 [ 10%] (Warmup)
#> Chain 3: Iteration: 1200 / 6000 [ 20%] (Warmup)
#> Chain 3: Iteration: 1800 / 6000 [ 30%] (Warmup)
#> Chain 3: Iteration: 2400 / 6000 [ 40%] (Warmup)
#> Chain 3: Iteration: 3000 / 6000 [ 50%] (Warmup)
#> Chain 3: Iteration: 3600 / 6000 [ 60%] (Warmup)
#> Chain 3: Iteration: 4200 / 6000 [ 70%] (Warmup)
#> Chain 3: Iteration: 4501 / 6000 [ 75%] (Sampling)
#> Chain 3: Iteration: 5100 / 6000 [ 85%] (Sampling)
#> Chain 3: Iteration: 5700 / 6000 [ 95%] (Sampling)
#> Chain 3: Iteration: 6000 / 6000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 1.537 seconds (Warm-up)
#> Chain 3: 0.615 seconds (Sampling)
#> Chain 3: 2.152 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'APmodel_SPT_virus' NOW (CHAIN 4).
#> Chain 4: Rejecting initial value:
#> Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
#> Chain 4: Stan can't start sampling from this initial value.
#> Chain 4:
#> Chain 4: Gradient evaluation took 3.7e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.37 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 6000 [ 0%] (Warmup)
#> Chain 4: Iteration: 600 / 6000 [ 10%] (Warmup)
#> Chain 4: Iteration: 1200 / 6000 [ 20%] (Warmup)
#> Chain 4: Iteration: 1800 / 6000 [ 30%] (Warmup)
#> Chain 4: Iteration: 2400 / 6000 [ 40%] (Warmup)
#> Chain 4: Iteration: 3000 / 6000 [ 50%] (Warmup)
#> Chain 4: Iteration: 3600 / 6000 [ 60%] (Warmup)
#> Chain 4: Iteration: 4200 / 6000 [ 70%] (Warmup)
#> Chain 4: Iteration: 4501 / 6000 [ 75%] (Sampling)
#> Chain 4: Iteration: 5100 / 6000 [ 85%] (Sampling)
#> Chain 4: Iteration: 5700 / 6000 [ 95%] (Sampling)
#> Chain 4: Iteration: 6000 / 6000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 1.381 seconds (Warm-up)
#> Chain 4: 0.408 seconds (Sampling)
#> Chain 4: 1.789 seconds (Total)
#> Chain 4:target=EVSPT_sim$array1 # fixing the first output array, contains the sample Markov chainsFirst of all, note the important message printed off by default to users of the estimate_virus_parameters_SPT and estimate_virus_parameters_PT functions. It is important that users recognise that they must evaluate the model fit as satisfactory before moving forward with and/or reporting parameter estimates.
In the above, the MCMC sampler sends updates to the console (e.g. chain iteration progress). Warnings relating to the initial value are not a problem - provided that the chain settles down to smoothing successfully.
Printing the dimensions and first parts of each of the seven output objects from estimate_virus_parameters_SPT:
# viewing output content
print(dim(EVSPT_sim$array0))#> [1] 1500 4 5print(head(EVSPT_sim$array0))#> # A draws_array: 6 iterations, 4 chains, and 5 variables
#> , , variable = c1[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.043 0.045 0.042 0.044
#> 2 0.046 0.051 0.041 0.040
#> 3 0.046 0.051 0.047 0.039
#> 4 0.052 0.041 0.042 0.041
#> 5 0.047 0.050 0.044 0.043
#>
#> , , variable = c3[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 1.6 1.5 2.0 2.2
#> 2 1.8 1.5 2.2 2.3
#> 3 1.9 1.5 1.7 2.3
#> 4 1.8 2.6 2.0 2.1
#> 5 1.9 1.3 2.0 1.8
#>
#> , , variable = c2[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 3.09 1.45 1.41 0.71
#> 2 1.14 0.99 1.47 1.51
#> 3 0.62 1.07 0.71 1.33
#> 4 0.92 1.37 1.00 1.09
#> 5 1.51 0.55 1.18 1.02
#>
#> , , variable = bD[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 107 83 59 335
#> 2 308 229 120 325
#> 3 350 168 239 276
#> 4 134 239 138 224
#> 5 106 139 33 329
#>
#> # ... with 1 more iterations, and 1 more variablesprint(dim(EVSPT_sim$array1))#> [1] 1500 4 3print(head(EVSPT_sim$array1))#> # A draws_array: 6 iterations, 4 chains, and 3 variables
#> , , variable = al[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 1.582 0.295 0.151 0.045
#> 2 0.125 0.124 0.148 0.150
#> 3 0.041 0.144 0.054 0.111
#> 4 0.091 0.108 0.077 0.086
#> 5 0.223 0.045 0.109 0.092
#>
#> , , variable = be[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.14 0.32 0.80 1.52
#> 2 0.75 0.59 0.92 0.93
#> 3 1.28 0.56 1.07 1.08
#> 4 0.94 1.38 1.13 1.07
#> 5 0.61 0.83 0.97 0.85
#>
#> , , variable = mu[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 1.50 1.15 1.25 0.66
#> 2 1.02 0.87 1.32 1.36
#> 3 0.58 0.93 0.65 1.22
#> 4 0.83 1.26 0.92 1.01
#> 5 1.29 0.50 1.07 0.93
#>
#> # ... with 1 more iterationsprint(dim(EVSPT_sim$array2))#> [1] 44 10print(head(EVSPT_sim$array2))#> # A tibble: 6 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 c2[1] 1.23e+0 1.15e+0 5.20e-1 4.81e-1 0.549 2.20e+0 1.00 3133.
#> 2 c3[1] 1.95e+0 1.92e+0 3.50e-1 3.43e-1 1.42 2.57e+0 1.00 3696.
#> 3 c1[1] 4.44e-2 4.41e-2 4.24e-3 3.82e-3 0.0382 5.16e-2 1.00 2784.
#> 4 bD[1] 1.79e+2 1.78e+2 1.04e+2 1.34e+2 17.8 3.42e+2 1.00 4328.
#> 5 binPam_succA[1] 5.42e-1 5.45e-1 4.16e-2 4.12e-2 0.469 6.06e-1 1.00 4180.
#> 6 binPam_succA[2] 5.73e-1 5.73e-1 3.10e-2 3.10e-2 0.521 6.22e-1 1.00 4588.
#> # ℹ 1 more variable: ess_tail <dbl>print(length(EVSPT_sim$array3))#> [1] 4print(head(EVSPT_sim$array3))#> $lenA
#> [1] 2.0 3.0 3.5 4.0 4.5 5.0 6.0 8.0
#>
#> $propA
#> [1] 0.5333333 0.6333333 0.6666667 0.5000000 0.6000000 0.5000000 0.5333333
#> [8] 0.5333333
#>
#> $simulLA
#> [1] 0.3333333 0.3666667 0.4000000 0.4000000 0.4000000 0.4000000 0.4000000
#> [8] 0.4000000
#>
#> $simulUA
#> [1] 0.7333333 0.7666667 0.7666667 0.7666667 0.7666667 0.7666667 0.7666667
#> [8] 0.7666667print(length(EVSPT_sim$array4))#> [1] 4print(head(EVSPT_sim$array4))#> $lenI
#> [1] 0.08333333 0.16666667 0.25000000 0.33333333 0.41666667 0.50000000 0.66666667
#> [8] 0.83333333 1.00000000
#>
#> $propI
#> [1] 0.1333333 0.2000000 0.3000000 0.4333333 0.4333333 0.2666667 0.4666667
#> [8] 0.4333333 0.7333333
#>
#> $simulLI
#> [1] 0.03333333 0.06666667 0.13333333 0.16666667 0.20000000 0.23333333 0.26666667
#> [8] 0.30000000 0.33333333
#>
#> $simulUI
#> [1] 0.2666667 0.3666667 0.4666667 0.5333333 0.5666667 0.6000000 0.6666667
#> [8] 0.7000000 0.7000000print(length(EVSPT_sim$array5))#> [1] 2print(head(EVSPT_sim$array5))#> $bayesR2_mn
#> [1] 0.4025417 0.5997492
#>
#> $bayesR2_sd
#> [1] 0.1060072 0.1250199print(length(EVSPT_sim$converge_results))#> [1] 2print(head(EVSPT_sim$converge_results))#> $divergent_transitions
#> [1] 0
#>
#> $max_treedepth_exceeded
#> [1] FALSEruns=dim(target)[1]
nChains=dim(target)[2]estimate_virus_parameters_SPT function output and model fit
The above 7 objects that are returned from the estimate_virus_parameters_SPT function are briefly as follows: array 0 and array 1 are the samples for estimated parameters. Array 0 includes the 3 fitted parameters of the model which are composite parameters and insect mortality as well as lp; array 1 is simply the viral parameters that we are chiefly interested in here (they are derived from the composite parameters). Array 2 contains summary statistics for a range of parameters. It is returned primarily for model diagnostics purposes: the columns rhat and ess_bullk can be checked if needed. Array 3 and array 4 contain forward simulation of the access period data (array 3 is the acquisition varying sub-assay, array 4 is the inoculation varying sub-assay). Array 5 contains Bayesian R2 values which users may wish to consult as measures of model fit (the first and second elements correspond to the acquisition and inoculation varying sub-assays respectively). Finally, converge_results provides values that confirm model convergence (there must be no divergent transitions, and no max tree depth exceeded instances).
NOTE that these are returned to provide flexibility for users to examine model fit - but the absence of warnings that are printed to the console is itself indicative of model convergence. The warnings themselves can be checked with the warnings() command. Again this is not really necessary as rstan convergence warnings will print to console, but rechecking of warnings is here demonstrated for convenience.
An important diagnostic step in Bayesian model fitting is inspecting pairs plots, which are particularly useful for visualising issues such as divergent transitions and correlations - these complement the explicit warnings of the stan sampler however which are the main port of call. These plots display relationships between model parameters and also include lp__, the log-posterior of the model, in addition bD represents the rate of insect mortality in the laboratory,# modelled here as a form of nuisance parameter (can theoretically influence the duration of the final IAP phase in assays - but not generally a focal model component). bD is assigned a uniform prior bounded by the minimum possible survival time (0.1h) and the maximum insect survival under natural conditions, as specified by the user (in days). Given this broad prior and the typically fixed durations of experiments, the posterior distribution of bD is expected to remain largely uninformative.
if (length(warnings()) > 0) {
print("Warnings occurred during fitting; check convergence diagnostics.")
}We use mcmc_pairs from the bayesplot package to perform pair plots of parameters. Not displaying here to reduce runtime/size but can be run by user to visualise the output.
Since we are satisfied with the model fits we proceed with the viral parameter estimates, which we now display alone as density plots. Not displaying here to reduce runtime/size but can be run by user to visualise the output.
We can also plot that part of the function output concerning forward predictions of the assay data (based upon the parameter estimates). One can plot the original assay data (also included in the output) as well as the 95% credible interval. Plots confirm that the 95% credible interval is very good fit for the data.
We now extract the parameter estimates for use in the calculation of local epidemic probability. The posterior distributions now have the form of matrices within a 3D array, third dimension is the parameter (dimensions: iterations x no. Markov-chains). We also prepare the local parameters.
# P EPIDEMIC inferences
# accessing the chains from estimate_virus_parameters
# VIRUS PARAMETERS
# convert from per hours to per day
al_fits=target[,,dimnames(target)$variable=="al[1]", drop = TRUE]*24
be_fits=target[,,dimnames(target)$variable=="be[1]", drop = TRUE]*24
mu_fits=target[,,dimnames(target)$variable=="mu[1]", drop = TRUE]*24
#confirm there are no negative values
print(sum((al_fits<0)))#> [1] 0print(sum((be_fits<0)))#> [1] 0print(sum((mu_fits<0)))#> [1] 0# LOCAL PARAMETERS
# set the local parameters rates (per day)
thet_external <- 0.45 # dispersal
r_external <- 1/28 # roguing
bf_external <- 1/14 # vector mortality
h_external <- 1/365 # harvesting rate
nu_pl_external <- 1/14 # plant latent progression rate (1/plant latent period)
localParams=c(thet_external, r_external, h_external, bf_external, nu_pl_external)
# generate p Epdemic results for given insect burden (per plant)
numInsects=3Typically if one wishes to convert the virus parameter estimates into epidemic probability estimates one would use all of the posterior samples. In this vignette for simplicity however we only look at the last 10 posterior samples (runsPrag=10). Users alternatively may already have a set of virus parameter values they wish to calculate epidemic probability from - in which case they would typically just run the function for a single parameter set.
runsPrag=10
# initial guess for expected time in hours until next event - the function will use this
# to work out efficient time step (optional, default = 0.1)
interval_ind=1/10 #
# initialising vectors for storing chain samples
result_vec_byPl <- array(rep(0, runsPrag*nChains), dim=c(runsPrag,nChains))
result_vec_byIns <- result_vec_byPl
if (identical(Sys.getenv("NOT_CRAN"), "true")) {
# loop through MCMC chains to sample posterior distributions of virus parameters and
# calculate epidemic probability for a given vector abundance
for (ggg in 1:nChains){
print(paste0('nChains',ggg))
for (ppp in 1:runsPrag) {
al_estim=al_fits[nrow(al_fits)+1-ppp,ggg]
be_estim=be_fits[nrow(al_fits)+1-ppp,ggg]
mu_estim=mu_fits[nrow(al_fits)+1-ppp,ggg]
virusParams=c(al_estim,be_estim,mu_estim)
# returns vector of epidemic probabilities for different inoculum states (see ms)
numVars <- ((numInsects+1)*3)-1 # insects per plant number no. inoculum states
qm_out=calculate_epidemic_probability(numInsects,interval_ind,localParams,virusParams)
result_vec_byPl[ppp,ggg]=qm_out[1] # storing epi prob for state: inoculum=single infected plant
result_vec_byIns[ppp,ggg]=qm_out[(numVars-(numInsects-1))] # storing for state: single infected insect
}
}
} else {
# Load precomputed object for CRAN
result_vec_byPl <- readRDS(system.file("extdata", "result_vec_byPl.rds", package = "EpiPvr"))
result_vec_byIns <- readRDS(system.file("extdata", "result_vec_byIns.rds", package = "EpiPvr"))
}#> [1] "nChains1"
#> [1] "nChains2"
#> [1] "nChains3"
#> [1] "nChains4"The output data is itself a posterior distribution - though in this vignette we have only looked at a small number of samples from the virus parameters posterior distributions so results generated should be seen as illustrative only. We can now use as_draws and summarize_draws from the posterior package to extract median values and credible intervals for the epidemic probabilities.
# data wrangling the set of epidemic probabilities from the chains
data_table_Pl=array(rep(0,3), dim=c(3, 1))
data_table_Ins=array(rep(0,3), dim=c(3, 1))
target_A1A2=array(rep(0,2*runsPrag*nChains), dim=c(2,runsPrag*nChains))
# for a given level of insect use 'summarise_draws' to calculate credible intervals
target_A1=array(rep(0,runsPrag*nChains), dim=c(runsPrag, nChains))
target_A2=target_A1
for (gg in 1:nChains){
target_A1[,gg]=t(result_vec_byPl[,gg])
target_A2[,gg]=t(result_vec_byIns[,gg])
}
target_A1A2=rbind(as.vector(target_A1),as.vector(target_A2))
# first the first few elements of epidemic probability from plant inocula
print(head(target_A1A2))#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.48539955 0.4864752 0.3902995 0.4248782 0.4567002 0.4734675 0.4966216
#> [2,] 0.05809464 0.0910801 0.2247828 0.0755253 0.1530937 0.0960898 0.1899132
#> [,8] [,9] [,10] [,11] [,12] [,13] [,14]
#> [1,] 0.3082258 0.4937234 0.3721803 0.4628061 0.4122478 0.3714351 0.49292920
#> [2,] 0.1478117 0.1257100 0.1553661 0.1558168 0.2272802 0.2157636 0.04712676
#> [,15] [,16] [,17] [,18] [,19] [,20] [,21]
#> [1,] 0.4214287 0.4394569 0.4345773 0.4717516 0.3788192 0.4110947 0.4607851
#> [2,] 0.3131646 0.3387703 0.3283832 0.3614309 0.2314139 0.1352918 0.1168496
#> [,22] [,23] [,24] [,25] [,26] [,27] [,28]
#> [1,] 0.4467450 0.4493908 0.3620150 0.4730346 0.4799682 0.4764021 0.40929680
#> [2,] 0.2200593 0.1785746 0.1115594 0.3666578 0.3416713 0.2383932 0.09331562
#> [,29] [,30] [,31] [,32] [,33] [,34] [,35]
#> [1,] 0.4003134 0.3820375 0.4622633 0.47799041 0.4704895 0.4376102 0.3933992
#> [2,] 0.1124377 0.1314315 0.1137368 0.08607046 0.1356931 0.2316856 0.1277673
#> [,36] [,37] [,38] [,39] [,40]
#> [1,] 0.4627149 0.3969122 0.4296629 0.4005608 0.4247246
#> [2,] 0.1837299 0.1785707 0.2067816 0.1479656 0.1462172# and from insect inocula
print(head(target_A1A2))#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.48539955 0.4864752 0.3902995 0.4248782 0.4567002 0.4734675 0.4966216
#> [2,] 0.05809464 0.0910801 0.2247828 0.0755253 0.1530937 0.0960898 0.1899132
#> [,8] [,9] [,10] [,11] [,12] [,13] [,14]
#> [1,] 0.3082258 0.4937234 0.3721803 0.4628061 0.4122478 0.3714351 0.49292920
#> [2,] 0.1478117 0.1257100 0.1553661 0.1558168 0.2272802 0.2157636 0.04712676
#> [,15] [,16] [,17] [,18] [,19] [,20] [,21]
#> [1,] 0.4214287 0.4394569 0.4345773 0.4717516 0.3788192 0.4110947 0.4607851
#> [2,] 0.3131646 0.3387703 0.3283832 0.3614309 0.2314139 0.1352918 0.1168496
#> [,22] [,23] [,24] [,25] [,26] [,27] [,28]
#> [1,] 0.4467450 0.4493908 0.3620150 0.4730346 0.4799682 0.4764021 0.40929680
#> [2,] 0.2200593 0.1785746 0.1115594 0.3666578 0.3416713 0.2383932 0.09331562
#> [,29] [,30] [,31] [,32] [,33] [,34] [,35]
#> [1,] 0.4003134 0.3820375 0.4622633 0.47799041 0.4704895 0.4376102 0.3933992
#> [2,] 0.1124377 0.1314315 0.1137368 0.08607046 0.1356931 0.2316856 0.1277673
#> [,36] [,37] [,38] [,39] [,40]
#> [1,] 0.4627149 0.3969122 0.4296629 0.4005608 0.4247246
#> [2,] 0.1837299 0.1785707 0.2067816 0.1479656 0.1462172Printing median values and and credible intervals for the epidemic probabilities. This is conditioned on a local environment with numInsects=3 insect vectors per host plant, and on the earlier defined local rates.
# ensure at least ~20 data points for each of 4 chains (ie cant summarise distribution if v small)
ts=posterior::summarize_draws(posterior::as_draws(t(target_A1A2))) #> Warning: The ESS has been capped to avoid unstable estimates.data_table_Pl=c(ts[1,]$median,ts[1,]$q5,ts[1,]$q95)
data_table_Ins=c(ts[2,]$median,ts[2,]$q5,ts[2,]$q95)And then printing the summary statistics of local epidemic probability (median and default 90% credible interval from summarize_draws() here for simplicity). Adjust code for e.g. 95% intervals as needed.
names(data_table_Pl)=c('50','5','95')
print(data_table_Pl)#> 50 5 95
#> 0.4385336 0.3709641 0.4929689names(data_table_Ins)=c('50','5','95')
print(data_table_Ins)#> 50 5 95
#> 0.15422988 0.07465377 0.34265926Repeating the entire process for a PT virus dataset.
Loading the example dataset for a PT virus. See package documentation entry for the ap_data_sim_PT dataset for explanation of the dataset columns which correspond to AP experiment structure.
virus_type='PT'
data("ap_data_sim_PT", package = "EpiPvr")
print(ap_data_sim_PT)#> $d_AAP
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> T_vec 1 1.5 1.75 2 2.25 2.5 3 4
#> R_vec 30 30.0 30.00 30 30.00 30.0 30 30
#> I_vec 11 20.0 20.00 23 23.00 25.0 28 30
#>
#> $d_LAP
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> T_vec 0.5 1 2 3 4 5 6 7 8
#> R_vec 30.0 30 30 30 30 30 30 30 30
#> I_vec 29.0 25 25 28 26 28 28 27 26
#>
#> $d_IAP
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> T_vec 0.08333333 0.1666667 0.25 0.3333333 0.4166667 0.5 0.6666667
#> R_vec 30.00000000 30.0000000 30.00 30.0000000 30.0000000 30.0 30.0000000
#> I_vec 4.00000000 12.0000000 9.00 19.0000000 21.0000000 13.0 22.0000000
#> [,8] [,9]
#> T_vec 0.8333333 1
#> R_vec 30.0000000 30
#> I_vec 22.0000000 26
#>
#> $d_durations
#> [,1] [,2] [,3]
#> AAPfixedComponent -1 0.5 1
#> LAPfixedComponent 2 -1.0 1
#> IAPfixedComponent 2 0.5 -1
#>
#> $d_vectorspp
#> [1] 20
#>
#> $d_virusType
#> [1] "PT"
#>
#> attr(,"alpha")
#> [1] 0.1
#> attr(,"beta")
#> [1] 1
#> attr(,"gamma")
#> [1] 0.5
#> attr(,"mu")
#> [1] 0.01Note: for runtime reasons, CRAN runs a precomputed fit. To rerun the full analysis, set Sys.setenv(NOT_CRAN = “true”) before rebuilding the vignette.
if (identical(Sys.getenv("NOT_CRAN"), "true")) {
EVPT_sim=estimate_virus_parameters_PT(ap_data_sim_PT,lsEst_in,D_numPtsPdin=1,mcmcOptions=c(1000,2000),numChainsIn=4,mc.parallel=0)
} else {
# Load precomputed object for CRAN
EVPT_sim <- readRDS(system.file("extdata", "EVPT_sim.rds", package = "EpiPvr"))
}#> numChains = 4#>
#> SAMPLING FOR MODEL 'APmodel_PT_virus' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000134 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.34 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 3.966 seconds (Warm-up)
#> Chain 1: 4.178 seconds (Sampling)
#> Chain 1: 8.144 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'APmodel_PT_virus' NOW (CHAIN 2).
#> Chain 2: Rejecting initial value:
#> Chain 2: Error evaluating the log probability at the initial value.
#> Chain 2: Exception: beta_lpdf: Random variable is 2.50873, but must be in the interval [0, 1] (in 'string', line 205, column 2 to column 25)
#> Chain 2:
#> Chain 2: Gradient evaluation took 9.5e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.95 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 4.293 seconds (Warm-up)
#> Chain 2: 3.259 seconds (Sampling)
#> Chain 2: 7.552 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'APmodel_PT_virus' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 9.8e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 4.156 seconds (Warm-up)
#> Chain 3: 4.921 seconds (Sampling)
#> Chain 3: 9.077 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'APmodel_PT_virus' NOW (CHAIN 4).
#> Chain 4: Rejecting initial value:
#> Chain 4: Error evaluating the log probability at the initial value.
#> Chain 4: Exception: beta_lpdf: Random variable is 6.99809, but must be in the interval [0, 1] (in 'string', line 205, column 2 to column 25)
#> Chain 4:
#> Chain 4: Gradient evaluation took 0.000114 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.14 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 4.464 seconds (Warm-up)
#> Chain 4: 3.471 seconds (Sampling)
#> Chain 4: 7.935 seconds (Total)
#> Chain 4:#> E-BFMI indicated no pathological behavior.#> MESSAGE: Check core convergence diagnostics and model fit!
#> - First: Review RStan sampling messages for warnings (do NOT suppress them).
#> - R-hat < 1.05 for all parameters - see array 2.
#> - Effective Sample Size (ESS) ~>= 100 per chain - see array 2.
#> - Posterior predictive fit: compare simulated vs. observed (arrays 3, 4, 5).
#> - No divergent transitions or treedepth exceedances - see converge_results.
#> More details:
#> Function help & vignette (includes advanced diagnostics)
#> Stan diagnostics guide: https://mc-stan.org/learn-stan/diagnostics-warnings.htmltarget=EVPT_sim$array1
runs=dim(target)[1]
nChains=dim(target)[2]Taking a quick look at estimate_virus_parameters_PT’s output:
print(dim(EVPT_sim$array0))#> [1] 1000 4 6print(head(EVPT_sim$array0))#> # A draws_array: 6 iterations, 4 chains, and 6 variables
#> , , variable = al[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.10 0.11 0.096 0.30
#> 2 0.14 0.09 0.109 0.24
#> 3 0.16 0.13 0.107 0.20
#> 4 0.13 0.13 0.187 0.23
#> 5 0.11 0.27 0.309 0.23
#>
#> , , variable = be[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 1.4 1.54 1.78 0.94
#> 2 1.7 1.37 1.64 0.33
#> 3 1.4 1.00 1.61 0.45
#> 4 1.1 0.82 0.62 0.41
#> 5 1.5 0.40 0.52 0.46
#>
#> , , variable = mu_lat[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.085 0.092 0.081 0.5528
#> 2 0.380 0.079 0.117 0.0015
#> 3 0.302 0.071 0.119 0.0043
#> 4 0.126 0.027 0.095 0.0575
#> 5 0.155 0.073 0.266 0.0291
#>
#> , , variable = lat[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.46 0.38 0.40 0.19
#> 2 0.24 0.52 0.36 0.80
#> 3 0.23 0.41 0.37 0.61
#> 4 0.37 0.43 0.45 0.70
#> 5 0.38 0.52 0.30 0.55
#>
#> # ... with 1 more iterations, and 2 more variablesprint(dim(EVPT_sim$array1))#> [1] 1000 4 4print(head(EVPT_sim$array1))#> # A draws_array: 6 iterations, 4 chains, and 4 variables
#> , , variable = al[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.10 0.11 0.096 0.30
#> 2 0.14 0.09 0.109 0.24
#> 3 0.16 0.13 0.107 0.20
#> 4 0.13 0.13 0.187 0.23
#> 5 0.11 0.27 0.309 0.23
#>
#> , , variable = be[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 1.4 1.54 1.78 0.94
#> 2 1.7 1.37 1.64 0.33
#> 3 1.4 1.00 1.61 0.45
#> 4 1.1 0.82 0.62 0.41
#> 5 1.5 0.40 0.52 0.46
#>
#> , , variable = mu[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.039 0.035 0.032 0.1048
#> 2 0.090 0.041 0.043 0.0012
#> 3 0.069 0.029 0.044 0.0026
#> 4 0.047 0.011 0.042 0.0404
#> 5 0.059 0.038 0.080 0.0161
#>
#> , , variable = lat[1]
#>
#> chain
#> iteration 1 2 3 4
#> 1 0.46 0.38 0.40 0.19
#> 2 0.24 0.52 0.36 0.80
#> 3 0.23 0.41 0.37 0.61
#> 4 0.37 0.43 0.45 0.70
#> 5 0.38 0.52 0.30 0.55
#>
#> # ... with 1 more iterationsprint(dim(EVPT_sim$array2))#> [1] 63 10print(head(EVPT_sim$array2))#> # A tibble: 6 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 lat[1] 0.406 0.381 0.143 0.129 0.226 0.670 1.00 1569.
#> 2 bD[1] 179. 178. 104. 132. 15.4 340. 1.00 2280.
#> 3 al[1] 0.135 0.120 0.0634 0.0362 0.0784 0.245 1.00 1115.
#> 4 be[1] 1.31 1.28 0.503 0.507 0.549 2.18 1.00 1513.
#> 5 mu_lat[1] 0.142 0.101 0.132 0.108 0.00667 0.415 1.00 1110.
#> 6 mu[1] 0.0446 0.0392 0.0311 0.0334 0.00359 0.102 1.00 1150.
#> # ℹ 1 more variable: ess_tail <dbl>print(length(EVPT_sim$array3))#> [1] 4print(head(EVPT_sim$array3))#> $lenA
#> [1] 1.00 1.50 1.75 2.00 2.25 2.50 3.00 4.00
#>
#> $propA
#> [1] 0.3666667 0.6666667 0.6666667 0.7666667 0.7666667 0.8333333 0.9333333
#> [8] 1.0000000
#>
#> $simulLA
#> [1] 0.3333333 0.5000000 0.5666667 0.6333333 0.7000000 0.7333333 0.8000000
#> [8] 0.9000000
#>
#> $simulUA
#> [1] 0.7000000 0.8333333 0.9000000 0.9333333 0.9666667 0.9666667 1.0000000
#> [8] 1.0000000print(length(EVPT_sim$array4))#> [1] 4print(head(EVPT_sim$array4))#> $lenL
#> [1] 0.5 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
#>
#> $propL
#> [1] 0.9666667 0.8333333 0.8333333 0.9333333 0.8666667 0.9333333 0.9333333
#> [8] 0.9000000 0.8666667
#>
#> $simulLL
#> [1] 0.6333333 0.7000000 0.7333333 0.7666667 0.7666667 0.8000000 0.7666667
#> [8] 0.7666667 0.7666667
#>
#> $simulUL
#> [1] 0.9333333 0.9666667 0.9666667 1.0000000 1.0000000 1.0000000 1.0000000
#> [8] 1.0000000 1.0000000print(length(EVPT_sim$array5))#> [1] 4print(head(EVPT_sim$array5))#> $lenI
#> [1] 0.08333333 0.16666667 0.25000000 0.33333333 0.41666667 0.50000000 0.66666667
#> [8] 0.83333333 1.00000000
#>
#> $propI
#> [1] 0.1333333 0.4000000 0.3000000 0.6333333 0.7000000 0.4333333 0.7333333
#> [8] 0.7333333 0.8666667
#>
#> $simulLI
#> [1] 0.03333333 0.13333333 0.23333333 0.30000000 0.36666667 0.40000000 0.50000000
#> [8] 0.60000000 0.63333333
#>
#> $simulUI
#> [1] 0.3000000 0.4666667 0.6000000 0.6666667 0.7333333 0.7666667 0.8666667
#> [8] 0.9000000 0.9333333print(length(EVPT_sim$array6))#> [1] 2print(head(EVPT_sim$array6))#> $bayesR2_mn
#> [1] 0.7638187 0.3892205 0.7316994
#>
#> $bayesR2_sd
#> [1] 0.13317879 0.08920427 0.09085339print(length(EVPT_sim$converge_results))#> [1] 2print(head(EVPT_sim$converge_results))#> $divergent_transitions
#> [1] 0
#>
#> $max_treedepth_exceeded
#> [1] FALSEAnd indeed seeing if there were any warnings - NB may not be essential persé - as rstan convergence warnings will print to console - rechecking of warnings, however, demonstrated here for reference are crucial and users must make sure that they have assessed convergence diagnostics and model fit.
if (length(warnings()) > 0) {
print("Warnings occurred during fitting; check convergence diagnostics.")
}Displaying pairs plots for the parameter estimates of a PT virus dataset. Not displaying here to reduce runtime/size but can be run by user to visualise the output.
Since we are satisfied with the model fits we proceed with the viral parameter estimates, which we now display alone as density plots. Not displaying here to reduce runtime/size but can be run by user to visualise the output.
Plotting forward simulated assay data based on the parameter estimates against the original assay data.
Extracting the Markov chain samples for the estimate virus parameters, and calcuating epidemic probability based on a small number of samples (*for illustration only).
# P EPIDEMIC inferences
# accessing the chains from estimate_virus_parameters
# VIRUS PARAMETERS
# convert from per min to per day
al_fits=target[,,dimnames(target)$variable=="al[1]", drop = TRUE]*24
be_fits=target[,,dimnames(target)$variable=="be[1]", drop = TRUE]*24
mu_fits=target[,,dimnames(target)$variable=="mu[1]", drop = TRUE]*24
print(sum((al_fits<0)))#> [1] 0print(sum((be_fits<0)))#> [1] 0print(sum((mu_fits<0)))#> [1] 0# LOCAL PARAMETERS
# set the local parameters (per day)
thet_external <- 0.45 # dispersal
r_external <- 1/28 # roguing
bf_external <- 1/14 # vector mortality
h_external <- 1/365 # harvesting rate
nu_pl_external <- 1/14 # plant latent progression rate (1/plant latent period)
localParams=c(thet_external, r_external, h_external, bf_external, nu_pl_external)
# generate p Epdemic results for various insect burdens (per plant)
interval_ind=(1/10) # this is an optional start setting relating to efficiency of time-step
result_vec_byPl_PT <- array(rep(0, runsPrag*nChains), dim=c(runsPrag,nChains))
result_vec_byIns_PT <- result_vec_byPl
if (identical(Sys.getenv("NOT_CRAN"), "true")) {
# loop through MCMC chains to sample posterior distributions of virus parameters
# and calc epi prob for given vector abundance
for (ggg in 1:nChains){
print(paste0('nChains',ggg))
for (ppp in 1:runsPrag) {
al_estim=al_fits[nrow(al_fits)+1-ppp,ggg]
be_estim=be_fits[nrow(al_fits)+1-ppp,ggg]
mu_estim=mu_fits[nrow(al_fits)+1-ppp,ggg]
virusParams=c(al_estim,be_estim,mu_estim)
# returns vector of epidemic probabilities for different inoculum states (see ms)
numVars <- ((numInsects+1)*3)-1 # insects per plant number no. inoculum states
qm_out=calculate_epidemic_probability(numInsects,interval_ind,localParams,virusParams)
result_vec_byPl_PT[ppp,ggg]=qm_out[1] # storing epi prob for state: inoculum=single infected plant
result_vec_byIns_PT[ppp,ggg]=qm_out[(numVars-(numInsects-1))] # storingfor state: single infected insect
}
}
} else {
# Load precomputed object for CRAN
result_vec_byPl_PT <- readRDS(system.file("extdata", "result_vec_byPl_PT.rds", package = "EpiPvr"))
result_vec_byIns_PT <- readRDS(system.file("extdata", "result_vec_byIns_PT.rds", package = "EpiPvr"))
}#> [1] "nChains1"
#> [1] "nChains2"
#> [1] "nChains3"
#> [1] "nChains4"# now run the epidemic probability calculator
# runtime for a single parameter set is fast
# runtime for a MCMC chains (potentially very many parameter sets) may take substantial time
# data wrangling the set of epidemic probabilities from the chains
data_table_Pl=array(rep(0,3), dim=c(3, 1))
data_table_Ins=array(rep(0,3), dim=c(3, 1))
target_A1A2=array(rep(0,2*runsPrag*nChains), dim=c(2,runsPrag*nChains))
target_A1=array(rep(0,runsPrag*nChains), dim=c(runsPrag, nChains))
target_A2=target_A1
for (gg in 1:nChains){
target_A1[,gg]=t(result_vec_byPl_PT[,gg])
target_A2[,gg]=t(result_vec_byIns_PT[,gg])
}
target_A1A2=rbind(as.vector(target_A1),as.vector(target_A2))
# ensure at least ~20 data points for each of 4 chains
ts=posterior::summarize_draws(posterior::as_draws(t(target_A1A2)))
data_table_Pl=c(ts[1,]$median,ts[1,]$q5,ts[1,]$q95)
data_table_Ins=c(ts[2,]$median,ts[2,]$q5,ts[2,]$q95)
print(head(target_A1A2[1,]))#> [1] 0.9922785 0.9903904 0.9924311 0.9924553 0.9899436 0.9908326print(head(target_A1A2[2,]))#> [1] 0.9169592 0.8627879 0.9209993 0.9219654 0.9388406 0.9462078# nb default summarize_draws 90% interval - please adjust for 95%
names(data_table_Pl)=c('50','5','95')
print(data_table_Pl)#> 50 5 95
#> 0.9908183 0.9853185 0.9937999names(data_table_Ins)=c('50','5','95')
print(data_table_Ins)#> 50 5 95
#> 0.9258639 0.8703072 0.9534720