Preparations

Load the necessary libraries

#library(tidyverse) #for data wrangling, includes dplyr and ggplot2
library(dplyr)
library(ggplot2)
library(pspearman)  #for spearman correlation
library(tidybayes)    #for data wrangling
library(performance) #for VIF estimation
library(GGally)    #for plotting graphs
library(ggpubr)    #for combining multiple plots 
library(DHARMa)   #for residuals and diagnostics
library(rstan)      #for interfacing with STAN
library(brms)       #for fitting models in STAN
library(ggmcmc)     #for MCMC diagnostics
library(mice)        #for imputing missing data
library(ppcor)      #for estimating partial correlations
library(ggnewscale)   #for manual plotting of Bayesian results
library(MCMCvis)     #for automatic plotting of Bayesian results
library(tidybayes)   #for tidymcmc Bayesian results
library(ggeffects)    #for model testing
library(ggridges)   #for plotting
library(stringr)     #for pars variables

Introduction

The aim of the analysis is to investigate whether a comprehensive operationalization of actor-specific adaptive capacity could usefully explain variations in observed responses among tourism operators, while controlling for contextual conditions related to the disturbance characteristics and governance settings where the tourism operators were based.

Description of method

As described in greater detail in our manuscript, we conducted surveys with representatives of reef tourism companies (operators) in the APAC Region. The surveys requested information on the actions that each operator took in response to a specific climate disturbance and indicators of the operator’s adaptive capacity before the disturbance. Our study design included a priori treatment and control groups of tourism operators based on whether their locations had been directly affected by the climate disturbance. We used a combination of descriptive and multivariate statistical approaches to explore the relationships between adaptive capacity, adaptive responses, and contextual variables that were relevant to the relationship between adaptive capacity and adaptive responses.

Read in the data

data = read.csv('HB_adaptive_capacity_raw.csv')
glimpse(data)
## Rows: 231
## Columns: 38
## $ operator           <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, …
## $ location           <chr> "id", "id", "id", "id", "mi", "id", "id", "fp", "id…
## $ gov_effect         <dbl> -0.04, -0.04, -0.04, -0.04, 0.44, -0.04, -0.04, 1.4…
## $ resp_region        <chr> "eur", "eur", "other", "nam", "nam", "nam", "nam", …
## $ resp_age           <chr> "45 - 54", "35 - 44", "35 - 44", "55 - 64", "55 - 6…
## $ resp_gender        <int> 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, …
## $ dist_type          <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ dist_severity      <dbl> 0.00, 0.75, 0.50, 0.00, 0.50, 0.25, 0.25, 0.00, 0.0…
## $ divbl              <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, …
## $ chgsites           <int> 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, …
## $ chgact             <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, …
## $ chgopmode          <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, …
## $ nrm                <int> 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, …
## $ insurance          <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ monitor            <int> 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, …
## $ relief             <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ support            <int> 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ co2                <int> 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, …
## $ educate            <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ none               <int> 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, …
## $ topresponse        <chr> "monitor", "monitor", "none", "none", "nrm", "chgsi…
## $ exclaccess         <int> 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ particmgmt_tourism <int> 1, 2, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 2, 2, 2, 2, …
## $ particmgmt_coral   <int> 1, 2, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 2, 1, 2, 2, …
## $ psgseats           <int> 1, 1, 1, 5, 1, 2, 2, 1, 2, 2, 1, 3, 1, 1, 2, 1, 2, …
## $ savings            <int> 2, 3, 2, 4, 1, 2, 3, 1, 3, 4, 1, 3, 2, 2, 3, 2, 4, …
## $ distmarket         <int> 90, 90, 60, 90, 17, 135, 135, 5, 50, 18, 7, 46, 4, …
## $ othservices        <int> 2, 3, 0, 0, 0, 2, 3, 0, 2, 2, 1, 7, 0, 0, 0, 2, 0, …
## $ accsites           <int> 4, 5, 3, 5, 2, 5, 5, 2, 3, 5, 5, 5, 1, 4, 4, 5, 2, …
## $ edu                <int> 3, 4, 3, 4, 1, 3, 3, 3, 3, 3, 3, 4, 4, 3, 4, 3, 1, …
## $ experience         <int> 3, 3, 1, 4, 4, 3, 3, 3, 2, 3, 4, 3, 2, 4, 3, 4, 3, …
## $ employ_involve     <int> 2, 3, 3, 2, 3, 3, 3, 0, 2, 2, 1, 2, 2, 0, 2, 2, 2, …
## $ ind_memberships    <int> 1, 1, 0, 1, 0, 0, 1, 2, 0, 1, 0, 0, 1, 2, 0, 2, 2, …
## $ research_rely      <int> 1, 2, 3, 0, 2, 0, 0, 0, 2, 1, 1, 2, 1, 0, 0, 0, 2, …
## $ govt_ties          <int> 2, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 2, …
## $ reef_attach        <int> 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 0, 2, 2, 1, 2, …
## $ comp_concerns      <int> 1, 0, 0, 3, 0, 1, 1, 0, 2, 2, 1, 2, 0, 0, 1, 1, 0, …
## $ adapt_conf         <int> 2, 2, 2, 3, 0, 3, 3, 2, 2, 2, 1, 2, 0, 2, 0, 1, 2, …

Data preparation

Predictors correlation (potential for collinearity in models)

We first evaluated the paired correlation between our predictors to see whether there are any highly correlated variables. The highest correlation was found, as we expected between our two indicators of management participation (in reef tourism mgmt., and in coral reef mgmt.). These predictors had a correlation of 74% using Spearman’s rank test. All other correlations were below 40%. We merged our two management participation indicators into a combined predictor that took the maximum of any of the two values.

For participation in management (x-scale): 0 for none, 1 for passive attendance (attended meetings but did not usually speak up), 2 for active participation (spoke up in meetings), and 3 for power (had a formal position / power over decision making)

