Skip to content

Commit 101ccec

Browse files
committed
Create phylolm.R
1 parent 7d50858 commit 101ccec

File tree

1 file changed

+42
-0
lines changed

1 file changed

+42
-0
lines changed
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
library(phylolm)
2+
3+
4+
5+
set.seed(123456)
6+
tre = rcoal(60)
7+
taxa = sort(tre$tip.label)
8+
b0=0; b1=1;
9+
x <- rTrait(n=1, phy=tre,model="BM",
10+
parameters=list(ancestral.state=0,sigma2=10))
11+
y <- b0 + b1*x +
12+
rTrait(n=1,phy=tre,model="lambda",parameters=list(
13+
ancestral.state=0,sigma2=1,lambda=0.5))
14+
dat = data.frame(trait=y[taxa],pred=x[taxa])
15+
fit = phylolm(trait~pred,data=dat,phy=tre,model="lambda")
16+
summary(fit)
17+
18+
# adding measurement errors and bootstrap
19+
z <- y + rnorm(60,0,1)
20+
dat = data.frame(trait=z[taxa],pred=x[taxa])
21+
fit = phylolm(trait~pred,data=dat,phy=tre,model="BM",measurement_error=TRUE,boot=100)
22+
summary(fit)
23+
24+
25+
26+
phyloglm
27+
28+
29+
30+
set.seed(123456)
31+
tre = rtree(50)
32+
x = rTrait(n=1,phy=tre)
33+
X = cbind(rep(1,50),x)
34+
y = rbinTrait(n=1,phy=tre, beta=c(-1,0.5), alpha=1 ,X=X)
35+
dat = data.frame(trait01 = y, predictor = x)
36+
fit = phyloglm(trait01~predictor,phy=tre,data=dat,boot=100)
37+
summary(fit)
38+
coef(fit)
39+
vcov(fit)
40+
41+
predict(fit)
42+
simulate(fit)

0 commit comments

Comments
 (0)