Statistical models fit to compositional data are often difficult to interpret due to the sum to one constraint on data variables. `DImodelsVis`

provides novel visualisations tools to aid with the interpretation for models where the predictor space is compositional in nature. All visualisations in the package are created using the `ggplot2`

plotting framework and can be extended like every other ggplot object.

You can install the released version of `DImodelsVis`

from CRAN by running:

You can install the development version of `DImodelsVis`

from GitHub with:

While sometimes it is of interest to model a compositional data response, there are times when the predictors of a response are compositional, rather than the response itself. Diversity-Interactions (DI) models (Kirwan et al., 2009, Connolly et al., 2013, Moral et al., 2023) are a regression based modelling technique for analysing and interpreting data from biodiversity experiments that explore the effects of species diversity on the different outputs (called ecosystem functions) produced by an ecosystem. Traditional techniques for analysing diversity experiments quantify species diversity in terms of species richness (i.e., the number of species present in a community). The DI method builds on top of this richness approach by taking the relative abundances of the species within in the community into account, thus the predictors in the model are compositional in nature. The DI approach can differentiate among different species identities as well as between communities with same set of species but with different relative proportions, thereby enabling us to better capture the relationship between diversity and ecosystem functions within an ecosystem. The `DImodels`

and `DImodelsMulti`

R packages are available to aid the user in fitting these models. The `DImodelsVis`

(DI models Visualisation) package is a complimentary package for visualising and interpreting the results from these models. However, the package is versatile and can be used with any standard statistical model object in R where the predictor space is compositional in nature.

The functions in the package can be categorised as functions for visualising model selection and validation or functions to aid with model interpretation. Here is a list of important visualisation functions present in the package along with a short description.

: Create diagnostics plots for a statistical model with the additional ability to overlay the points with pie-glyphs showing the proportions of the compositional predictor variables.`model_diagnostics`

: Show a visual comparison of selection criteria of different models. Can also show the split of an information criteria into deviance and penalty components to visualise why a parsimonious model would be preferable over a complex one.`model_selection`

: The predicted response for observations is visualised as a stacked bar-chart showing the contributions of each term in the regression model.`prediction_contributions`

: The predicted response for specific observations are shown using pie-glyphs along with the average change in the predicted response over the richness or evenness diversity gradients.`gradient_change`

: Assuming we have`conditional_ternary`

`n`

compositional variables, fix`n-3`

variables to have specific values and visualise the change in the predicted response across the remaining three variables as a contour plot in a ternary diagram.: Visualise the effect of increasing or decreasing a predictor variable (from a set of compositional predictor variables) on the predicted response whilst keeping the ratio of the other`visualise_effects`

`n-1`

compositional predictor variables constant.: Visualise the change in the predicted response along a straight line between two points in the simplex space.`simplex_path`

: A utility function to add prediction and associated uncertainty to data using a statistical model object or raw model coefficients.`add_prediction`

: Utility function to create all possible combinations of equi-proportional communities at a given level of richness from a set of n compositional variables.`get_equi_comms`

: A handy wrapper around the dplyr`custom_filter`

`filter()`

function enabling the user to filter rows which satisfy specific conditions for compositional data like all equi-proportional communities, or communities with a given value of richness without having to make any changes to the data or adding any additional columns.and`prop_to_tern_proj`

: Helper functions for converting between 3-d compositional data and their 2-d projection.`tern_to_prop_proj`

and`ternary_data`

: Visualise the change in the predicted response across a set of three compositional predictor variables as a contour map within a ternary diagram.`ternary_plot`

This dataset originates from a grassland biodiversity experiment conducted in Switzerland as part of the “Agrodiversity Experiment” Kirwan et al 2014. In this study, 68 grassland plots consisting of 1 to 4 species were established across a gradient of species diversity. The proportions of four species were varied across the plots: there were plots with 100% of a single species (called the monoculture of a species), and 2- and 4-species mixtures with varying proportions (e.g., (0.5, 0.5, 0, 0) and (0.7, 0.1, 0.1, 0.1)). Nitrogen fertilizer (at 50 or 150 kg/ha/yr) and seeding density (low or high) treatments were also manipulated across the plots. The total annual yield per plot was recorded for the first year after establishment. The data is available in the `DImodels`

R package. An analysis of the this dataset can be found in Kirwan et al 2009. For our example we only consider the plots that received the 150 kg nitrogen treatment. The four species proportions form our compositional predictors while the annual yield is our continuous response.

