Biotic interactions: modelling species dependencies with sequential SDMs

Overview

Species distribution models (SDMs) typically rely on abiotic predictors — climate, topography, remote sensing — but many species are shaped just as strongly by who they live with. Obligate specialists, in particular, cannot persist where their host is absent, no matter how suitable the climate. A sequential modelling approach captures this by first modelling the host species and then feeding those predictions as a biotic predictor into the dependent species’ model.

The interaction dataset is designed for exactly this workflow. It contains occurrence records for two Sierra Nevada (SE Spain) species: the host plant Androsace vitaliana (host_plant) and the obligate specialist butterfly Agriades zullichi (butterfly), together with random background points.

The dataset includes 3 response variables and 10 environmental predictors derived from remote sensing, climate, and topographic data at 100 m resolution.

Setup

The code below loads the R libraries and example data required to run this tutorial.

if (!requireNamespace("pak", quietly = TRUE)) {
  install.packages("pak")
}
pak::pkg_install(
  c(
    "dplyr",
    "terra",
    "sf",
    "stars",
    "mapview",
    "collinear",
    "spatialRF"
  ),
  ask = FALSE
)
#> ℹ Loading metadata database
#> ✔ Loading metadata database ... done
#> 
#> 
#> → Will install 1 system package:
#> + texlive  - huxtable
#> ℹ No downloads are needed
#> ℹ Installing system requirements
#> ℹ Executing `sudo sh -c apt-get -y update`
#> Get:1 file:/etc/apt/apt-mirrors.txt Mirrorlist [144 B]
#> Hit:2 http://azure.archive.ubuntu.com/ubuntu noble InRelease
#> Hit:3 http://azure.archive.ubuntu.com/ubuntu noble-updates InRelease
#> Hit:4 http://azure.archive.ubuntu.com/ubuntu noble-backports InRelease
#> Hit:5 http://azure.archive.ubuntu.com/ubuntu noble-security InRelease
#> Hit:6 https://packages.microsoft.com/repos/azure-cli noble InRelease
#> Hit:7 https://packages.microsoft.com/ubuntu/24.04/prod noble InRelease
#> Reading package lists...
#> ℹ Executing `sudo sh -c apt-get -y install libgdal-dev gdal-bin libgeos-dev libproj-dev libsqlite3-dev libabsl-dev cmake libssl-dev libudunits2-dev make libuv1-dev zlib1g-dev pandoc libpng-dev libicu-dev libfontconfig1-dev libfreetype6-dev libfribidi-dev libharfbuzz-dev texlive libxml2-dev`
#> Reading package lists...
#> Building dependency tree...
#> Reading state information...
#> libgdal-dev is already the newest version (3.8.4+dfsg-3ubuntu3).
#> gdal-bin is already the newest version (3.8.4+dfsg-3ubuntu3).
#> libgeos-dev is already the newest version (3.12.1-3build1).
#> libproj-dev is already the newest version (9.4.0-1build2).
#> libsqlite3-dev is already the newest version (3.45.1-1ubuntu2.5).
#> libabsl-dev is already the newest version (20220623.1-3.1ubuntu3.2).
#> cmake is already the newest version (3.28.3-1build7).
#> libssl-dev is already the newest version (3.0.13-0ubuntu3.7).
#> libudunits2-dev is already the newest version (2.2.28-7build1).
#> make is already the newest version (4.3-4.1build2).
#> libuv1-dev is already the newest version (1.48.0-1.1build1).
#> zlib1g-dev is already the newest version (1:1.3.dfsg-3.1ubuntu2.1).
#> pandoc is already the newest version (3.1.3+ds-2).
#> libpng-dev is already the newest version (1.6.43-5ubuntu0.5).
#> libicu-dev is already the newest version (74.2-1ubuntu3.1).
#> libfontconfig1-dev is already the newest version (2.15.0-1.1ubuntu2).
#> libfreetype-dev is already the newest version (2.13.2+dfsg-1ubuntu0.1).
#> libfribidi-dev is already the newest version (1.0.13-3build1).
#> libharfbuzz-dev is already the newest version (8.3.0-2build2).
#> libxml2-dev is already the newest version (2.9.14+dfsg-1.3ubuntu3.7).
#> The following additional packages will be installed:
#>   dvisvgm fonts-lmodern fonts-texgyre fonts-texgyre-math libgumbo2
#> libkpathsea6 libmujs3 libpotrace0 libptexenc1 libsynctex2 libteckit0
#>   libtexlua53-5 libwoff1 libzzip-0-13t64 lmodern mupdf-tools tex-gyre
#>   texlive-base texlive-binaries texlive-fonts-recommended texlive-latex-base
#>   texlive-latex-recommended tipa
#> Suggested packages:
#>   perl-tk xpdf | pdf-viewer xzdec texlive-binaries-sse2 hintview
#>   texlive-fonts-recommended-doc texlive-latex-base-doc wp2latex
#>   texlive-latex-recommended-doc texlive-luatex texlive-pstricks tipa-doc
#> The following NEW packages will be installed:
#>   dvisvgm fonts-lmodern fonts-texgyre fonts-texgyre-math libgumbo2
#> libkpathsea6 libmujs3 libpotrace0 libptexenc1 libsynctex2 libteckit0
#>   libtexlua53-5 libwoff1 libzzip-0-13t64 lmodern mupdf-tools tex-gyre texlive
#>   texlive-base texlive-binaries texlive-fonts-recommended texlive-latex-base
#>   texlive-latex-recommended tipa
#> 0 upgraded, 24 newly installed, 0 to remove and 69 not upgraded.
#> Need to get 131 MB of archives.
#> After this operation, 322 MB of additional disk space will be used.
#> Get:1 file:/etc/apt/apt-mirrors.txt Mirrorlist [144 B]
#> Get:2 http://azure.archive.ubuntu.com/ubuntu noble/main amd64 libkpathsea6 amd64 2023.20230311.66589-9build3 [63.0 kB]
#> Get:3 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 libpotrace0 amd64 1.16-2build1 [17.7 kB]
#> Get:4 http://azure.archive.ubuntu.com/ubuntu noble/main amd64 libwoff1 amd64 1.0.2-2build1 [45.3 kB]
#> Get:5 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 dvisvgm amd64 3.2.1+ds-1build1 [1042 kB]
#> Get:6 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 fonts-lmodern all 2.005-1 [4799 kB]
#> Get:7 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 fonts-texgyre all 20180621-6 [8350 kB]
#> Get:8 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 fonts-texgyre-math all 20180621-6 [2373 kB]
#> Get:9 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 libgumbo2 amd64 0.12.0+dfsg-2build1 [126 kB]
#> Get:10 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 libmujs3 amd64 1.3.3-3build2 [134 kB]
#> Get:11 http://azure.archive.ubuntu.com/ubuntu noble/main amd64 libptexenc1 amd64 2023.20230311.66589-9build3 [40.4 kB]
#> Get:12 http://azure.archive.ubuntu.com/ubuntu noble/main amd64 libsynctex2 amd64 2023.20230311.66589-9build3 [59.6 kB]
#> Get:13 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 libteckit0 amd64 2.5.12+ds1-1 [411 kB]
#> Get:14 http://azure.archive.ubuntu.com/ubuntu noble/main amd64 libtexlua53-5 amd64 2023.20230311.66589-9build3 [123 kB]
#> Get:15 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 libzzip-0-13t64 amd64 0.13.72+dfsg.1-1.2build1 [28.1 kB]
#> Get:16 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 lmodern all 2.005-1 [9542 kB]
#> Get:17 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 mupdf-tools amd64 1.23.10+ds1-1build3 [49.3 MB]
#> Get:18 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 tex-gyre all 20180621-6 [6396 kB]
#> Get:19 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 texlive-binaries amd64 2023.20230311.66589-9build3 [8529 kB]
#> Get:20 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 texlive-base all 2023.20240207-1 [21.7 MB]
#> Get:21 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 texlive-fonts-recommended all 2023.20240207-1 [4973 kB]
#> Get:22 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 texlive-latex-base all 2023.20240207-1 [1238 kB]
#> Get:23 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 texlive-latex-recommended all 2023.20240207-1 [8826 kB]
#> Get:24 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 texlive all 2023.20240207-1 [14.0 kB]
#> Get:25 http://azure.archive.ubuntu.com/ubuntu noble/universe amd64 tipa all 2:1.3-21 [2967 kB]
#> Preconfiguring packages ...
#> Fetched 131 MB in 2s (54.5 MB/s)
#> Selecting previously unselected package libkpathsea6:amd64.
#> (Reading database ...
#> (Reading database ... 5%(Reading database ... 10%(Reading database ... 15%(Reading database ... 20%(Reading database ... 25%(Reading database ... 30%(Reading database ... 35%(Reading database ... 40%(Reading database ... 45%(Reading database ... 50%(Reading database ... 55%
#> (Reading database ... 60%
#> (Reading database ... 65%
#> (Reading database ... 70%
#> (Reading database ... 75%
#> (Reading database ... 80%
#> (Reading database ... 85%
#> (Reading database ... 90%
#> (Reading database ... 95%
#> (Reading database ... 100%(Reading database ... 257223 files and directories currently installed.)
#> Preparing to unpack .../00-libkpathsea6_2023.20230311.66589-9build3_amd64.deb ...
#> Unpacking libkpathsea6:amd64 (2023.20230311.66589-9build3) ...
#> Selecting previously unselected package libpotrace0:amd64.
#> Preparing to unpack .../01-libpotrace0_1.16-2build1_amd64.deb ...
#> Unpacking libpotrace0:amd64 (1.16-2build1) ...
#> Selecting previously unselected package libwoff1:amd64.
#> Preparing to unpack .../02-libwoff1_1.0.2-2build1_amd64.deb ...
#> Unpacking libwoff1:amd64 (1.0.2-2build1) ...
#> Selecting previously unselected package dvisvgm.
#> Preparing to unpack .../03-dvisvgm_3.2.1+ds-1build1_amd64.deb ...
#> Unpacking dvisvgm (3.2.1+ds-1build1) ...
#> Selecting previously unselected package fonts-lmodern.
#> Preparing to unpack .../04-fonts-lmodern_2.005-1_all.deb ...
#> Unpacking fonts-lmodern (2.005-1) ...
#> Selecting previously unselected package fonts-texgyre.
#> Preparing to unpack .../05-fonts-texgyre_20180621-6_all.deb ...
#> Unpacking fonts-texgyre (20180621-6) ...
#> Selecting previously unselected package fonts-texgyre-math.
#> Preparing to unpack .../06-fonts-texgyre-math_20180621-6_all.deb ...
#> Unpacking fonts-texgyre-math (20180621-6) ...
#> Selecting previously unselected package libgumbo2:amd64.
#> Preparing to unpack .../07-libgumbo2_0.12.0+dfsg-2build1_amd64.deb ...
#> Unpacking libgumbo2:amd64 (0.12.0+dfsg-2build1) ...
#> Selecting previously unselected package libmujs3:amd64.
#> Preparing to unpack .../08-libmujs3_1.3.3-3build2_amd64.deb ...
#> Unpacking libmujs3:amd64 (1.3.3-3build2) ...
#> Selecting previously unselected package libptexenc1:amd64.
#> Preparing to unpack .../09-libptexenc1_2023.20230311.66589-9build3_amd64.deb ...
#> Unpacking libptexenc1:amd64 (2023.20230311.66589-9build3) ...
#> Selecting previously unselected package libsynctex2:amd64.
#> Preparing to unpack .../10-libsynctex2_2023.20230311.66589-9build3_amd64.deb ...
#> Unpacking libsynctex2:amd64 (2023.20230311.66589-9build3) ...
#> Selecting previously unselected package libteckit0:amd64.
#> Preparing to unpack .../11-libteckit0_2.5.12+ds1-1_amd64.deb ...
#> Unpacking libteckit0:amd64 (2.5.12+ds1-1) ...
#> Selecting previously unselected package libtexlua53-5:amd64.
#> Preparing to unpack .../12-libtexlua53-5_2023.20230311.66589-9build3_amd64.deb ...
#> Unpacking libtexlua53-5:amd64 (2023.20230311.66589-9build3) ...
#> Selecting previously unselected package libzzip-0-13t64:amd64.
#> Preparing to unpack .../13-libzzip-0-13t64_0.13.72+dfsg.1-1.2build1_amd64.deb ...
#> Unpacking libzzip-0-13t64:amd64 (0.13.72+dfsg.1-1.2build1) ...
#> Selecting previously unselected package lmodern.
#> Preparing to unpack .../14-lmodern_2.005-1_all.deb ...
#> Unpacking lmodern (2.005-1) ...
#> Selecting previously unselected package mupdf-tools.
#> Preparing to unpack .../15-mupdf-tools_1.23.10+ds1-1build3_amd64.deb ...
#> Unpacking mupdf-tools (1.23.10+ds1-1build3) ...
#> Selecting previously unselected package tex-gyre.
#> Preparing to unpack .../16-tex-gyre_20180621-6_all.deb ...
#> Unpacking tex-gyre (20180621-6) ...
#> Selecting previously unselected package texlive-binaries.
#> Preparing to unpack .../17-texlive-binaries_2023.20230311.66589-9build3_amd64.deb ...
#> Unpacking texlive-binaries (2023.20230311.66589-9build3) ...
#> Selecting previously unselected package texlive-base.
#> Preparing to unpack .../18-texlive-base_2023.20240207-1_all.deb ...
#> Unpacking texlive-base (2023.20240207-1) ...
#> Selecting previously unselected package texlive-fonts-recommended.
#> Preparing to unpack .../19-texlive-fonts-recommended_2023.20240207-1_all.deb ...
#> Unpacking texlive-fonts-recommended (2023.20240207-1) ...
#> Selecting previously unselected package texlive-latex-base.
#> Preparing to unpack .../20-texlive-latex-base_2023.20240207-1_all.deb ...
#> Unpacking texlive-latex-base (2023.20240207-1) ...
#> Selecting previously unselected package texlive-latex-recommended.
#> Preparing to unpack .../21-texlive-latex-recommended_2023.20240207-1_all.deb ...
#> Unpacking texlive-latex-recommended (2023.20240207-1) ...
#> Selecting previously unselected package texlive.
#> Preparing to unpack .../22-texlive_2023.20240207-1_all.deb ...
#> Unpacking texlive (2023.20240207-1) ...
#> Selecting previously unselected package tipa.
#> Preparing to unpack .../23-tipa_2%3a1.3-21_all.deb ...
#> Unpacking tipa (2:1.3-21) ...
#> Setting up fonts-texgyre-math (20180621-6) ...
#> Setting up libwoff1:amd64 (1.0.2-2build1) ...
#> Setting up libmujs3:amd64 (1.3.3-3build2) ...
#> Setting up libzzip-0-13t64:amd64 (0.13.72+dfsg.1-1.2build1) ...
#> Setting up libgumbo2:amd64 (0.12.0+dfsg-2build1) ...
#> Setting up libteckit0:amd64 (2.5.12+ds1-1) ...
#> Setting up libtexlua53-5:amd64 (2023.20230311.66589-9build3) ...
#> Setting up fonts-texgyre (20180621-6) ...
#> Setting up libkpathsea6:amd64 (2023.20230311.66589-9build3) ...
#> Setting up fonts-lmodern (2.005-1) ...
#> Setting up tex-gyre (20180621-6) ...
#> Setting up libsynctex2:amd64 (2023.20230311.66589-9build3) ...
#> Setting up libpotrace0:amd64 (1.16-2build1) ...
#> Setting up dvisvgm (3.2.1+ds-1build1) ...
#> Setting up mupdf-tools (1.23.10+ds1-1build3) ...
#> Setting up libptexenc1:amd64 (2023.20230311.66589-9build3) ...
#> Setting up texlive-binaries (2023.20230311.66589-9build3) ...
#> update-alternatives: using /usr/bin/xdvi-xaw to provide /usr/bin/xdvi.bin (xdvi.bin) in auto mode
#> update-alternatives: using /usr/bin/bibtex.original to provide /usr/bin/bibtex (bibtex) in auto mode
#> Setting up lmodern (2.005-1) ...
#> Setting up texlive-base (2023.20240207-1) ...
#> tl-paper: setting paper size for dvips to a4: /var/lib/texmf/dvips/config/config-paper.ps
#> tl-paper: setting paper size for dvipdfmx to a4: /var/lib/texmf/dvipdfmx/dvipdfmx-paper.cfg
#> tl-paper: setting paper size for xdvi to a4: /var/lib/texmf/xdvi/XDvi-paper
#> tl-paper: setting paper size for pdftex to a4: /var/lib/texmf/tex/generic/tex-ini-files/pdftexconfig.tex
#> Setting up texlive-latex-base (2023.20240207-1) ...
#> Setting up texlive-latex-recommended (2023.20240207-1) ...
#> Setting up texlive-fonts-recommended (2023.20240207-1) ...
#> Setting up tipa (2:1.3-21) ...
#> Setting up texlive (2023.20240207-1) ...
#> Processing triggers for man-db (2.12.0-4build2) ...
#> Not building database; man-db/auto-update is not 'true'.
#> Processing triggers for tex-common (6.18) ...
#> Running mktexlsr. This may take some time...
#> done.
#> Running updmap-sys. This may take some time...
#> done.
#> Running mktexlsr /var/lib/texmf ...
#> done.
#> Building format(s) --all.
#>  This may take some time...
#> done.
#> Processing triggers for install-info (7.1-3build2) ...
#> Processing triggers for fontconfig (2.15.0-1.1ubuntu2) ...
#> Processing triggers for libc-bin (2.39-0ubuntu8.7) ...
#> Running kernel seems to be up-to-date.
#> 
#> Restarting services...
#> Service restarts being deferred:
#>  systemctl restart getty@tty1.service
#>  systemctl restart serial-getty@ttyS0.service
#>  systemctl restart systemd-logind.service
#> No containers need to be restarted.
#> 
#> User sessions running outdated binaries:
#>  runner @ user manager service: systemd[1117]
#> No VM guests are running outdated hypervisor (qemu) binaries on this host.
#> ✔ 7 pkgs + 126 deps: kept 131 [42.6s]

