Dear r-sig-phylo,
In the ouch package, the behavior of the simulate function been puzzling me
for a bit. I think I have pinned down the trouble in the package, would
appreciate if anyone could confirm this. I believe the problem stems from
the fact that the value of alpha in the hansen model as reported by
mo...@alpha seems to be the square root of the actual value of alpha used by
the model. Perhaps this is a somewhat subtle bug or I have misunderstood
how some of this works.
For instance, consider the example provided in the package documentation:
data(bimac)
tree <- with(bimac,ouchtree(node,ancestor,time/max(time),species))
print(h2 <- hansen(log(bimac['size']),tree,bimac['OU.1'],alpha=1,sigma=1))
The print command says alpha = 0.1921554
The call h...@alpha says alpha = .4383554
which is the sqrt of the other. Certainly these should all agree. The
summary command agrees with the print command:
smry = summary(h2)
smry$alpha
A little exploration seems to show that the print command is the correct
one. If we want to simulate a model with a particular alpha value, we first
have to generate a fitted model as above, then update the alpha value and
call simulate on the fitted model. I believe ouch is accidentally squaring
the value of alpha we give it before running the simulation. For example:
> h...@alpha = 5
> replicates <- as.data.frame(simulate(h2, 1000))
I find the expected variance by averaging the variance across the
replicates,
> mean(sapply(1:1000, function(i){var(replicates[i], na.rm=T)} ) )
[1] 0.0009774488
So did it use 5 or 25 as the alpha value? Since either value is large
compared to the length of the tree, the tip variance should be stationary
independent of tree structure so we can check it against the analytic
predicted variance (sigma^2/2 alpha):
> (h...@sigma)^2/(2...@alpha)
[1] 0.004836469
> (h...@sigma)^2/(2*(h...@alpha)^2)
[1] 0.0009672938
As the second equation agrees with the simulation, it believe that ouch has
erroneously squared the value of alpha in the simulation.
Have I done this correctly? Thanks for the check.
Cheers,
Carl
--
Carl Boettiger
Population Biology, UC Davis
http://two.ucdavis.edu/~cboettig
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
