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