Data Structure

The dataset interaction is an sf data frame with 1000 rows and 14 columns.

dplyr::glimpse(interaction)
#> Rows: 1,000
#> Columns: 14
#> $ butterfly               <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ host_plant              <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
#> $ class                   <fct> butterfly, butterfly, butterfly, butterfly, bu…
#> $ landsat_ndvi            <dbl> -0.3300000, -0.3299874, -0.3500096, -0.3599888…
#> $ landsat_pca_bands_123   <dbl> 82.99309, 71.99628, 81.99959, 81.99986, 85.995…
#> $ landsat_pca_bands_457   <dbl> 184.0043, 191.0012, 183.0027, 184.9953, 165.00…
#> $ rainfall_annual         <dbl> 833.9741, 802.0032, 1009.0140, 1020.0392, 934.…
#> $ solar_radiation         <dbl> 6408.431, 5598.753, 5580.747, 5335.484, 6403.8…
#> $ temperature_annual_mean <dbl> 8.000202, 8.299868, 7.199986, 7.200083, 7.8995…
#> $ temperature_summer_max  <dbl> 25.40020, 25.70001, 24.39999, 24.50018, 24.599…
#> $ temperature_winter_min  <dbl> -4.899726, -4.500000, -5.299973, -5.300069, -4…
#> $ topographic_complexity  <dbl> 269.9580, 211.9756, 235.9967, 241.9926, 267.03…
#> $ topographic_position    <dbl> -0.997110546, 2.998544216, 6.007016182, 14.000…
#> $ geometry                <POINT [m]> POINT (463352.8 4096288), POINT (463332.…

