Preparations

library(caret) # load this first as it loads plyr
library(gbm)
library(car)
library(grid)
library(gridExtra)
library(randomForest)
library(xtable)
library(kableExtra)
library(abt)
library(tidyverse)
source('FISH_functions.R')
Raw fish data 
#fish = read.csv('data/Fish benthic physical sitelevel 2019.csv', strip.white=TRUE)
fish = read.csv('data/Selected fish benthic physical sitelevel 2019.csv', strip.white=TRUE)
fish %>% glimpse
## Rows: 619
## Columns: 302
## $ YEAR                        <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, …
## $ REGION                      <chr> "Keppel", "Keppel", "Keppel", "Keppel", "K…
## $ EXPOSURE                    <chr> "Semi-Exposed", "Exposed", "Semi-Exposed",…
## $ NTR                         <chr> "Fished", "Fished", "Fished", "Fished", "F…
## $ NTR.Pooled                  <chr> "Fished", "Fished", "Fished", "Fished", "F…
## $ SITE                        <chr> "GK1", "GK2", "GK3", "GK4", "GK7", "GK8", …
## $ lut.arge                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.carp                    <dbl> 13.3320, 11.3322, 4.6662, 8.6658, 16.6650,…
## $ lut.flma                    <dbl> 0.0000, 3.3330, 0.0000, 0.0000, 7.3326, 0.…
## $ lut.fulv                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.kas                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.lemn                    <dbl> 0.0000, 0.6666, 0.0000, 0.0000, 0.0000, 0.…
## $ lut.lutj                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.mono                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.quin                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ lut.rivu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.russ                    <dbl> 67.3266, 4.6662, 0.0000, 5.3328, 116.6550,…
## $ lut.seba                    <dbl> 1.9998, 0.6666, 9.9990, 0.0000, 0.0000, 0.…
## $ lut.vitt                    <dbl> 377.9622, 10.6656, 16.6650, 9.3324, 0.0000…
## $ mac.macu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pagrus                      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sym.nema                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.barb                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ par.boides                  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.bifa                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.indi                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ par.multi                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.cili                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ upe.trag                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ platax                      <dbl> 0.0000, 0.6666, 0.0000, 0.6666, 0.6666, 0.…
## $ Platycephalus.spp.          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.atki                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ let.hara                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.lati                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ let.lent                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ let.mini                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.nebu                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ let.oliv                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.obso                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.orna                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ gym.spp                     <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 5.…
## $ mon.gran                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sco.bili                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sco.marg                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ sco.mono                    <dbl> 4.6662, 1.9998, 3.9996, 5.3328, 3.9996, 3.…
## $ sco.line                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ any.leuc                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ath.roga                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cep.argu                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cep.boen                    <dbl> 21.9978, 5.9994, 7.3326, 8.6658, 3.3330, 0…
## $ cep.cyan                    <dbl> 1.3332, 0.0000, 1.9998, 1.3332, 0.0000, 0.…
## $ cep.micr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cro.alti                    <dbl> 0.6666, 0.0000, 0.6666, 0.6666, 0.0000, 0.…
## $ epi.caer                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.coio                    <dbl> 0.0000, 0.0000, 0.0000, 1.3332, 0.0000, 1.…
## $ epi.cora                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.fasc                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 0.0000, 0.…
## $ epi.fusc                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.howl                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.lanc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ epi.hexa                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.merr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.ongu                    <dbl> 0.6666, 0.0000, 0.6666, 0.0000, 0.0000, 0.…
## $ epi.quoy                    <dbl> 6.6660, 2.6664, 12.6654, 1.9998, 1.9998, 5…
## $ epi.sexf                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pms.laev                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pms.leop                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.6666, 0.…
## $ pms.macu                    <dbl> 5.9994, 3.9996, 17.3316, 7.3326, 9.9990, 3…
## $ pte.ante                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pte.voli                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ dip.bifa                    <dbl> 1.9998, 1.9998, 7.9992, 1.9998, 1.9998, 0.…
## $ dia.pict                    <dbl> 23.3310, 3.3330, 4.6662, 1.9998, 9.9990, 0…
## $ ple.chae                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ ple.flav                    <dbl> 1.9998, 1.3332, 2.6664, 0.6666, 0.0000, 0.…
## $ ple.gibb                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.less                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.line                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.picu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.unic                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ bodianus                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ oxy.diag                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ che.chlo                    <dbl> 1.3332, 0.6666, 5.9994, 0.6666, 1.3332, 1.…
## $ che.fasc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ che.tril                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ che.undu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cho.anch                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cho.ceph                    <dbl> 0.6666, 0.0000, 0.6666, 1.3332, 0.6666, 0.…
## $ cho.cya                     <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ cho.fasc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ cho.grap                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 1.3332, 0.…
## $ cho.scho                    <dbl> 1.3332, 0.0000, 0.6666, 0.6666, 0.0000, 0.…
## $ cho.vitt                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ epb.insi                    <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 1.3332, 0.…
## $ gom.vari                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ hem.melt                    <dbl> 17.3316, 9.3324, 25.9974, 13.9986, 25.3308…
## $ hem.fasc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ psa.wai                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.bloc                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 4.6662, 0.…
## $ aca.duss                    <dbl> 2.6664, 0.0000, 0.0000, 0.0000, 9.9990, 0.…
## $ aca.gram                    <dbl> 0.0000, 0.0000, 0.6666, 1.3332, 4.6662, 0.…
## $ aca.line                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.nans                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.ncus                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.nuda                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.xant                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ aca.thom                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.spp                     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cte.bino                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cte.stri                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.annu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.brach                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.brevi                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.litu                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.tong                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.unic                    <dbl> 0.0000, 1.9998, 0.0000, 0.0000, 0.0000, 0.…
## $ pri.spp                     <dbl> 0.0000, 0.0000, 0.0000, 3.9996, 18.6648, 0…
## $ kyp.spp                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ zeb.scop                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ zeb.veli                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ zan.corn                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ bol.muri                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cal.caro                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cet.bico                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.blee                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.micr                    <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ chs.sord                    <dbl> 10.6656, 2.6664, 0.0000, 0.0000, 0.0000, 0…
## $ hip.long                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.alti                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.cham                    <dbl> 6.6660, 0.6666, 0.0000, 0.0000, 0.0000, 0.…
## $ sca.dimi                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.flav                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.fors                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.fren                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ sca.ghob                    <dbl> 7.9992, 3.9996, 1.9998, 0.6666, 13.3320, 2…
## $ sca.glob                    <dbl> 0.6666, 7.9992, 0.0000, 0.0000, 0.6666, 0.…
## $ sca.nigr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.ovic                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.psit                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.rivu                    <dbl> 37.3296, 25.3308, 19.3314, 7.9992, 41.3292…
## $ sca.rubr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.schl                    <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ sca.spin                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.tric                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.spp                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.arge                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.cana                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.cora                    <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.0000, 0.…
## $ sig.doli                    <dbl> 5.3328, 3.9996, 1.3332, 2.6664, 6.6660, 6.…
## $ sig.fusc                    <dbl> 0.0000, 0.0000, 11.9988, 6.6660, 0.0000, 1…
## $ sig.javu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.line                    <dbl> 0.0000, 0.0000, 0.0000, 1.3332, 0.0000, 0.…
## $ sig.puel                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 0.0000, 0.…
## $ sig.ptus                    <dbl> 0.0000, 0.0000, 0.0000, 1.9998, 0.0000, 0.…
## $ sig..ptiss                  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.stell                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.vulp                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ mon.arge                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.auri                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ cha.afas                    <dbl> 63.3270, 89.9910, 39.3294, 57.3276, 51.328…
## $ cha.baro                    <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.0000, 0.…
## $ cha.benn                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.ephi                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ cha.citr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.flav                    <dbl> 0.0000, 0.0000, 0.0000, 0.6666, 0.0000, 0.…
## $ cha.line                    <dbl> 5.9994, 5.3328, 0.0000, 1.3332, 5.3328, 0.…
## $ cha.lunu                    <dbl> 1.3332, 0.0000, 0.6666, 1.3332, 0.0000, 0.…
## $ cha.ttus                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.mela                    <dbl> 4.6662, 1.3332, 1.9998, 0.0000, 0.6666, 0.…
## $ cha.ocel                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.orna                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.pleb                    <dbl> 0.0000, 0.0000, 0.6666, 0.6666, 0.0000, 0.…
## $ cha.rain                    <dbl> 3.3330, 13.9986, 1.3332, 3.3330, 3.9996, 0…
## $ cha.raff                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.seme                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.spec                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ cha.tlis                    <dbl> 1.3332, 2.6664, 0.0000, 0.0000, 0.0000, 0.…
## $ cha.ulie                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ cha.vaga                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chm.rost                    <dbl> 7.3326, 3.3330, 1.9998, 1.3332, 6.6660, 1.…
## $ cor.alti                    <dbl> 3.9996, 1.9998, 0.0000, 0.6666, 1.9998, 0.…
## $ cor.chry                    <dbl> 0.0000, 1.3332, 0.0000, 0.6666, 1.9998, 0.…
## $ hen.acum                    <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.0000, 0.…
## $ hen.mono                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ hen.vari                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ mct.stri                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.…
## $ prc.ocel                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.bico                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.bisp                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.nox                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.tibi                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.vrol                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chd.doub                    <dbl> 5.9994, 3.9996, 0.6666, 1.3332, 7.9992, 0.…
## $ chd.mere                    <dbl> 5.3328, 5.3328, 1.9998, 2.6664, 4.6662, 0.…
## $ pmc.impe                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pmc.semi                    <dbl> 0.0000, 3.3330, 0.0000, 0.0000, 0.6666, 0.…
## $ pmc.sext                    <dbl> 1.3332, 1.3332, 0.0000, 0.6666, 0.6666, 0.…
## $ pmc.xant                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pyg.diac                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ abu.spp                     <int> 406, 48, 0, 0, 238, 0, 30, 112, 38, 12, 4,…
## $ acn.poly                    <int> 90, 48, 74, 58, 124, 22, 152, 78, 62, 80, …
## $ amb.aure                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amb.cura                    <int> 0, 0, 4, 2, 0, 0, 0, 0, 0, 2, 0, 0, 6, 0, …
## $ amb.leuc                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amp.spp                     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amp.peri                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pre.spp                     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.labi                    <int> 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.ambo                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.alis                    <int> 60, 0, 12, 0, 0, 0, 138, 34, 16, 64, 0, 0,…
## $ chr.apes                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.marg                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.niti                    <int> 8920, 9680, 7340, 15840, 12210, 240, 45860…
## $ chr.retr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.tern                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.webe                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.xant                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.roll                    <int> 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.talb                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.flav                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ das.spp                     <int> 12, 0, 0, 0, 0, 0, 0, 0, 0, 24, 4, 0, 0, 4…
## $ das.reti                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ das.trim                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ dis.spp                     <int> 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ hgy.plag                    <int> 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0,…
## $ neg.mela                    <int> 14, 6, 2, 4, 12, 0, 4, 0, 8, 0, 0, 2, 0, 0…
## $ parma                       <int> 10, 14, 0, 0, 4, 2, 2, 2, 4, 0, 4, 0, 0, 0…
## $ neg.nigr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pgy.dick                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pgy.lacr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.adel                    <int> 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.ambo                    <int> 0, 0, 8, 0, 0, 0, 18, 6, 4, 32, 0, 0, 0, 0…
## $ pom.aust                    <int> 58, 140, 14, 54, 28, 0, 216, 156, 210, 360…
## $ pom.bank                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.brac                    <int> 84, 0, 4, 6, 52, 0, 90, 14, 10, 0, 6, 0, 1…
## $ pom.chry                    <int> 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.coel                    <int> 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.gram                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.lepi                    <int> 8, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.molu                    <int> 602, 294, 426, 634, 542, 190, 564, 198, 22…
## $ pom.naga                    <int> 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.phil                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.reid                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.trip                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.vaiu                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.ward                    <int> 258, 224, 206, 338, 200, 456, 136, 156, 18…
## $ ste.apic                    <int> 2, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, …
## $ ste.fasc                    <int> 8, 8, 4, 10, 8, 12, 0, 4, 8, 0, 4, 6, 0, 0…
## $ ste.livi                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ste.nigr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ oxy.long                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Anampses                    <int> 0, 6, 0, 2, 8, 0, 0, 0, 0, 0, 6, 2, 8, 0, …
## $ Coris                       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Halichoeres                 <int> 14, 10, 20, 30, 8, 8, 26, 6, 10, 48, 16, 4…
## $ Labroides                   <int> 18, 10, 10, 12, 16, 4, 10, 10, 4, 2, 4, 6,…
## $ Labropsis                   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Psuedolabrus.guentheri      <int> 20, 58, 20, 24, 26, 104, 16, 22, 20, 86, 3…
## $ Stethojoulis                <int> 12, 4, 16, 28, 2, 30, 16, 2, 0, 84, 6, 0, …
## $ Thalasoma                   <int> 24, 8, 12, 12, 24, 2, 30, 10, 16, 18, 12, …
## $ Macropharyngodon            <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Labrichthys                 <int> 8, 0, 0, 6, 0, 2, 0, 2, 4, 4, 4, 0, 4, 0, …
## $ Total.fish.density          <dbl> 11381.9272, 10805.9752, 8401.3122, 17230.6…
## $ Total.fish.species.richness <dbl> 35.2, 29.4, 27.6, 26.4, 29.2, 17.8, 33.2, …
## $ Prey.density                <int> 10654, 10558, 8190, 17060, 13502, 1076, 47…
## $ Prey.biomass                <dbl> 44.760041, 30.942239, 23.789262, 46.075893…
## $ Plectropomus.total.density  <dbl> 5.9994, 3.9996, 17.3316, 7.3326, 10.6656, …
## $ Plectropomus.total.biomass  <dbl> 3.980542, 5.371889, 9.546584, 8.926352, 6.…
## $ Plectropomus.legal.density  <dbl> 1.9998, 2.6664, 4.6662, 3.9996, 1.9998, 0.…
## $ Plectropomus.legal.biomass  <dbl> 2.4194697, 4.5827655, 5.0227904, 7.1055744…
## $ BE                          <dbl> 121.9966, 119.9976, 114.6630, 129.3312, 12…
## $ GRAZ                        <dbl> 72.6594, 47.9952, 36.6630, 26.6640, 99.990…
## $ GRAZ2                       <dbl> 7.9992, 7.3326, 15.3318, 17.9982, 44.6622,…
## $ Parrot                      <dbl> 64.6602, 40.6626, 21.3312, 8.6658, 55.3278…
## $ COR                         <dbl> 103.3246, 119.3214, 45.9954, 71.9934, 67.9…
## $ OM                          <dbl> 499.9996, 99.9996, 74.0000, 59.9998, 366.6…
## $ PL                          <dbl> 9000, 9680, 7366, 15842, 12210, 240, 45998…
## $ CA                          <dbl> 521.9478, 45.9954, 65.9934, 43.9956, 159.9…
## $ PI                          <dbl> 7.9992, 6.6660, 25.9974, 10.6656, 12.6654,…
## $ FA                          <int> 290, 246, 218, 348, 212, 474, 138, 164, 19…
## $ slope                       <dbl> 2.60, 2.64, 2.00, 2.40, 1.96, 2.20, 2.20, …
## $ rugosity                    <dbl> 2.68, 3.08, 3.00, 3.00, 3.00, 2.68, 3.00, …
## $ SCI                         <dbl> 7.040, 8.088, 6.000, 7.200, 5.880, 5.912, …
## $ LCC_.                       <dbl> 60.0, 78.8, 77.6, 77.6, 57.2, 13.6, 57.2, …
## $ LHC_.                       <dbl> 56.4, 68.8, 77.6, 75.6, 55.6, 13.6, 57.2, …
## $ SC_.                        <dbl> 3.6, 10.0, 0.0, 2.0, 1.6, 0.0, 0.0, 4.0, 2…
## $ MAC_.                       <dbl> 2.4, 0.4, 9.6, 4.4, 3.6, 4.0, 29.2, 0.8, 6…
## $ SPO_.                       <dbl> 1.6, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.…
## $ Turf_.                      <dbl> 26.0, 17.6, 4.4, 11.6, 28.0, 71.2, 6.4, 8.…
## $ Unconsolidated_.            <dbl> 9.6, 3.2, 8.4, 6.0, 11.2, 11.2, 7.2, 10.8,…
## $ Benthic.richness            <dbl> 6.6, 5.8, 2.2, 6.0, 3.8, 2.0, 3.4, 5.0, 6.…
## $ Coral_Morph.Diversity       <dbl> 5.4, 5.6, 1.2, 5.0, 3.4, 1.4, 2.4, 4.8, 5.…
## $ ChlA                        <dbl> 0.45221510, 0.45221510, 0.41940367, 0.4194…
## $ kd490                       <dbl> 0.06519040, 0.06519040, 0.06261333, 0.0626…
## $ SSTmean                     <chr> "23.80386083", "23.80357108", "23.79124967…
## $ SSTanom                     <dbl> -0.3130472, -0.3140394, -0.3159858, -0.316…
## $ wave.exposure.index         <dbl> 0.007256172, 0.160425811, 0.000990070, 0.0…
## $ Corrected.depth             <dbl> 4.03, 6.39, 5.12, 5.47, 3.73, 3.02, 5.28, …
## $ maxDHW                      <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
## $ deltaLHC                    <dbl> 14.0, 7.6, 6.4, 16.8, 8.0, 2.8, 4.8, 4.0, …
## $ deltaMA                     <dbl> -6.8, -1.2, -3.6, 0.4, -6.4, -41.2, 12.0, …
## $ Cyclone                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ LT.Fprimary                 <dbl> 0.087, 0.087, 0.087, 0.031, 0.025, 0.098, …
## $ Exposure.to.primary.weeks   <int> 2, 2, 2, 2, 3, 2, 3, 1, 1, 0, 3, 2, 2, 2, …
pco <- read.csv('data/Fish 2007-2019 PCO scores_allRegions.csv',  strip.white=TRUE)
pco %>% glimpse
## Rows: 99
## Columns: 112
## $ Island.group     <chr> "Keppels", "Keppels", "Keppels", "Keppels", "Keppels"…
## $ Site             <chr> "GK1", "GK2", "GK3", "GK4", "GK7", "GK8", "H3", "MI1"…
## $ PCO1             <dbl> 19.77982, 28.47255, 27.97316, 27.58896, 28.28861, 36.…
## $ PCO2             <dbl> 10.288750, 7.923383, 12.571452, 10.753819, 10.836750,…
## $ PCO3             <dbl> 5.31666080, 10.34580534, 3.21687927, 3.05770960, 7.32…
## $ PCO4             <dbl> 10.0053714, 2.9547334, 4.1084158, 5.7423832, 3.406272…
## $ PCO5             <dbl> 12.3337640, 6.5326895, -5.7922550, 2.5072552, 4.88353…
## $ PCO6             <dbl> 1.9316235, -11.8438838, 5.5197674, 4.5707637, -5.0634…
## $ PCO7             <dbl> 8.91222617, -1.34000339, 0.23260096, -0.46350746, 2.3…
## $ PCO8             <dbl> -0.5368321, -1.7728024, 1.7431962, 3.8489396, -0.1890…
## $ PCO9             <dbl> -1.48115610, 4.23348641, 0.04026016, 0.58709933, 3.86…
## $ PCO10            <dbl> -3.2738567, -6.7162154, 4.9418862, 0.7882569, 0.17374…
## $ PCO11            <dbl> 2.04980026, -0.59878395, -6.22277208, -9.97847909, -0…
## $ PCO12            <dbl> -0.3160799, -3.2929525, 1.7593295, -1.0537705, 1.4487…
## $ PCO13            <dbl> -1.8280337, -7.2551437, 0.4875337, -0.7849804, -5.143…
## $ PCO14            <dbl> -6.475337290, 3.881576192, -0.007853996, -1.823631756…
## $ PCO15            <dbl> 4.24610215, -3.08374848, -4.53975405, -0.26258604, 4.…
## $ PCO16            <dbl> 0.24046658, -0.48257590, 1.36665847, 3.79776412, -4.2…
## $ PCO17            <dbl> -1.23120787, -1.18489406, -7.00721119, 0.97583163, -0…
## $ PCO18            <dbl> 0.05306787, 2.41681891, -1.31546008, 1.85294338, 3.06…
## $ PCO19            <dbl> 0.3952058, 0.9122742, 0.1451317, 3.0856467, 0.6180767…
## $ PCO20            <dbl> -5.0232258, -3.6597008, 0.5088019, -0.7955180, 2.8493…
## $ PCO21            <dbl> -3.4147369, 3.9898593, -1.9365749, 0.2374281, -1.6786…
## $ PCO22            <dbl> 0.07697152, -2.15068026, 5.57267153, -1.06704713, 3.1…
## $ PCO23            <dbl> -4.4766509, -1.7161279, 0.8543020, -2.9526822, 1.9117…
## $ PCO24            <dbl> -1.4285919, -2.7013037, -1.7825005, -2.6839986, 0.812…
## $ PCO25            <dbl> -1.88682848, -2.11097870, -0.02824356, -2.89245494, 6…
## $ PCO26            <dbl> -0.36028133, 1.03192967, 3.91673251, -2.21306156, -4.…
## $ PCO27            <dbl> -0.6826462, -1.9152682, -2.8978772, -3.8511504, -0.39…
## $ PCO28            <dbl> 0.37826361, 5.65860553, -0.50546339, -0.80182794, 1.7…
## $ PCO29            <dbl> -0.21514192, -2.77971241, 0.47265410, -0.42124574, 0.…
## $ PCO30            <dbl> 2.63644764, -1.47104898, -0.06829036, 1.40988309, -0.…
## $ PCO31            <dbl> -1.53276514, -0.06413772, 1.54428936, -0.16163066, 3.…
## $ PCO32            <dbl> 0.2978719, -1.1760474, 2.4927713, 0.2973524, 3.646646…
## $ PCO33            <dbl> 1.4275302, -0.6607154, -2.8813920, -1.5245984, 1.3078…
## $ PCO34            <dbl> -1.657503078, 0.387159574, -1.678408256, 2.719817186,…
## $ PCO35            <dbl> 0.41938817, -2.61129210, -2.09536540, -0.90144359, 0.…
## $ PCO36            <dbl> 1.007770787, 2.337730545, 2.172037794, 0.760260969, -…
## $ PCO37            <dbl> 0.68321810, 0.29520404, 3.57885599, 0.38061761, -1.77…
## $ PCO38            <dbl> -0.40439268, -0.02722549, 0.40983402, 1.98480750, -0.…
## $ PCO39            <dbl> -0.94896136, 0.57403476, -0.10123012, 0.20601555, -0.…
## $ PCO40            <dbl> -1.31087655, -1.21234119, -0.18311826, -1.01916576, -…
## $ PCO41            <dbl> -1.19013371, -1.27908822, 1.61292094, 2.49380211, 0.8…
## $ PCO42            <dbl> -2.27630024, 0.39360650, -0.12971877, -0.01380194, -0…
## $ PCO43            <dbl> -2.10880483, 1.09943545, -0.60439525, 1.04520522, -0.…
## $ PCO44            <dbl> 1.19819517, 0.04760478, -0.92610622, -0.46389300, -3.…
## $ PCO45            <dbl> -0.44358713, 0.17489519, 1.31403038, -0.63687136, 1.7…
## $ PCO46            <dbl> -0.52213104, -1.34637687, -3.03213935, 3.47848838, 2.…
## $ PCO47            <dbl> 0.929193523, -0.637830194, -1.641008245, -0.628368720…
## $ PCO48            <dbl> -0.2195081, -1.1872956, 1.5649802, 0.5725595, 0.30987…
## $ PCO49            <dbl> 2.17285861, 0.17531771, 0.04014705, 0.96170327, -0.96…
## $ PCO50            <dbl> -1.2101279, 1.0437858, -0.7014482, 1.1492174, 0.23397…
## $ PCO51            <dbl> -0.02183201, -0.44816488, -0.41723680, -1.52186153, -…
## $ PCO52            <dbl> 0.37589539, -0.55249421, -1.34486386, 0.16085620, -0.…
## $ PCO53            <dbl> 0.36317140, 0.04484013, -1.73529095, 1.19508879, -0.6…
## $ PCO54            <dbl> 0.680657807, -0.359419722, 1.393585526, -1.492765585,…
## $ PCO55            <dbl> 0.39786511, 0.43774699, 0.21828214, 0.31756448, 0.894…
## $ PCO56            <dbl> -0.16257420, 0.11273957, 0.05564049, 0.16413144, -0.0…
## $ PCO57            <dbl> -0.24615059, 0.58323689, -0.32474593, 0.13805566, -0.…
## $ PCO58            <dbl> 0.53382272, 0.30987656, -0.40117793, 0.20374165, -0.0…
## $ PCO59            <dbl> 0.33093067, 0.04815707, -0.55687138, 0.97433711, -0.1…
## $ PCO60            <dbl> 0.65272370, -0.21681139, 0.28477075, -0.93189150, -0.…
## $ PCO61            <dbl> -0.122594033, 0.131776794, -0.270263415, 0.011724829,…
## $ PCO62            <dbl> 0.048876707, 0.455989497, 0.504482056, -0.319259661, …
## $ PCO63            <dbl> -0.275563372, -0.091476773, -1.059562161, 0.486282546…
## $ PCO64            <dbl> -0.10874177, -0.10252986, 0.05732430, 0.55672239, -0.…
## $ PCO65            <dbl> 0.091634733, -0.128971133, 0.298156225, 0.178421809, …
## $ PCO66            <dbl> -0.070124813, -0.003057020, 0.001705034, -0.012147311…
## $ PCO67            <dbl> -0.043879914, 0.107276489, -0.319738267, 0.251953599,…
## $ PCO68            <dbl> 0.46436640, 0.42175906, -0.06030514, -0.83895821, -0.…
## $ PCO69            <dbl> 0.542360094, 0.122058874, 0.567197694, 0.435691301, 0…
## $ PCO70            <dbl> -0.25734378, -0.59217435, 0.11926417, 0.23113661, -0.…
## $ PCO71            <dbl> 0.10455216, -0.14220872, -0.48193160, 1.47205836, -0.…
## $ PCO72            <dbl> 0.08311740, 0.74091599, 0.43239781, -0.98687017, -0.9…
## $ PCO73            <dbl> -0.34703304, 0.53767106, -0.79308552, -0.28986033, 0.…
## $ PCO74            <dbl> -0.260147153, -0.456894482, -0.007457675, -0.20545764…
## $ PCO75            <dbl> -0.38161772, -0.79414377, 0.37232686, -0.57257861, -0…
## $ PCO76            <dbl> 1.372693406, -0.005522137, -0.888201644, -0.044413108…
## $ PCO77            <dbl> 0.1780861, -0.4630424, 0.3992732, -0.1283059, -1.1427…
## $ PCO78            <dbl> 1.123559320, -0.289170212, -0.787158873, -0.309267932…
## $ PCO79            <dbl> 0.14439711, -0.95740350, 0.39057447, 0.70904360, 0.03…
## $ PCO80            <dbl> -1.06705696, 0.01990699, 1.41810230, 0.02475525, 0.61…
## $ PCO81            <dbl> -0.16758932, -0.87675058, -0.45808418, -0.78127057, -…
## $ PCO82            <dbl> -2.1588321, -0.5958911, -0.1271015, 1.7839509, 1.3585…
## $ PCO83            <dbl> -0.21820529, -2.46432548, 0.36256361, 1.57605993, 0.4…
## $ PCO84            <dbl> 1.76742994, -0.72915719, 0.41142592, 0.38923040, 2.36…
## $ PCO85            <dbl> -2.08497364, -0.15052212, -2.30609348, 0.39161981, -0…
## $ PCO86            <dbl> -2.71615266, -1.56824993, 0.93681997, 0.76935080, -0.…
## $ PCO87            <dbl> 1.97756830, -2.88665147, -0.87843317, -1.15067481, -1…
## $ PCO88            <dbl> 2.83684940, 2.30755220, -0.30444015, 0.92541879, -1.1…
## $ PCO89            <dbl> 0.60758383, -0.40037440, 1.29134214, 0.04331432, -3.9…
## $ PCO90            <dbl> 1.8751664, -1.6297199, 1.2760172, -0.8513889, 1.38680…
## $ PCO91            <dbl> 1.37735108, -0.61324110, 1.67290430, -0.96739365, 1.5…
## $ PCO92            <dbl> -1.90054564, 2.69654825, 0.40242463, 1.34181660, -2.2…
## $ PCO93            <dbl> 0.581197290, -0.666370647, -2.872479867, 1.750333951,…
## $ PCO94            <dbl> -1.618284495, 1.881134414, 1.532679818, 3.750732560, …
## $ PCO95            <dbl> 2.3666030, -3.8313315, 0.4800660, 3.0730310, -1.88522…
## $ PCO96            <dbl> -2.798925686, -2.457247581, 1.119049270, 2.946038020,…
## $ PCO97            <dbl> 1.9004347, -3.6946658, -2.0603858, 1.3464237, -1.0177…
## $ PCO98            <dbl> 5.2305380, -2.5855596, 1.1161859, 3.3560485, -0.54454…
## $ X                <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ Site.1           <chr> "GK1", "GK2", "GK3", "GK4", "GK7", "GK8", "H3", "MI1"…
## $ Year             <chr> "Undefined!", "Undefined!", "Undefined!", "Undefined!…
## $ Region           <chr> "Keppel", "Keppel", "Keppel", "Keppel", "Keppel", "Ke…
## $ Viz              <chr> "Undefined!", "Undefined!", "Undefined!", "Undefined!…
## $ Depth            <chr> "Undefined!", "Undefined!", "Undefined!", "Undefined!…
## $ Exposure         <chr> "Semi-Exposed", "Exposed", "Semi-Exposed", "Exposed",…
## $ NTR              <chr> "Fished", "Fished", "Fished", "Fished", "Fished", "Fi…
## $ NTR.Pooled       <chr> "Fished", "Fished", "Fished", "Fished", "Fished", "Fi…
## $ Transect         <chr> "Undefined!", "Undefined!", "Undefined!", "Undefined!…
## $ YearNTR.Pooled   <chr> "Undefined!", "Undefined!", "Undefined!", "Undefined!…
## $ RegionNTR.Pooled <chr> "KeppelFished", "KeppelFished", "KeppelFished", "Kepp…
fish = fish %>%
  left_join(pco %>% dplyr::select(SITE=Site, PCO1, PCO2) %>% distinct)