cor.test(data$particmgmt_tourism, data$particmgmt_coral,
         alternative = c("two.sided"),
         method = c("spearman"))
## Warning in cor.test.default(data$particmgmt_tourism, data$particmgmt_coral, :
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  data$particmgmt_tourism and data$particmgmt_coral
## S = 536050, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.7390672
data <- data %>% 
  mutate(particmgmt = pmax(particmgmt_coral,particmgmt_tourism)
  )

ggplot(data, aes(x=particmgmt)) +
  geom_bar()

Imputation of missing data

We used the ‘mice’ package (Buuren and Groothuis-Oudshoorn, 2011; Zhang, 2016) in R modeling software to impute missing values in our dataset. For the 230 reef tourism operator surveys in our dataset, we had missing data for savings (15), competing concerns (7), adaptation confidence (5), education (1), accessible dive / snorkel sites (1), and reef attachment (1). Without imputing missing data, we would have had to exclude 20 operators from our sample and this would have led to the exclusion of other relevant data for these operators. We imputed the missing data based on the rest of the dataset but we did not use the disturbance type, severity, topresponse as predictors for missing values. We only used the combined mgmt. participation indicator.

imp0 <- mice(data, maxit = 0) #First, do a set-up run of mice(), without any iterations (maxit = 0). This is to get all the objects that we may need to adjust in order to customize the imputation to our data.
## Warning: Number of logged events: 4
pred <- imp0$predictorMatrix #in the code below we tell the model which variables not to use for predicting missing values

pred[c("savings", "accsites", "edu", "reef_attach", "comp_concerns", "adapt_conf"), "dist_type"] <- 0
pred[c("savings", "accsites", "edu", "reef_attach", "comp_concerns", "adapt_conf"), "dist_severity"] <- 0
pred[c("savings", "accsites", "edu", "reef_attach", "comp_concerns", "adapt_conf"), "operator"] <- 0
pred[c("savings", "accsites", "edu", "reef_attach", "comp_concerns", "adapt_conf"), "topresponse"] <- 0
pred[c("savings", "accsites", "edu", "reef_attach", "comp_concerns", "adapt_conf"), "particmgmt_coral"] <- 0
pred[c("savings", "accsites", "edu", "reef_attach", "comp_concerns", "adapt_conf"), "particmgmt_tourism"] <- 0

imp <- mice(data, predictorMatrix = pred, m = 1, print = FALSE,meth='pmm', seed=123)  #impute data using predictive mean matching. We set a seed value to offset the random number generator (so that our imputed results have reproducibility, without seed there would be variations in the predictions each time we run the model)

densityplot(imp)   #compare the distribution of the imputed values against the distribution of the observed values.

imp$imp$savings    #check the imputed values for 'savings' indicator
imp$imp$comp_concerns
imp$imp$adapt_conf
imp$imp$edu
imp$imp$accsites
imp$imp$reef_attach
imputed_data <- complete(imp,1)   #The missing values have been replaced with the imputed values and we now have a complete dataset to work with in our consequent analysis. 

We collected data on the indicators in our surveys using multiple-choice categories to provide a consistent and directly comparable level of detail in the answers. Based on the distribution of the collected data, we transformed each indicator into a binary variable. While transforming some of our indicators into continuous or ordinal variables would likely have provided slightly more nuance, the nature of our research warranted the comparison of indicators at similar levels of complexity. We asked respondents about each indicator of adaptive capacity as they judged it to be before they were affected by the particular climate disturbance. Below we detail how we decided to cut-off ordinal and continuous into binary variables. These graphs do not include missing data point, which will be handled later.

For government effectiveness (‘gov_effect’): -0.37 for Fiji, -0.04 for Indonesia, 0.44 for Mariana Islands, 1.42 for the Hawaiian Islands, 1.42 for French Polynesia, 1.61 for Australia, and 1.81 for JP. For exclusivity of reef access (‘exclaccess’): 0 for open access (accessed same reef sites as everyone else), 1 for limited access (accessed sites that were shared with only a few other operators), and 2 for exclusive access (accessed sites where no one else was able to go). For participation in management (‘particmgmt_coral’): 0 for none, 1 for passive attendance (attended meetings but did not usually speak up), 2 for active participation (spoke up in meetings), and 3 for power (had a formal position / power over decision making). For savings, 1 for 0 – 3 months, 2 for 3 – 6 months, 3 for 6 – 12 months, and 4 for >12 months.

For boats (‘psgseats’): 1 for 0 – 10 seats, 2 for 10 – 20 seats, 3 for 20 – 50 seats, 4 for 50 – 100 seats, 5 for 100 – 200 seats, 6 for 200 – 300 seats, 7 for 300 – 400 seats, 8 for 400 – 500 seats, and 9 for >500 seats. For other services (‘othservices’): The number reflects the other types of products / services the operator was selling besides tourism activities (e.g. accommodation, first-aid training). For accessible dive / snorkel sites (‘accsites’): 1 for 0 – 2 sites, 2 for 2 – 5 sites, 3 for 5 – 10 sites, 4 for 10 – 20 sites, and 5 for >20 sites. For (GM) education (‘edu_ordinal’): 1 for secondary school, 2 for trade certification, 3 for graduate degree, and 4 for postgraduate degree.

For (GM) experience (‘experience’): 1 for 0 – 2 years, 2 for 2 – 5 years, 3 for 5 – 10 years, and 4 for >10 years. For employee involvement (‘employ_involve’): 0 for not involved, 1 for passively involved (they were involved prior to change), 2 for moderately involved (they were asked for, and voiced their opinions, and these were taken into account), and 3 for highly involved (they actively participated in the decision making process). For government ties (‘govt_ties’): 1 for weak ties (did not interact much with government stakeholder), 2 for moderate ties (had some interaction with government stakeholders, mainly through formal events), and 3 for strong ties (management or staff had informal ties with key government stakeholders). For reliance on research in decision making (‘research_rely’): 0 for not at all, 1 for low reliance, 2 for moderate reliance, and 3 for high reliance.

