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.
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)))