Charles Pouchon
2018-03-20 16:50:40 UTC
Hi everybody,
I'm sending this email because I did (or tried to do) a QuaSSE analysis
from the diversitree package to understand the impact of my traits values
into speciation rate. However, I am not very familar with such analysis and
with the tutorial script provided, and I would like to know if somebody
could help me understand some functions and results.
My trait vector corresponds to the positions of species into the first
component of a PCA made on morphological traits:
*>trait1 [1] 1.53635718 0.11435966 0.40744255 -0.16969599 0.83903277
0.27982589 0.34246181 1.40108911 0.92971423 0.58037380 1.01094382
0.48161104[13] 0.91112199 1.19083790 1.21212025 0.73398302 0.93707017
0.01386014 1.16228321 0.21442138 1.02680456 1.07615452 -1.13881116
1.10192106[25] 1.13232218 1.38341060 1.35456782 -3.81104255 -3.51734609
-3.73474404 -2.75423989 -3.96242786 -3.61512213 -3.98707154 -3.26697792
1.06852034[37] -3.69987363 -3.91920521 1.10428403 -0.25527186 1.25529456
1.41681721 1.27082040 -4.76293180 -1.49813425 -0.01049742 -1.49607314
-0.07286070[49] -7.63586519 -8.65295750*
*>range(trait1)*
*[1] -8.652957 1.536357*
*>f.c <- make.quasse(tree, trait1, states.sd <http://states.sd/>=sd(trait1)
constant.x, constant.x, control=control) >p <- starting.point.quasse(tree,
trait1)>f.c.nodrift <- constrain(f.c, drift~0)>mle.c <-
find.mle(f.c.nodrift, p, control=control, lower=0, verbose=0)*
Next, I would like to make linear, sigmoidal and hump-shaped functions with
and without drift. In order to make this, I used this script:
*>xr <- range(trait1)*
*>p.c <- mle.c$par>p.l <- c(p.c[1], l.m=0, p.c[2:3])>p.s <- p.h <-
c(p.c[1], p.c[1], mean(xr), 1, p.c[2:3])*
*>linear.x <- make.linear.x(xr[1],xr[2])>f.l <- make.quasse(tree, trait1,
states.sd <http://states.sd/>=Axis1M.sd, linear.x, constant.x,
control=control)>f.l.nodrift <- constrain(f.l, drift~0)>mle.l <-
find.mle(f.l.nodrift, p.l, control=control, lower=0, verbose=0)>mle.l.drift
<- find.mle(f.l, coef(mle.l, TRUE), control=control, verbose=0) >f.s <-
make.quasse(tree, trait1, states.sd <http://states.sd/>=Axis1M.sd,
sigmoid.x, constant.x, control=control)>f.s.nodrift <- constrain(f.s,
drift~0)>mle.s <- find.mle(f.s.nodrift, p.s, control=control,
verbose=0)>mle.s.drift <- find.mle(f.s, coef(mle.s, TRUE), control=control,
verbose=0) >f.h <- make.quasse(tree, trait1, states.sd
<http://states.sd/>=Axis1M.sd, noroptimal.x, constant.x, control=control)*
*>f.h.nodrift <- constrain(f.h, drift~0)*
*>mle.h <- find.mle(f.h.nodrift, p.h, control=control, verbose=0)>mle.d.h
<- find.mle(f.h, coef(mle.h, TRUE), control=control, verbose=0)*
Finally, I compared these models into ANOVA:
*>anova(mle.c, linear=mle.l, sigmoid=mle.s,hump=mle.h,
drift.linear=mle.l.drift, drift.sigmoid=mle.s.drift,*
*drift.hump=mle.d.h)*
Df lnLik AIC ChiSq Pr(>|Chi|)
minimal 3 -138.48 282.96
linear 4 -138.13 284.27 0.697 0.4039
sigmoid 6 -137.12 286.23 2.733 0.4347
hump 6 -137.24 286.48 2.489 0.4772
drift.linear 5 -121.61 253.21 33.753 4.684e-08 ***
drift.sigmoid 7 -119.73 253.47 37.495 1.424e-07 ***
drift.hump 7 -121.48 256.97 33.996 7.467e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
If I check this result, considering AIC values and Df, the linear+drift
model seems to be the best one. I now want to plot the speciation rate
according to my trait1 values assuming this model. But I don't really
understand the parameters given by this model and how to code the drift
into a model to plot:
*>mle.l.drift$par* l.c l.m m.c
drift diffusion
-7.450302e+02 4.868530e+02 2.805533e-06 -1.125048e+01 4.334433e-01
I’m guessing that l.c and l.m refer to lambda estimates and m.c to mu, but
could you tell me what the .c and .m attributes are? Moreover, in order to
plot this function in the form y=ax+b, do you know which lambda estimates I
have to use? And how could I consider drift and diffusion into this model?
I am really sorry,
thank you
best regards
Charles Pouchon
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/
I'm sending this email because I did (or tried to do) a QuaSSE analysis
from the diversitree package to understand the impact of my traits values
into speciation rate. However, I am not very familar with such analysis and
with the tutorial script provided, and I would like to know if somebody
could help me understand some functions and results.
My trait vector corresponds to the positions of species into the first
component of a PCA made on morphological traits:
*>trait1 [1] 1.53635718 0.11435966 0.40744255 -0.16969599 0.83903277
0.27982589 0.34246181 1.40108911 0.92971423 0.58037380 1.01094382
0.48161104[13] 0.91112199 1.19083790 1.21212025 0.73398302 0.93707017
0.01386014 1.16228321 0.21442138 1.02680456 1.07615452 -1.13881116
1.10192106[25] 1.13232218 1.38341060 1.35456782 -3.81104255 -3.51734609
-3.73474404 -2.75423989 -3.96242786 -3.61512213 -3.98707154 -3.26697792
1.06852034[37] -3.69987363 -3.91920521 1.10428403 -0.25527186 1.25529456
1.41681721 1.27082040 -4.76293180 -1.49813425 -0.01049742 -1.49607314
-0.07286070[49] -7.63586519 -8.65295750*
*>range(trait1)*
*[1] -8.652957 1.536357*
control <- list(method='fftC')
I made a constant speciation and extinction model:*>f.c <- make.quasse(tree, trait1, states.sd <http://states.sd/>=sd(trait1)
constant.x, constant.x, control=control) >p <- starting.point.quasse(tree,
trait1)>f.c.nodrift <- constrain(f.c, drift~0)>mle.c <-
find.mle(f.c.nodrift, p, control=control, lower=0, verbose=0)*
Next, I would like to make linear, sigmoidal and hump-shaped functions with
and without drift. In order to make this, I used this script:
*>xr <- range(trait1)*
*>p.c <- mle.c$par>p.l <- c(p.c[1], l.m=0, p.c[2:3])>p.s <- p.h <-
c(p.c[1], p.c[1], mean(xr), 1, p.c[2:3])*
*>linear.x <- make.linear.x(xr[1],xr[2])>f.l <- make.quasse(tree, trait1,
states.sd <http://states.sd/>=Axis1M.sd, linear.x, constant.x,
control=control)>f.l.nodrift <- constrain(f.l, drift~0)>mle.l <-
find.mle(f.l.nodrift, p.l, control=control, lower=0, verbose=0)>mle.l.drift
<- find.mle(f.l, coef(mle.l, TRUE), control=control, verbose=0) >f.s <-
make.quasse(tree, trait1, states.sd <http://states.sd/>=Axis1M.sd,
sigmoid.x, constant.x, control=control)>f.s.nodrift <- constrain(f.s,
drift~0)>mle.s <- find.mle(f.s.nodrift, p.s, control=control,
verbose=0)>mle.s.drift <- find.mle(f.s, coef(mle.s, TRUE), control=control,
verbose=0) >f.h <- make.quasse(tree, trait1, states.sd
<http://states.sd/>=Axis1M.sd, noroptimal.x, constant.x, control=control)*
*>f.h.nodrift <- constrain(f.h, drift~0)*
*>mle.h <- find.mle(f.h.nodrift, p.h, control=control, verbose=0)>mle.d.h
<- find.mle(f.h, coef(mle.h, TRUE), control=control, verbose=0)*
Finally, I compared these models into ANOVA:
*>anova(mle.c, linear=mle.l, sigmoid=mle.s,hump=mle.h,
drift.linear=mle.l.drift, drift.sigmoid=mle.s.drift,*
*drift.hump=mle.d.h)*
Df lnLik AIC ChiSq Pr(>|Chi|)
minimal 3 -138.48 282.96
linear 4 -138.13 284.27 0.697 0.4039
sigmoid 6 -137.12 286.23 2.733 0.4347
hump 6 -137.24 286.48 2.489 0.4772
drift.linear 5 -121.61 253.21 33.753 4.684e-08 ***
drift.sigmoid 7 -119.73 253.47 37.495 1.424e-07 ***
drift.hump 7 -121.48 256.97 33.996 7.467e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
If I check this result, considering AIC values and Df, the linear+drift
model seems to be the best one. I now want to plot the speciation rate
according to my trait1 values assuming this model. But I don't really
understand the parameters given by this model and how to code the drift
into a model to plot:
*>mle.l.drift$par* l.c l.m m.c
drift diffusion
-7.450302e+02 4.868530e+02 2.805533e-06 -1.125048e+01 4.334433e-01
I’m guessing that l.c and l.m refer to lambda estimates and m.c to mu, but
could you tell me what the .c and .m attributes are? Moreover, in order to
plot this function in the form y=ax+b, do you know which lambda estimates I
have to use? And how could I consider drift and diffusion into this model?
I am really sorry,
thank you
best regards
Charles Pouchon
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/