parametric.bootstrap  R Documentation 
The function performs a parametric bootstrap for confidence intervals for estimates of the evolutionary model. The user may specify what parameters are to have their confidence intervals returned. The user is recommended to install the suggested package PCMBaseCpp which significantly speeds up the calculations (see Details).
parametric.bootstrap(estimated.model, phyltree,
values.to.bootstrap = NULL, regimes = NULL,
root.regime = NULL, M.error = NULL, predictors = NULL,
kY = NULL, numboot = 100, Atype = NULL, Syytype = NULL,
diagA = NULL, parameter_signs = NULL, start_point_for_optim = NULL,
parscale = NULL, min_bl = 0.0003, maxiter = c(10,50,100), estimateBmethod="ML")
estimated.model 
An estimated by evolutionary model. It can be e.g. the output of 
phyltree 
The phylogeny in 
values.to.bootstrap 
A vector of parameter/composite statistic names that the user is interested in. They are extracted from the bootstrapped elements for easy access. 
regimes 
A vector or list of regimes. If vector then each entry corresponds to each of 
root.regime 
The regime at the root of the tree. If not given, then it is taken as the regime that is present
on the root's daughter lineages and is the most frequent one in the 
M.error 
An optional measurement error covariance structure. The measurement errors between species are assumed independent. The program tries to recognize the structure of the passed matrix and accepts the following possibilities :
From version 2.0.0 of mvSLOUCH it is impossible to pass a single joint measurement error matrix for all the species and traits. 
predictors 
A vector giving the numbers of the columns from the original data which are to be considered predictor ones, i.e. conditioned on in the program output. If not provided then the "X" variables are treated as predictors, but this only for the OUBM models (for the others in this case none are treated as predictors). 
kY 
Number of "Y" (response) variables, for the OUBM models.
The first 
numboot 
The number of bootstraps to perform. 
Atype 
The class of the 
Syytype 
The class of the Syy matrix, ignored if 
diagA 
Should the diagonal of 
parameter_signs 
WARNING: ONLY use this option if you understand what you are doing! This option
is still in an experimental stage so some setups might not work (please report).
A list allowing the user to control whether specific entries for each model parameter
should be positive, negative, zero or set to a specific (other) value. The entries
of the list have to be named, the admissible names are 
start_point_for_optim 
A named list with starting parameters for of the parameters for be optimized by start_point_for_optim=list(A=rbind(c(2,0),(0,4)), Syy=rbind(c(1,0.5),c(0,2))). This starting point is always jittered in each bootstrap replicate as the employed

