Converted Simple Phylogenetic model example from BRMS to INLA

[Sam Passmore]

Hi,

I am trying to implement the simple phylogenetic model described here in BRMS into INLA.

Effectively, this model uses the IntraClass correlation to determine the effect of phylogeny in the model.

Although I think I have implemented the model in INLA correctly, I can't seem to figure out how to translate the precision estimates to a form that would allow my to calculate the intraclass correlation and replicate the example.

Could anyone offer any advice as to why I get somewhat similar results between these two approaches?

I have stored the example code in this gist & pasted it below if that is useful.

https://gist.github.com/SamPassmore/05db876b22d699fd5bd29f7696104a2e

Thanks,

Sam

library(geiger)
library(INLA)
library(dplyr)
library(brms)

# data
phylo <- ape::read.nexus("https://paul-buerkner.github.io/data/phylo.nex")
data_simple <- read.table(
  "https://paul-buerkner.github.io/data/data_simple.txt",
  header = TRUE
)

#### Brms model ####
A <- ape::vcv.phylo(phylo)
model_simple <- brm(
  phen ~ cofactor + (1|gr(phylo, cov = A)),
  data = data_simple,
  family = gaussian(),
  data2 = list(A = A),
  prior = c(
    prior(normal(0, 10), "b"),
    prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"),
    prior(student_t(3, 0, 20), "sigma")
  ),
  chains = 1, file = "brms_continuousphylo_test"
)

#### INLA model ####
# precision matrix
phylo_covar_mat <- ape::vcv(phylo)
phylo_covar_mat <- phylo_covar_mat/max(phylo_covar_mat)
phylo_covar_mat <- phylo_covar_mat / exp(determinant(phylo_covar_mat)$modulus[1] / nrow(phylo_covar_mat))
phylo_prec_mat <- solve(phylo_covar_mat)

#
pcprior_phy = list(prec = list(
  prior="pc.prec",
  param = c(0.5, 0.5))
)

data_simple$phylo_id = 1:nrow(data_simple)

inla_model = inla(phen ~ cofactor +
                    f(phylo_id,
                      model = "generic0",
                      Cmatrix = phylo_prec_mat,
                      constr = TRUE,
                      hyper = pcprior_phy),
                  data = data_simple)


#### Calculate Lambda ####
# BRMS Lambda formula
brms_model = summary(model_simple)
# from example
hyp <- "sd_phylo__Intercept^2 / (sd_phylo__Intercept^2 + sigma^2) = 0"
(hyp <- hypothesis(model_simple, hyp, class = NULL))

# Alternative brms calculation
# (sd_phylogeny^2)^2 / (sd_phylogeny^2 + sigma^2)
phylo_brms = brms_model$random$phylo$Estimate
sigma_brms = brms_model$spec_pars["sigma","Estimate"]
phylo_eff^2 / (phylo_eff^2 + sigma^2)

# INLA lambda formula
inla_summary = summary(inla_model)
phylo_inla = inla_summary$hyperpar["Precision for phylo_id","mean"]
sigma_inla =
  inla_summary$hyperpar["Precision for the Gaussian observations","mean"]

# INLA returns precision estimates but to copy BRMS we need SD
# according to https://www.flutterbys.com.au/stats/tut/tut12.10.html
# tau = 1 / sigma^2 so sigma^2 = 1 / tau
# But according to https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frspb.2019.2817&file=rspb20192817supp2.pdf
# sigma = 1 / phi

# So the calculation should either be:
(1 / (phylo_inla))^2 / ((1 / (phylo_inla))^2 + (1 / (sigma_inla))^2)
# or
(1 / (phylo_inla)) / ((1 / (phylo_inla)) + (1 / (sigma_inla)))
# But neither are correct

# 0.7 - Is this a coincidence? Not sure how this would be justified.
(1 / sqrt(sigma_inla)) / ((1 / sqrt(phylo_inla)) + (1 / sqrt(sigma_inla)))

[Helpdesk]

you can do inla.hyperpar.sample(...) and just compute the ratio by
sampling ?

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/b99bc103-25b9-41ae-9379-dd7153bfd48an%40googlegroups.com.

--
Håvard Rue
he...@r-inla.org

[Sam Passmore]

Dear Håvard,

Thank-you for your reply. Are you suggesting something along the lines of:

inla_hyperpar = inla.hyperpar.sample(1000, inla_model, improve.marginals = TRUE)

lambda = inla_hyperpar[,"Precision for phylo_id"] /
  (inla_hyperpar[,"Precision for phylo_id"] +
     inla_hyperpar[,"Precision for the Gaussian observations"])

This gives an estimate of ~0.84 where BRMS gives an estimate of 0.7

summary(lambda)
#   Min.      1st Qu.  Median   Mean   3rd Qu.   Max.
 # 0.6297  0.8087  0.8487  0.8426  0.8829  0.9591

So, if this calculation is correct, perhaps I need to consider the differences in priors between BRMS and INLA?

Thanks again,

Sam

[Helpdesk]

if you want to compare results, then yes, the priors must also be the
same.

[r.di...@gmail.com]

Hi Sam,

Try this:

phylo_covar_mat <- ape::vcv(phylo)
phylo_prec_mat <- solve(phylo_covar_mat)

pcprior_phy = list(prec = list(
  prior="pc.prec",

  param = c(20, 0.1))

)

data_simple$phylo_id = 1:nrow(data_simple)

inla_model = inla(phen ~ cofactor +
                    f(phylo_id,
                      model = "generic0",
                      Cmatrix = phylo_prec_mat,
                      constr = TRUE,
                      hyper = pcprior_phy),

                  data = data_simple,
                  control.family = list(hyper = pcprior_phy))

This should get you very close. To get even closer, try using exactly the same priors as in brms. There were two issues with your original code.

1) For the brms model you used an unscaled phylogenetic variance-covariance matrix, but for the INLA model you were scaling the matrix so that it had a generalized variance of 1, so this changes the scale of any estimated variance parameters. For the calculation you are doing, I would not scale the covariance matrix.

2) Your original prior was very unrealistic as it was essentially saying only put 50% of the prior probability density above 0.5 for the phylogenetic SD. The actual estimated value in brms is somewhere around 13, so this prior was shrinking your estimate considerably. Giving 10% probability density over 20 worked much better (I kind of arbitrarily chose those values but it seemed like those numbers would make sense). I also applied the same prior to the residual error as well, as was done in brms.

Using the above code I get:

> (1 / (phylo_inla)) / ((1 / (phylo_inla)) + (1 / (sigma_inla))) ## this is the correct version (1 / precision = variance)

[1] 0.7058824

So not exactly the same but pretty close.