For industry memberships (‘ind_memberships’): The number reflects the number of tourism industry associations the operator was a member of. For reef attachment (‘reef_attach’): 0 for coral reef sites were not at all important in company’s identity, 2 for coral reef sites were somewhat important in company’s identity, and 3 for coral reef sites were major part of company’s identity. For competing concerns (‘comp_concerns’): 0 for no competing concerns, 1 for 1 – 2 competing concerns, 2 for 3 – 4 competing concerns, and 3 for >5 competing concerns. For adaptation confidence in ability to adapt to climate disturbances (‘adapt_conf’): 0 for not at all confident, 1 for low confidence, 2 for moderately confident, and 3 for highly confident.

We now apply these modifications to our data using the code below. We also create z-scored variables for the continuous predictors in our models: disturbance severity and distance to markets

data_modified <- imputed_data %>% 
  mutate(operator = factor(operator),
         resp_gender = factor(resp_gender),
         resp_region = factor (resp_region),
         resp_age = factor(resp_age),
         location = factor(location),
         gov_effect_binary = ifelse(gov_effect > 0.5, 1, 0),
         exclaccess_binary = ifelse(exclaccess > 0, 1, 0),
         particmgmt_binary = ifelse(particmgmt > 1, 1, 0),
         psgseats_binary = ifelse(psgseats > 2, 1, 0),
         savings_binary = ifelse(savings > 2, 1, 0),
         othservices_binary = ifelse(othservices > 0, 1, 0),
         accsites_binary = ifelse(accsites > 3, 1, 0),
         edu_binary = ifelse(edu > 2, 1, 0),
         experience_binary = ifelse(experience > 3, 1, 0),
         employ_involve_binary = ifelse(employ_involve > 1, 1, 0),
         ind_memberships_binary = ifelse(ind_memberships > 0, 1, 0),
         research_rely_binary = ifelse(research_rely > 0, 1, 0),
         govt_ties_binary = ifelse(govt_ties > 1, 1, 0),
         reef_attach_binary = ifelse(reef_attach > 1, 1, 0),
         comp_concerns_binary = ifelse(comp_concerns > 0, 1, 0),
         adapt_conf_binary = ifelse(adapt_conf > 1, 1, 0),
         z.dist_severity =   (dist_severity-mean(dist_severity))/(2*sd(dist_severity)),
         z.distmarket = (distmarket-mean(distmarket))/(2*sd(distmarket))
                                  )

Setting up the multinomial response model

We now calculate the partial correlations between the adaptive responses that operators adopted in response to climate disturbances. These partial correlations reflect whether particular responses were more frequently implemented together than others. We used Spearman’s Rank correlation because our responses are measures on a binary scale.

response_cor <- data_modified[, c("divbl","relief","support","insurance","chgopmode","chgact", "educate","chgsites","monitor","nrm","co2")]