#save(fish, file='data/fish.RData')

## the following are the PCO scores for calculated purely WITHIN a region
## these will be added as PCO1r and PCO2r
## the 'r' stands for region and for region specific analyses,  the formulae will be ammended to add the 'r'
pco.k <- read.csv('data/Fish 2007-2019 PCO scores_Keppel.csv',  strip.white=TRUE) %>% dplyr::select(SITE=Site, PCO1r=PC1, PCO2r=PC2) 
pco.m <- read.csv('data/Fish 2007-2019 PCO scores_Magnetic.csv',  strip.white=TRUE) %>% dplyr::select(SITE=Site, PCO1r=PC1, PCO2r=PC2)
pco.p <- read.csv('data/Fish 2007-2019 PCO scores_Palm.csv',  strip.white=TRUE) %>% dplyr::select(SITE=Site, PCO1r=PC1, PCO2r=PC2)
pco.w <- read.csv('data/Fish 2007-2019 PCO scores_Whitsunday.csv',  strip.white=TRUE) %>% dplyr::select(SITE=Site, PCO1r=PC1, PCO2r=PC2)
pco.r <- bind_rows(pco.k, pco.m, pco.p, pco.w)
pco.r %>% glimpse
## Rows: 99
## Columns: 3
## $ SITE  <chr> "GK1", "GK2", "GK3", "GK4", "GK7", "GK8", "H3", "MI1", "MI2", "N…
## $ PCO1r <dbl> -15.6612542, -11.8407823, 3.6681599, -2.6515472, -7.7527675, 4.3…
## $ PCO2r <dbl> 2.0982480, -9.4099666, 5.7880978, 0.6778477, -3.1508521, -12.666…
fish=fish %>%
  left_join(pco.r)