The three response variables are:

• class: factor with levels "butterfly", "host_plant", and "background", identifying the observation type of each row. You can interpret background points as the presence of “anything that is not butterfly or host plant”.

table(interaction$class)
#> 
#> background  butterfly host_plant 
#>        628        197        175

The map below shows a spatial representation of this column.

mapview::mapview(interaction, zcol = "class")

• butterfly: integer with values 1 and 0 representing butterfly and background records, and NA for the host plant records (to simplify subsetting).

table(na.omit(interaction$butterfly))
#> 
#>   0   1 
#> 628 197

• host_plant: integer with values 1 and 0 indicating host plant presence and background records, and NA for the butterfly records.

table(na.omit(interaction$host_plant))
#> 
#>   0   1 
#> 628 175

The dataset includes 10 environmental predictors covering remote sensing, climate, and topography at 100m resolution.

interaction_predictors
#>  [1] "landsat_ndvi"            "landsat_pca_bands_123"  
#>  [3] "landsat_pca_bands_457"   "rainfall_annual"        
#>  [5] "solar_radiation"         "temperature_annual_mean"
#>  [7] "temperature_summer_max"  "temperature_winter_min" 
#>  [9] "topographic_complexity"  "topographic_position"

Example Usage

In this example we train two models:

• Host plant model: Random forest model of the host plant using abiotic predictors.