pcor(response_cor,method="spearman")
## $estimate
##                  divbl       relief     support   insurance   chgopmode
## divbl      1.000000000 -0.020357860  0.07702295  0.01024531  0.20864800
## relief    -0.020357860  1.000000000  0.19952106 -0.02632262  0.05180185
## support    0.077022948  0.199521055  1.00000000 -0.01375645  0.02830850
## insurance  0.010245306 -0.026322622 -0.01375645  1.00000000 -0.07616094
## chgopmode  0.208648005  0.051801846  0.02830850 -0.07616094  1.00000000
## chgact     0.037796189  0.003162198  0.06836009 -0.08821421  0.35515039
## educate   -0.127592314 -0.114291373 -0.03133336  0.20047507  0.13208329
## chgsites   0.042165317  0.171738350 -0.06275669  0.18695861  0.18979091
## monitor    0.031665550 -0.107723152  0.07054858 -0.10563299 -0.03138290
## nrm       -0.000335073  0.046177719  0.28126118 -0.00538762 -0.10650199
## co2       -0.012193354  0.054519949 -0.03944470  0.05810986  0.10597687
##                 chgact     educate    chgsites     monitor          nrm
## divbl      0.037796189 -0.12759231  0.04216532  0.03166555 -0.000335073
## relief     0.003162198 -0.11429137  0.17173835 -0.10772315  0.046177719
## support    0.068360085 -0.03133336 -0.06275669  0.07054858  0.281261181
## insurance -0.088214211  0.20047507  0.18695861 -0.10563299 -0.005387620
## chgopmode  0.355150391  0.13208329  0.18979091 -0.03138290 -0.106501985
## chgact     1.000000000  0.09800120  0.24784670 -0.07271320  0.106511089
## educate    0.098001196  1.00000000 -0.02280847  0.19694840 -0.012888401
## chgsites   0.247846702 -0.02280847  1.00000000  0.19071162  0.099265614
## monitor   -0.072713196  0.19694840  0.19071162  1.00000000  0.253184790
## nrm        0.106511089 -0.01288840  0.09926561  0.25318479  1.000000000
## co2        0.006726194 -0.09624868 -0.05649930  0.24329687  0.326590348
##                    co2
## divbl     -0.012193354
## relief     0.054519949
## support   -0.039444704
## insurance  0.058109855
## chgopmode  0.105976873
## chgact     0.006726194
## educate   -0.096248682
## chgsites  -0.056499300
## monitor    0.243296870
## nrm        0.326590348
## co2        1.000000000
## 
## $p.value
##                 divbl      relief      support   insurance    chgopmode
## divbl     0.000000000 0.762923308 2.531067e-01 0.879349540 1.773768e-03
## relief    0.762923308 0.000000000 2.825634e-03 0.696498264 4.424939e-01
## support   0.253106676 0.002825634 0.000000e+00 0.838494957 6.748579e-01
## insurance 0.879349540 0.696498264 8.384950e-01 0.000000000 2.584727e-01
## chgopmode 0.001773768 0.442493875 6.748579e-01 0.258472688 0.000000e+00
## chgact    0.575362391 0.962632918 3.105928e-01 0.190366218 5.323994e-08
## educate   0.057680843 0.089345380 6.424078e-01 0.002693854 4.935562e-02
## chgsites  0.531983880 0.010363087 3.520125e-01 0.005196749 4.544065e-03
## monitor   0.638883591 0.109456883 2.953207e-01 0.116553592 6.418817e-01
## nrm       0.996039085 0.493650869 2.107227e-05 0.936379511 1.135607e-01
## co2       0.856636039 0.418892263 5.588015e-01 0.388875769 1.153620e-01
##                 chgact     educate     chgsites      monitor          nrm
## divbl     5.753624e-01 0.057680843 0.5319838798 0.6388835912 9.960391e-01
## relief    9.626329e-01 0.089345380 0.0103630867 0.1094568834 4.936509e-01
## support   3.105928e-01 0.642407764 0.3520125277 0.2953207303 2.107227e-05
## insurance 1.903662e-01 0.002693854 0.0051967492 0.1165535916 9.363795e-01
## chgopmode 5.323994e-08 0.049355623 0.0045440655 0.6418816509 1.135607e-01
## chgact    0.000000e+00 0.145545398 0.0001912091 0.2807149326 1.135297e-01
## educate   1.455454e-01 0.000000000 0.7353902986 0.0032106414 8.485595e-01
## chgsites  1.912091e-04 0.735390299 0.0000000000 0.0043483617 1.403958e-01
## monitor   2.807149e-01 0.003210641 0.0043483617 0.0000000000 1.370727e-04
## nrm       1.135297e-01 0.848559511 0.1403958382 0.0001370727 0.000000e+00
## co2       9.206195e-01 0.152918672 0.4021780239 0.0002524891 6.500265e-07
##                    co2
## divbl     8.566360e-01
## relief    4.188923e-01
## support   5.588015e-01
## insurance 3.888758e-01
## chgopmode 1.153620e-01
## chgact    9.206195e-01
## educate   1.529187e-01
## chgsites  4.021780e-01
## monitor   2.524891e-04
## nrm       6.500265e-07
## co2       0.000000e+00
## 
## $statistic
##                  divbl      relief    support   insurance  chgopmode
## divbl      0.000000000 -0.30201846  1.1458389  0.15197042  3.1643958
## relief    -0.302018459  0.00000000  3.0200990 -0.39056291  0.7693785
## support    1.145838866  3.02009900  0.0000000 -0.20406038  0.4200513
## insurance  0.151970418 -0.39056291 -0.2040604  0.00000000 -1.1329399
## chgopmode  3.164395843  0.76937853  0.4200513 -1.13293988  0.0000000
## chgact     0.561008940  0.04690321  1.0163214 -1.31354903  5.6350876
## educate   -1.908095283 -1.70639656 -0.4649771  3.03514292  1.9764280
## chgsites   0.625969422  2.58570829 -0.9326706  2.82281677  2.8671662
## monitor    0.469911653 -1.60714464  1.0490183 -1.57560570 -0.4657130
## nrm       -0.004969936  0.68565769  4.3472714 -0.07991248 -1.5887156
## co2       -0.180870115  0.80986605 -0.5855152  0.86336736  1.5807931
##                chgact    educate   chgsites    monitor          nrm         co2
## divbl      0.56100894 -1.9080953  0.6259694  0.4699117 -0.004969936 -0.18087012
## relief     0.04690321 -1.7063966  2.5857083 -1.6071446  0.685657692  0.80986605
## support    1.01632139 -0.4649771 -0.9326706  1.0490183  4.347271371 -0.58551518
## insurance -1.31354903  3.0351429  2.8228168 -1.5756057 -0.079912476  0.86336736
## chgopmode  5.63508758  1.9764280  2.8671662 -0.4657130 -1.588715557  1.58079314
## chgact     0.00000000  1.4606236  3.7945535 -1.0813735  1.588852923  0.09976783
## educate    1.46062365  0.0000000 -0.3383923  2.9795752 -0.191181764 -1.43425746
## chgsites   3.79455355 -0.3383923  0.0000000  2.8815991  1.479655056 -0.83936081
## monitor   -1.08137349  2.9795752  2.8815991  0.0000000  3.881814847  3.72046889
## nrm        1.58885292 -0.1911818  1.4796551  3.8818148  0.000000000  5.12515025
## co2        0.09976783 -1.4342575 -0.8393608  3.7204689  5.125150248  0.00000000
## 
## $n
## [1] 231
## 
## $gp
## [1] 9
## 
## $method
## [1] "spearman"

Response clusters definitions

data_modified <- data_modified %>% 
  mutate(top_mn = case_when(
    topresponse == 'none' | topresponse == 'relief' ~ 'cope',
    topresponse == 'chgsites' | topresponse == 'chgact' | topresponse == 'chgopmode' |    topresponse == 'divbl' | topresponse == 'insurance' | topresponse == 'monitor' | topresponse == 'educate' ~   'adapt',
    topresponse == 'nrm' | topresponse == 'support'| topresponse == 'co2' ~ 'transform'),
    top_mn = factor(top_mn,levels = c('cope', 'adapt', 'transform'))
    )

ggplot(data_modified, aes(x=top_mn)) +
  geom_bar()

table(data_modified$top_mn)
## 
##      cope     adapt transform 
##        63        92        76

Model formula:

The relationship between adaptive capacity and the primary responses to climate disturbances were explored using multinomial logistic regression in a Bayesian framework (Gelman et al., 2013). The Bayesian approach for binary regression extends to multinomial models. We used a baseline-category logit model in which the primary response of ‘coping measures (e.g. doing nothing or looking for relief)’ was used as the baseline category. The likelihood of other response clusters (‘adaptive’ and ‘transformative’ measures) were then estimated against the baseline category. As such with sample size ni and αi representing the probability of an operator either adopting a coping (αi1), adaptive (αi2), or transformative (αi3) primary response to a climate disturbance. The multinomial generalized linear model is parameterized in terms of the logarithm of the ratio of the probability of each category relative to that of our baseline category (coping: αi1), where β(j) is a vector of parameters for each category (Gelman et al., 2013, p. 426). We used weakly-informative normal priors on all population effects. Although we had some indications of the effect of our predictor variables on adaptive responses (Bartelet et al., 2022a), we decided to not use informative priors in our Bayesian models because most of the prior knowledge providing the foundation for our study is based on a different context, i.e., smallholder farming in Asia and Africa.

\[ \begin{align} y_i &\sim{} \mathcal{Multin}(n_i;α_i1, α_i2, α_i3)\\ with\ \alpha_i1 + \alpha_i2 + \alpha_i3 = 1\\ \log(\alpha_i2/\alpha_i1) = X_i \beta^j \\ \log(\alpha_i3/\alpha_i1) = X_i \beta^j \\ \beta_0 &\sim{} \mathcal{N}(0,3)\\ \beta_{1-19} &\sim{} \mathcal{N}(0,3)\\ \end{align} \] ## Test for collinearity

We evaluated whether there was potential for collinearity issues in our model by running a frequentist statistical model (logistic regression), in which we predict the likelihood of transformative action as primary responses versus doing something else. We found that all variance inflation factors (VIFs) are below the critical threshold of ‘3’.

data_modified <- data_modified %>% 
  mutate(top_transform = ifelse(top_mn == 'transform', 1, 0)
  )

model.glm <- glm(top_transform ~ gov_effect_binary + z.distmarket + 
    dist_type + z.dist_severity +
    exclaccess_binary + particmgmt_binary + 
    psgseats_binary + savings_binary + 
    othservices_binary + accsites_binary + 
    edu_binary + experience_binary + employ_involve_binary +
    ind_memberships_binary + research_rely_binary + govt_ties_binary +
    reef_attach_binary + comp_concerns_binary + adapt_conf_binary,
                 data=data_modified, family=binomial(link='logit'))

check_collinearity(model.glm)

Fit the Bayesian multinomial logistic regression model

We now fit the Bayesian model using the ‘categorical’ data distribution family reflecting the multinomial nature of our response data. We set adapt_delta to 0.99 (from the default of 0.95) to reduce the step size in our models which will lead to lower likelihood of divergences. This will usually slow down the sampling but lead to more robust samples.

Setting one prior on all effects makes the computation much easier. We use iteratively-determined priors by running the model first only with priors and seeing whether the outcomes are non-informative and covering the full ribbon (0 to 100%). We use zero slope coefficients to not prime the model on either positive or negative effects. We also use a reasonably wide variance (3) that is reasonable for logistic regression.

priors <- prior(normal(0, 3), class = 'Intercept', dpar = 'muadapt') +
          prior(normal(0, 3), class = 'Intercept', dpar = 'mutransform') +
          prior(normal(0,3), class = 'b', dpar = 'muadapt') +
          prior(normal(0,3), class = 'b', dpar = 'mutransform')

We first run the model with priors and find that our selected priors are uninformative.

mod_top_multinomial = brms::brm(  
  top_mn ~ gov_effect_binary + z.distmarket + 
    dist_type + z.dist_severity +
    exclaccess_binary + particmgmt_binary + 
    psgseats_binary + savings_binary + 
    othservices_binary + accsites_binary + 
    edu_binary + experience_binary + employ_involve_binary +
    ind_memberships_binary + research_rely_binary + govt_ties_binary +
    reef_attach_binary + comp_concerns_binary + adapt_conf_binary,
  data = data_modified,
  prior = priors,
  sample_prior='only',
  family = 'categorical',
  iter = 5000, warmup = 1000,
  chains = 3, cores = 3,
  thin = 5,
  refresh = 0,
  control = list(adapt_delta = 0.99),
  seed = 123)
## Compiling Stan program...
## Start sampling
mod_top_multinomial %>% ggpredict() %>% plot(add.data=TRUE, jitter=FALSE)  
## $gov_effect_binary

## 
## $z.distmarket

## 
## $dist_type

## 
## $z.dist_severity

## 
## $exclaccess_binary

## 
## $particmgmt_binary

## 
## $psgseats_binary

## 
## $savings_binary

## 
## $othservices_binary

## 
## $accsites_binary

## 
## $edu_binary

## 
## $experience_binary

## 
## $employ_involve_binary

## 
## $ind_memberships_binary

## 
## $research_rely_binary

## 
## $govt_ties_binary

## 
## $reef_attach_binary

## 
## $comp_concerns_binary

## 
## $adapt_conf_binary

After checking the priors, run the Bayesian model using both prior and data information.

mod_top_multinomial <- mod_top_multinomial %>% 
  update(sample_prior='yes', refresh = 0)
## The desired updates require recompiling the model
## Compiling Stan program...
## Start sampling

Model validation

preds <- posterior_predict(mod_top_multinomial, ndraws=250,summary=FALSE)

fert.resids <- createDHARMa(simulatedResponse=t(preds),
                            observedResponse = as.numeric(as.factor(data_modified$top_mn)), 
                            fittedPredictedResponse=apply(preds,2,median), 
                            integerResponse=FALSE)
testResiduals (fert.resids)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.040303, p-value = 0.8472
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0138, p-value = 0.872
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 231, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0001095949 0.0238823445
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                            0.004329004
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.040303, p-value = 0.8472
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0138, p-value = 0.872
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 231, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0001095949 0.0238823445
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                            0.004329004
plot(fert.resids)

We tested model convergence and diagnostics (traceplots, autocorrelation, convergence, and effective sample size). Tests reveal no problems with our model convergence.

pars <- mod_top_multinomial$fit %>% get_variables()

pars_adapt <- str_extract(pars, '^b_muadapt_.*') %>% na.omit()

