February 09, 2018 Code R lhreg phylogeny detectability

It all started with this paper in *Methods in Ecol. Evol.* where we looked at
detectability of many species. So we wanted to use life history
traits to validate our results. But we had to cut the manuscript,
and there was this leftover with some neat patterns, but without much focus.
It took a few years, and the most positive peer-review experience ever,
and the paper is now early view in *Ecography*. This post is a quick summary of the goodies stuffed inside the **lhreg** R package that makes the whole analysis reproducible, and provides some functions for similar PGLMM models.

The R package is hosted on GitHub
(no CRAN version yet),
please submit any issues here.
The package is also archived on Zenodo with DOI 10.5281/zenodo.596410.
To install the package, use
`devtools::install_github("borealbirds/lhreg")`

.

Here, I am going to skim the implementation based on the more
complete supporting information of the paper which has all the
reproducible code (try `vignette(topic = "lhreg", package = "lhreg")`

after
installing and loading the package).
Here is the rendered html version.

The most important function is `lhreg`

which takes the following main arguments:

`Y`

: response vector,`X`

: model matrix for the mean.`SE`

: standard error estimate (observation error) for the response,`V`

: correlation matrix,

and fits a Multivariate Normal model to the observed `Y`

vector
with phylogenetically based (or any other known) correlations
and optionally with observation error (`SE`

), and covariate effects (`X`

).
The function is pretty bare-bones (i.e. no formula interface,
the design matrix `X`

needs to be properly specified through
e.g. `model.matrix()`

). The `lambda`

argument
is a non-negative number modifying the strength of phylogenetic effects.
`lambda = 0`

is equivalent to `lm`

with
`weights = 1/(SE^2)`

, `lambda = 1`

implies Brownian motion evolution,
`lambda = NA`

lets the function estimate it based on the data.

In terms of optimization, besides the algorithms from `stats::optim`

,
we also have differential evolution algorithm based on the
**DEoptim** package (a bit time consuming but very reliable).
The output object class has some methods defined (like `logLik`

and `summary`

)
and as a result AIC/BIC will work out of the box. The vignette also
describes a few techniques which are pretty nice to have in
a multivariate setting (i.e. profile likelihood, parametric bootstrap)
to support advanced hypothesis testing and model selection.

We used leave one out cross-validation to see how well we could predict the
values based on data from the other species, traits and phylogeny.
The conditional distribution we used for that is described in the paper which
made this exercise very straightforward.
Maybe it is just ignorance, but I couldn’t find another paper
that would have described it in a nice and useful manner,
however, if one wishes to make trait/phylogeny based
predictions for detectability, this formula is going to be
very useful (look inside the `loo2`

function for implementation).

At the end of the vignette, there is a hack based on `phytools::contMap`

function to produce *non-rainbow* colors.
(It was surprisingly *non-straightforward* to hack the code —
modular code please!)
The following figure shows the two input data vectors mirrored side-by-side:

I realize this is not a very detailed post, but the paper and the vignette should satisfy your curiosity. If you still have unanswered questions, feel free to ask them below!

**UPDATE**

- 2018-09-26: post on Ecography blog and video abstract.
- 2018-04-13: a related blog post in Hungarian.

I moved to Canada in 2008 to start a postdoctoral fellowship with Prof. Subhash Lele at the stats department of the University of Alberta. Subhash at the time just published a paper about a statistical technique called data cloning. Data cloning is a way to use Bayesian MCMC algorithms to do frequentist inference. Yes, you read that right.

- How many birds are out there?
- Fitting removal models with the detect R package
- Shiny slider examples with the intrval R package
- PVA: Publication Viability Analysis, round 3
- The progress bar just got a lot cheaper

ABMI (7) ARU (1) Alberta (1) BAM (1) C (1) CRAN (1) Hungary (2) JOSM (2) MCMC (1) PVA (2) PVAClone (1) QPAD (3) R (20) R packages (1) abundance (1) bioacoustics (1) biodiversity (1) birds (2) course (2) data (1) data cloning (4) datacloning (1) dclone (3) density (1) dependencies (1) detect (3) detectability (3) footprint (3) forecasting (1) functions (3) intrval (4) lhreg (1) mefa4 (1) monitoring (2) pbapply (5) phylogeny (1) plyr (1) poster (2) processing time (2) progress bar (4) publications (2) report (1) sector effects (1) shiny (1) single visit (1) site (1) slider (1) slides (2) special (3) species (1) trend (1) tutorials (2) video (4) workshop (1)