0

I aim to calculate total heterogeneity and heterogeneity arising from each random factor from models generated with brm() function. Example of the model:

model <- brm(formula = r | se(sqrt(Vr)) ~ 1 + ( 1|STUDY_ID.1) + (1|INDEX) + (1|TREE_T) + (1|gr(animal, cov=A)),
                                data=Data_intra,
                                data2=list(A = forcedC_intra),
                                seed=1,
                                cores=4, chains=4, iter=4000,
                                control = list(adapt_delta = 0.9999, max_treedepth = 15))

This is a script of what I aim to do but for a model obtained with MCMCglmm:

  ## I2 Observation ID variance - Mode and HPD
  Index2_null_models <- which_model$VCV[, "Index"]/(which_model$VCV[, "animal"] +  which_model$VCV[, "Study_ID"] + which_model$VCV[, "Index"] + which_model$VCV[, "units"])
  Mode_Index2 <- posterior.mode(Index2_null_models)
  HPD_Index2 <- HPDinterval(Index2_null_models)
  ALL_Index2 <- as.data.frame(cbind(Mode_Index2, HPD_Index2))

The vcov() function for the brm models only works for the population-level factors, not the random ones. The ranef() seemed to be good, but this has a different length for each random factor, because it's not a cov matrix, it's the variance for each term (I think).

I tried the following:

BRMS_heterogeneity_intra <- function(which_model, model_name) {
  # Extract the posterior samples of the random effects
  random_effects <- ranef(which_model)
  
  # Calculate the variance for each random effect and normalize them
  index_variance <- var(random_effects$INDEX) / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T))
  study_variance <- var(random_effects$STUDY_ID.1) / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T))
  phylogeny_variance <- var(random_effects$animal) / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T))
  species_variance <- var(random_effects$TREE_T) / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T))

# Calculate the Highest Posterior Density (HPD) intervals for each normalized variance component
  index_hpd <- hdi(random_effects$INDEX / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T)))
  study_hpd <- hdi(random_effects$STUDY_ID.1 / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T)))
  phylogeny_hpd <- hdi(random_effects$animal / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T)))
  species_hpd <- hdi(random_effects$TREE_T / (var(random_effects$INDEX) + var(random_effects$STUDY_ID.1) + var(random_effects$animal) + var(random_effects$TREE_T)))

 # Create a data frame with the results
  results_BRMS_heterogeneity_intra <- data.frame(
    Index_Variance = index_variance,
    Index_HPD_Lower = index_hpd[1],
    Index_HPD_Upper = index_hpd[2],     
    Study_Variance = study_variance,
    Study_HPD_Lower = study_hpd[1],
    Study_HPD_Upper = study_hpd[2],
    Phylogeny_Variance = phylogeny_variance,
    Phylogeny_HPD_Lower = phylogeny_hpd[1],
    Phylogeny_HPD_Upper = phylogeny_hpd[2],
    Species_Variance = species_variance,
    Species_HPD_Lower = species_hpd[1],
    Species_HPD_Upper = species_hpd[2]
  )
  results_BRMS_heterogeneity_intra <- data.frame(Model = model_name, results_BRMS_heterogeneity_intra)
  return(results_BRMS_heterogeneity_intra)
}

But this is wrong because I can only calculate the hpd from the posterior.

Do you have any guesses on how could I do it?

bfrt
  • 61
  • 1

0 Answers0