pars_transform <- str_extract(pars, '^b_mutransform_.*') %>% na.omit()

mod_top_multinomial$fit %>% stan_trace(pars=pars_adapt, inc_warmup=TRUE)

mod_top_multinomial$fit %>% stan_trace(pars=pars_transform, inc_warmup=TRUE)

mod_top_multinomial$fit %>% stan_ac(pars=pars_adapt) #autocorrelation test (you only want autocorrelation at beginning of plot. If evidence of autocorrelation, use more thinning. If lag 2 is greater than 0.2 you might have autocorrelation, adjust tinning to get beyond 0.2)

mod_top_multinomial$fit %>% stan_ac(pars=pars_transform)

mod_top_multinomial$fit %>% stan_rhat()            #measure of convergence between chains. We want values to be less than 1.03 or 1.01. Values greater than 1.05 shows that chains were not consistent, they did not converge. There is a bar for every parameter. If problem with converge, run it longer and potentially narrow the priors. 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mod_top_multinomial$fit %>% stan_ess()    #effective sample size. If chains don’t converge a lot of samples will be thrown away. This tests the efficiency of the sampler. If too many samples are rejected you don’t them contributing. Above 50% effective sample size is high enough. Less than 50% means sampler has drifted into unsupported areas. You need to tighten the priors.            
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

test whether prior and posterior match for different predictors and responses. For example:

#mod_top_multinomial %>% get_variables()

mod_top_multinomial %>% hypothesis('muadapt_savings_binary=0') %>% plot()

mod_top_multinomial %>% hypothesis('mutransform_comp_concerns_binary=0') %>% plot()

Analysis of results

We use median highest probability density intervals and present credibility intervals at both 80% (weak evidence) and 95% (strong evidence).

Strong evidence: log odds scale

mod_top_multinomial %>% tidy_draws() %>% dplyr::select(starts_with("b_")) %>%
  summarise_draws(median,~HDInterval::hdi(.x, credMass = 0.95),rhat,ess_bulk,ess_tail) %>% knitr::kable()
variable median lower upper rhat ess_bulk ess_tail
b_muadapt_Intercept 0.7896507 -0.5707917 2.0941099 1.0008616 2571.411 2388.281
b_mutransform_Intercept 0.0670525 -1.2767953 1.5644713 1.0002720 2411.670 1700.096
b_muadapt_gov_effect_binary -0.5225526 -1.3643373 0.3846161 0.9996307 2421.821 2402.461
b_muadapt_z.distmarket -0.6829090 -1.5891268 0.1761063 1.0013428 2207.214 2446.568
b_muadapt_dist_type -0.8526734 -2.0980241 0.2723958 1.0001378 2283.985 2403.939
b_muadapt_z.dist_severity 1.8699248 1.0372663 2.8064329 1.0002079 1998.903 2232.692
b_muadapt_exclaccess_binary 0.1784635 -0.7582171 1.0884664 1.0000032 2526.078 2491.499
b_muadapt_particmgmt_binary 0.1266520 -0.7636110 1.0242901 1.0024361 2573.241 2028.750
b_muadapt_psgseats_binary -0.3566502 -1.1862244 0.5644566 0.9998594 2515.490 2361.012
b_muadapt_savings_binary 0.2700554 -0.4824477 1.0701633 0.9999330 2304.356 2028.181
b_muadapt_othservices_binary 0.3527455 -0.4330590 1.1631295 0.9995843 2376.174 2321.149
b_muadapt_accsites_binary -0.3448160 -1.1143163 0.4192817 0.9997169 2161.249 2021.134
b_muadapt_edu_binary -0.1493597 -0.9888886 0.6771764 1.0003925 2498.604 2489.500
b_muadapt_experience_binary -0.6657694 -1.4866936 0.1902605 0.9997352 2341.851 2209.791
b_muadapt_employ_involve_binary -0.1325657 -0.8930583 0.7040542 1.0002235 2288.784 2282.263
b_muadapt_ind_memberships_binary 0.8142634 -0.0118199 1.6387235 1.0002858 2289.751 2251.340
b_muadapt_research_rely_binary 0.5983230 -0.2101715 1.5074950 0.9997970 2449.534 2403.939
b_muadapt_govt_ties_binary -0.3869888 -1.3596592 0.5281761 1.0002337 2317.363 2301.522
b_muadapt_reef_attach_binary -0.2293992 -0.9913047 0.5620005 1.0002022 2290.691 2289.754
b_muadapt_comp_concerns_binary 0.0224119 -0.7269714 0.8204324 1.0006043 2255.706 2281.905
b_muadapt_adapt_conf_binary 0.0550231 -0.8848561 0.9381379 0.9996399 2478.925 2026.598
b_mutransform_gov_effect_binary -1.1892413 -2.2207860 -0.2254530 1.0004771 2474.323 2185.953
b_mutransform_z.distmarket 0.4316454 -0.4663618 1.3228894 0.9996506 2395.812 2397.392
b_mutransform_dist_type -1.1525783 -2.4497045 0.2416389 0.9998151 2378.624 2245.774
b_mutransform_z.dist_severity 0.9272731 0.0018739 1.8896830 1.0002159 2255.098 2304.976
b_mutransform_exclaccess_binary -0.2385664 -1.2819340 0.7744782 0.9999840 2483.451 2125.298
b_mutransform_particmgmt_binary 0.5052244 -0.4361672 1.4344391 1.0004089 2445.127 2283.747
b_mutransform_psgseats_binary 0.6958288 -0.1800360 1.6172336 0.9999440 2400.006 2208.597
b_mutransform_savings_binary 0.0103069 -0.8387588 0.8880051 0.9996355 2234.431 2242.833
b_mutransform_othservices_binary -0.0692409 -0.8909545 0.7944451 0.9997495 2473.133 2471.211
b_mutransform_accsites_binary 0.0485982 -0.8703754 0.9321499 1.0001868 2092.162 2051.294
b_mutransform_edu_binary -0.3284547 -1.1369170 0.6146089 1.0007582 2557.956 2248.297
b_mutransform_experience_binary -0.2103759 -1.1678801 0.6649130 1.0011946 2248.228 2020.691
b_mutransform_employ_involve_binary 0.7677624 -0.0591760 1.7498694 1.0001015 2339.526 2263.377
b_mutransform_ind_memberships_binary 1.0782039 0.2372243 1.9834708 0.9998260 2314.764 2362.630
b_mutransform_research_rely_binary 1.0827900 0.1683884 1.9822058 0.9999343 2484.064 2211.670
b_mutransform_govt_ties_binary 0.7790557 -0.1886954 1.7402248 0.9996178 2452.987 2282.263
b_mutransform_reef_attach_binary 0.1680752 -0.6719244 1.0598517 0.9997563 2251.871 2099.048
b_mutransform_comp_concerns_binary -0.6072639 -1.5106822 0.2136596 1.0001265 2327.044 2234.775
b_mutransform_adapt_conf_binary -1.2904208 -2.2719004 -0.3514563 0.9997993 2388.340 2362.428

