Description Usage Arguments Details Value Author(s) References See Also Examples
View source: R/expoTree.optim.R
Perform optimization to find a maximum likelihood estimate.
1 2 3 |
forest |
List of trees in two-column format. Column 1 are the event times (branching or sampling) and column 2 is the event type (0 = sampling, 1 = branching). |
fix |
Logical vector specifying which parameters to keep constant. |
fix.val |
Values for fixed parameters. Also specify dummy values for variable parameters. These are ignored. |
pars |
Starting values for the parameters. |
lo |
Lower bound for parameters. |
hi |
Upper bound for parameters. |
vflag |
Set verbosity level. |
survival |
Condition on the likelihood of observing the tree. |
method |
Choose optimization method. Passed on to 'optim'. See 'optim' for details. |
control |
Control parameters for the optimization method. See 'optim' for details. |
Special case: It is possible to fix the ratio between mu and psi. If mu is fixed at a negtive value, it is seen as a fixed ratio. Example: r <- 0.5 fix <- c(F,F,T,F,T) fix.val <- c(0,0,-r,0,0) At every evaluation, mu is calculated accordingly: psi/(psi+mu) = r ==> mu = psi*(1/r-1)
Output from the optimization. List with the following entries: par : parameter estimates value : negative log-likelihood For the other returned values, see help(optim) respectively
Gabriel E. Leventhal
Leventhal, Guenthard, Bonhoeffer & Stadler, 2013
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 | # simulate trees
N <- 15
beta <- 1
mu <- 0.1
psi <- 0.1
rho <- 0
nsamp <- 20
epi <- sim.epi(N,beta,mu,psi,nsamp)
tree <- cbind(epi$times,epi$ttypes)
extant <- sum(2*tree[,2]-1)
lineages <- sum(2*tree[,2]-1)+cumsum(1-2*tree[,2])
max.lineages <- max(lineages)
# calculate likelihood for the forest
lik <- runExpoTree(pars=c(N,beta,mu,psi,rho),times=tree[,1],ttypes=tree[,2],
survival=TRUE,return.full=FALSE)
cat("Likelihood = ",lik,"\n")
if (! any(is.nan(lik))) {
expoTree.optim(list(tree),pars=c(N,beta,psi),
lo=c(max(max.lineages),0,0),hi=c(50,3,1),
fix=c(FALSE,FALSE,TRUE,FALSE,TRUE),
fix.val=c(0,0,-psi/(psi+mu),0,0))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.