parscale 
A vector to calculate the 
min_bl 
Value to which PCMBase's 
maxiter 
The maximum number of iterations for different components of the estimation
algorithm. A vector of three integers. The first is the number of iterations for phylogenetic
GLS evaluations, i.e. conditional on the other parameters, the regime optima, perhaps 
estimateBmethod 
Only relevant for OUBM models, should 
The likelihood calculations are done by the PCMBase package. However, there is a C++ backend, PCMBaseCpp. If it is not available, then the likelihood is calculated slower using pure R. However, with the calculations in C++ up to a 100fold increase in speed is possible (more realistically 1020 times). The PCMBaseCpp package is available from https://github.com/venelin/PCMBaseCpp.
The setting Atype="Any"
means that one assumes the matrix A
is eigendecomposable.
If the estimation algorithm hits a defective A
, then it sets the loglikelihood at
the minimum value and will try to get out of this dip.
A list with all the bootstrap simulations is returned. The elements of the list are the following.
paramatric.bootstrap.estimation.replicates 
A list of length
equalling 
bootstrapped.parameters 
If 
The estimation can take a long time and hence many bootstrap replicates will take even more time.The code can produce (a lot of) warnings and errors during the search procedure, this is nothing to worry about.
The ouch package implements a parametric bootstrap and reading about it could be helpful.
Krzysztof Bartoszek
Bartoszek, K. and Pienaar, J. and Mostad. P. and Andersson, S. and Hansen, T. F. (2012) A phylogenetic comparative method for studying multivariate adaptation. Journal of Theoretical Biology 314:204215.
Butler, M.A. and A.A. King (2004) Phylogenetic comparative analysis: a modeling approach for adaptive evolution. American Naturalist 164:683695.
BrownianMotionModel
, estimate.evolutionary.model
,
mvslouchModel
, ouchModel
, bootstrap
,
optim
RNGversion(min(as.character(getRversion()),"3.6.1"))
set.seed(12345, kind = "MersenneTwister", normal.kind = "Inversion")
### We will first simulate a small phylogenetic tree using functions from ape.
### For simulating the tree one could also use alternative functions, e.g. sim.bd.taxa
### from the TreeSim package
phyltree<ape::rtree(5)
## The line below is not necessary but advisable for speed
phyltree<phyltree_paths(phyltree)
BMparameters<list(vX0=matrix(0,nrow=3,ncol=1),
Sxx=rbind(c(1,0,0),c(0.2,1,0),c(0.3,0.25,1)))
### Now simulate the data.
BMdata<simulBMProcPhylTree(phyltree,X0=BMparameters$vX0,Sigma=BMparameters$Sxx)
BMdata<BMdata[phyltree$tip.label,,drop=FALSE]
### Recover the parameters of the Brownian motion.
BMestim<BrownianMotionModel(phyltree,BMdata)
### And finally obtain bootstrap confidence intervals for some parameters
BMbootstrap<parametric.bootstrap(estimated.model=BMestim,phyltree=phyltree,
values.to.bootstrap=c("vX0","StS"),M.error=NULL,numboot=2)
RNGversion(as.character(getRversion()))
## Not run: ##It takes too long to run this
### Define a vector of regimes.
regimes<c("small","small","large","small","small","large","large","large")
### Define SDE parameters to be able to simulate data under the mvOUBM model.
OUBMparameters<list(vY0=matrix(c(1,1),ncol=1,nrow=2),A=rbind(c(9,0),c(0,5)),
B=matrix(c(2,2),ncol=1,nrow=2),mPsi=cbind("small"=c(1,1),"large"=c(1,1)),
Syy=rbind(c(1,0.25),c(0,1)),vX0=matrix(0,1,1),Sxx=matrix(1,1,1),
Syx=matrix(0,ncol=1,nrow=2),Sxy=matrix(0,ncol=2,nrow=1))
### Now simulate the data.
OUBMdata<simulMVSLOUCHProcPhylTree(phyltree,OUBMparameters,regimes,NULL)
OUBMdata<OUBMdata[phyltree$tip.label,,drop=FALSE]
### Try to recover the parameters of the mvOUBM model.
OUBMestim<mvslouchModel(phyltree,OUBMdata,2,regimes,Atype="DecomposablePositive",
Syytype="UpperTri",diagA="Positive",maxiter=c(10,50,100))
### And finally bootstrap with particular interest in the evolutionary and optimal
### regressions
OUBMbootstrap<parametric.bootstrap(estimated.model=OUBMestim,phyltree=phyltree,
values.to.bootstrap=c("evolutionary.regression","optimal.regression"),
regimes=regimes,root.regime="small",M.error=NULL,predictors=c(3),kY=2,
numboot=5,Atype="DecomposablePositive",Syytype="UpperTri",diagA="Positive",
maxiter=c(10,50,100))
### We now demonstrate an alternative setup
### Define SDE parameters to be able to simulate data under the OUOU model.
OUOUparameters<list(vY0=matrix(c(1,1,0.5),nrow=3,ncol=1),
A=rbind(c(9,0,0),c(0,5,0),c(0,0,1)),mPsi=cbind("small"=c(1,1,0.5),"large"=c(1,1,0.5)),
Syy=rbind(c(1,0.25,0.3),c(0,1,0.2),c(0,0,1)))
### Now simulate the data.
OUOUdata<simulOUCHProcPhylTree(phyltree,OUOUparameters,regimes,NULL)
OUOUdata<OUOUdata[phyltree$tip.label,,drop=FALSE]
### Try to recover the parameters of the OUOU model.
estimResults<estimate.evolutionary.model(phyltree,OUOUdata,regimes=regimes,
root.regime="small",M.error=NULL,repeats=3,model.setups=NULL,predictors=c(3),kY=2,
doPrint=TRUE,pESS=NULL,maxiter=c(10,50,100))
### And finally bootstrap with particular interest in the evolutionary regression
OUOUbootstrap<parametric.bootstrap(estimated.model=estimResults,phyltree=phyltree,
values.to.bootstrap=c("evolutionary.regression"),
regimes=regimes,root.regime="small",M.error=NULL,predictors=c(3),kY=NULL,
numboot=5,Atype=NULL,Syytype=NULL,diagA=NULL)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.