Strong evidence: odds scale

mod_top_multinomial %>% tidy_draws() %>% exp() %>% dplyr::select(starts_with("b_")) %>%
  summarise_draws(median,~HDInterval::hdi(.x, credMass = 0.95),rhat,ess_bulk,ess_tail) %>% knitr::kable()
variable median lower upper rhat ess_bulk ess_tail
b_muadapt_Intercept 2.2026269 0.3015927 7.0769799 1.0013365 2571.411 2388.281
b_mutransform_Intercept 1.0693516 0.1339271 3.4688070 1.0003494 2411.670 1700.096
b_muadapt_gov_effect_binary 0.5930049 0.1543910 1.2597059 0.9996492 2421.821 2402.461
b_muadapt_z.distmarket 0.5051455 0.1547902 1.0848722 1.0010931 2207.214 2446.568
b_muadapt_dist_type 0.4262739 0.0629299 1.1863878 1.0005549 2283.985 2403.939
b_muadapt_z.dist_severity 6.4878151 1.8870845 14.1827026 1.0002079 1998.903 2232.692
b_muadapt_exclaccess_binary 1.1953792 0.3277001 2.5570613 1.0000032 2526.078 2491.499
b_muadapt_particmgmt_binary 1.1350221 0.3403974 2.4776616 1.0023723 2573.241 2028.750
b_muadapt_psgseats_binary 0.7000173 0.2142925 1.4396566 0.9998594 2515.490 2361.012
b_muadapt_savings_binary 1.3100375 0.4859544 2.6177932 0.9999330 2304.356 2028.181
b_muadapt_othservices_binary 1.4229689 0.5265716 2.8821483 0.9995842 2376.174 2321.149
b_muadapt_accsites_binary 0.7083507 0.2606377 1.4033812 0.9997169 2161.249 2021.134
b_muadapt_edu_binary 0.8612592 0.3092182 1.8009614 1.0002441 2498.604 2489.500
b_muadapt_experience_binary 0.5138780 0.1774303 1.0255894 0.9997352 2341.851 2209.791
b_muadapt_employ_involve_binary 0.8758454 0.3236181 1.8358319 1.0002879 2288.784 2282.263
b_muadapt_ind_memberships_binary 2.2575123 0.7707048 4.6246852 1.0003575 2289.751 2251.340
b_muadapt_research_rely_binary 1.8190658 0.6053523 3.8275378 0.9997970 2449.534 2403.939
b_muadapt_govt_ties_binary 0.6790987 0.1658744 1.5293200 1.0002576 2317.363 2301.522
b_muadapt_reef_attach_binary 0.7950112 0.2853164 1.5526568 1.0002867 2290.691 2289.754
b_muadapt_comp_concerns_binary 1.0226649 0.4306755 2.0995293 1.0006932 2255.706 2281.905
b_muadapt_adapt_conf_binary 1.0565652 0.2554295 2.2907905 0.9996399 2478.925 2026.598
b_mutransform_gov_effect_binary 0.3044522 0.0746546 0.7009630 1.0005021 2474.323 2185.953
b_mutransform_z.distmarket 1.5397890 0.4553948 3.2913945 0.9996506 2395.812 2397.392
b_mutransform_dist_type 0.3158214 0.0533050 0.9960727 0.9996653 2378.624 2245.774
b_mutransform_z.dist_severity 2.5276075 0.7277176 5.7413034 1.0002159 2255.098 2304.976
b_mutransform_exclaccess_binary 0.7877564 0.2160790 1.9336998 0.9999840 2483.451 2125.298
b_mutransform_particmgmt_binary 1.6573575 0.4917394 3.7234518 1.0004440 2445.127 2283.747
b_mutransform_psgseats_binary 2.0053704 0.5681610 4.2707586 0.9999606 2400.006 2208.597
b_mutransform_savings_binary 1.0103602 0.3254636 2.0607464 0.9996147 2234.431 2242.833
b_mutransform_othservices_binary 0.9331019 0.3145282 1.9634308 0.9997495 2473.133 2471.211
b_mutransform_accsites_binary 1.0497986 0.3207284 2.3042775 1.0000160 2092.162 2051.294
b_mutransform_edu_binary 0.7200356 0.1866646 1.5573561 1.0007582 2557.956 2248.297
b_mutransform_experience_binary 0.8102796 0.2202343 1.7301554 1.0012981 2248.228 2020.691
b_mutransform_employ_involve_binary 2.1549389 0.6031091 4.7117546 1.0001244 2339.526 2263.377
b_mutransform_ind_memberships_binary 2.9393961 0.8296821 6.3490269 0.9998909 2314.764 2362.630
b_mutransform_research_rely_binary 2.9529067 0.9164544 6.4130129 0.9997787 2484.064 2211.670
b_mutransform_govt_ties_binary 2.1794135 0.6798297 5.2131803 0.9996115 2452.987 2282.263
b_mutransform_reef_attach_binary 1.1830256 0.3851240 2.5246507 0.9997563 2251.871 2099.048
b_mutransform_comp_concerns_binary 0.5448396 0.1892597 1.1427649 1.0000649 2327.044 2234.775
b_mutransform_adapt_conf_binary 0.2751550 0.0687799 0.6035765 0.9997993 2388.340 2362.428