• Butterfly model: Random forest model of the butterfly using the prediction of the host plant model as a biotic predictor along with the abiotic ones.

Host plant model

In this section we train and analyze a habitat suitability model for the host plant based on the abiotic predictors.

Data Preparation

To train the host plant model we first prepare a training data frame with the plant and background records, and a data frame with case coordinates to be used in spatial cross-validation.

host_plant_df <- interaction |> 
  dplyr::select(
    -butterfly
  ) |> 
  na.omit()

host_plant_xy <- host_plant_df |> 
    sf::st_coordinates() |> 
    as.data.frame() |> 
    dplyr::rename(
      x = X,
      y = Y
    )

Multicollinearity Filtering

The function collinear::collinear() selects a non-redundant subset of the most relevant predictors. The function ranks all predictors by the AUC of univariate binomial random forest models with case weights, and then applies a multicollinearity filtering that preserves the most important ones.

host_plant_selection <- collinear::collinear(
  df = host_plant_df,
  responses = "host_plant",
  predictors = interaction_predictors,
  f = collinear::f_binomial_rf,
  quiet = TRUE
)

Here are the selected predictors for the host plant:

host_plant_selection$host_plant$selection
#> [1] "temperature_annual_mean" "landsat_ndvi"           
#> [3] "landsat_pca_bands_457"   "topographic_position"   
#> [5] "topographic_complexity"  "solar_radiation"        
#> attr(,"validated")
#> [1] TRUE

