leave_one_out() and orchard_leave1out() with phylogenetic matrix

library(metafor)
library(devtools)
library(rotl)
library(ape)
library(corrplot)
load_all() # Load dev package

Example with subset from Pottier data

Use only a subset to avoid long runtimes in leave-one-out

data("pottier")
pottier_subset <- pottier[1:20, ]

Get the phylo tree using rotl

unique_names <- unique(pottier_subset$search_string)
taxa <- tnrs_match_names(unique_names)

tree <- tol_induced_subtree(taxa$ott_id)
tree$tip.label <- strip_ott_ids(tree$tip.label)

plot(tree)

tree <- ape::compute.brlen(tree)
phylo_matrix <- ape::vcv(tree, corr = TRUE)

phylo_matrix

##                    Amalosia_lesueurii Tor_putitora Labeo_rohita Danio_rerio
## Amalosia_lesueurii               1.00         0.25         0.25        0.25
## Tor_putitora                     0.25         1.00         0.50        0.50
## Labeo_rohita                     0.25         0.50         1.00        0.75
## Danio_rerio                      0.25         0.50         0.75        1.00
## Propylea_japonica                0.00         0.00         0.00        0.00
##                    Propylea_japonica
## Amalosia_lesueurii                 0
## Tor_putitora                       0
## Labeo_rohita                       0
## Danio_rerio                        0
## Propylea_japonica                  1

# Simplified model just for the example
res <- rma.mv(yi = dARR,
              V = Var_dARR,
              test = "t",
              random = list(~1 | species_ID,
                            ~1 | phylogeny),
              R = list(phylogeny = phylo_matrix),
              data = pottier_subset)
res

## 
## Multivariate Meta-Analysis Model (k = 20; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor    R 
## sigma^2.1  0.0921  0.3035      5     no  species_ID   no 
## sigma^2.2  0.0000  0.0002      5     no   phylogeny  yes 
## 
## Test for Heterogeneity:
## Q(df = 19) = 3805.7554, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval    ci.lb   ci.ub    
##   0.0718  0.1365  0.5263  19  0.6047  -0.2138  0.3575    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

loo_test <- leave_one_out(res,
                          group = "phylogeny",
                          phylo_args = list(tree = tree, species_colname = "phylogeny"))

loo_test

##                 name    estimate    lowerCL   upperCL    lowerPR   upperPR
## 1 Amalosia_lesueurii  0.07920414 -0.2994704 0.4578787 -0.7641136 0.9225219
## 2       Tor_putitora -0.02777859 -0.2885881 0.2330309 -0.6064177 0.5508605
## 3       Labeo_rohita  0.13611075 -0.1684389 0.4406604 -0.4594178 0.7316393
## 4        Danio_rerio  0.03422792 -0.3252445 0.3937004 -0.7660659 0.8345217
## 5  Propylea_japonica  0.09647882 -0.2707291 0.4636867 -0.7221414 0.9150990

orchard_leave1out(res,
                  loo_output = loo_test,
                  group = "phylogeny",
                  xlab = "dARR",
                  ylab = "Species left out",
                  trunk.size = 1.2,
                  branch.size = 1.5,
                  alpha = 0.1,
                  legend.pos = "top.out")

Just for testing: remove Danio_rerio by hand

subset2 <- subset(pottier_subset, phylogeny != "Danio_rerio")

The function .prune_tree() is a simplified version of Dan’s code. It is used internally

kept_spp <- unique(subset2$phylogeny)
tree2 <- .prune_tree(tree, kept_spp)
plot(tree2)

tree2 <- ape::compute.brlen(tree2)
phylo_matrix2 <- ape::vcv(tree2, corr = TRUE)
phylo_matrix2

##                    Amalosia_lesueurii Tor_putitora Labeo_rohita
## Amalosia_lesueurii          1.0000000    0.3333333    0.3333333
## Tor_putitora                0.3333333    1.0000000    0.6666667
## Labeo_rohita                0.3333333    0.6666667    1.0000000
## Propylea_japonica           0.0000000    0.0000000    0.0000000
##                    Propylea_japonica
## Amalosia_lesueurii                 0
## Tor_putitora                       0
## Labeo_rohita                       0
## Propylea_japonica                  1

res2 <- rma.mv(yi = dARR,
               V = Var_dARR,
               test = "t",
               random = list(~1 | species_ID,
                             ~1 | phylogeny),
               R = list(phylogeny = phylo_matrix2),
               data = subset2)

By-hand result:

mod_results(res2, group = "phylogeny")

##      name   estimate    lowerCL   upperCL    lowerPR   upperPR
## 1 Intrcpt 0.03422792 -0.3252445 0.3937004 -0.7660659 0.8345217

leave-one-out result:

subset(loo_test$mod_table, name == "Danio_rerio")

##          name   estimate    lowerCL   upperCL    lowerPR   upperPR
## 4 Danio_rerio 0.03422792 -0.3252445 0.3937004 -0.7660659 0.8345217

Notes:

• To-do: Write good documentation and examples

• To-do: Add error messages for models that use phylogeny but don’t pass phylo_args (same for vcalc and robust)

• It doesn’t work with moderators. For example, I got an error if I use the model from the vignette. Don’t know if it is worth the effort to make this work.