Uses Rstan and the No U-turn sampler to approximate the growth rate using the likelihood from Stadler 2009 "On incomplete sampling under birth–death models and connections to the sampling-based coalescent"
Usage
birthDeathMCMC(
tree,
maxGrowthRate = 4,
alpha = 0.05,
verbose = TRUE,
nChains = 4,
nCores = 1,
chainLength = 2000
)Arguments
- tree
An ultrametric tree subset to include only the clone of interest. Alternatively, a list with several such trees.
- maxGrowthRate
Sets upper bound on birth rate. Default is 4 but this will depend on the nature of the data
- alpha
Used for calculation of confidence intervals. 1-alpha confidence intervals used with default of alpha = 0.05 (95 percent confidence intervals)
- verbose
TRUE or FALSE, should the Rstan MCMC intermediate output and progress be printed?
- nChains
Number of chains to run in MCMC. Default is 4
- nCores
Number of cores to perform MCMC. Default is 1, but chains can be run in parallel
- chainLength
Number of iterations for each chain in MCMC. Default is 2000 (1000 warm-up + 1000 sampling), increase if stan tells you to
Value
A dataframe including the net growth rate estimate, confidence intervals, and other important details (clone age estimate, runtime, n, etc.)
See also
internalLengths maxLikelihood which use alternative methods for growth rate estimation from an ultrametric tree.
Examples
# \donttest{
df <- birthDeathMCMC(cloneRate::exampleUltraTrees[[1]])
#>
#> CHECKING DATA AND PREPROCESSING FOR MODEL 'bdSampler' NOW.
#>
#> COMPILING MODEL 'bdSampler' NOW.
#>
#> STARTING SAMPLER FOR MODEL 'bdSampler' NOW.
#>
#> SAMPLING FOR MODEL 'bdSampler' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000132 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.32 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: 14.666 seconds (Warm-up)
#> Chain 1: 15.323 seconds (Sampling)
#> Chain 1: 29.989 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'bdSampler' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 0.00011 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.1 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: 16.962 seconds (Warm-up)
#> Chain 2: 16.485 seconds (Sampling)
#> Chain 2: 33.447 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'bdSampler' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 9e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.9 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: 15.03 seconds (Warm-up)
#> Chain 3: 17.164 seconds (Sampling)
#> Chain 3: 32.194 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'bdSampler' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 8.5e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.85 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: 14.397 seconds (Warm-up)
#> Chain 4: 15.744 seconds (Sampling)
#> Chain 4: 30.141 seconds (Total)
#> Chain 4:
# }