Model Training

Here we train a random forest model with spatialRF::rf(), which computes case weights internally to control for class imbalance.

host_plant_model <- spatialRF::rf(
  data = sf::st_drop_geometry(host_plant_df),
  dependent.variable.name = "host_plant",
  predictor.variable.names = host_plant_selection$host_plant$selection,
  verbose = FALSE
)

To understand how well the model works on unseen data, we perform spatial cross-validation via spatialRF::rf_evaluate(), which partitions the data into spatially disjunct training and testing folds, providing an honest assessment of the model’s ability to predict habitat suitability for the host plant.

host_plant_model <- spatialRF::rf_evaluate(
  model = host_plant_model,
  xy = host_plant_xy,
  metrics = "auc",
  verbose = FALSE
)

host_plant_model$evaluation$aggregated |> 
  dplyr::filter(model == "Testing") |> 
  dplyr::glimpse()
#> Rows: 1
#> Columns: 11
#> $ model                     <chr> "Testing"
#> $ metric                    <chr> "auc"
#> $ median                    <dbl> 0.962
#> $ median_absolute_deviation <dbl> 0.0133434
#> $ q1                        <dbl> 0.95125
#> $ q3                        <dbl> 0.969
#> $ mean                      <dbl> 0.9618667
#> $ se                        <dbl> 0.002
#> $ sd                        <dbl> 0.01338999
#> $ min                       <dbl> 0.937
#> $ max                       <dbl> 0.994

