phylostep: Stepwise model selection for Phylogenetic Linear Model

In phylolm: Phylogenetic Linear Regression

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)

Example output

Loading required package: ape
----------
Starting model: trait ~ 1 + pred1 + pred2 + pred3
Direction: both
AIC(k=2): -8.4779982199083
	Proposed: trait ~ 1 + pred2 + pred3
	AIC(k=2): 137.180491864307
	Proposed: trait ~ 1 + pred1 + pred3
	AIC(k=2): -10.4039496699306
	Proposed: trait ~ 1 + pred1 + pred2
	AIC(k=2): -9.0442487442954
----------
Step 1
Current model: trait ~ 1 + pred1 + pred3
AIC(k=2): -10.4039496699306
	Proposed: trait ~ 1 + pred3
	AIC(k=2): 136.420559959007
	Proposed: trait ~ 1 + pred1 + pred2 + pred3
	AIC(k=2): -8.4779982199083
	Proposed: trait ~ 1 + pred1
	AIC(k=2): -10.9498137284556
----------
Step 2
Current model: trait ~ 1 + pred1
AIC(k=2): -10.9498137284556
	Proposed: trait ~ 1
	AIC(k=2): 134.477417907675
	Proposed: trait ~ 1 + pred1 + pred2
	AIC(k=2): -9.0442487442954
	Proposed: trait ~ 1 + pred1 + pred3
	AIC(k=2): -10.4039496699306
----------
Step 3
Current model: trait ~ 1 + pred1
AIC(k=2): -10.9498137284556
---END

Call:
phylolm(formula = create.formula(plm), data = data, phy = phy, 
    model = model, lower.bound = lower.bound, upper.bound = upper.bound, 
    starting.value = starting.value)

    AIC  logLik 
-10.950   8.475 

Raw residuals:
    Min      1Q  Median      3Q     Max 
-1.6702  0.4548  0.7314  1.4221  2.3152 

Mean tip height: 1.650307
Parameter estimate(s) using ML:
sigma2: 0.7141216 

Coefficients:
            Estimate   StdErr t.value p.value    
(Intercept) 0.403778 0.552382   0.731  0.4677    
pred1       0.903001 0.036297  24.878  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

phylolm documentation built on July 2, 2020, 3:44 a.m.