Fits the phylogenetic logistic regression described in Ives and Garland (2010) and the Poisson regression described in Paradis and Claude (2002). The computation uses an algorithm that is linear in the number of tips in the tree.
phyloglm(formula, data, phy, method = c("logistic_MPLE","logistic_IG10", "logistic_MLE",
"poisson_GEE"), btol = 10, log.alpha.bound = 4,
start.beta=NULL, start.alpha=NULL,
boot = 0, full.matrix = TRUE, save = FALSE)
the named vector of coefficients.
(logistic regression only) the phylogenetic correlation parameter.
(Poisson regression only) the scale parameter which estimates the overdispersion.
standard deviation for the regression coefficients.
covariance matrix for the regression coefficients.
(logistic regression only) log likelihood.
(logistic regression only) AIC.
(logistic regression only) penalized log likelihood, using Firth's penalty for coefficients.
response.
number of observations (tips in the tree).
number of dependent variables.
the model formula.
the original call to the function.
the estimation method.
An integer code with '0' for successful optimization. With logistic_MPLE, this is the convergence code from the optim routine.
(logistic regression only) An interger code with '0' for the estimate of alpha is not near the lower and upper bounds, code with '1' for the estimate of alpha near the lower bound, code with '2' for the estimate of alpha near the upper bound.
a design matrix with one row for each tip in the phylogenetic tree.
(boot > 0 only) bootstrap means of the parameters estimated.
(boot > 0 only) bootstrap standard deviations of the estimated parameters.
(boot > 0 only) bootstrap 95% confidence interval.
(boot > 0 only) bootstrap mean of the logs of the estimated alphas.
(boot > 0 only) bootstrap standard deviation of the logs of the estimated alphas.
(boot > 0 only) number of independent bootstrap replicates for which
phyloglm failed. These failures may be due to the bootstrap data having too few 0's or too few 1's.
(boot > 0 and full.matrix = TRUE only) matrix of all bootstrap estimates.
(boot > 0 and save = TRUE only) matrix of all bootstrap data (each column is a dataset).
a model formula.
a data frame containing variables in the model. If not
found in data, the variables are taken from the current environment.
a phylogenetic tree of type phylo with branch lengths.
The "logistic_IG10" method optimizes a GEE approximation to the penalized likelihood of the logistic regression. "logistic_MPLE" maximizes the penalized likelihood of the logistic regression. In both cases, the penalty is Firth's correction. The "poisson_GEE" method solves the generalized estimating equations (GEE) for Poisson regression.
(logistic regression only) bound on the linear predictor to bound the searching space.
(logistic regression only) bound for the log of the parameter alpha.
starting values for beta coefficients.
(logistic regression only) starting values for alpha (phylogenetic correlation).
number of independent bootstrap replicates, 0 means no bootstrap.
if TRUE, the full matrix of bootstrap estimates (coefficients and alpha) will be returned.
if TRUE, the simulated bootstrap data will be returned.
Lam Si Tung Ho, Robert Lachlan, Rachel Feldman and Cecile Ane
This function uses an algorithm that is linear in the number of tips in the tree.
Bootstrapping can be parallelized using the future package on any future
compatible back-end. For example, run library(future); plan(multiprocess)),
after which bootstrapping will automatically occur in parallel. See
plan for options.
Ho, L. S. T. and Ane, C. 2014. "A linear-time algorithm for Gaussian and non-Gaussian trait evolution models". Systematic Biology 63(3):397-408.
Ives, A. R. and T. Garland, Jr. 2010. "Phylogenetic logistic regression for binary dependent variables". Systematic Biology 59:9-26.
Paradis E. and Claude J. 2002. "Analysis of Comparative Data Using Generalized Estimating Equations". Journal of Theoretical Biology 218:175-185.
compar.gee.
set.seed(123456)
tre = rtree(50)
x = rTrait(n=1,phy=tre)
X = cbind(rep(1,50),x)
y = rbinTrait(n=1,phy=tre, beta=c(-1,0.5), alpha=1 ,X=X)
dat = data.frame(trait01 = y, predictor = x)
fit = phyloglm(trait01~predictor,phy=tre,data=dat,boot=100)
summary(fit)
coef(fit)
vcov(fit)
Run the code above in your browser using DataLab