The dataframe shows the stats of the AUC scores on 30 cross-validation repetitions. The median score indicates that the model has no issues discriminating between host plant and background records.

Finally, we predict the host plant suitability for all rows in the dataframe interaction, and write the result to the new column host_plant_suitability.

interaction$host_plant_suitability <- predict(
  host_plant_model,
  data = sf::st_drop_geometry(interaction)
)$predictions

Butterfly Model

Now that we have a trained and evaluated host plant model, let’s use its predictions as a biotic predictor for the butterfly.

Data Preparation

We prepare butterfly_df and butterfly_xy in the same way we did for the host plant. Additionally, we add the new biotic predictor host_plant_suitability to the vector interaction_predictors.

butterfly_df <- interaction |> 
  dplyr::select(
    -host_plant
  ) |> 
  na.omit()

butterfly_xy <- butterfly_df |> 
    sf::st_coordinates() |> 
    as.data.frame() |> 
    dplyr::rename(
      x = X,
      y = Y
    )

butterfly_predictors <- c(
  interaction_predictors, 
  "host_plant_suitability"
  )

Multicollinearity Filtering

Here we use collinear::collinear() again to select the best subset of non-collinear predictors for the butterfly.

butterfly_selection <- collinear::collinear(
  df = butterfly_df,
  responses = "butterfly",
  predictors = butterfly_predictors,
  f = collinear::f_binomial_rf,
  quiet = TRUE
)

