Spatial cross-validation

In this tutorial, we will TK

using SpeciesDistributionToolkit
using CairoMakie
using PrettyTables
using Statistics

Getting the data

records = OccurrencesInterface.__demodata()
landmass = getpolygon(PolygonData(OpenStreetMap, Places); place = "Oregon")
records = Occurrences(mask(records, landmass))
spatial_extent = SpeciesDistributionToolkit.boundingbox(landmass)

(left = -124.70354461669922, right = -116.46315002441406, bottom = 41.99179458618164, top = 46.292396545410156)

out

provider = RasterData(CHELSA2, BioClim)
L = SDMLayer{Float32}[
    SDMLayer(
        provider;
        layer = x,
        spatial_extent...,
    ) for x in eachindex(layers(provider))
];

mask!(L, landmass)

19-element Vector{SDMLayer{Float32}}:
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)
 🗺️  A 517 × 990 layer (411408 Float32 cells)

We will generate pseudo-absences at random in a radius around each known observations. The distance between an observation and a pseudo-absence is between 40 and 130 km.

presencelayer = mask(first(L), records)
absencemask =
    pseudoabsencemask(BetweenRadius, presencelayer; closer = 40.0, further = 120.0)
absencelayer = backgroundpoints(absencemask, 2sum(presencelayer))

🗺️  A 517 × 990 layer with 411408 Bool cells
   Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs

Training on model v1

model = SDM(PCATransform, Logistic, L, presencelayer, absencelayer);
hyperparameters!(classifier(model), :interactions, :self);

variogram for temperature

Lᵢ = L[1]
x, y, n = variogram(Lᵢ; width=15., shift=3.)