var.lookup = rbind(
    data.frame(pretty.name='Total density', Field.name='Total.fish.density', Abbreviation='TFD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Total species richness', Field.name='Total.fish.species.richness', Abbreviation='TFSR', Family='gaussian', Type='Response', Transform='I', Groupby=''),
    data.frame(pretty.name='Benthic invertivores', Field.name='BE', Abbreviation='BE', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Grazers', Field.name='GRAZ', Abbreviation='GRAZ', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Grazers2', Field.name='GRAZ2', Abbreviation='GRAZ2', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Parrot', Field.name='Parrot', Abbreviation='PA', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Corallivores', Field.name='COR', Abbreviation='COR', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Omnivores', Field.name='OM', Abbreviation='OM', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Planktivores', Field.name='PL', Abbreviation='PL', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Carnivores', Field.name='CA', Abbreviation='CA', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Piscivores', Field.name='PI', Abbreviation='PI', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Farmers', Field.name='FA', Abbreviation='FA', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus total density', Field.name='Plectropomus.total.density', Abbreviation='PTD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus total biomass', Field.name='Plectropomus.total.biomass', Abbreviation='PTB', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus legal density', Field.name='Plectropomus.legal.density', Abbreviation='PLD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus legal biomass', Field.name='Plectropomus.legal.biomass', Abbreviation='PLB', Family='gaussian', Type='Response', Transform='log', Groupby=''),
  data.frame(pretty.name='PCO1',  Field.name='PCO1', Abbreviation='PCO1', Family='gaussian', Type='Response', Transform='I', Groupby=''), 

  data.frame(pretty.name='Region', Field.name='REGION', Abbreviation='REGION', Family=NA, Type='Predictor', Transform='I', Groupby=''),
    data.frame(pretty.name='Zone', Field.name='NTR.Pooled', Abbreviation='NTR.Pooled', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    ## data.frame(pretty.name='LHC % (live hard coral cover)', Field.name='LHC_.', Abbreviation='LHC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='% hard coral', Field.name='LHC_.', Abbreviation='LHC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    ## data.frame(pretty.name='SC % (soft coral)', Field.name='SC_.', Abbreviation='SC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='% soft coral', Field.name='SC_.', Abbreviation='SC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    ## data.frame(pretty.name='MA % (macroalgal cover)', Field.name='MAC_.', Abbreviation='MA', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='% macroalgae', Field.name='MAC_.', Abbreviation='MA', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    ## data.frame(pretty.name='Turf %', Field.name='Turf_.', Abbreviation='TURF', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='% turf', Field.name='Turf_.', Abbreviation='TURF', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    ## data.frame(pretty.name='Unconsolidated %', Field.name='Unconsolidated_.', Abbreviation='UC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='% unconsolidated', Field.name='Unconsolidated_.', Abbreviation='UC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Benthic richness', Field.name='Benthic.richness', Abbreviation='BR', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Coral morphological diversity', Field.name='Coral_Morph.Diversity', Abbreviation='CMD', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Slope', Field.name='slope', Abbreviation='SLOPE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Rugosity', Field.name='rugosity', Abbreviation='RUG', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='SCI (structural complexity)', Field.name='SCI', Abbreviation='SCI', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Chla', Field.name='ChlA', Abbreviation='CHL', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Kd490', Field.name='kd490', Abbreviation='KD490', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='SST mean', Field.name='SSTmean', Abbreviation='SSTMEAN', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='SST anom', Field.name='SSTanom', Abbreviation='SSTANOM', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Wave exposure', Field.name='wave.exposure.index', Abbreviation='WAVE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Corrected depth', Field.name='Corrected.depth', Abbreviation='DEPTH', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Max DHW', Field.name='maxDHW', Abbreviation='DHW', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Cyclone', Field.name='Cyclone', Abbreviation='CYCLONE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='High turbidity exposure', Field.name='Exposure.to.primary.weeks', Abbreviation='EXP', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Prey density', Field.name='Prey.density', Abbreviation='PREY.DENSITY', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
  data.frame(pretty.name='Prey biomass', Field.name='Prey.biomass', Abbreviation='PREY.BIOMASS', Family=NA, Type='Predictor', Transform='I', Groupby='Region')
)
save(var.lookup, file='data/var.lookup.RData') 
names = with(var.lookup, setNames(as.character(Field.name), Abbreviation))
## exclude those whose names equal their values otherwise there will be duplicate fields created in the fish data
names = names[names(names)!=names]
## Create duplicates of the fields that are not already abbreviated 
fish = fish %>%
    mutate(SSTmean=ifelse(as.character(SSTmean)=='#N/A',NA,as.numeric(as.character(SSTmean)))) %>%
    mutate_at(as.character(var.lookup$Field.name), list(A=~I(.))) %>%
    rename(!!! gsub('(.*)','\\1_A',names)) %>%
    mutate(REGION=factor(REGION, levels=c('Palm','Magnetic','Whitsunday','Keppel')))
save(fish, file='data/fish.RData')

Response Variables:

Variable Data field Abbreviation Distribution
Total density Total.fish.density TFD Gaussian (log)
Total species richness Total.fish.species.richness TFSR Gaussian
Benthic invertivores BE BE Gaussian (log)
Grazers GRAZ GRAZ Gaussian (log)
Corallivores COR COR Gaussian (log)
Omnivores OM OM Gaussian (log)
Planktivores PL PL Gaussian (log)
Carnivores CA CA Gaussian (log)
Piscivores PI PI Gaussian (log)
Farmers FA FA Gaussian (log)
Plectropomus total density Plectropomus.total.density PTD Gaussian (log)
Plectropomus total biomass Plectropomus.total.biomass PTB Gaussian (log)
Plectropomus legal density Plectropomus.legal.density PLD Gaussian (log)
Plectropomus legal biomass Plectropomus.legal.biomass PLB Gaussian (log)
PCO1 PCO1 PCO1 Identity

Predictor Variables:

Variable Data field Abbreviation
Region REGION REGION
NTR Pooled NTR.Pooled NTR.Pooled
% hard cover LHC_. LHC
% soft coral SC_. SC
% macroalgae MAC_. MA
% turf Turf_. TURF
% unconsolidated Unconsolidated_. UC
Benthic richness Benthic.richness BR
Coral morphological diversity Coral_Morph.Diversity CMD
Slope slope SLOPE
Rugosity rugosity RUG
SCI (structural complexity) SCI SCI
Chla ChlA CHL
Kd490 kd490 KD490
SST mean SSTmean SSTMEAN
SST anom SSTanom SSTANOM
Wave exposure wave.exposure.index WAVE
Corrected depth Corrected.depth DEPTH
Max DHW maxDHW DHW
Cyclone Cyclone CYCLONE
Exposure to primary weeks Exposure.to.primary.weeks EXP

FOR Plectropomus response variables ONLY, add:

Variable Data field Abbreviation
Prey density (for Plectropomus density and legal density) Prey.density PREY.DENSITY
Prey biomass (for Plectropomus biomass and legal biomass) Prey.biomass PREY.BIOMASS

Species-specific analyses: list of species in the attached file. There are two lists there: the single-column list is for the initial analysis between regions, and the table is for the species lists to use for each region individually.

Raw selections data 
selections = read.csv('data/Species selection.csv', strip.white=TRUE)
selections %>% glimpse
## Rows: 248
## Columns: 6
## $ Palms.model       <chr> "sca.rivu", "sig.doli", "sca.flav", "chs.sord", "sca…
## $ Magnetic.model    <chr> "sca.rivu", "sca.ghob", "aca.duss", "aca.bloc", "sig…
## $ Whitsunday.model  <chr> "sca.rivu", "sig.doli", "sca.flav", "sca.nigr", "chs…
## $ Keppel.model      <chr> "sca.rivu", "sca.ghob", "sig.arge", "sig.fusc", "sig…
## $ X                 <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ All.regions.model <chr> "aca.line", "aca.ncus", "nas.brach", "nas.litu", "ca…

Results - full analyses

When running the full analyses, the substantially influential predictors for any single analysis are identified. These predictors are then expressed as their spatial and temporal components by centering them against their respective temporal (for each location) and spatial (for each year) means. The analyses are then repeated using just the important (spatial and temporal) versions of the influential predictors. These later analyses are used to identify which predictors should feature in spatially focussed (and later temporally focussed) analyses. Finally, the spatial and temporal versions of the influential predictors are combined with all other predictors and the analyses are re-run as an alternative way to identify predictors should feature in spatially focussed (and later temporally focussed) analyses.

formulas = list(
    all = ~REGION + NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    all1 = ~NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Palm = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Magnetic = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Whitsunday = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Keppel = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP
)
load('data/fish.RData')
load('data/var.lookup.RData')
resp.lookup = var.lookup %>% filter(Type=='Response') %>% droplevels
pred.lookup = var.lookup %>% filter(Type=='Predictor') %>% droplevels

groupings.all = vector('list', nrow(pred.lookup))
groupings.all1 = vector('list', nrow(pred.lookup))
names(groupings.all) = pred.lookup$Abbreviation
groupings.region = vector('list', nrow(pred.lookup))
names(groupings.region) = pred.lookup$Abbreviation
for (i in 1:nrow(pred.lookup)) {
    pred = as.character(pred.lookup[i,'Abbreviation'])
    groupings.all[[pred]] = ifelse(pred=='REGION',NA,'REGION')
    groupings.all1[[pred]] = NA
    groupings.region[[pred]] = ifelse(pred=='NTR.Pooled',NA,'NTR.Pooled')
}
groupings.all = do.call('c',groupings.all)
groupings.all1 = do.call('c',groupings.all1)
groupings.region = do.call('c',groupings.region)

gfun = function(f, form) {
    #print(f)
    #print(form[[f]])
    form = attr(terms(form[[f]]),'term.labels')
    if (f=='all') return(groupings.all[form])
    if (f=='all1') return(groupings.all1[form])
    else return(groupings.region[form])
}

analyses = vector('list',nrow(resp.lookup))
names(analyses) = resp.lookup$Abbreviation
for (i in 1:nrow(resp.lookup)) {
    resp = as.character(resp.lookup[i,'Abbreviation'])
    fun = as.character(resp.lookup[i,'Transform'])
    fam = as.character(resp.lookup[i,'Family'])
    
    analyses[[resp]] = list('formulas' = lapply(formulas, function(f) update(f, paste0(fun,'(',resp,') ~.'))),
                            'family' = fam)
    if (resp %in% c('PI','PTD','PLD')) analyses[[resp]][['formulas']] = lapply(analyses[[resp]][['formulas']], function(f) f=update(f, .~.+PREY.DENSITY))
    if (resp %in% c('PTB','PLB')) analyses[[resp]][['formulas']] = lapply(analyses[[resp]][['formulas']], function(f) f=update(f, .~.+PREY.BIOMASS))
    analyses[[resp]][['groups']] = sapply(c('all','all1','Palm','Magnetic','Whitsunday','Keppel'), gfun, form=analyses[[resp]][['formulas']], USE.NAMES = TRUE,simplify=FALSE)
    if (resp %in% c('PCO1', 'PCO2')) {
      for (ff in 3:6) {  # used to be 2:5 (as in the regions only)
        analyses[[resp]]$formulas[[ff]] = update(analyses[[resp]]$formulas[[ff]],  paste0(fun, '(', resp,'r', ') ~.'))
      }
    }
    analyses[[resp]][['groups']] = sapply(c('all','all1','Palm','Magnetic','Whitsunday','Keppel'), gfun, form=analyses[[resp]]$formulas, USE.NAMES = TRUE,simplify=FALSE)
}
save(analyses, file='data/analyses.RData')
for (a in 1:length(analyses)) {
    resp=names(analyses)[a]
    print(paste('Response =',resp))
    ## for (f in 1:length(analyses[[a]]$formulas)) {
    for (f in 2) {  # temporary just so that we can run the new all1 set
        mod.name = names(analyses[[a]]$formulas)[f]
        print(paste('Model =',mod.name))
        MONOTONE = assignMonotone(fish, analyses[[a]]$formulas[[f]])
        if (mod.name %in% c('all','all1')) {fish.sub=fish
        } else {fish.sub = fish %>% filter(REGION==mod.name)}
        if (any(fish.sub[,as.character(get_response(analyses[[a]]$formulas[[f]]))]==0)) {
            val = fish.sub[, as.character(get_response(analyses[[a]]$formulas[[f]]))]
            val=min(val[val>0], na.rm=TRUE)
            fish.sub[, as.character(get_response(analyses[[a]]$formulas[[f]]))] = fish.sub[,as.character(get_response(analyses[[a]]$formulas[[f]]))] + val
        }
        set.seed(123)
        fish.sub = fish.sub %>% mutate_if(is.character,  as.factor)
        mod = abt(analyses[[a]]$formulas[[f]], data=fish.sub, distribution=analyses[[a]]$family,
                  cv.folds=10,interaction.depth=10,n.trees=10000, shrinkage=0.001, n.minobsinnode=2,
                  var.monotone=as.vector(MONOTONE))
        if (f>2) {  # only applies to the regional models
          var.lookup1 = var.lookup %>%
            mutate(Field.name=ifelse(Field.name %in% c('PCO1', 'PCO2'), paste0(Field.name, 'r'), Field.name),
                   Abbreviation=ifelse(Abbreviation %in% c('PCO1', 'PCO2'), paste0(Abbreviation, 'r'), Abbreviation))
        } else {
          var.lookup1 = var.lookup
        }
        gr <- na.omit(unique(analyses[[a]]$groups[[f]]))
        if (is.logical(gr)) gr=NULL
        p=plot.abts(mod, var.lookup1, center=FALSE, type='response', return.grid=TRUE,
                    groupby=gr,pt.size=14)#
        thresholds=p$thresholds
        save(thresholds, file=paste0('data/thresholds_',mod.name,'_',resp,'.RData'))
        ps = p[['ps']]
        p = p[['p']]
        if ((length(levels(fish$REGION))>1 || length(levels(fish$NTR.Pooled))>1) & mod.name!='all1') {
            p = common_legend(p)
            ps = common_legend(ps)
        }
        ## version for the supplimentary
                                        #do.call('grid.arrange', p) ## version for the supplimentary
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT.png'), do.call('grid.arrange', p), width=15, height=10, dpi=300)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT.pdf'), do.call('grid.arrange', p), width=15, height=10)

        pred.1=stats.abt(mod, fitMethod=1, analysis = analyses[[a]]$groups[[f]])
        
        save(pred.1, file=paste0('data/pred.1_',mod.name,'_',resp,'.RData'))
        
        optim=summarize_values(pred.1$optim, type='optim') %>%
            dplyr::rename_at(vars(-contains('Var')), paste0, '.optim')
        
        rel.inf = summarize_values(pred.1$rel.imp, type='Rel.inf') %>%
            dplyr::rename_at(vars(-contains('Var')), paste0, '.rel.inf')
        
        R2=summarize_values(pred.1$R2.value, 'R2') %>%
            dplyr::rename_at(vars(-contains('Var')), paste0, '.R2')
        
        stats = optim %>% full_join(R2) %>% full_join(rel.inf) %>%
            left_join(var.lookup %>% dplyr::select(Var=Abbreviation, pretty.name))
        
        save(stats, file=paste0('data/stats_',mod.name,'_',resp,'.RData'))
        write.csv(stats %>% as.data.frame, file=paste0('output/data/stats_',mod.name,'_',resp,'.csv'), quote=FALSE, row.names=FALSE)

        ## Version full model with just the relative importance and the substantial panels
        ps1 = arrangeGrob(ps[[1]],ps[[length(ps)]], widths=c(3,1))
        numberOfPlots = length(ps)-2 #minus 1 since one of the items is the legend and minus one for the relative importance plot
        numberOfPlotRows = ceiling(numberOfPlots / 2)
        if (length(ps)<3) {
          g=ps1 ##arrangeGrob(ps1, ps[2])    #do.call('arrangeGrob', c(ps[c(-1, -length(ps))], list(ncol=2))), heights=c(2,numberOfPlotRows))
        } else {
          g=arrangeGrob(ps1, do.call('arrangeGrob', c(ps[c(-1, -length(ps))], list(ncol=2))), heights=c(2,numberOfPlotRows))
        }
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short.png'), g, width=10, height=2+(1.7*numberOfPlotRows), dpi=300)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short.pdf'), g, width=10, height=2+(1.7*numberOfPlotRows))

        ## Version with a grid of partials in the lower right corner of the influence figure
        ## The location and size of the figure in figure will be determined by the scale of relative influence
        ## ideally, we want the subfigure to be 3/4 of the width and 3/4 of the height
        ymax=max(pretty(as.numeric(as.character(stats$upper.rel.inf))))
        ymin=ymax*1/4
        xmax=(length(p)-2)*3/4
        if (ymin <= 100/(length(p)-2)) {  # if the left side of the subfigure is too close to the dashed vertical line
            ymax = ymax + 5
            ymin=ymax*1/4
            ps[[1]] = ps[[1]] + scale_y_continuous('Relative Importance', limits=c(0,ymax))
        }
        ps[2:(length(ps)-1)]=lapply(ps[2:(length(ps)-1)], function(f) f+theme(axis.title.y=element_blank())) 
        ps2=do.call('arrangeGrob', c(ps[c(length(ps),2:(length(ps)-1))], list(ncol=2)))
                                        #g= ps[[1]] + annotation_custom(grob=ps2, xmin=1, xmax=15, ymin=10, ymax=Inf)
        g= ps[[1]] + annotation_custom(grob=ps2, xmin=1, xmax=xmax, ymin=ymin, ymax=Inf)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short_in.png'), g, width=10, height=8, dpi=300)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short_in.pdf'), g, width=10, height=8)
        
        ## Version for vertical table of regional relative importance and substantial panels
        if (length(ps)<2) {
          g = ps[[1]]
        } else {
          g=do.call('arrangeGrob', c(ps[1:4], list(ncol=1, heights=c(1.5, rep(1,3)))))
        }
                                        #save(g, file=paste0('data/g.abt_',mod_num,'_',resp,'.RData'))
        
                                        #save(data.all.abt, file=paste0('data/data.all.abt_',mod_num,'_',resp,'.RData'))
        
        
                                        #save(rel.importance, file=paste0('data/rel.importance_',mod_num,'_',resp,'.RData'))
                                        #save(thresholds, file=paste0('data/thresholds_', mod_num,'_',resp,'.RData'))
                                        #write.table(bind_rows(thresholds, .id='Var'), file=paste0('data/thresholds_',mod_num,'_',resp,'.csv'), quote=FALSE, row.names=FALSE)
        if (1==2) { # all this is to attempt to address co-authors comments
            ## Many of the predictors vary over space and time.  As is, the models just indicate whether there is a relationship between
            ## the response and the various predictors.  Although fitting the models separately for the different Regions does partly tease
            ## apart the spatial from temporal, it does not completely do this (as there is still spatiality within Regions.
            ## One way we might be able to address this is to take the 'important predictors' from a given model and then refit with
            ## specifically centred versions of the predictors.  For example, if we centre a predictor within a site, then we effectively remove
            ## the spatial component of this variable.  Similarly, if we centre on time (year), then we remove the temporal element.

            ## 1. start by identifying the important predictors
            wch=rel.inf$Var[which(as.numeric(rel.inf$Mean.rel.inf)> (100/nrow(rel.inf)))]
            wch.n = wch[unlist(lapply(wch, function(x) is.numeric(fish.sub[,x])))]
            wch.c = wch[unlist(lapply(wch, function(x) is.factor(fish.sub[,x])))]

            ## Remove spatial element
            fish.sub1 = fish.sub %>% group_by_at(vars(one_of(c(wch.c,'SITE')))) %>%
                mutate_at(vars(wch.n), function(x) x-mean(x)) %>%
                dplyr::rename_at(vars(wch.n), ~paste0(wch.n,'.t')) %>%
                dplyr::select(!!c('TFD','YEAR',paste0(wch.n,'.t')))
            fish.sub1.means = fish.sub %>% group_by_at(vars(one_of(c(wch.c,'SITE')))) %>%
                summarize_at(vars(wch.n), function(x) mean(x)) 

            ## #fish.sub1 %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC, color=SITE)) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC, color=factor(YEAR))) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_line(aes(y=log(TFD),x=LHC, color=SITE)) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_line(aes(y=log(TFD),x=LHC, color=factor(YEAR))) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_line(aes(y=log(TFD),x=YEAR, color=factor(SITE))) + facet_wrap(~REGION)
            ## fish.sub1 %>% ggplot() + geom_line(aes(y=LHC.t,x=YEAR, color=factor(SITE))) + facet_wrap(~REGION)
            ## #fish.sub %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC)) + facet_wrap(~REGION)
            ## fish.sub %>% ggplot() + geom_line(aes(y=LHC,x=YEAR, color=factor(SITE))) + facet_wrap(~REGION)
            
            ## Remove temporal element
            fish.sub2 = fish.sub %>% group_by_at(vars(one_of(c(wch.c,'YEAR')))) %>%
                mutate_at(vars(wch.n), function(x) x-mean(x)) %>%
                dplyr::rename_at(vars(wch.n), ~paste0(wch.n,'.s')) %>%
                dplyr::select(!!c('TFD','SITE',paste0(wch.n,'.s')))

            fish.sub3 = fish.sub1 %>% full_join(fish.sub2)
            
            fish.sub2.means = fish.sub %>% group_by_at(vars(one_of(c(wch.c,'YEAR')))) %>%
                summarize_at(vars(wch.n), function(x) mean(x)) 

            ## fish.sub1 %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC.t, color=factor(YEAR))) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_line(aes(y=log(TFD),x=LHC, color=factor(SITE))) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_line(aes(y=log(TFD),x=LHC, color=factor(YEAR))) + facet_wrap(~REGION)
            ## fish.sub2 %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC.s, color=factor(YEAR))) + facet_wrap(~REGION)
            ## fish.sub %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC)) + facet_wrap(~REGION)

            ## fish.sub1 %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC.s, color=factor(SITE))) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_line(aes(y=log(TFD),x=LHC, color=factor(SITE))) + facet_wrap(~REGION)
            ## #fish.sub1 %>% ggplot() + geom_line(aes(y=log(TFD),x=LHC, color=factor(YEAR))) + facet_wrap(~REGION)
            ## fish.sub2 %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC, color=factor(SITE))) + facet_wrap(~REGION)
            ## fish.sub %>% ggplot() + geom_point(aes(y=log(TFD),x=LHC)) + facet_wrap(~REGION)

            
            ## fish.sub2 %>% ggplot() + geom_line(aes(y=LHC,x=YEAR, color=factor(SITE))) + facet_wrap(~REGION)
            ## fish.sub2 %>% ggplot() + geom_line(aes(y=LHC,x=as.numeric(SITE), color=factor(YEAR))) + facet_wrap(~REGION)

            p.names=colnames(fish.sub3)[!colnames(fish.sub3) %in% c('SITE','TFD','YEAR')]
            ff = update(analyses[[a]]$formulas[[f]], . ~ 1)
            ff=as.formula(c(gsub('1','',deparse(ff)), paste(p.names, collapse='+')))
            MONOTONE = assignMonotone(fish.sub3, ff)
            set.seed(123)
            mod = abt(ff, data=fish.sub3, distribution=analyses[[a]]$family,
                      cv.folds=10,interaction.depth=10,n.trees=10000, shrinkage=0.001, n.minobsinnode=2,
                      var.monotone=as.vector(MONOTONE))
            summary(mod)
            ## pdf(file='TFD.W.test.pdf', width=10, height=30)
            ## par(mfrow=c(5,2))
            ## #for (i in c(2,7,3,8,4,9,5,10,6,11,1)) plot(mod,i, ylim=c(-1.5,1.5))
            ## for (i in c(1,6,2,7,3,8,4,9,5,10)) plot(mod,i, ylim=c(-1,1))
            ## dev.off()

            ## And now including the other predictors as well
            fish.sub1 = fish.sub %>% group_by_at(vars(one_of(c(wch.c,'SITE')))) %>%
                mutate_at(vars(wch.n), function(x) x-mean(x)) %>%
                dplyr::rename_at(vars(wch.n), ~paste0(wch.n,'.t')) %>%
                dplyr::select(!!c('TFD','YEAR',paste0(wch.n,'.t')))
            fish.sub2 = fish.sub %>% group_by_at(vars(one_of(c(wch.c,'YEAR')))) %>%
                mutate_at(vars(wch.n), function(x) x-mean(x)) %>%
                dplyr::rename_at(vars(wch.n), ~paste0(wch.n,'.s')) %>%
                dplyr::select(!!c('TFD','SITE',paste0(wch.n,'.s')))

            fish.sub3 = fish.sub1 %>% full_join(fish.sub2) %>% full_join(fish.sub)
            
            ff = deparse(analyses[[a]]$formulas[[f]])
            for (w in wch.n) ff = gsub(w,paste0(w,'.t + ',w,'.s'), ff)
            MONOTONE = assignMonotone(fish.sub3, ff)
            set.seed(123)
            mod = abt(ff, data=fish.sub3, distribution=analyses[[a]]$family,
                      cv.folds=10,interaction.depth=10,n.trees=10000, shrinkage=0.001, n.minobsinnode=2,
                      var.monotone=as.vector(MONOTONE))
            summary(mod)
            ## Capture the list of influential spatial and temporal predictors...

            
            ## pdf(file='TFD.W.test.pdf', width=10, height=30)
            ## pdf(file='TFD.Region.test.pdf', width=10, height=30)
            ## par(mfrow=c(5,2))
            ##                             #for (i in c(2,7,3,8,4,9,5,10,6,11,1)) plot(mod,i, ylim=c(-1.5,1.5))
            ##                             #for (i in c(5,6,7,8,2,3,19,20,11,22)) plot(mod,i, ylim=c(-1,1))
            ## for (i in c(3,4,8,9,20,21,16,1,6,7)) plot(mod,i, ylim=c(-1,1))
            ## dev.off()

        }
        
        rm(mod,pred.1,p,ps,ps1,g,stats,optim,rel.inf,R2,fit)
        gc()
    }
}
#summary(mod)
#plot(mod,c(15,1))
load('data/var.lookup.RData')
resp.lookup = var.lookup %>% filter(Type=='Response') %>% droplevels
mods = names(formulas)
for (a in 1:length(analyses)) {
    resp=names(analyses)[a]
    cat(paste('## ',resp.lookup$pretty.name[resp.lookup$Abbreviation==resp],' {.tabset .tabset-pills} \n\n'))
    for (f in 1:length(analyses[[a]]$formulas)) {
        mod.name = names(analyses[[a]]$formulas)[f]
        cat(paste('### ',mod.name,'\n\n'))
        cat(paste0('![](output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short_in.png)\n\n'))
        ## Now for a tables
        load(file=paste0('data/stats_',mod.name,'_',resp,'.RData'))
        pretty.stats(stats, var.lookup) %>% rmarkdown:::paged_table_html() %>% cat
        load(file=paste0('data/thresholds_',mod.name,'_',resp,'.RData'))
        pretty.thresholds(thresholds, stats, var.lookup) %>% rmarkdown:::paged_table_html() %>% cat
    }
}