The selected non-collinear predictors include the host plant’s suitability.

butterfly_selection$butterfly$selection
#> [1] "host_plant_suitability" "rainfall_annual"        "landsat_pca_bands_123" 
#> [4] "topographic_complexity" "topographic_position"   "solar_radiation"       
#> attr(,"validated")
#> [1] TRUE

Model Training

We use the same modelling pattern we used with the host plant to model the suitability of the butterfly and evaluate the model via cross-validation.

butterfly_model <- spatialRF::rf(
  data = sf::st_drop_geometry(butterfly_df),
  dependent.variable.name = "butterfly",
  predictor.variable.names = butterfly_selection$butterfly$selection,
  verbose = FALSE
)

butterfly_model <- spatialRF::rf_evaluate(
  model = butterfly_model,
  xy = butterfly_xy,
  metrics = "auc",
  verbose = FALSE
)

butterfly_model$evaluation$aggregated |> 
  dplyr::filter(model == "Testing") |> 
  dplyr::glimpse()
#> Rows: 1
#> Columns: 11
#> $ model                     <chr> "Testing"
#> $ metric                    <chr> "auc"
#> $ median                    <dbl> 0.993
#> $ median_absolute_deviation <dbl> 0.0044478
#> $ q1                        <dbl> 0.98825
#> $ q3                        <dbl> 0.994
#> $ mean                      <dbl> 0.9914
#> $ se                        <dbl> 0.001
#> $ sd                        <dbl> 0.004140631
#> $ min                       <dbl> 0.984
#> $ max                       <dbl> 0.998

The cross-validation shows that the model has an almost perfect ability to discriminate between butterfly and background records.

spatialRF::plot_importance(butterfly_model)

The permutation importance scores plotted above indicate that host_plant_suitability contributes most to model error when permuted.

spatialRF::plot_response_curves(
  model = butterfly_model,
  variables = "host_plant_suitability"
  )

The response curves show a clear sigmoidal increase in butterfly suitability as host plant suitability increases, with little variation across the quantiles of the remaining predictors — consistent with an obligate dependency.

