Last chance! 50% off unlimited learning
Sale ends in
Model selection with Approximate Bayesian Computation (ABC).
postpr(target, index, sumstat, tol, subset = NULL, method, corr=TRUE,
kernel="epanechnikov", numnet = 10, sizenet = 5, lambda =
c(0.0001,0.001,0.01), trace = TRUE, maxit = 500, ...)An object of class
"postpr", containing the following
components:
a vector of model probabilities when method is
"mnlogistic" or "neuralnet".
the vector of model indices in the accepted region using the rejection method.
vector of regression weights when method is
"mnlogistic" or "neuralnet".
summary statistics in the accepted region.
the original call.
a logical vector indicating the elements or rows that
were excluded, including both NA/NaN's and
elements/rows selected by subset
a character string indicating the method used, i.e.
"rejection", "mnlogistic" or "neuralnet".
logical, if TRUE the posterior model probabilities
are corrected for the number of simulations performed for each
model.
the number of simulations performed for each model a priori.
a character vector of model names (a priori).
the number of summary statistics used.
a list of two elements: model contains the model names,
and statistics.names the names of the summary statistics.
a vector of the observed summary statistics.
a vector of model indices. It can be character or numeric and
will be coerced to factor. It must have the same length as
sumstat to indicate which row of sumstat belong to
which model.
a vector, matrix or data frame of the simulated summary statistics.
numeric, the required proportion of points nearest the target values (tolerance), or a vector of the desired tolerance values. If a vector is given
a logical expression indicating elements or rows to keep. Missing
values in index and/or sumstat are taken as FALSE.
a character string indicating the type of simulation required.
Possible values are "rejection", "mnlogistic",
"neuralnet". See Details.
logical, if TRUE (default) posterior model probabilities are corrected
for the number of simulations performed for each model. If equal
number of simulations are available for all models, corr has
no effect.
a character string specifying the kernel to be used when method is
"mnlogistic" or "neuralnet". Defaults to
"epanechnikov". See density for details.
the number of neural networks when method is "neuralnet". It
corresponds to the number of times the function nnet is called.
the number of units in the hidden layer. Can be zero if there are no skip-layer units.
a numeric vector or a single value indicating the weight decay when
method is "neuralnet". By default, 0.0001, 0.001, or 0.01 is
randomly chosen for the each of the networks. See nnet
for more details.
logical, TRUE switches on tracing the
optimization of nnet (applies when method is
"neuralnet").
numeric, the maximum number of iterations. Defaults to 500. See also
nnet.
other arguments passed on from nnet.
Katalin Csillery, Olivier Francois and Michael Blum with some initial code from Mark Beaumont.
The function computes the posterior model probabilities. Simulations
have to be performed with at least two distinct models. When method is
"rejection", the posterior probability of a given model is
approximated by the proportion of accepted simulations given this
model. This approximation holds when the different models are a priori
equally likely, and the same number of simulations is performed for
each model. When method is "mnlogistic" the posterior model
probabilities are estimated using a multinomial logistic regression as
implemented in the function multinom from the package
nnet. When method is "neuralnet", neural networks
are used to predict the probabilities of models based on the observed
statistics using nnet. This method can be useful if many
summary statistics are used.
Names for the summary statistics are strongly recommended. Names can
be supplied as colnames to sumstat (and target). If no names are
supplied S1, S2, ... to summary statistics will be assigned to
parameters and the user will be warned.
Beaumont, M.A. (2008) Joint determination of topology, divergence time, and immigration in population trees. In Simulation, Genetics, and Human Prehistory (Matsumura, S., Forster, P. and Renfrew, C., eds) McDonald Institute for Archaeological Research
summary.postpr
## an artifical example
ss <- cbind(runif(1000),rt(1000,df=20))
postpr(target=c(3), index=c(rep("norm",500),rep("t",500)),
sumstat=ss[,1], tol=.1, method="rejection")
## human demographic history
require(abc.data)
data(human)
## five R objects are loaded. See ?human and vignette("abc") for details.
## illustrate the informativeness of two summary statistics: mean and
## variance of Tajima's D
par(mfcol = c(1,3))
boxplot(stat.3pops.sim[,"pi"]~models, main="Mean nucleotide diversity")
boxplot(stat.3pops.sim[,"TajD.m"]~models, main="Mean Tajima's D")
boxplot(stat.3pops.sim[,"TajD.v"]~models, main="Var in Tajima's D")
## model selection with ABC for the European population
modsel.it <- postpr(stat.voight["italian",], models, stat.3pops.sim, tol=.05, method="mnlogistic")
summary(modsel.it)
## In Europe, the most supported model
## is a population bottleneck
## Check that in Africa, population expansion is the most supported model, while in
## Asia, it is a population bottleneck
##modsel.ha <- postpr(stat.voight["hausa",], models, stat.3pops.sim,
##tol=.05, method="mnlogistic")
##modsel.ch <- postpr(stat.voight["chinese",], models, stat.3pops.sim,
## tol=.05, method="mnlogistic")
Run the code above in your browser using DataLab