Results - spatial focus

load('data/var.lookup_spatial.RData')
#load(file='data/new_analyses.RData')
load(file='data/analyses.RData')
new_analyses=analyses
resp.lookup = var.lookup %>% filter(Type=='Response') %>% droplevels
mods = names(formulas)
for (a in 1:length(new_analyses)) {
    resp=names(new_analyses)[a]
    cat(paste('## ',resp.lookup$pretty.name[resp.lookup$Abbreviation==resp],' {.tabset .tabset-pills} \n\n'))
    for (f in 1:length(new_analyses[[a]]$formulas)) {
        mod.name = names(new_analyses[[a]]$formulas)[f]
        cat(paste('### ',mod.name,'\n\n'))
        cat(paste0('![](output/figures/data.all.abt_spatial.',mod.name,'_',resp,'_ABT_short_in.png)\n\n'))
        ## Now for a tables
        load(file=paste0('data/stats_spatial_',mod.name,'_',resp,'.RData'))
        pretty.stats(stats, var.lookup) %>% rmarkdown:::paged_table_html() %>% cat
        load(file=paste0('data/thresholds_spatial_',mod.name,'_',resp,'.RData'))
        pretty.thresholds(thresholds, stats, var.lookup) %>% rmarkdown:::paged_table_html() %>% cat 
    }
}