Habitat Suitability Maps

With both models trained, we can now project habitat suitability across the Sierra Nevada. The package spatialData provides the raster dataset sierra_nevada_env for exactly this purpose.

sierra_nevada_env <- spatialData::interaction_extra()
#> spatialData::interaction_extra(): Downloading file 'sierra_nevada_env.tif'.
class(sierra_nevada_env)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"
names(sierra_nevada_env)
#>  [1] "landsat_ndvi"            "landsat_pca_bands_123"  
#>  [3] "landsat_pca_bands_457"   "rainfall_annual"        
#>  [5] "solar_radiation"         "temperature_annual_mean"
#>  [7] "temperature_summer_max"  "temperature_winter_min" 
#>  [9] "topographic_complexity"  "topographic_position"

We can project the model onto this raster dataset using terra::predict().

host_plant_prediction <- terra::predict(
  object = sierra_nevada_env,
  model = host_plant_model,
  na.rm = TRUE
)
names(host_plant_prediction) <- "host_plant_suitability"

mapview(host_plant_prediction, layer.name = "Host plant suitability")

To predict the habitat suitability for the butterfly, we first need to add the raster prediction of the host plant to sierra_nevada_env, and then apply terra::predict().

sierra_nevada_env <- c(
  sierra_nevada_env, 
  host_plant_prediction
  )

butterfly_prediction <- terra::predict(
  object = sierra_nevada_env,
  model = butterfly_model,
  na.rm = TRUE
)

names(butterfly_prediction) <- "butterfly_suitability"

mapview(butterfly_prediction, layer.name = "Butterfly suitability")

The habitat suitability maps of both species look very similar, but there are actual differences between them. The code below generates a map of differences between predicted suitabilities.

#raster of prediction differences
difference_predictions <- butterfly_prediction - host_plant_prediction

#compute maximum absolute value
difference_predictions_max <- terra::global(
  x = difference_predictions, 
  fun = range, 
  na.rm = TRUE
  ) |> 
  abs() |> 
  max()

#compute breaks for color palette
color_breaks <- seq(
  from = -difference_predictions_max, 
  to = difference_predictions_max, 
  length.out = 9
  )

#color palette
color_palette <- grDevices::colorRampPalette(c("#2166AC", "white", "#B2182B"))(length(color_breaks))

#differences map
mapview(
  difference_predictions, 
  col.regions = color_palette, 
  at = color_breaks, 
  layer.name = "Suitability difference"
  )

The divergent palette shows three distinct zones:

• Blue (negative values): indicating areas where the habitat is more suitable for the plant than the butterfly.
• White (values near zero): areas where both species show similar predicted suitability.
• Red (positive values): areas where butterfly suitability exceeds host plant suitability, which may reflect the butterfly’s mobility, or just be a model extrapolation artefact.

Conclusion

The interaction dataset illustrates a two-stage species distribution modelling workflow in which the output of one model becomes a predictor for another. By first modelling the host plant Androsace vitaliana from abiotic predictors and then using those predictions as a biotic input for the butterfly Agriades zullichi, we showed how obligate ecological dependencies can be encoded in a sequential modelling framework.

The resulting suitability maps and difference raster make interaction well suited for tutorials on biotic interactions, predictor engineering, spatial cross-validation, and divergent palette visualisation in R.

References

Barea-Azcón, J.M., Benito, B.M., Olivares, F.J., Ruiz, H., Martín, J., García, A.L., & López, R. (2014). Distribution and conservation of the relict interaction between the butterfly Agriades zullichi and its larval foodplant (Androsace vitaliana nevadensis). Biodiversity and Conservation, 23(4), 927–944. https://doi.org/10.1007/s10531-014-0643-4

Benito, B., Lorite, J., & Peñas, J. (2011). Simulating potential effects of climatic warming on altitudinal patterns of key species in Mediterranean-alpine ecosystems. Climatic Change, 108, 471–483. https://doi.org/10.1007/s10584-010-0015-3