```
data(Switzerland)
my_data <- Switzerland[Switzerland$nitrogen == 150, ]
head(my_data)
#> plot nitrogen density p1 p2 p3 p4 yield
#> 1 1 150 high 0.70 0.10 0.10 0.10 13.51823
#> 2 2 150 high 0.10 0.70 0.10 0.10 13.16549
#> 3 3 150 high 0.10 0.10 0.70 0.10 19.95682
#> 4 4 150 high 0.10 0.10 0.10 0.70 17.93976
#> 5 5 150 high 0.25 0.25 0.25 0.25 13.74719
#> 6 6 150 high 0.40 0.40 0.10 0.10 15.11899
```

We fit different models with different interaction structures as described in Moral et al 2023.

```
mod_ID <- DI(y = "yield", prop = 4:7,
DImodel = "ID", data = my_data)
#> Fitted model: Species identity 'ID' DImodel
mod_AV <- DI(y = "yield", prop = 4:7,
DImodel = "AV", data = my_data)
#> Fitted model: Average interactions 'AV' DImodel
mod_FG <- DI(y = "yield", prop = 4:7,
DImodel = "FG", data = my_data,
FG = c("G", "G", "H", "H"))
#> Fitted model: Functional group effects 'FG' DImodel
mod_ADD <- DI(y = "yield", prop = 4:7,
DImodel = "ADD", data = my_data)
#> Fitted model: Additive species contributions to interactions 'ADD' DImodel
mod_FULL <- DI(y = "yield", prop = 4:7,
DImodel = "FULL", data = my_data)
#> Fitted model: Separate pairwise interactions 'FULL' DImodel
```

We can visualise model selection by passing our models as a list to the `model_selection`

function and visualising the best performing metric across different information criteria. Run `?model_selection`

or see the associated vignette for more information on customising the plot.

```
mods = list("ID" = mod_ID, "AV" = mod_AV, "FG" = mod_FG,
"ADD" = mod_ADD, "FULL" = mod_FULL)
model_selection(models = mods, metric = c("AIC", "AICc", "BIC", "BICc"))
```

The model `mod_FG`

(labelled as “FG”) is the best model as it has the lowest value for all the information criteria. We proceed with this model.

The coefficients are as follows

```
summary(mod_FG)
#>
#> Call:
#> glm(formula = fmla, family = family, data = data)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> p1_ID 8.5406 0.7627 11.198 2.50e-14 ***
#> p2_ID 8.7926 0.7627 11.528 9.70e-15 ***
#> p3_ID 16.0825 0.7627 21.086 < 2e-16 ***
#> p4_ID 11.9263 0.7627 15.637 < 2e-16 ***
#> FG_bfg_G_H 17.3817 2.1713 8.005 4.66e-10 ***
#> FG_wfg_G 7.6604 4.4234 1.732 0.0905 .
#> FG_wfg_H 1.0119 4.4234 0.229 0.8201
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 2.370592)
#>
#> Null deviance: 10290.75 on 50 degrees of freedom
#> Residual deviance: 101.94 on 43 degrees of freedom
#> AIC: 193.51
#>
#> Number of Fisher Scoring iterations: 2
```

After choosing a model we can create diagnostics plot where the points are replaced by pie-glyphs showing the proportions of the compositional variables. Run `?model_diagnostics`

or see the associated vignette for more information on customising the plot.

Replacing the points with pie-glyphs could help us to quickly identify any problematic observations in the model. For example, we can see here that the diagnostics plots look fine and no assumptions seem to be violated. However, we can quickly spot that the all monocultures (communities with only 1 species) and certain communities with 2 species have high leverage values compared to all other communities in the data.

We visualise the predicted response for specific observations as a stacked bar-chart showing the contribution (predictor coefficient * predictor value) of each term in the model.

```
prediction_contributions(model = mod_FG,
data = my_data[c(1:5, 12:15),])
#> ✔ Finished data preparation.
#> ✔ Created plot.
```

The coloured bars show the contributions of the different terms in the model. The contribution is defined as the product of the coefficient and value for each predictor variable. Thus, the contribution for a term would be zero if it’s value in an observation is zero regardless of it’s coefficient value (e.g. prediction bars for the monocultures at the right of the graph).

This plot would aid in understanding why certain observations have higher predictions. For example, we can see that higher predictions are primarily driven by the `p3_ID`

and `p4_ID`

terms and hence the `p1`

and `p2`

monocultures have low predictions as the all the other terms have a value of zero here. Similarly, we can also see that mixtures dominated by `p3`

perform the best. Run `?prediction_contributions`

or see the associated vignette for more information on creating and customising the plot.

Next we show a scatterplot of the predicted response across all equi-proportional mixtures at each level of richness (number of species in a community). The points are replaced with pie-glyphs to show the proportions of the different species while the dashed black line shows the average change in response over the richness gradient.