Total density

all

all1

Palm

Magnetic

Whitsunday

Keppel

Total species richness

all

all1

Palm

Magnetic

Whitsunday

Keppel

Benthic invertivores

all

all1

Palm

Magnetic

Whitsunday

Keppel

Grazers

all

all1

Palm

Magnetic

Whitsunday

Keppel

Grazers2

all

all1

Palm

Magnetic

Whitsunday

Keppel

Parrot

all

all1

Palm

Magnetic

Whitsunday

Keppel

Corallivores

all

all1

Palm

Magnetic

Whitsunday

Keppel

Omnivores

all

all1

Palm

Magnetic

Whitsunday

Keppel

Planktivores

all

all1

Palm

Magnetic

Whitsunday

Keppel

Carnivores

all

all1

Palm

Magnetic

Whitsunday

Keppel

Piscivores

all

all1

Palm

Magnetic

Whitsunday

Keppel

Farmers

all

all1

Palm

Magnetic

Whitsunday

Keppel

Plectropomus total density

all

all1

Palm

Magnetic

Whitsunday

Keppel

Plectropomus total biomass

all

all1

Palm

Magnetic

Whitsunday

Keppel

PCO1

all

all1

Palm

Magnetic

Whitsunday

Keppel

Heatmaps

all

all1

Palm

Magnetic

Whitsunday

Keppel