I am using Phylogenetic Least Squares models, fitted with the phylolm package, to understand the effect of continuous and categorical variables. I used the MuMIn package to gel model-averaged coefficients, and my response was log-transformed.
How do I visualize slopes for multiple factor levels, with confidence bands, using model-averaged coefficients? My aim is a plot that combines the raw data, and the slopes for each factor level (in other words, to visualize two effects in a single plot).
I would normally rely on predict(se.fit = TRUE) or the ggpredict package, but unfortunately, these methods are not compatible with model-averaged phylolm objects.
Therefore, I am thinking on relying on model-averaged coefficients, standard errors / confidence intervals. In particular, my model has 3 levels, and a common slope (the reprex below has significant differences among factor levels, but it serves to illustrate my doubts).
From previous SO questions (here, here), I've got an idea about how to plot slopes when there is a transformation: inverse.function(intercept + slope * variable).
However, I am stumped when it comes to visualizing the uncertainty around the slopes for multiple factor levels. For instance, if using the method of +/- 1.96 * Standard error, how should I take into account the uncertainty around the factor level and the slope? In the reprex below I only used the standard error of the intercept and each factor level. A similar doubt would apply if relying on confidence intervals.
Note that in the reprex below, phylolm is loaded before MuMIn, this is needed for compatibility of logLik method (https://github.com/lamho86/phylolm/issues/65).
library(caper)
#> Warning: package 'caper' was built under R version 4.1.3
#> Loading required package: ape
#> Warning: package 'ape' was built under R version 4.1.3
#> Loading required package: MASS
#> Loading required package: mvtnorm
#> Warning: package 'mvtnorm' was built under R version 4.1.1
library(phylolm)
#> Warning: package 'phylolm' was built under R version 4.1.3
library(MuMIn)
#> Warning: package 'MuMIn' was built under R version 4.1.3
#> Registered S3 methods overwritten by 'MuMIn':
#> method from
#> nobs.pgls caper
#> nobs.phylolm phylolm
#> logLik.phylolm phylolm
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 4.1.3
#> Warning: package 'ggplot2' was built under R version 4.1.3
#> Warning: package 'tibble' was built under R version 4.1.3
#> Warning: package 'tidyr' was built under R version 4.1.3
#> Warning: package 'readr' was built under R version 4.1.3
#> Warning: package 'purrr' was built under R version 4.1.3
#> Warning: package 'dplyr' was built under R version 4.1.3
#> Warning: package 'stringr' was built under R version 4.1.3
#> Warning: package 'forcats' was built under R version 4.1.3
## Create data
data(shorebird)
shorebird.data$M.Mass.std <- (shorebird.data$M.Mass - mean(shorebird.data$M.Mass))/sd(shorebird.data$M.Mass)
str(shorebird.data)
#> 'data.frame': 71 obs. of 7 variables:
#> $ Species : chr "Actophilornis_africanus" "Aphriza_virgata" "Arenaria_interpres" "Arenaria_melanocephala" ...
#> $ M.Mass : num 143 186 108 114 151 ...
#> $ F.Mass : num 261 216 113 124 164 ...
#> $ Egg.Mass : num 8.6 22.4 17.9 17.3 23.5 11.2 10.7 9.6 19.3 13.3 ...
#> $ Cl.size : num 4 4 3.5 4 3.99 3.9 3.9 4 3.7 3.9 ...
#> $ Mat.syst : Factor w/ 3 levels "MO","PA","PG": 2 1 1 1 1 2 1 1 1 1 ...
#> $ M.Mass.std: num -0.1658 0.0666 -0.3557 -0.3255 -0.1238 ...
str(shorebird.tree)
#> List of 4
#> $ edge : num [1:132, 1:2] 72 73 74 75 76 77 78 79 79 80 ...
#> $ edge.length: num [1:132] 18.1 13.4 13.4 2.6 2.3 ...
#> $ tip.label : chr [1:71] "Catoptrophorus_semipalmatus" "Tringa_ochropus" "Tringa_stagnatilis" "Tringa_flavipes" ...
#> $ Nnode : num 62
#> - attr(*, "class")= chr "phylo"
#> - attr(*, "origin")= chr "/Users/dorme/Scripts/CAIC/caic/pkg/data/shorebird.nex"
## Fit model
mod <- phylolm::phylolm(log(Egg.Mass) ~ M.Mass.std + Mat.syst, data = shorebird.data, phy = shorebird.tree, model = "lambda")
## Model average
options(na.action = "na.fail")
mod.d <- MuMIn::dredge(mod, rank = "AICc")
#> Fixed term is "(Intercept)"
mod.avg <- MuMIn::model.avg(mod.d, revised.var = TRUE, fit = TRUE)
# Coefficients and model-averaged confidence intervalss
summary(mod.avg)$coefmat.full
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.9477111 0.24188761 12.186284 0.00000000
#> M.Mass.std 0.4704710 0.05739476 8.197107 0.00000000
#> Mat.systPA -0.2154969 0.09765281 2.206766 0.02733044
#> Mat.systPG 0.2837279 0.17353949 1.634947 0.10206008
# Keep coefficients
intercept <- summary(mod.avg)$coefmat.full[1, 1]
slope <- summary(mod.avg)$coefmat.full[2, 1]
level.pa <- summary(mod.avg)$coefmat.full[3, 1]
level.pg <- summary(mod.avg)$coefmat.full[4, 1]
# Keep standard errors
intercept.std <- summary(mod.avg)$coefmat.full[1, 2]
slope.std <- summary(mod.avg)$coefmat.full[2, 2]
level.pa.std <- summary(mod.avg)$coefmat.full[3, 2]
level.pg.std <- summary(mod.avg)$coefmat.full[4, 2]
# Slopes at the response scale, from model-averaged coefficients
seq.x <- seq(min(shorebird.data$M.Mass.std), max(shorebird.data$M.Mass.std), by = 0.05)
slope.mo <- exp(intercept + slope * seq.x)
slope.pa <- exp(intercept + level.pa + slope * seq.x)
slope.pg <- exp(intercept + level.pg + slope * seq.x)
# Uncertainty around slope, from model-averaged standard errors
uncert.mo.low <- exp(intercept + slope * seq.x - 1.96 * intercept.std)
uncert.mo.up <- exp(intercept + slope * seq.x + 1.96 * intercept.std)
uncert.pa.low <- exp(intercept + level.pa + slope * seq.x - 1.96 * level.pa.std)
uncert.pa.up <- exp(intercept + level.pa + slope * seq.x + 1.96 * level.pa.std)
uncert.pg.low <- exp(intercept + level.pg + slope * seq.x - 1.96 * level.pg.std)
uncert.pg.up <- exp(intercept + level.pg + slope * seq.x + 1.96 * level.pg.std)
ggplot() +
geom_jitter(data = shorebird.data, aes(x = M.Mass.std, y = Egg.Mass, col = Mat.syst)) +
geom_line(aes(x = seq.x, y = slope.mo), col = "red") +
geom_ribbon(aes(x = seq.x, y = slope.mo, ymin = uncert.mo.low, ymax = uncert.mo.up), fill = "red", alpha = 0.5) +
geom_line(aes(x = seq.x, y = slope.pa), col = "green") +
geom_ribbon(aes(x = seq.x, y = slope.pa, ymin = uncert.pa.low, ymax = uncert.pa.up), fill = "green", alpha = 0.5) +
geom_line(aes(x = seq.x, y = slope.pg), col = "blue") +
geom_ribbon(aes(x = seq.x, y = slope.pg, ymin = uncert.pg.low, ymax = uncert.pg.up), fill = "blue", alpha = 0.5) +
theme_bw()
confint(mod.avg, full = TRUE)
#> 2.5 % 97.5 %
#> (Intercept) 2.47362007 3.42180208
#> M.Mass.std 0.35797932 0.58296266
#> Mat.systPA -0.40689283 -0.02410087
#> Mat.systPG -0.05640321 0.62385908
sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252
#> [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
#> [5] LC_TIME=Spanish_Spain.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] forcats_0.5.1 stringr_1.5.0 dplyr_1.1.1 purrr_1.0.1
#> [5] readr_2.1.2 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.2
#> [9] tidyverse_1.3.1 MuMIn_1.46.0 phylolm_2.6.2 caper_1.0.1
#> [13] mvtnorm_1.1-3 MASS_7.3-54 ape_5.6-2
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.8.3 lubridate_1.8.0 lattice_0.20-44 listenv_0.8.0
#> [5] assertthat_0.2.1 digest_0.6.29 utf8_1.2.2 parallelly_1.31.1
#> [9] cellranger_1.1.0 R6_2.5.1 backports_1.4.1 reprex_2.0.1
#> [13] stats4_4.1.0 evaluate_0.15 httr_1.4.2 highr_0.9
#> [17] pillar_1.9.0 rlang_1.1.0 readxl_1.4.0 rstudioapi_0.13
#> [21] Matrix_1.3-3 rmarkdown_2.13 labeling_0.4.2 munsell_0.5.0
#> [25] broom_0.8.0 modelr_0.1.8 compiler_4.1.0 xfun_0.30
#> [29] pkgconfig_2.0.3 globals_0.15.0 htmltools_0.5.2 tidyselect_1.2.0
#> [33] codetools_0.2-18 fansi_1.0.3 future_1.25.0 crayon_1.5.1
#> [37] tzdb_0.3.0 dbplyr_2.1.1 withr_2.5.0 grid_4.1.0
#> [41] jsonlite_1.8.0 nlme_3.1-152 gtable_0.3.0 lifecycle_1.0.3
#> [45] DBI_1.1.2 magrittr_2.0.3 scales_1.2.0 future.apply_1.9.0
#> [49] cli_3.6.0 stringi_1.7.6 farver_2.1.0 fs_1.5.2
#> [53] xml2_1.3.3 ellipsis_0.3.2 generics_0.1.2 vctrs_0.6.1
#> [57] tools_4.1.0 glue_1.6.2 hms_1.1.1 parallel_4.1.0
#> [61] fastmap_1.1.0 yaml_2.3.5 colorspace_2.0-3 rvest_1.0.2
#> [65] knitr_1.38 haven_2.5.0
Created on 2023-11-06 by the reprex package (v2.0.1)