Weak evidence: log odds scale

mod_top_multinomial %>% tidy_draws() %>% dplyr::select(starts_with("b_")) %>%
  summarise_draws(median,~HDInterval::hdi(.x, credMass = 0.80),rhat,ess_bulk,ess_tail) %>% knitr::kable()
variable median lower upper rhat ess_bulk ess_tail
b_muadapt_Intercept 0.7896507 -0.0380937 1.7356287 1.0008616 2571.411 2388.281
b_mutransform_Intercept 0.0670525 -0.8337771 0.9951409 1.0002720 2411.670 1700.096
b_muadapt_gov_effect_binary -0.5225526 -1.1205329 0.0350963 0.9996307 2421.821 2402.461
b_muadapt_z.distmarket -0.6829090 -1.1981734 -0.0546751 1.0013428 2207.214 2446.568
b_muadapt_dist_type -0.8526734 -1.6025678 -0.0187438 1.0001378 2283.985 2403.939
b_muadapt_z.dist_severity 1.8699248 1.2915806 2.4347703 1.0002079 1998.903 2232.692
b_muadapt_exclaccess_binary 0.1784635 -0.4386495 0.7553614 1.0000032 2526.078 2491.499
b_muadapt_particmgmt_binary 0.1266520 -0.4600502 0.6838187 1.0024361 2573.241 2028.750
b_muadapt_psgseats_binary -0.3566502 -0.8711731 0.2383807 0.9998594 2515.490 2361.012
b_muadapt_savings_binary 0.2700554 -0.2215155 0.7723556 0.9999330 2304.356 2028.181
b_muadapt_othservices_binary 0.3527455 -0.1768710 0.8669532 0.9995843 2376.174 2321.149
b_muadapt_accsites_binary -0.3448160 -0.8640243 0.1637284 0.9997169 2161.249 2021.134
b_muadapt_edu_binary -0.1493597 -0.6583141 0.4133722 1.0003925 2498.604 2489.500
b_muadapt_experience_binary -0.6657694 -1.1474395 -0.0813802 0.9997352 2341.851 2209.791
b_muadapt_employ_involve_binary -0.1325657 -0.6602729 0.3924227 1.0002235 2288.784 2282.263
b_muadapt_ind_memberships_binary 0.8142634 0.2848200 1.3499991 1.0002858 2289.751 2251.340
b_muadapt_research_rely_binary 0.5983230 0.0317533 1.1710166 0.9997970 2449.534 2403.939
b_muadapt_govt_ties_binary -0.3869888 -0.9452191 0.2472671 1.0002337 2317.363 2301.522
b_muadapt_reef_attach_binary -0.2293992 -0.7550111 0.2277267 1.0002022 2290.691 2289.754
b_muadapt_comp_concerns_binary 0.0224119 -0.5096420 0.5048475 1.0006043 2255.706 2281.905
b_muadapt_adapt_conf_binary 0.0550231 -0.4734947 0.7062899 0.9996399 2478.925 2026.598
b_mutransform_gov_effect_binary -1.1892413 -1.7976446 -0.4952732 1.0004771 2474.323 2185.953
b_mutransform_z.distmarket 0.4316454 -0.1568487 1.0036763 0.9996506 2395.812 2397.392
b_mutransform_dist_type -1.1525783 -2.0418327 -0.3024858 0.9998151 2378.624 2245.774
b_mutransform_z.dist_severity 0.9272731 0.3217585 1.5276140 1.0002159 2255.098 2304.976
b_mutransform_exclaccess_binary -0.2385664 -0.8939594 0.4510119 0.9999840 2483.451 2125.298
b_mutransform_particmgmt_binary 0.5052244 -0.1062662 1.1126471 1.0004089 2445.127 2283.747
b_mutransform_psgseats_binary 0.6958288 0.1060918 1.2609721 0.9999440 2400.006 2208.597
b_mutransform_savings_binary 0.0103069 -0.5187061 0.6019658 0.9996355 2234.431 2242.833
b_mutransform_othservices_binary -0.0692409 -0.5846318 0.5337172 0.9997495 2473.133 2471.211
b_mutransform_accsites_binary 0.0485982 -0.5921316 0.6011409 1.0001868 2092.162 2051.294
b_mutransform_edu_binary -0.3284547 -0.9333058 0.2339310 1.0007582 2557.956 2248.297
b_mutransform_experience_binary -0.2103759 -0.8034390 0.3830542 1.0011946 2248.228 2020.691
b_mutransform_employ_involve_binary 0.7677624 0.1737649 1.3233784 1.0001015 2339.526 2263.377
b_mutransform_ind_memberships_binary 1.0782039 0.4694498 1.6282849 0.9998260 2314.764 2362.630
b_mutransform_research_rely_binary 1.0827900 0.4311505 1.5979607 0.9999343 2484.064 2211.670
b_mutransform_govt_ties_binary 0.7790557 0.1154874 1.4110606 0.9996178 2452.987 2282.263
b_mutransform_reef_attach_binary 0.1680752 -0.4217511 0.7080049 0.9997563 2251.871 2099.048
b_mutransform_comp_concerns_binary -0.6072639 -1.1572139 -0.0521454 1.0001265 2327.044 2234.775
b_mutransform_adapt_conf_binary -1.2904208 -1.9138298 -0.6752319 0.9997993 2388.340 2362.428

Conditional effects

mod_top_multinomial %>% conditional_effects(categorical = TRUE) %>% plot(points=TRUE)