([16.83556311451665, 16.83556311451665, 23.103782688601665, 25.753389435159313, 28.298867191669384, 30.657022913912577, 32.042241346580234, 36.18293768530572, 39.46223456155963, 41.90342885771339, 45.107900624058374, 46.52381014461397, 49.70005576227489, 54.24387545925999, 56.11604329122996, 60.49043747545854, 61.26933286755387, 64.52196415855728, 67.91589943191376, 72.30545444364446, 76.73737261568482, 79.09598697039293, 81.44237740177397, 83.46404158764715, 86.11157571161974, 89.59391123865103, 92.70727213604165, 95.16044491310574, 98.55411190118235, 101.87565718292366, 105.39766959897783, 108.38374250779165, 110.72073962745966, 112.60285576329443, 115.94970658684183, 119.13032562470215, 122.76457842600783, 125.85079452432353, 127.2919906066322, 131.2284841926793, 134.63159724549587, 137.68456898252785, 140.96964498293562, 143.78770506742262, 146.85351968095247, 150.55520602056322, 152.33058577231373, 154.59737290558706, 158.43964681271254, 161.34549715833595, 164.6447895800366, 166.92667993085445, 168.85183646080665, 173.67823790072296, 177.02007474290096, 179.86113946207746, 182.39287035423123, 185.99029153114267, 189.56585066798834, 192.19426797963453, 194.7957206741648, 197.0883049662358, 199.58000589162927, 202.7424724806915, 205.36777749536822, 210.30658690823404, 213.92915253423223, 216.2185136090813, 218.76633075782172, 220.25192662990827, 224.29580352071295, 227.93866692550117, 231.00587161265187, 233.29421313625897, 235.40173809103658, 238.6667182190976, 242.3180300654441, 246.35807670434954, 249.30848272728159, 251.46427346961306, 254.54272188590946, 257.3735730149165, 260.4936997956129, 263.39783416706126, 266.0697914324131, 268.94233314425776, 272.78783458125514, 275.4212309213428, 278.1600449260087, 281.34040052572885, 284.14034053696116, 287.9891968899236, 291.2937949233772, 293.6369647765982, 297.035815103732, 300.05437008051416, 301.73877248074274, 304.9078672145201, 307.30989579831464, 311.0059483902359, 315.91282755040334, 317.7215792767945, 319.91913540874486, 323.47625743131124, 325.9153388161379, 329.40351389244853, 332.7856661226507, 335.1126518604, 338.62614005622606, 341.682691227933, 344.26007359194597, 347.1929719302386, 350.5075694110729, 353.9137640459961, 357.6734935075512, 359.80097649711945, 362.63984192473396, 364.9184615772338, 367.390090253687, 370.32596781092184, 373.6580333110306, 378.07779476022483, 380.1248456003553, 383.92988735690244, 386.6404042981029, 389.55315587799, 393.03174743384614, 395.3570687288378, 397.7314048998272, 401.82321795486706, 404.0332626054033, 407.71057787102535, 410.8532826778768, 413.2171912664151, 416.7607317825822, 419.06052107869925, 422.86073370963356, 425.56050725282137, 428.7289416032822, 431.5868951700548, 434.1638202990803, 436.65716705553865, 439.50293173686146, 442.7347796873833, 446.8874007476432, 450.2403780621976, 453.29198720792925, 456.0342974460868, 458.6153037837133, 461.56380177647026, 464.20368693901656, 466.691966383324, 468.8079649404163, 470.79128780439, 474.71462438964727, 480.82951675071496, 484.71160080693284, 486.3215064902958, 488.4197692356378, 491.17013007637155, 493.7507269683708, 497.34286854916746, 500.98519709660104, 504.03977938824954, 507.32635817465786, 507.88656427042935, 510.3501966220603, 515.529828721657, 518.4690906143065, 522.5469982253443, 525.9251599750156, 529.0499904487481, 531.2864505203889, 531.2864505203889, 533.6091546781496, 537.0521007260379, 542.9583993555486, 547.9147594698899, 549.3299646609445, 551.9453014974921, 552.7488398237151, 556.3564169971266, 560.404652390278, 564.3962059348734, 566.6363629326391, 568.2952225456394, 573.419827775658, 575.2943236050885, 577.9137045909403, 579.5219994947479, 582.6703551290922, 585.811132164887, 587.002196452197, 588.7117097061853, 597.1677915531716, 603.616630778402, 605.8892249214309, 605.8892249214309, 609.1046460364297, 611.8486536654434, 615.0679077437967, 620.194133245999, 621.8600314133139, 623.968212462276, 626.0923250498265, 627.6209614329342, 632.6670022394472, 635.5716404015416, 635.5716404015416, 638.6786179906055, 638.6786179906055, 655.41120125659, 655.41120125659, 655.41120125659, 655.41120125659, 655.41120125659, 695.7250727971176, 695.7250727971176, 695.7250727971176, 695.7250727971176, 695.7250727971176, 709.3739142789392, 709.3739142789392], [0.19874987125399457, 0.19874987125399457, 0.9521426963806496, 0.8006249094009661, 1.6907692520435376, 1.7642308301192402, 2.350312603712112, 2.2414706746269997, 2.748181938691609, 2.7780556133058667, 2.688913087430256, 2.341666673478677, 2.841500094890665, 2.2847223758698254, 3.3002382923308677, 3.343421368849979, 3.795000177684666, 4.087631888640622, 5.191764917934707, 4.0015001735687745, 4.032083484729189, 4.266818365039299, 3.990277844005211, 3.6926250724792964, 3.3643750548363185, 4.71762500810627, 4.445694509612224, 4.6560295264861455, 5.974305756886833, 7.1307895416962905, 5.953421160296434, 6.503255742317069, 7.247659445011915, 7.101829061275621, 6.000238108407933, 7.288690188952912, 6.816749813079897, 6.311710226661352, 5.820768896494224, 5.768749594688475, 4.7344998519421155, 4.979729584101137, 4.531875023841906, 4.049387782544556, 4.173478220234798, 4.5065383408620265, 4.373301750579018, 4.107550868793461, 4.328979446060852, 5.25413440484273, 4.529545350508318, 3.840581291331858, 3.9733332320622567, 3.6045236612502114, 3.028026209630501, 3.1547434584300094, 3.4106249117851632, 4.508301715311057, 4.449326716203003, 5.298017125951875, 5.233571282038658, 5.739782429847211, 5.684090696681675, 6.069907282016947, 5.335760672818117, 6.209230559972687, 6.1410575383443735, 6.028770524869261, 6.068749951124243, 6.530634935999724, 5.9370689873038085, 5.880999986171768, 5.3287498784065725, 5.275172360847789, 4.80827577722489, 6.261636389125531, 6.944999930063949, 9.580094021311416, 8.93279631808659, 8.784416421413482, 8.474296696037108, 7.653880537659424, 5.450396934766658, 6.830000030469133, 6.047063415542528, 5.045727090575534, 5.633727015495352, 5.640188473035761, 5.670095804654702, 6.296914586513674, 6.4793133381302, 6.369693621518745, 6.363272473508743, 5.756481468677588, 4.742803092436429, 4.427361328469448, 3.96246053286967, 4.924914170626995, 4.436885598448513, 5.482041082187621, 5.464130484539589, 5.316603758110245, 4.289791411956204, 4.947346590587126, 4.808673295196283, 4.787558070116315, 5.472976217042891, 5.477500085416108, 4.874418763227302, 6.245777626885365, 7.084090676090992, 7.598604113778466, 7.273295027559509, 7.045760571853031, 5.70110001506811, 5.344509843377554, 4.630454867102929, 5.130686536676789, 6.326196023070343, 6.625139205190913, 7.041805818875679, 8.062575722752234, 8.348499809742025, 6.470908668836002, 6.449078555734586, 5.323242851592455, 4.387999736785946, 3.541249947547962, 3.5860000212533643, 3.6169445180893356, 3.969375182092225, 4.5057353827532935, 5.475625167191069, 5.2980883034538335, 6.442413902940479, 6.222187623083645, 6.338666921933543, 6.241551741896054, 6.456833218574559, 5.3751559740305375, 7.252902918784901, 6.638147801646538, 7.264814521648317, 7.851666663487819, 7.5741307723004585, 5.50934817314156, 4.700178804397654, 3.485345097245911, 5.772343835532724, 5.425302892453809, 5.760142698832976, 6.4788332195282665, 7.825384631156997, 4.836750075340359, 5.104062861800273, 5.515333927154651, 7.567778353161391, 6.241389353010498, 6.722143275851382, 6.338478587606751, 5.522222394413474, 3.3517857960292594, 3.8855884310778825, 3.388437778353733, 4.637917126814567, 4.819091470891821, 4.493636822267046, 6.3188889694213985, 6.5474999650319505, 4.869999938011176, 7.556785623005486, 6.604705870572301, 5.575333408991544, 5.575333408991544, 6.049642844200191, 2.0933334891002082, 3.225999893188532, 3.9893747758865885, 5.003999863624619, 5.384444553587223, 5.9350000226498025, 6.548500012397824, 7.614285830089086, 5.428571496691104, 4.230714124952117, 4.971666545867985, 2.514444289737336, 2.1099999067519093, 2.0831248831749463, 3.290624678134975, 2.6307140159607507, 4.351249771118205, 5.7999996821085915, 8.449999523162887, 3.632500391006488, 0.045000057220477174, 0.9249996662140347, 0.9249996662140347, 0.7516663297017961, 1.1049994659424556, 1.5433328374227433, 1.0016664886475155, 4.396250638961874, 5.726667633056725, 8.646250865459507, 10.721667960484902, 14.343334252039654, 14.224999833107006, 14.224999833107006, 11.04499910354616, 11.04499910354616, 11.519998626709025, 11.519998626709025, 11.519998626709025, 11.519998626709025, 11.519998626709025, 7.220000724792499, 7.220000724792499, 7.220000724792499, 7.220000724792499, 7.220000724792499, 13.005001945495678, 13.005001945495678], [4.0, 4.0, 7.0, 8.0, 13.0, 13.0, 16.0, 17.0, 22.0, 18.0, 23.0, 21.0, 20.0, 18.0, 21.0, 19.0, 19.0, 19.0, 17.0, 20.0, 24.0, 33.0, 36.0, 40.0, 40.0, 40.0, 36.0, 34.0, 36.0, 38.0, 38.0, 43.0, 47.0, 41.0, 42.0, 42.0, 40.0, 38.0, 39.0, 40.0, 40.0, 37.0, 40.0, 49.0, 46.0, 52.0, 53.0, 49.0, 49.0, 52.0, 44.0, 43.0, 42.0, 42.0, 38.0, 39.0, 40.0, 53.0, 52.0, 58.0, 63.0, 69.0, 55.0, 54.0, 46.0, 52.0, 52.0, 61.0, 60.0, 63.0, 58.0, 60.0, 60.0, 58.0, 58.0, 55.0, 51.0, 53.0, 59.0, 60.0, 64.0, 67.0, 63.0, 59.0, 63.0, 55.0, 55.0, 53.0, 52.0, 47.0, 51.0, 49.0, 55.0, 54.0, 66.0, 72.0, 63.0, 58.0, 61.0, 49.0, 46.0, 53.0, 48.0, 49.0, 49.0, 43.0, 42.0, 46.0, 43.0, 45.0, 44.0, 43.0, 44.0, 46.0, 50.0, 51.0, 55.0, 51.0, 46.0, 36.0, 36.0, 33.0, 30.0, 33.0, 38.0, 37.0, 35.0, 40.0, 35.0, 36.0, 32.0, 34.0, 32.0, 34.0, 29.0, 32.0, 30.0, 29.0, 30.0, 32.0, 31.0, 27.0, 27.0, 24.0, 23.0, 23.0, 28.0, 29.0, 32.0, 33.0, 35.0, 30.0, 26.0, 20.0, 16.0, 15.0, 18.0, 18.0, 21.0, 23.0, 18.0, 14.0, 17.0, 16.0, 12.0, 11.0, 11.0, 9.0, 6.0, 8.0, 14.0, 17.0, 15.0, 15.0, 14.0, 9.0, 5.0, 8.0, 10.0, 9.0, 8.0, 10.0, 7.0, 7.0, 7.0, 9.0, 9.0, 9.0, 8.0, 8.0, 7.0, 4.0, 3.0, 2.0, 2.0, 1.0, 2.0, 2.0, 3.0, 2.0, 3.0, 3.0, 4.0, 3.0, 4.0, 3.0, 3.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

Code for the figure

f = Figure()
ax = Axis(f[1, 1]; xlabel = "Distance (km)", ylabel = "Semivariogram")
hlines!(ax, [var(Lᵢ)]; color = :grey40, linestyle = :dash)
scatter!(
    ax,
    x,
    y;
    markersize = n ./ maximum(n) .* 8 .+ 2,
    color = :grey55,
)
xlims!(ax; low = 0.0)
ylims!(ax, 0.0, quantile(y, 0.95))

note that with an unexported function

vario = SpeciesDistributionToolkit.fitvariogram(x, y, n; family = :gaussian)

(sill = 5.658900323073887, nugget = 0.6785190841053046, range = 110.44888993154616, error = 1.3313725982256805, model = SpeciesDistributionToolkit.var"#__mod#__variogram_gaussian##0"{Float64, Float64, Float64}(5.658900323073887, 0.6785190841053046, 110.44888993154616))

Code for the figure

vx = LinRange(extrema(x)..., 100)
lines!(ax, vx, vario.model.(vx); color = :orange, linewidth = 2, alpha = 0.6)

this gives a range of approx (in km)

vario.range

110.44888993154616

we can also check the sill

vario.sill

5.658900323073887

and compare to layer variance

var(Lᵢ)

5.451709f0

and check the base amount of autocorrelation for close points

vario.nugget

0.6785190841053046

we tile the model observations - ideally would be done for a surface corresponding to the range in the variogram but here, too large

tiles = tessellate(model, 20.0; tile = :hexagons, pointy = true)

FeatureCollection with 109 features, each with 4 properties

Code for the figure

f = Figure()
ax = Axis(f[1, 1]; aspect = DataAspect())
poly!(ax, landmass; color = :grey98)
poly!(ax, intersect(tiles, landmass); color = :darkgreen, alpha = 0.1)
lines!(ax, intersect(tiles, landmass); color = :grey60)
scatter!(ax, presences(model); color = :darkgreen)
scatter!(ax, absences(model); color = :grey20, markersize = 4)
lines!(ax, landmass; color = :black)

we now assign the folds for cross validation

assignfolds!(tiles; n = 6, order = :balanced)

FeatureCollection with 109 features, each with 6 properties

check what the folds look like

Code for the figure

f = Figure()
ax = Axis(f[1, 1]; aspect = DataAspect())
colpal = Makie.wong_colors()[1:maximum(uniqueproperties(tiles)["__fold"])]
for k in uniqueproperties(tiles)["__fold"]
    poly!(ax, tiles["__fold" => k]; color = colpal[k], alpha = 0.3)
    for f in tiles["__fold" => k].features
        text!(
            ax,
            f.properties["__centroid"];
            text = string(f.properties["__fold"]),
            align = (:center, :center),
            color = colpal[k],
        )
    end
end
lines!(ax, tiles; color = :grey40)

this turns into folds

blockfolder = spatialfold(tiles)
folds = blockfolder(model)

6-element Vector{Tuple{Vector{Int64}, Vector{Int64}}}:
 ([1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 13, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 39, 40, 41, 42, 43, 44, 45, 46, 47, 49, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 63, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 77, 78, 79, 80, 82, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 102, 104, 105, 106, 110, 111, 114, 115, 116, 117, 118, 120, 121, 122, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 202, 203, 204, 205, 206, 207, 210, 211, 212, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255], [9, 12, 14, 35, 36, 37, 38, 48, 50, 62, 64, 75, 76, 81, 83, 84, 101, 103, 107, 108, 109, 112, 113, 119, 123, 124, 141, 152, 201, 208, 209, 213, 214, 215])
 ([1, 2, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 42, 43, 44, 45, 46, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 59, 60, 61, 62, 63, 64, 65, 66, 67, 69, 70, 71, 72, 73, 75, 76, 78, 79, 81, 83, 84, 86, 87, 88, 89, 90, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 130, 131, 132, 133, 134, 136, 137, 138, 140, 141, 142, 143, 146, 147, 149, 150, 152, 153, 154, 155, 156, 157, 159, 160, 161, 162, 163, 165, 167, 168, 169, 170, 171, 173, 174, 176, 177, 178, 179, 181, 182, 184, 186, 187, 188, 190, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 207, 208, 209, 210, 213, 214, 215, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255], [3, 4, 39, 40, 41, 47, 58, 68, 74, 77, 80, 82, 85, 91, 129, 135, 139, 144, 145, 148, 151, 158, 164, 166, 172, 175, 180, 183, 185, 189, 191, 206, 211, 212, 216, 217, 218, 219, 236])
 ([1, 2, 3, 4, 5, 6, 7, 9, 11, 12, 13, 14, 15, 16, 17, 20, 21, 23, 25, 26, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 41, 42, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 58, 59, 60, 61, 62, 64, 68, 72, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 93, 95, 97, 98, 100, 101, 102, 103, 104, 105, 107, 108, 109, 112, 113, 117, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 141, 142, 144, 145, 146, 148, 151, 152, 153, 155, 156, 157, 158, 160, 162, 163, 164, 166, 169, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 200, 201, 202, 203, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255], [8, 10, 18, 19, 22, 24, 27, 28, 29, 33, 43, 44, 56, 57, 63, 65, 66, 67, 69, 70, 71, 73, 92, 94, 96, 99, 106, 110, 111, 114, 115, 116, 118, 140, 143, 147, 149, 150, 154, 159, 161, 165, 167, 168, 170, 198, 199, 204, 205])
 ([1, 2, 3, 4, 5, 8, 9, 10, 12, 14, 18, 19, 21, 22, 23, 24, 25, 26, 27, 28, 29, 31, 33, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 88, 89, 91, 92, 93, 94, 95, 96, 97, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 147, 148, 149, 150, 151, 152, 154, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 180, 182, 183, 184, 185, 187, 189, 191, 192, 193, 196, 198, 199, 201, 204, 205, 206, 208, 209, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 222, 223, 224, 225, 226, 229, 231, 232, 233, 235, 236, 237, 239, 240, 245, 249, 250, 251, 252, 253, 254, 255], [6, 7, 11, 13, 15, 16, 17, 20, 30, 32, 34, 52, 53, 54, 86, 87, 90, 98, 146, 153, 155, 179, 181, 186, 188, 190, 194, 195, 197, 200, 202, 203, 207, 210, 221, 227, 228, 230, 234, 238, 241, 242, 243, 244, 246, 247, 248])
 ([3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 60, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 80, 81, 82, 83, 84, 85, 86, 87, 90, 91, 92, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 118, 119, 120, 121, 122, 123, 124, 125, 128, 129, 132, 134, 135, 136, 137, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 158, 159, 161, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 179, 180, 181, 182, 183, 185, 186, 188, 189, 190, 191, 194, 195, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 232, 233, 234, 235, 236, 237, 238, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253], [1, 2, 21, 59, 61, 72, 78, 79, 88, 89, 95, 100, 117, 126, 127, 130, 131, 133, 138, 156, 157, 160, 162, 163, 177, 178, 184, 187, 192, 193, 196, 220, 231, 239, 254, 255])
 ([1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 27, 28, 29, 30, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43, 44, 47, 48, 50, 52, 53, 54, 56, 57, 58, 59, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 94, 95, 96, 98, 99, 100, 101, 103, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 123, 124, 126, 127, 129, 130, 131, 133, 135, 138, 139, 140, 141, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 170, 172, 175, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 227, 228, 230, 231, 234, 236, 238, 239, 241, 242, 243, 244, 246, 247, 248, 254, 255], [5, 23, 25, 26, 31, 42, 45, 46, 49, 51, 55, 60, 93, 97, 102, 104, 105, 120, 121, 122, 125, 128, 132, 134, 136, 137, 142, 169, 171, 173, 174, 176, 182, 222, 223, 224, 225, 226, 229, 232, 233, 235, 237, 240, 245, 249, 250, 251, 252, 253])

now we can cross-validate

spatial_cv = crossvalidate(model, folds)

(validation = ConfusionMatrix[(tp: 14, fp: 2; fn: 2, tn: 16), (tp: 9, fp: 1; fn: 4, tn: 25), (tp: 19, fp: 4; fn: 3, tn: 23), (tp: 10, fp: 5; fn: 4, tn: 28), (tp: 4, fp: 3; fn: 4, tn: 25), (tp: 9, fp: 1; fn: 3, tn: 37)], training = ConfusionMatrix[(tp: 60, fp: 13; fn: 9, tn: 139), (tp: 64, fp: 12; fn: 8, tn: 132), (tp: 53, fp: 10; fn: 10, tn: 133), (tp: 60, fp: 5; fn: 11, tn: 132), (tp: 67, fp: 10; fn: 10, tn: 132), (tp: 63, fp: 12; fn: 10, tn: 120)])

we get the value of the first variable by fold

fold_mean = byzone(mean, Lᵢ, tiles, "__fold")

Dict{Int64, Float32} with 6 entries:
  5 => 8.34475
  4 => 8.86887
  6 => 8.0143
  2 => 8.17922
  3 => 8.75059
  1 => 9.34716

look at the results

Code for the figure

f = Figure()
ax = Axis(f[1, 1]; ylabel = "MCC", xlabel = layers(provider)[1])
ax2 = Axis(f[1, 2])
_valstyle = (color = :white, strokewidth = 1, strokecolor = :darkgreen, markersize = 18)
_trnstyle = (
    color = :grey90,
    strokewidth = 1,
    strokecolor = :orange,
    markersize = 14,
    marker = :rect,
)
X = [fold_mean[i] for i in eachindex(folds)]
scatter!(ax, X, mcc.(spatial_cv.training); _trnstyle..., label = "Training")
scatter!(ax, X, mcc.(spatial_cv.validation); _valstyle..., label = "Validation")
errorbars!(
    ax2,
    [1],
    [mcc(spatial_cv.validation)],
    [0.5 * ci(spatial_cv.validation)];
    whiskerwidth = 12,
    color = :darkgreen,
)
scatter!(ax2, [1], [mcc(spatial_cv.validation)]; _valstyle...)
errorbars!(
    ax2,
    [2],
    [mcc(spatial_cv.training)],
    [0.5 * ci(spatial_cv.training)];
    whiskerwidth = 12,
    color = :orange,
)
scatter!(ax2, [2], [mcc(spatial_cv.training)]; _trnstyle...)
xlims!(ax2, 0, 3)
hidexdecorations!(ax2)
ylims!(ax, 0.5, 1)
ylims!(ax2, 0.5, 1)
colsize!(f.layout, 2, Relative(0.3))
axislegend(ax; position = :lt)

now we can look at output

now we can check variables with these folds

variables!(model, ForwardSelection, folds; included = [1])

☑️  PCATransform → Logistic → P(x) ≥ 0.234 🗺️

beter model?

selection_cv = crossvalidate(model, folds)

(validation = ConfusionMatrix[(tp: 16, fp: 3; fn: 0, tn: 15), (tp: 10, fp: 1; fn: 3, tn: 25), (tp: 19, fp: 4; fn: 3, tn: 23), (tp: 14, fp: 5; fn: 0, tn: 28), (tp: 7, fp: 4; fn: 1, tn: 24), (tp: 12, fp: 3; fn: 0, tn: 35)], training = ConfusionMatrix[(tp: 63, fp: 16; fn: 6, tn: 136), (tp: 68, fp: 16; fn: 4, tn: 128), (tp: 61, fp: 15; fn: 2, tn: 128), (tp: 65, fp: 12; fn: 6, tn: 125), (tp: 72, fp: 13; fn: 5, tn: 129), (tp: 68, fp: 15; fn: 5, tn: 117)])

Code for the figure

f = Figure()
ax = Axis(f[1, 1]; ylabel = "MCC", xlabel = layers(provider)[1])
ax2 = Axis(f[1, 2])
_valstyle = (color = :white, strokewidth = 1, strokecolor = :darkgreen, markersize = 18)
_trnstyle = (
    color = :grey90,
    strokewidth = 1,
    strokecolor = :orange,
    markersize = 14,
    marker = :rect,
)
X = [fold_mean[i] for i in eachindex(folds)]
scatter!(ax, X, mcc.(selection_cv.training); _trnstyle..., label = "Training")
scatter!(ax, X, mcc.(selection_cv.validation); _valstyle..., label = "Validation")
errorbars!(
    ax2,
    [1],
    [mcc(selection_cv.validation)],
    [0.5 * ci(selection_cv.validation)];
    whiskerwidth = 12,
    color = :darkgreen,
)
scatter!(ax2, [1], [mcc(selection_cv.validation)]; _valstyle...)
errorbars!(
    ax2,
    [2],
    [mcc(selection_cv.training)],
    [0.5 * ci(selection_cv.training)];
    whiskerwidth = 12,
    color = :orange,
)
scatter!(ax2, [2], [mcc(selection_cv.training)]; _trnstyle...)
xlims!(ax2, 0, 3)
hidexdecorations!(ax2)
ylims!(ax, 0.5, 1)
ylims!(ax2, 0.5, 1)
colsize!(f.layout, 2, Relative(0.3))
axislegend(ax; position = :lt)

table

todo turn into PrettyTables

cv_res = [
    [
        i,
        mcc(spatial_cv.training[i]),
        mcc(spatial_cv.validation[i]),
        mcc(selection_cv.training[i]),
        mcc(selection_cv.validation[i]),
        mcc(selection_cv.validation[i]) >= mcc(spatial_cv.validation[i]) ? "↑" : "↓",
    ] for i in 1:length(spatial_cv.validation)
];

this is the result as a table - see that variable selection makes cross-validation better

pretty_table(
    permutedims(hcat(cv_res...));
    alignment = [:l, :c, :c, :c, :c, :c],
    backend = :markdown,
    column_labels = ["Fold", "Train.", "Val.", "Train. + sel.", "Val. + sel.", "Trend"],
    formatters = [fmt__printf("%5.3f", [2, 3, 4, 5])],
)

Fold	Train.	Val.	Train. + sel.	Val. + sel.	Trend
1	0.772	0.764	0.781	0.838	↑
2	0.795	0.706	0.806	0.766	↑
3	0.771	0.713	0.824	0.713	↑
4	0.827	0.552	0.813	0.791	↑
5	0.800	0.413	0.826	0.661	↑
6	0.768	0.773	0.798	0.858	↑

make spatial prediction

P = predict(model, L; threshold = false)

🗺️  A 517 × 990 layer with 411408 Float64 cells
   Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs

Code for the figure

f = Figure()
ax = Axis(f[1, 1]; aspect = DataAspect())
heatmap!(
    ax,
    predict(model, L; threshold = false);
    colormap = Reverse(:navia),
    colorrange = (0, 1),
)
scatter!(ax, presences(model); color = :transparent, strokecolor = :orange, strokewidth = 2)

tile by any shape we want, for example this is prediction aggregated in triangular cells

Code for the figure

f = Figure()
ax = Axis(f[1, 1]; aspect = DataAspect())
newtiles = tessellate(P, 20.0; tile = :triangles)
heatmap!(
    ax,
    mosaic(mean, P, newtiles, "__centroid");
    colormap = Reverse(:navia),
    colorrange = (0, 1),
)
lines!(ax, landmass; color = :black)