phylostep: Stepwise model selection for Phylogenetic Linear Model

In lamho86/phylolm: Phylogenetic Linear Regression

Stepwise model selection for Phylogenetic Linear Model

Description

Performs stepwise model selection for phylogenetic linear models, using the criterion -2*log-likelihood + k*npar, where npar is the number of estimated parameters and k=2 for the usual AIC.

Usage

phylostep(formula, starting.formula = NULL, keeping.formula = NULL, data = list(), 
       phy, model = c("BM", "OUrandomRoot","OUfixedRoot", 
       "lambda", "kappa", "delta", "EB", "trend"),
       direction = c("both", "backward", "forward"), trace = 2,
       lower.bound = NULL, upper.bound = NULL, 
       starting.value = NULL, k=2, ...)

Arguments

formula	formula of the full model.
starting.formula	optional formula of the starting model.
keeping.formula	optional formula of the keeping model. Covariates of the keeping model are always included in the model.
data	a data frame containing variables in the model. If not found in data, the variables are taken from current environment.
phy	a phylogenetic tree of type phylo with branch lengths.
model	a model for the phylogenetic covariance of residuals.
direction	direction for stepwise search, can be both, forward, and backward.
trace	if positive, information on each searching step is printed. Larger values may give more detailed information.
lower.bound	optional lower bound for the optimization of the phylogenetic model parameter.
upper.bound	optional upper bound for the optimization of the phylogenetic model parameter.
starting.value	optional starting value for the optimization of the phylogenetic model parameter.
k	optional weight for the penalty.
...	further arguments to be passed to the function optim.

Details

The default k=2 corresponds to the usual AIC penalty. Use k=\log(n) for the usual BIC, although it is unclear how BIC should be defined for phylogenetic regression.

See phylolm for details on the possible phylogenetic models for the error term, for default bounds on the phylogenetic signal parameters, or for matching tip labels between the tree and the data.

Value

A phylolm object correponding to the best model is returned.

Author(s)

Lam Si Tung Ho and Cecile Ane

See Also

phylolm.

Examples

set.seed(123456)
tre = rcoal(60)
taxa = sort(tre$tip.label)
b0=0; b1=1;
x1 = rTrait(phy=tre,model="BM",
           parameters=list(ancestral.state=0,sigma2=10))
x2 = rTrait(phy=tre,model="BM",
            parameters=list(ancestral.state=0,sigma2=10))
x3 = rTrait(phy=tre,model="BM",
            parameters=list(ancestral.state=0,sigma2=10))
y <- b0 + b1*x1 + 
     rTrait(n=1,phy=tre,model="BM",parameters=list(
              ancestral.state=0,sigma2=1))
dat = data.frame(trait=y[taxa],pred1=x1[taxa],pred2=x2[taxa],pred3=x3[taxa])
fit = phylostep(trait~pred1+pred2+pred3,data=dat,phy=tre,model="BM",direction="both")
summary(fit)

lamho86/phylolm documentation built on Nov. 11, 2023, 6:17 a.m.