phylostep {phylolm}

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, 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.

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)