```
# Create data including all equi-proportional communities at
# each level of richness
plot_data <- get_equi_comms(nvars = 4, variables = paste0("p", 1:4))
# Show the average change over richness
gradient_change(mod_FG, data = plot_data)
#> ✔ Finished data preparation
#> ✔ Created plot.
```

This shows that on average the predicted response increases as richness increases but at a saturating rate. Run `?gradient_change`

or see the associated vignette for more information on creating and customising the plot.

Ternary diagrams are a great tool for visualising the change in a continuous response, however they can only be created for examples with three compositional variables. If we have more than three compositional variables we create conditional ternary diagrams where fix n-3 compositional variables to have specific values and visualise the change in the predicted response across the remaining three variables as a contour plot in a ternary diagram.

For this example, since we have four species we can condition one of the species (say `p2`

) to have a specific values 0.2, 0.5, and 0.8 and see how the response is affected as we change the proportions of the other three species whilst ensuring that the sum of the four species proportions is 1.

```
conditional_ternary(model = mod_FG,
tern_vars = c("p1", "p3", "p4"),
conditional = data.frame("p2" = c(0.2, 0.5, 0.8)))
#> ✔ Finished data preparation.
#> ✔ Created plot.
```

This figure shows that the predicted response decreases as we increase the proportion of `p2`

and is maximised as we increase the proportion of `p3`

. Run `?conditional_ternary`

or see the associated vignette for more information on creating and customising the plot.

Effects plots are great for visualising the average effect of a predictor in a model. However, if the predictors are compositional in nature, then standard effects plots are not very useful because of the sum to 1 constraint. The `visualise_effects`

function creates effects plot for the compositional predictors in a model by ensuring that the sum to one constraint is respected as we increase or decrease the proportion of a particular variable.

In this example we specify few communities using the `data`

argument and see the change in the predicted response as we increase the proportions of the species `p2`

and `p3`

in each community whilst keeping the ratio of the other species constant.

```
visualise_effects(model = mod_FG, data = my_data[1:15, ],
var_interest = c("p2", "p3"))
#> ✔ Finished data preparation.
#> ✔ Created plot.
```

The grey lines show the effect (on the predicted response) of increasing the species of interest within a particular community while the solid black line shows the average effect of increasing the proportion of a species on the predicted response. It can be seen that for all communities increasing `p2`

results in a decrease in the predicted response while increasing `p3`

has a positive effect on the predicted response. Run `?visualise_effects`

or see the associated vignette for more information on creating and customising the plot.

The concept used in `visualise_effects`

can be extended to look at the change in the predicted response as we move in a straight line between any two points within the simplex space. The interpolation constant (shown on the X-axis) is a number between 0 and 1 identifying points along the straight line between the start and end points. We can even traverse a path comprising of multiple points within the simplex and see the change in the predicted response. In this example we show the change in the response as we move from the centroid mixture to the monoculture of each of the four species.

```
simplex_path(model = mod_FG,
starts = my_data[5, ],
ends = my_data[12:15, ])
#> ✔ Finished data preparation.
#> ✔ Created plot.
```

We can see that moving from the centroid community to `p1`

, `p2`

, and `p4`

decreases the predicted response, while moving towards a monoculture of `p3`

increases the response. Run `?simplex_path`

or see the associated vignette for more information on creating and customising the plot.

- DI models website: https://dimodels.com
- Package website: https://rishvish.github.io/DImodelsVis/
- Github repo: https://github.com/rishvish/DImodelsVis
- Report bugs: https://github.com/rishvish/DImodelsVis/issues

- Connolly J, T Bell, T Bolger, C Brophy, T Carnus, JA Finn, L Kirwan, F Isbell, J Levine, A Lüscher, V Picasso, C Roscher, MT Sebastia, M Suter and A Weigelt (2013) An improved model to predict the effects of changing biodiversity levels on ecosystem function. Journal of Ecology, 101, 344-355.
- Moral, R.A., Vishwakarma, R., Connolly, J., Byrne, L., Hurley, C., Finn, J.A. and Brophy, C., 2023. Going beyond richness: Modelling the BEF relationship using species identity, evenness, richness and species interactions via the DImodels R package. Methods in Ecology and Evolution, 14(9), pp.2250-2258
- Kirwan L., Connolly J., Finn J.A., Brophy C., Lüscher A., Nyfeler D. and Sebastia M.T. 2009. Diversity-interaction modelling - estimating contributions of species identities and interactions to ecosystem function. Ecology, 90, 2032-2038.
- Kirwan, L., Connolly, J., Brophy, C., Baadshaug, O.H., Belanger, G., Black, A., Carnus, T., Collins, R.P., Čop, J., Delgado, I. and De Vliegher, A., 2014. The Agrodiversity Experiment: three years of data from a multisite study in intensively managed grasslands.