Socio-economic trends in MPAs

1 Load libraries

Code

library(here)
library(glmmTMB)
library(tidyr)
library(tidyverse)
library(corrplot)
library(Hmisc)
library(ggeffects)
library(tinytable)
library(emmeans)
library(ggpp)
library(gtools)
library(DescTools)

2 Load data

Code

MPA_data <- read.csv(here("Data/MPA.csv"))
names(MPA_data)

 [1] "fid"                                 
 [2] "WDPAID"                              
 [3] "DESIG_TYPE"                          
 [4] "IUCN_CAT"                            
 [5] "NO_TK_AREA"                          
 [6] "GOV_TYPE"                            
 [7] "OWN_TYPE"                            
 [8] "ISO3"                                
 [9] "ADM0_A3"                             
[10] "INCOME_GRP"                          
[11] "ECONOMY"                             
[12] "CONTINENT"                           
[13] "SUBREGION"                           
[14] "REGION_WB"                           
[15] "removal_of_marine_life_is_prohibited"
[16] "category_name"                       
[17] "protection_focus"                    
[18] "designation"                         
[19] "lat"                                 
[20] "lon"                                 
[21] "pop_power_transform"                 
[22] "pop_raw"                             
[23] "mig_power_transform"                 
[24] "mig_raw"                             
[25] "gdp_power_transform"                 
[26] "gdp_raw"                             
[27] "age_power_transform"                 
[28] "age_raw"                             
[29] "fh_power_transform"                  
[30] "fh_raw"                              
[31] "GIS_M_AREA_power_transform"          
[32] "GIS_M_AREA_raw"                      
[33] "directional_balance_index"           
[34] "signed_magnitude_index"

Code

MPA_data <- MPA_data %>%
  mutate(Population = pop_power_transform,
         Migration = mig_power_transform,
         GDP = gdp_power_transform,
         Age = age_power_transform,
         Fishing_hours = fh_power_transform,
         MPA_area = GIS_M_AREA_raw,
         SMI = signed_magnitude_index,
         DBI = directional_balance_index)

Little data descriptor:

IUCN_CAT: I = more restrictive VI = less restrictive

DBI: +1 = change in the positive direction -1 ~ change in the negative direction

SMI: +1 = change is positive (e.g., values increased in magnitude) -1 = change is negative (e.g., values decreased in magnitude)

Code

MPA_data <- MPA_data %>%
  mutate_if(is.character, as.factor)

categorical_variables <- sapply(MPA_data,class)
categorical_variables <- categorical_variables[which(categorical_variables=="factor")]
categorical_variables <- names(categorical_variables)

lapply(categorical_variables, function(x){
  MPA_data %>%
    group_by(!!sym(x)) %>%
    tally()
})

[[1]]
# A tibble: 4 × 2
  DESIG_TYPE         n
  <fct>          <int>
1 International     78
2 National        3654
3 Not Applicable    15
4 Regional        1055

[[2]]
# A tibble: 10 × 2
   IUCN_CAT           n
   <fct>          <int>
 1 Ia               487
 2 Ib                63
 3 II               198
 4 III              189
 5 IV              1310
 6 Not Applicable    20
 7 Not Assigned     353
 8 Not Reported    1419
 9 V                399
10 VI               364

[[3]]
# A tibble: 10 × 2
   GOV_TYPE                                   n
   <fct>                                  <int>
 1 Collaborative governance                 282
 2 Federal or national ministry or agency  2921
 3 Government-delegated management           14
 4 Indigenous peoples                        34
 5 Individual landowners                     10
 6 Joint governance                          39
 7 Local communities                         29
 8 Non-profit organisations                 143
 9 Not Reported                             302
10 Sub-national ministry or agency         1028

[[4]]
# A tibble: 7 × 2
  OWN_TYPE                     n
  <fct>                    <int>
1 Communal                    11
2 Individual landowners       44
3 Joint ownership             19
4 Multiple ownership          43
5 Non-profit organisations   130
6 Not Reported              3473
7 State                     1082

[[5]]
# A tibble: 134 × 2
   ISO3      n
   <fct> <int>
 1 AGO       1
 2 AIA       2
 3 ALA       6
 4 ALB       2
 5 ARE       7
 6 ARG      10
 7 ASM       6
 8 ATG       7
 9 AUS     292
10 AZE       1
# ℹ 124 more rows

[[6]]
# A tibble: 125 × 2
   ADM0_A3     n
   <fct>   <int>
 1 ""         71
 2 "AGO"       1
 3 "AIA"       2
 4 "ALB"       2
 5 "ARE"       7
 6 "ARG"      10
 7 "ASM"       6
 8 "ATG"       7
 9 "AUS"     292
10 "AZE"       1
# ℹ 115 more rows

[[7]]
# A tibble: 6 × 2
  INCOME_GRP                    n
  <fct>                     <int>
1 ""                           71
2 "1. High income: OECD"     3777
3 "2. High income: nonOECD"   189
4 "3. Upper middle income"    472
5 "4. Lower middle income"    234
6 "5. Low income"              59

[[8]]
# A tibble: 8 × 2
  ECONOMY                          n
  <fct>                        <int>
1 ""                              71
2 "1. Developed region: G7"     1540
3 "2. Developed region: nonG7"  2241
4 "3. Emerging region: BRIC"      46
5 "4. Emerging region: MIKT"     240
6 "5. Emerging region: G20"      158
7 "6. Developing region"         414
8 "7. Least developed region"     92

[[9]]
# A tibble: 8 × 2
  CONTINENT                     n
  <fct>                     <int>
1 ""                           71
2 "Africa"                    157
3 "Asia"                      474
4 "Europe"                   2375
5 "North America"            1016
6 "Oceania"                   606
7 "Seven seas (open ocean)"    29
8 "South America"              74

[[10]]
# A tibble: 23 × 2
   SUBREGION                       n
   <fct>                       <int>
 1 ""                              1
 2 "Australia and New Zealand"   479
 3 "Caribbean"                   189
 4 "Central America"              66
 5 "Central Asia"                  4
 6 "Eastern Africa"               99
 7 "Eastern Asia"                273
 8 "Eastern Europe"               42
 9 "Melanesia"                    66
10 "Micronesia"                   33
# ℹ 13 more rows

[[11]]
# A tibble: 8 × 2
  REGION_WB                        n
  <fct>                        <int>
1 ""                              71
2 "East Asia & Pacific"         1037
3 "Europe & Central Asia"       2357
4 "Latin America & Caribbean"    290
5 "Middle East & North Africa"    70
6 "North America"                800
7 "South Asia"                     9
8 "Sub-Saharan Africa"           168

[[12]]
# A tibble: 9 × 2
  category_name                         n
  <fct>                             <int>
1 ""                                   14
2 "Fisheries Management Area"        1692
3 "IUCN MPA"                         1487
4 "Jurisdictional Authority Area"    1158
5 "OECM"                               20
6 "Other"                             346
7 "Recreational Area"                  28
8 "Vessel Restricted Area"             45
9 "Water Quality/Human Health Area"    12

[[13]]
# A tibble: 12 × 2
   protection_focus                      n
   <fct>                             <int>
 1 ""                                   21
 2 "Biodiversity Protection"            21
 3 "Cultural Heritage"                  11
 4 "Ecosystem"                        3750
 5 "Ecosystem "                          9
 6 "Ecosystem and Cultural Heritage"     1
 7 "Ecosystem, focal species"            3
 8 "Ecosysytem"                         16
 9 "Ecosytem"                            2
10 "Focal species"                     191
11 "Focal Species"                     751
12 "Sovereignty"                        26

[[14]]
# A tibble: 310 × 2
   designation                                                                 n
   <fct>                                                                   <int>
 1 ""                                                                         14
 2 "Absolute Nature Reserve"                                                   1
 3 "Aquatic National Park"                                                     1
 4 "Aquatic Nature Sanctuary"                                                  5
 5 "Aquatic Nature Sanctuary - Sustainable Fisheries Zone"                     1
 6 "Aquatic Park"                                                              3
 7 "Aquatic Reserve"                                                          12
 8 "Aquatic Tourism Park"                                                      7
 9 "Area Marina Protetta"                                                      1
10 "Area Marina Protetta; Riserva Naturale Marina e Aree Naturali Marine …     1
# ℹ 300 more rows

Why there are so few MPAs in Brazil?

Code

# classify IUCN_CAT
table(MPA_data$IUCN_CAT)


            Ia             Ib             II            III             IV 
           487             63            198            189           1310 
Not Applicable   Not Assigned   Not Reported              V             VI 
            20            353           1419            399            364

Code

MPA_data$IUCN_CAT <- factor(MPA_data$IUCN_CAT, levels = c("Ia","Ib","II","III","IV","V","VI"))

ggplot(MPA_data %>% 
         group_by(IUCN_CAT) %>%
         tally(),
       aes(x = IUCN_CAT, y = n)) +
  geom_bar(stat = "identity") +
  labs(x = "IUCN category of MPA", y = "Count")+
  theme_classic()

Code

# classify income group
MPA_data %>% 
  group_by(INCOME_GRP) %>%
  tally()

# A tibble: 6 × 2
  INCOME_GRP                    n
  <fct>                     <int>
1 ""                           71
2 "1. High income: OECD"     3777
3 "2. High income: nonOECD"   189
4 "3. Upper middle income"    472
5 "4. Lower middle income"    234
6 "5. Low income"              59

Code

INCOME_GRP <- sapply(as.character(MPA_data$INCOME_GRP), function(x){
  tmp <- strsplit(x,":",fixed = TRUE)[[1]][1]
  strsplit(tmp,". ",fixed = TRUE)[[1]][2]
})

MPA_data$INCOME_GRP <- factor(INCOME_GRP, levels = c("High income","Upper middle income","Lower middle income","Low income"))
table(MPA_data$INCOME_GRP)


        High income Upper middle income Lower middle income          Low income 
               3966                 472                 234                  59

Code

ggplot(MPA_data %>% 
         group_by(INCOME_GRP) %>%
         tally(),
       aes(x = INCOME_GRP, y = n)) +
  geom_bar(stat = "identity") +
  labs(x = "IUCN category of MPA", y = "Count")+
  theme_classic()

Why there are so many MPAs in the high income group?

Code

# make continuous protection category
MPA_data <- MPA_data %>%
  mutate(IUCN_CONT = case_when(
    IUCN_CAT == "Ia" ~ 1,
    IUCN_CAT == "Ib" ~ 1,
    IUCN_CAT == "II" ~ 2,
    IUCN_CAT == "III" ~ 3,
    IUCN_CAT == "IV" ~ 4,
    IUCN_CAT == "VI" ~ 5,
  ))

Code

hist(MPA_data$GIS_M_AREA_raw)

Code

hist(MPA_data$GIS_M_AREA_raw)
hist(log(MPA_data$GIS_M_AREA_raw))

Code

# log area
MPA_data <- MPA_data %>%
  mutate(MPA_area_log = log(GIS_M_AREA_raw))

3 IUCN protected area categories:

IUCN Protected Area Categories and Their Main Objectives
Category	Name	Objective
Ia	Strict Nature Reserve	Strict protection for biodiversity and geological/geomorphological features.
Ib	Wilderness Area	Large unmodified areas with no permanent human habitation; preserved in natural state.
II	National Park	Ecosystem protection and recreation; may allow visitor use and education.
III	Natural Monument or Feature	Protection of specific natural monuments (e.g., landforms, seamounts, geological features).
IV	Habitat/Species Management Area	Active management for conservation through intervention (e.g., habitat restoration).
V	Protected Landscape/Seascape	Protection of areas where interaction of people and nature has produced distinct character.
VI	Protected Area with Sustainable Use of Natural Resources	Conservation with sustainable use of natural resources; may support local livelihoods.

4 Correlation plot

Code

subset_data <- data.frame(MPA_data[,
                                   c("Population", "Migration", "GDP", "Age", "Fishing_hours", "MPA_area", "SMI","DBI","IUCN_CONT")])
str(subset_data)

'data.frame':   4802 obs. of  9 variables:
 $ Population   : num  0.129 0.114 -0.107 -0.229 -0.2 ...
 $ Migration    : num  -0.0849 0.2312 -0.282 0.4405 0.4576 ...
 $ GDP          : num  0.114 0.231 -0.105 0.624 0.624 ...
 $ Age          : num  -0.482 -0.326 -0.492 -0.492 -0.492 ...
 $ Fishing_hours: num  0 0 0 0 0 0 0 0 0 0 ...
 $ MPA_area     : num  0.04504 0.00788 0.00388 0.00102 0.02278 ...
 $ SMI          : num  -0.0723 0.1821 -0.471 0.0683 0.0794 ...
 $ DBI          : num  0 0.314 -0.632 0 0 ...
 $ IUCN_CONT    : num  3 1 3 3 3 3 3 3 3 3 ...

Code

c_df <- cor(subset_data, method = 'spearman', use = "pairwise.complete.obs")

corrplot::corrplot(c_df,
                   method = 'color', 
                   addCoef.col = 'black',
                   insig='blank',
                   diag = FALSE,
                   number.cex = 0.5,
                   type='lower')

This is very similar to the values on the presentation that Florin shared with me!

Based on correlations, SMI and DBI primarily reflect patterns of population and migration, but to a lesser extent also capture variation related to GDP and age.

5 Run models

5.1 Test for autocorrelation among variables

Start of with a model including IUCN category, MPA area and Income group

Code

subdata <- MPA_data |>
  select("DBI", "SMI", "IUCN_CAT", "MPA_area_log", "INCOME_GRP", "ADM0_A3", "REGION_WB") |>
  na.omit()

# DBI
mod_DBI_1 <- glm(DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP, 
                 data = subdata)
car::vif(mod_DBI_1)

                 GVIF Df GVIF^(1/(2*Df))
IUCN_CAT     1.360731  6        1.026001
MPA_area_log 1.308204  1        1.143768
INCOME_GRP   1.202883  3        1.031266

Code

# SMI
mod_SMI_1 <- glm(SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP, 
                 data = subdata)
car::vif(mod_SMI_1)

                 GVIF Df GVIF^(1/(2*Df))
IUCN_CAT     1.360731  6        1.026001
MPA_area_log 1.308204  1        1.143768
INCOME_GRP   1.202883  3        1.031266

VIF values are acceptable (VIF < 3) in both models. So we keep all variables.

5.2 Find the best random structure

Here we test what is the best random structure: country or region. The random intercept on country controls for regional variability (economy, climate and potentially other stuff I have to think more about). In contrast, the random intercept on region controls for effects at a coarser scales, potentially representing socioeconomic related grouping.

5.2.1 DBI

Code

# country
mod_DBI_c <- glmmTMB(DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1|ADM0_A3), 
                     data = subdata)
summary(mod_DBI_c)

 Family: gaussian  ( identity )
Formula:          DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1 | ADM0_A3)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   611.3    689.2   -292.6    585.3     2937 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 ADM0_A3  (Intercept) 0.04736  0.2176  
 Residual             0.06702  0.2589  
Number of obs: 2950, groups:  ADM0_A3, 100

Dispersion estimate for gaussian family (sigma^2): 0.067 

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)   
(Intercept)                   -0.078211   0.041291  -1.894  0.05821 . 
IUCN_CATIb                     0.009268   0.038063   0.244  0.80762   
IUCN_CATII                     0.045144   0.026658   1.694  0.09037 . 
IUCN_CATIII                    0.039867   0.027011   1.476  0.13997   
IUCN_CATIV                     0.048163   0.017315   2.782  0.00541 **
IUCN_CATV                     -0.017335   0.020746  -0.836  0.40340   
IUCN_CATVI                     0.028867   0.022028   1.310  0.19003   
MPA_area_log                  -0.001754   0.001465  -1.197  0.23133   
INCOME_GRPUpper middle income  0.075380   0.057440   1.312  0.18941   
INCOME_GRPLower middle income  0.064236   0.073190   0.878  0.38013   
INCOME_GRPLow income           0.259120   0.121437   2.134  0.03286 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# region
mod_DBI_r <- glmmTMB(DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1|REGION_WB), 
                     data = subdata)
summary(mod_DBI_r)

 Family: gaussian  ( identity )
Formula:          DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1 | REGION_WB)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  1363.5   1441.4   -668.8   1337.5     2937 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 REGION_WB (Intercept) 0.02317  0.1522  
 Residual              0.09133  0.3022  
Number of obs: 2950, groups:  REGION_WB, 7

Dispersion estimate for gaussian family (sigma^2): 0.0913 

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    0.006146   0.062485   0.098 0.921651    
IUCN_CATIb                    -0.108699   0.040794  -2.665 0.007708 ** 
IUCN_CATII                     0.043102   0.027458   1.570 0.116473    
IUCN_CATIII                    0.003695   0.026751   0.138 0.890132    
IUCN_CATIV                    -0.050782   0.016206  -3.134 0.001727 ** 
IUCN_CATV                     -0.072326   0.021200  -3.412 0.000646 ***
IUCN_CATVI                    -0.037369   0.022700  -1.646 0.099716 .  
MPA_area_log                  -0.002018   0.001509  -1.338 0.181026    
INCOME_GRPUpper middle income  0.034743   0.023585   1.473 0.140724    
INCOME_GRPLower middle income  0.242807   0.028405   8.548  < 2e-16 ***
INCOME_GRPLow income           0.332976   0.096031   3.467 0.000526 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# compare
AIC(mod_DBI_c, mod_DBI_r)

          df      AIC
mod_DBI_c 13  611.287
mod_DBI_r 13 1363.537

5.2.2 SMI

Code

# country
mod_SMI_c <- glmmTMB(SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1|ADM0_A3), 
                     data = subdata)
summary(mod_SMI_c)

 Family: gaussian  ( identity )
Formula:          SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1 | ADM0_A3)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  -391.4   -313.5    208.7   -417.4     2937 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 ADM0_A3  (Intercept) 0.03948  0.1987  
 Residual             0.04751  0.2180  
Number of obs: 2950, groups:  ADM0_A3, 100

Dispersion estimate for gaussian family (sigma^2): 0.0475 

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -0.092710   0.037026  -2.504 0.012282 *  
IUCN_CATIb                     0.026205   0.032096   0.816 0.414241    
IUCN_CATII                     0.042046   0.022486   1.870 0.061499 .  
IUCN_CATIII                    0.037579   0.022770   1.650 0.098863 .  
IUCN_CATIV                     0.051836   0.014594   3.552 0.000382 ***
IUCN_CATV                      0.008025   0.017487   0.459 0.646306    
IUCN_CATVI                     0.030582   0.018569   1.647 0.099567 .  
MPA_area_log                  -0.002865   0.001236  -2.318 0.020453 *  
INCOME_GRPUpper middle income  0.044639   0.051676   0.864 0.387693    
INCOME_GRPLower middle income  0.031153   0.065553   0.475 0.634627    
INCOME_GRPLow income           0.248028   0.107955   2.298 0.021590 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

MuMIn::r.squaredGLMM(mod_SMI_c)

            R2m       R2c
[1,] 0.01121088 0.4599492

Code

# region
mod_SMI_r <- glmmTMB(SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1|REGION_WB), 
                     data = subdata)
summary(mod_SMI_r)

 Family: gaussian  ( identity )
Formula:          SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + (1 | REGION_WB)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   523.5    601.3   -248.7    497.5     2937 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 REGION_WB (Intercept) 0.01761  0.1327  
 Residual              0.06870  0.2621  
Number of obs: 2950, groups:  REGION_WB, 7

Dispersion estimate for gaussian family (sigma^2): 0.0687 

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -0.021553   0.054416  -0.396   0.6920    
IUCN_CATIb                    -0.070782   0.035380  -2.001   0.0454 *  
IUCN_CATII                     0.046088   0.023815   1.935   0.0530 .  
IUCN_CATIII                   -0.006727   0.023200  -0.290   0.7718    
IUCN_CATIV                    -0.031106   0.014055  -2.213   0.0269 *  
IUCN_CATV                     -0.025835   0.018386  -1.405   0.1600    
IUCN_CATVI                    -0.026295   0.019688  -1.336   0.1817    
MPA_area_log                  -0.003032   0.001308  -2.317   0.0205 *  
INCOME_GRPUpper middle income  0.019194   0.020467   0.938   0.3483    
INCOME_GRPLower middle income  0.226968   0.024637   9.212   <2e-16 ***
INCOME_GRPLow income           0.335104   0.083512   4.013    6e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

MuMIn::r.squaredGLMM(mod_SMI_r)

            R2m       R2c
[1,] 0.03650178 0.2330538

Code

# compare
AIC(mod_SMI_c, mod_SMI_r)

          df       AIC
mod_SMI_c 13 -391.3503
mod_SMI_r 13  523.4757

In terms of AIC and BIC, for both for DBI and SMI, the models including country in the random intercept are better than the models including region.

Note that: 1) Country captures a lot of variance as provided by conditional R2. 2) However, R2 conditional is higher for models including region as the random intercept. Therefore, even though models including country in as range intercept are better in terms of AIC, the variance explained the fixed effects is greater in the models including region in as random intercept.

6 Visualize effects

6.1 Country models

6.1.1 DBI

Code

# IUCN cat
emmeans_DBI <- emmeans(mod_DBI_c,
                       type = "response",
                       specs = "IUCN_CAT") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Income
emmeans_DBI <- emmeans(mod_DBI_c,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Area
preds <- ggeffects::ggeffect(
  mod_DBI_c, 
  terms = "MPA_area_log")

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- summary(mod_DBI_c)
sig_tmp <- sig_tmp$coefficients$cond
sig_tmp <- sig_tmp["MPA_area_log","Pr(>|z|)"]
preds$sig <- ifelse(sig_tmp <= 0.05,"significant","non-significant")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "DBI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()

No significant effects of DBI using country as random effect…

6.1.2 SMI

Code

# IUCN cat
emmeans_SMI <- emmeans(mod_SMI_c,
                       type = "response",
                       specs = "IUCN_CAT") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Income
emmeans_SMI <- emmeans(mod_SMI_c,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Area
preds <- ggeffects::ggeffect(
  mod_SMI_c, 
  terms = "MPA_area_log")

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- summary(mod_SMI_c)
sig_tmp <- sig_tmp$coefficients$cond
sig_tmp <- sig_tmp["MPA_area_log","Pr(>|z|)"]
preds$sig <- ifelse(sig_tmp <= 0.05,"significant","non-significant")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "SMI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()

No significant effects of SMI using country as random effect…

6.2 Region models

6.2.1 DBI

Code

# IUCN cat
emmeans_DBI <- emmeans(mod_DBI_r,
                       type = "response",
                       specs = "IUCN_CAT") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Income
emmeans_DBI <- emmeans(mod_DBI_r,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Area
preds <- ggeffects::ggeffect(
  mod_DBI_r, 
  terms = "MPA_area_log")

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- summary(mod_DBI_r)
sig_tmp <- sig_tmp$coefficients$cond
sig_tmp <- sig_tmp["MPA_area_log","Pr(>|z|)"]
preds$sig <- ifelse(sig_tmp <= 0.05,"significant","non-significant")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "DBI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()

Significant effects using region as random intercept…

Protected areas in categories more strict protection models (Ia, II and III) show positive DBIs.

Protected areas in middle to low and low income countries show positive DBIs.

6.2.2 SMI

Code

# IUCN cat
emmeans_SMI <- emmeans(mod_SMI_r,
                       type = "response",
                       specs = "IUCN_CAT") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Income
emmeans_SMI <- emmeans(mod_SMI_r,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

# Area
preds <- ggeffects::ggeffect(
  mod_SMI_r, 
  terms = "MPA_area_log")

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- summary(mod_SMI_r)
sig_tmp <- sig_tmp$coefficients$cond
sig_tmp <- sig_tmp["MPA_area_log","Pr(>|z|)"]
preds$sig <- ifelse(sig_tmp <= 0.05,"significant","non-significant")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "SMI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()

Same results when using SMI… Protected areas in categories more strict protection models (Ia, II and III) show positive SMIs.

Protected areas in middle-low and low income countries show positive SMIs.

7 Interpretation:

The models with country never provide significant relationships for DBI or SMI, even though they perform better in terms of AIC. However, the models with region as random intercept do provide some significant results. This can be because there is too much noise at the country level that the models cannot make sense of (variability in DBI/SMI values is not structured at the country level). However, this variability is better represented by regions. A model can have a better AIC if it captures more of the total variance, even if individual predictors are not statistically significant, especially if the added random effects increase uncertainty around fixed effect estimates.

I discussed the results with Mona, and this is what we think:

1) MPAs in more strictly protected categories tend to show more positive DBI and SMI values. Stricter MPAs often reflect more political or institutional commitment, which may be linked to areas with stronger planning capacity or more involvement from national or international organizations. These places may be better at attracting people, investments, or offering economic opportunities. That could explain why we see higher population growth, migration, and GDP trends (positive DBI and SMI) in these areas. In other words, stricter MPAs might be located in more dynamic or developing places, at least socio-economically.

2) MPAs located in lower and lower-middle income countries also show more positive DBI and SMI.
This might sound surprising at first, but many MPAs in lower-income countries are located in fast-changing regions, where population is growing quickly, migration is high, and economies are expanding. Even if the country is classified as lower-income overall, these coastal areas might be socio-economic hotspots. So DBI and SMI could be picking up high local growth rates, especially in terms of demographics (young and growing populations) and urban or economic expansion around MPAs.

3) Large MPAs show more negative DBI and SMI.
It is possible that large MPAs are often established in remote oceanic regions, far from dense population centers. These areas may have demographic or economic decline or imbalance relative to small MPAs.

8 Interaction model

The results are interesting so far. We found MPAs in low income countries show more positive DBI/SMI, but they do so across all protection categories? Also, MPAs in low income countries perform well across all sizes? We can test that adding interaction terms to the model.

8.1 DBI

Code

# country
mod_DBI_c_i <- glmmTMB(DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + 
                         IUCN_CAT : INCOME_GRP +
                         MPA_area_log : INCOME_GRP + (1|ADM0_A3), 
                       data = subdata)

dropping columns from rank-deficient conditional model: IUCN_CATIb:INCOME_GRPLow income, IUCN_CATIII:INCOME_GRPLow income, IUCN_CATVI:INCOME_GRPLow income

Code

summary(mod_DBI_c_i)

 Family: gaussian  ( identity )
Formula:          
DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + IUCN_CAT:INCOME_GRP +  
    MPA_area_log:INCOME_GRP + (1 | ADM0_A3)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   625.7    811.4   -281.8    563.7     2919 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 ADM0_A3  (Intercept) 0.04621  0.215   
 Residual             0.06656  0.258   
Number of obs: 2950, groups:  ADM0_A3, 100

Dispersion estimate for gaussian family (sigma^2): 0.0666 

Conditional model:
                                            Estimate Std. Error z value
(Intercept)                                -0.082182   0.041083  -2.000
IUCN_CATIb                                 -0.009069   0.039739  -0.228
IUCN_CATII                                  0.075520   0.031575   2.392
IUCN_CATIII                                 0.058448   0.029033   2.013
IUCN_CATIV                                  0.053500   0.018068   2.961
IUCN_CATV                                  -0.018134   0.021796  -0.832
IUCN_CATVI                                  0.026065   0.024511   1.063
MPA_area_log                               -0.002327   0.001564  -1.489
INCOME_GRPUpper middle income               0.131150   0.083716   1.567
INCOME_GRPLower middle income               0.072565   0.119754   0.606
INCOME_GRPLow income                        0.104599   0.266682   0.392
IUCN_CATIb:INCOME_GRPUpper middle income    0.263981   0.171972   1.535
IUCN_CATII:INCOME_GRPUpper middle income   -0.154017   0.079923  -1.927
IUCN_CATIII:INCOME_GRPUpper middle income  -0.181180   0.096074  -1.886
IUCN_CATIV:INCOME_GRPUpper middle income   -0.061791   0.072636  -0.851
IUCN_CATV:INCOME_GRPUpper middle income    -0.048877   0.084278  -0.580
IUCN_CATVI:INCOME_GRPUpper middle income   -0.061487   0.075509  -0.814
IUCN_CATIb:INCOME_GRPLower middle income    0.132952   0.212091   0.627
IUCN_CATII:INCOME_GRPLower middle income   -0.046468   0.133812  -0.347
IUCN_CATIII:INCOME_GRPLower middle income  -0.019915   0.142447  -0.140
IUCN_CATIV:INCOME_GRPLower middle income   -0.078187   0.115606  -0.676
IUCN_CATV:INCOME_GRPLower middle income    -0.003182   0.133660  -0.024
IUCN_CATVI:INCOME_GRPLower middle income    0.073199   0.110653   0.661
IUCN_CATIb:INCOME_GRPLow income                   NA         NA      NA
IUCN_CATII:INCOME_GRPLow income             0.168828   0.164063   1.029
IUCN_CATIII:INCOME_GRPLow income                  NA         NA      NA
IUCN_CATIV:INCOME_GRPLow income             0.278923   0.335550   0.831
IUCN_CATV:INCOME_GRPLow income              0.375985   0.232553   1.617
IUCN_CATVI:INCOME_GRPLow income                   NA         NA      NA
MPA_area_log:INCOME_GRPUpper middle income  0.004791   0.005466   0.877
MPA_area_log:INCOME_GRPLower middle income  0.001004   0.008061   0.125
MPA_area_log:INCOME_GRPLow income           0.002007   0.041558   0.048
                                           Pr(>|z|)   
(Intercept)                                 0.04546 * 
IUCN_CATIb                                  0.81949   
IUCN_CATII                                  0.01677 * 
IUCN_CATIII                                 0.04410 * 
IUCN_CATIV                                  0.00307 **
IUCN_CATV                                   0.40541   
IUCN_CATVI                                  0.28760   
MPA_area_log                                0.13660   
INCOME_GRPUpper middle income               0.11721   
INCOME_GRPLower middle income               0.54455   
INCOME_GRPLow income                        0.69489   
IUCN_CATIb:INCOME_GRPUpper middle income    0.12478   
IUCN_CATII:INCOME_GRPUpper middle income    0.05397 . 
IUCN_CATIII:INCOME_GRPUpper middle income   0.05932 . 
IUCN_CATIV:INCOME_GRPUpper middle income    0.39494   
IUCN_CATV:INCOME_GRPUpper middle income     0.56195   
IUCN_CATVI:INCOME_GRPUpper middle income    0.41547   
IUCN_CATIb:INCOME_GRPLower middle income    0.53075   
IUCN_CATII:INCOME_GRPLower middle income    0.72839   
IUCN_CATIII:INCOME_GRPLower middle income   0.88881   
IUCN_CATIV:INCOME_GRPLower middle income    0.49883   
IUCN_CATV:INCOME_GRPLower middle income     0.98100   
IUCN_CATVI:INCOME_GRPLower middle income    0.50828   
IUCN_CATIb:INCOME_GRPLow income                  NA   
IUCN_CATII:INCOME_GRPLow income             0.30346   
IUCN_CATIII:INCOME_GRPLow income                 NA   
IUCN_CATIV:INCOME_GRPLow income             0.40584   
IUCN_CATV:INCOME_GRPLow income              0.10593   
IUCN_CATVI:INCOME_GRPLow income                  NA   
MPA_area_log:INCOME_GRPUpper middle income  0.38069   
MPA_area_log:INCOME_GRPLower middle income  0.90086   
MPA_area_log:INCOME_GRPLow income           0.96148   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# region
mod_DBI_r_i <- glmmTMB(DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + 
                         IUCN_CAT : INCOME_GRP +
                         MPA_area_log : INCOME_GRP + (1|REGION_WB), 
                       data = subdata)

dropping columns from rank-deficient conditional model: IUCN_CATIb:INCOME_GRPLow income, IUCN_CATIII:INCOME_GRPLow income, IUCN_CATVI:INCOME_GRPLow income

Code

summary(mod_DBI_r_i)

 Family: gaussian  ( identity )
Formula:          
DBI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + IUCN_CAT:INCOME_GRP +  
    MPA_area_log:INCOME_GRP + (1 | REGION_WB)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  1348.2   1533.9   -643.1   1286.2     2919 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 REGION_WB (Intercept) 0.02081  0.1443  
 Residual              0.08978  0.2996  
Number of obs: 2950, groups:  REGION_WB, 7

Dispersion estimate for gaussian family (sigma^2): 0.0898 

Conditional model:
                                             Estimate Std. Error z value
(Intercept)                                -0.0002873  0.0598752  -0.005
IUCN_CATIb                                 -0.1202639  0.0427973  -2.810
IUCN_CATII                                  0.0983339  0.0339609   2.896
IUCN_CATIII                                 0.0329437  0.0288124   1.143
IUCN_CATIV                                 -0.0466687  0.0166843  -2.797
IUCN_CATV                                  -0.0726841  0.0222881  -3.261
IUCN_CATVI                                 -0.0660255  0.0250589  -2.635
MPA_area_log                               -0.0033272  0.0016344  -2.036
INCOME_GRPUpper middle income               0.0430461  0.0686671   0.627
INCOME_GRPLower middle income               0.1738389  0.1082448   1.606
INCOME_GRPLow income                        0.1894255  0.2130960   0.889
IUCN_CATIb:INCOME_GRPUpper middle income    0.2946576  0.1714974   1.718
IUCN_CATII:INCOME_GRPUpper middle income   -0.1414139  0.0799385  -1.769
IUCN_CATIII:INCOME_GRPUpper middle income  -0.2133204  0.0955891  -2.232
IUCN_CATIV:INCOME_GRPUpper middle income   -0.0379870  0.0737030  -0.515
IUCN_CATV:INCOME_GRPUpper middle income    -0.0134026  0.0831148  -0.161
IUCN_CATVI:INCOME_GRPUpper middle income    0.1167096  0.0806892   1.446
IUCN_CATIb:INCOME_GRPLower middle income   -0.0616567  0.2111392  -0.292
IUCN_CATII:INCOME_GRPLower middle income   -0.0858502  0.1461623  -0.587
IUCN_CATIII:INCOME_GRPLower middle income  -0.0319881  0.1540440  -0.208
IUCN_CATIV:INCOME_GRPLower middle income    0.0405864  0.1184397   0.343
IUCN_CATV:INCOME_GRPLower middle income     0.0154672  0.1325288   0.117
IUCN_CATVI:INCOME_GRPLower middle income    0.1693519  0.1211960   1.397
IUCN_CATIb:INCOME_GRPLow income                    NA         NA      NA
IUCN_CATII:INCOME_GRPLow income             0.0509339  0.1769980   0.288
IUCN_CATIII:INCOME_GRPLow income                   NA         NA      NA
IUCN_CATIV:INCOME_GRPLow income             0.2402776  0.3347080   0.718
IUCN_CATV:INCOME_GRPLow income              0.3526911  0.2308645   1.528
IUCN_CATVI:INCOME_GRPLow income                    NA         NA      NA
MPA_area_log:INCOME_GRPUpper middle income  0.0133263  0.0049818   2.675
MPA_area_log:INCOME_GRPLower middle income  0.0043720  0.0076394   0.572
MPA_area_log:INCOME_GRPLow income           0.0059557  0.0379492   0.157
                                           Pr(>|z|)   
(Intercept)                                 0.99617   
IUCN_CATIb                                  0.00495 **
IUCN_CATII                                  0.00379 **
IUCN_CATIII                                 0.25288   
IUCN_CATIV                                  0.00516 **
IUCN_CATV                                   0.00111 **
IUCN_CATVI                                  0.00842 **
MPA_area_log                                0.04178 * 
INCOME_GRPUpper middle income               0.53074   
INCOME_GRPLower middle income               0.10828   
INCOME_GRPLow income                        0.37405   
IUCN_CATIb:INCOME_GRPUpper middle income    0.08577 . 
IUCN_CATII:INCOME_GRPUpper middle income    0.07689 . 
IUCN_CATIII:INCOME_GRPUpper middle income   0.02564 * 
IUCN_CATIV:INCOME_GRPUpper middle income    0.60627   
IUCN_CATV:INCOME_GRPUpper middle income     0.87189   
IUCN_CATVI:INCOME_GRPUpper middle income    0.14806   
IUCN_CATIb:INCOME_GRPLower middle income    0.77027   
IUCN_CATII:INCOME_GRPLower middle income    0.55696   
IUCN_CATIII:INCOME_GRPLower middle income   0.83550   
IUCN_CATIV:INCOME_GRPLower middle income    0.73184   
IUCN_CATV:INCOME_GRPLower middle income     0.90709   
IUCN_CATVI:INCOME_GRPLower middle income    0.16231   
IUCN_CATIb:INCOME_GRPLow income                  NA   
IUCN_CATII:INCOME_GRPLow income             0.77353   
IUCN_CATIII:INCOME_GRPLow income                 NA   
IUCN_CATIV:INCOME_GRPLow income             0.47284   
IUCN_CATV:INCOME_GRPLow income              0.12659   
IUCN_CATVI:INCOME_GRPLow income                  NA   
MPA_area_log:INCOME_GRPUpper middle income  0.00747 **
MPA_area_log:INCOME_GRPLower middle income  0.56712   
MPA_area_log:INCOME_GRPLow income           0.87529   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

AIC(mod_DBI_c,
    mod_DBI_c_i,
    mod_DBI_r,
    mod_DBI_r_i)

            df      AIC
mod_DBI_c   13  611.287
mod_DBI_c_i 31  625.700
mod_DBI_r   13 1363.537
mod_DBI_r_i 31 1348.238

In terms of AIC, including interactions in the models doesn’t improve the model fit as much relative to the ones without interactions. However, note that the regional model show a slight improvement.

8.2 SMI

Code

# country
mod_SMI_c_i <- glmmTMB(SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + 
                         IUCN_CAT : INCOME_GRP +
                         MPA_area_log : INCOME_GRP + (1|ADM0_A3), 
                       data = subdata)

dropping columns from rank-deficient conditional model: IUCN_CATIb:INCOME_GRPLow income, IUCN_CATIII:INCOME_GRPLow income, IUCN_CATVI:INCOME_GRPLow income

Code

summary(mod_SMI_c_i)

 Family: gaussian  ( identity )
Formula:          
SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + IUCN_CAT:INCOME_GRP +  
    MPA_area_log:INCOME_GRP + (1 | ADM0_A3)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  -381.5   -195.8    221.7   -443.5     2919 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 ADM0_A3  (Intercept) 0.03852  0.1963  
 Residual             0.04711  0.2170  
Number of obs: 2950, groups:  ADM0_A3, 100

Dispersion estimate for gaussian family (sigma^2): 0.0471 

Conditional model:
                                            Estimate Std. Error z value
(Intercept)                                -0.097023   0.036805  -2.636
IUCN_CATIb                                  0.010787   0.033452   0.322
IUCN_CATII                                  0.069930   0.026587   2.630
IUCN_CATIII                                 0.053136   0.024436   2.174
IUCN_CATIV                                  0.056725   0.015211   3.729
IUCN_CATV                                   0.011759   0.018344   0.641
IUCN_CATVI                                  0.025032   0.020635   1.213
MPA_area_log                               -0.002952   0.001316  -2.243
INCOME_GRPUpper middle income               0.111139   0.072884   1.525
INCOME_GRPLower middle income               0.051884   0.103300   0.502
INCOME_GRPLow income                        0.136853   0.233435   0.586
IUCN_CATIb:INCOME_GRPUpper middle income    0.187430   0.146655   1.278
IUCN_CATII:INCOME_GRPUpper middle income   -0.142903   0.067512  -2.117
IUCN_CATIII:INCOME_GRPUpper middle income  -0.171559   0.081239  -2.112
IUCN_CATIV:INCOME_GRPUpper middle income   -0.059142   0.061407  -0.963
IUCN_CATV:INCOME_GRPUpper middle income    -0.080967   0.071326  -1.135
IUCN_CATVI:INCOME_GRPUpper middle income   -0.064129   0.063704  -1.007
IUCN_CATIb:INCOME_GRPLower middle income    0.161860   0.179528   0.902
IUCN_CATII:INCOME_GRPLower middle income   -0.070650   0.113030  -0.625
IUCN_CATIII:INCOME_GRPLower middle income  -0.009574   0.120512  -0.079
IUCN_CATIV:INCOME_GRPLower middle income   -0.080710   0.097853  -0.825
IUCN_CATV:INCOME_GRPLower middle income    -0.028320   0.113145  -0.250
IUCN_CATVI:INCOME_GRPLower middle income    0.088638   0.093442   0.949
IUCN_CATIb:INCOME_GRPLow income                   NA         NA      NA
IUCN_CATII:INCOME_GRPLow income             0.151378   0.139383   1.086
IUCN_CATIII:INCOME_GRPLow income                  NA         NA      NA
IUCN_CATIV:INCOME_GRPLow income             0.088465   0.286368   0.309
IUCN_CATV:INCOME_GRPLow income              0.337278   0.197473   1.708
IUCN_CATVI:INCOME_GRPLow income                   NA         NA      NA
MPA_area_log:INCOME_GRPUpper middle income  0.002142   0.004641   0.462
MPA_area_log:INCOME_GRPLower middle income -0.002924   0.006842  -0.427
MPA_area_log:INCOME_GRPLow income          -0.001030   0.035876  -0.029
                                           Pr(>|z|)    
(Intercept)                                0.008387 ** 
IUCN_CATIb                                 0.747103    
IUCN_CATII                                 0.008532 ** 
IUCN_CATIII                                0.029669 *  
IUCN_CATIV                                 0.000192 ***
IUCN_CATV                                  0.521502    
IUCN_CATVI                                 0.225100    
MPA_area_log                               0.024916 *  
INCOME_GRPUpper middle income              0.127290    
INCOME_GRPLower middle income              0.615481    
INCOME_GRPLow income                       0.557703    
IUCN_CATIb:INCOME_GRPUpper middle income   0.201237    
IUCN_CATII:INCOME_GRPUpper middle income   0.034284 *  
IUCN_CATIII:INCOME_GRPUpper middle income  0.034705 *  
IUCN_CATIV:INCOME_GRPUpper middle income   0.335486    
IUCN_CATV:INCOME_GRPUpper middle income    0.256307    
IUCN_CATVI:INCOME_GRPUpper middle income   0.314096    
IUCN_CATIb:INCOME_GRPLower middle income   0.367276    
IUCN_CATII:INCOME_GRPLower middle income   0.531934    
IUCN_CATIII:INCOME_GRPLower middle income  0.936679    
IUCN_CATIV:INCOME_GRPLower middle income   0.409477    
IUCN_CATV:INCOME_GRPLower middle income    0.802358    
IUCN_CATVI:INCOME_GRPLower middle income   0.342830    
IUCN_CATIb:INCOME_GRPLow income                  NA    
IUCN_CATII:INCOME_GRPLow income            0.277451    
IUCN_CATIII:INCOME_GRPLow income                 NA    
IUCN_CATIV:INCOME_GRPLow income            0.757382    
IUCN_CATV:INCOME_GRPLow income             0.087642 .  
IUCN_CATVI:INCOME_GRPLow income                  NA    
MPA_area_log:INCOME_GRPUpper middle income 0.644417    
MPA_area_log:INCOME_GRPLower middle income 0.669160    
MPA_area_log:INCOME_GRPLow income          0.977094    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# region
mod_SMI_r_i <- glmmTMB(SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + 
                         IUCN_CAT : INCOME_GRP +
                         MPA_area_log : INCOME_GRP + (1|REGION_WB), 
                       data = subdata)

dropping columns from rank-deficient conditional model: IUCN_CATIb:INCOME_GRPLow income, IUCN_CATIII:INCOME_GRPLow income, IUCN_CATVI:INCOME_GRPLow income

Code

summary(mod_SMI_r_i)

 Family: gaussian  ( identity )
Formula:          
SMI ~ IUCN_CAT + MPA_area_log + INCOME_GRP + IUCN_CAT:INCOME_GRP +  
    MPA_area_log:INCOME_GRP + (1 | REGION_WB)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   502.4    688.1   -220.2    440.4     2919 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 REGION_WB (Intercept) 0.01549  0.1244  
 Residual              0.06740  0.2596  
Number of obs: 2950, groups:  REGION_WB, 7

Dispersion estimate for gaussian family (sigma^2): 0.0674 

Conditional model:
                                             Estimate Std. Error z value
(Intercept)                                -0.0290754  0.0516658  -0.563
IUCN_CATIb                                 -0.0832589  0.0370817  -2.245
IUCN_CATII                                  0.0942689  0.0294260   3.204
IUCN_CATIII                                 0.0188204  0.0249652   0.754
IUCN_CATIV                                 -0.0248210  0.0144559  -1.717
IUCN_CATV                                  -0.0224417  0.0193111  -1.162
IUCN_CATVI                                 -0.0590832  0.0217128  -2.721
MPA_area_log                               -0.0036223  0.0014161  -2.558
INCOME_GRPUpper middle income               0.0559720  0.0595165   0.940
INCOME_GRPLower middle income               0.1856455  0.0937998   1.979
INCOME_GRPLow income                        0.2607575  0.1847848   1.411
IUCN_CATIb:INCOME_GRPUpper middle income    0.2576340  0.1485665   1.734
IUCN_CATII:INCOME_GRPUpper middle income   -0.1407196  0.0692604  -2.032
IUCN_CATIII:INCOME_GRPUpper middle income  -0.2141296  0.0828381  -2.585
IUCN_CATIV:INCOME_GRPUpper middle income   -0.0692035  0.0638614  -1.084
IUCN_CATV:INCOME_GRPUpper middle income    -0.0409085  0.0720146  -0.568
IUCN_CATVI:INCOME_GRPUpper middle income    0.0911260  0.0699095   1.303
IUCN_CATIb:INCOME_GRPLower middle income   -0.0370602  0.1829408  -0.203
IUCN_CATII:INCOME_GRPLower middle income   -0.1216633  0.1266221  -0.961
IUCN_CATIII:INCOME_GRPLower middle income  -0.0091558  0.1334788  -0.069
IUCN_CATIV:INCOME_GRPLower middle income   -0.0085441  0.1026039  -0.083
IUCN_CATV:INCOME_GRPLower middle income    -0.0305323  0.1148295  -0.266
IUCN_CATVI:INCOME_GRPLower middle income    0.1793807  0.1050126   1.708
IUCN_CATIb:INCOME_GRPLow income                    NA         NA      NA
IUCN_CATII:INCOME_GRPLow income             0.0459240  0.1533638   0.299
IUCN_CATIII:INCOME_GRPLow income                   NA         NA      NA
IUCN_CATIV:INCOME_GRPLow income             0.0655205  0.2900108   0.226
IUCN_CATV:INCOME_GRPLow income              0.2923805  0.2000361   1.462
IUCN_CATVI:INCOME_GRPLow income                    NA         NA      NA
MPA_area_log:INCOME_GRPUpper middle income  0.0082916  0.0043194   1.920
MPA_area_log:INCOME_GRPLower middle income  0.0013268  0.0066177   0.200
MPA_area_log:INCOME_GRPLow income          -0.0005159  0.0328810  -0.016
                                           Pr(>|z|)   
(Intercept)                                 0.57360   
IUCN_CATIb                                  0.02475 * 
IUCN_CATII                                  0.00136 **
IUCN_CATIII                                 0.45093   
IUCN_CATIV                                  0.08598 . 
IUCN_CATV                                   0.24519   
IUCN_CATVI                                  0.00651 **
MPA_area_log                                0.01053 * 
INCOME_GRPUpper middle income               0.34699   
INCOME_GRPLower middle income               0.04780 * 
INCOME_GRPLow income                        0.15820   
IUCN_CATIb:INCOME_GRPUpper middle income    0.08289 . 
IUCN_CATII:INCOME_GRPUpper middle income    0.04218 * 
IUCN_CATIII:INCOME_GRPUpper middle income   0.00974 **
IUCN_CATIV:INCOME_GRPUpper middle income    0.27852   
IUCN_CATV:INCOME_GRPUpper middle income     0.57000   
IUCN_CATVI:INCOME_GRPUpper middle income    0.19241   
IUCN_CATIb:INCOME_GRPLower middle income    0.83946   
IUCN_CATII:INCOME_GRPLower middle income    0.33663   
IUCN_CATIII:INCOME_GRPLower middle income   0.94531   
IUCN_CATIV:INCOME_GRPLower middle income    0.93363   
IUCN_CATV:INCOME_GRPLower middle income     0.79032   
IUCN_CATVI:INCOME_GRPLower middle income    0.08760 . 
IUCN_CATIb:INCOME_GRPLow income                  NA   
IUCN_CATII:INCOME_GRPLow income             0.76460   
IUCN_CATIII:INCOME_GRPLow income                 NA   
IUCN_CATIV:INCOME_GRPLow income             0.82126   
IUCN_CATV:INCOME_GRPLow income              0.14384   
IUCN_CATVI:INCOME_GRPLow income                  NA   
MPA_area_log:INCOME_GRPUpper middle income  0.05491 . 
MPA_area_log:INCOME_GRPLower middle income  0.84109   
MPA_area_log:INCOME_GRPLow income           0.98748   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

AIC(mod_SMI_c,
    mod_SMI_c_i,
    mod_SMI_r,
    mod_SMI_r_i)

            df       AIC
mod_SMI_c   13 -391.3503
mod_SMI_c_i 31 -381.4574
mod_SMI_r   13  523.4757
mod_SMI_r_i 31  502.3934

Same here… In terms of AIC, including interactions in the models doesn’t improve the model fit as much relative to the ones without interactions. However, note that the regional model show a slight improvement.

9 Visualize effects

9.1 Country models

9.1.1 DBI

Code

# IUCN cat x INCOME_GRP
emmeans_DBI <- emmeans(mod_DBI_c_i,
                       type = "response",
                       specs = c("IUCN_CAT","INCOME_GRP")) %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~INCOME_GRP)

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Code

# Income vs Area
preds <- ggeffects::ggeffect(
  mod_DBI_c, 
  terms = c("MPA_area_log","INCOME_GRP"))
preds <- data.frame(preds)

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- data.frame(emtrends(object = mod_DBI_c_i, 
                               specs = "INCOME_GRP",
                               var = "MPA_area_log"))

sig_tmp$sig <- ifelse(sign(sig_tmp$lower.CL)==sign(sig_tmp$upper.CL),
                      "significant","non-significant")

preds <- merge(preds, sig_tmp, by.x = "group", by.y = "INCOME_GRP")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "DBI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~group, nrow = 1)

9.1.2 SMI

Code

# IUCN cat x INCOME_GRP
emmeans_SMI <- emmeans(mod_SMI_c_i,
                       type = "response",
                       specs = c("IUCN_CAT","INCOME_GRP")) %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~INCOME_GRP)

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Code

# Income vs Area
preds <- ggeffects::ggeffect(
  mod_DBI_c, 
  terms = c("MPA_area_log","INCOME_GRP"))
preds <- data.frame(preds)

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- data.frame(emtrends(object = mod_SMI_c_i, 
                               specs = "INCOME_GRP",
                               var = "MPA_area_log"))

sig_tmp$sig <- ifelse(sign(sig_tmp$lower.CL)==sign(sig_tmp$upper.CL),
                      "significant","non-significant")

preds <- merge(preds, sig_tmp, by.x = "group", by.y = "INCOME_GRP")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "SMI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~group, nrow = 1)

No significant effects of SMI using country as random effect…

9.2 Region models

9.2.1 DBI

Code

# IUCN cat x INCOME_GRP
emmeans_DBI <- emmeans(mod_DBI_r_i,
                       type = "response",
                       specs = c("IUCN_CAT","INCOME_GRP")) %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~INCOME_GRP)

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Code

# Income vs Area
preds <- ggeffects::ggeffect(
  mod_DBI_c, 
  terms = c("MPA_area_log","INCOME_GRP"))
preds <- data.frame(preds)

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- data.frame(emtrends(object = mod_DBI_r_i, 
                               specs = "INCOME_GRP",
                               var = "MPA_area_log"))

sig_tmp$sig <- ifelse(sign(sig_tmp$lower.CL)==sign(sig_tmp$upper.CL),
                      "significant","non-significant")

preds <- merge(preds, sig_tmp, by.x = "group", by.y = "INCOME_GRP")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "DBI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~group, nrow = 1)

9.2.2 SMI

Code

# IUCN cat x INCOME_GRP
emmeans_SMI <- emmeans(mod_SMI_r_i,
                       type = "response",
                       specs = c("IUCN_CAT","INCOME_GRP")) %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = IUCN_CAT, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~INCOME_GRP)

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Code

# Income vs Area
preds <- ggeffects::ggeffect(
  mod_DBI_c, 
  terms = c("MPA_area_log","INCOME_GRP"))
preds <- data.frame(preds)

# reverse the log
# preds$x <- exp(preds$x)
sig_tmp <- data.frame(emtrends(object = mod_SMI_r_i, 
                               specs = "INCOME_GRP",
                               var = "MPA_area_log"))

sig_tmp$sig <- ifelse(sign(sig_tmp$lower.CL)==sign(sig_tmp$upper.CL),
                      "significant","non-significant")

preds <- merge(preds, sig_tmp, by.x = "group", by.y = "INCOME_GRP")

ggplot(preds, aes(x = x, y = predicted))+ 
  geom_line(aes(linetype = as.factor(sig)), linewidth = 1) +
  scale_linetype_manual(values = c("non-significant"="dashed",
                                   "significant"="solid")) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = .3) +
  labs(x = "MPA area (log)", 
       y  = "SMI",
       linetype = "Test")+
  geom_hline(yintercept = 0)+
  theme_classic()+
  facet_wrap(.~group, nrow = 1)

My interpretation is that the interactions we test doesn’t help explain the patterns on income group.

10 Supplementary

Now run similar models using a simplified IUCN classification scheme by grouping some categories.

Code

# create a more balanced 
# Merge categories Ia - III (high protection)
# Keep category IV as middle term protection
# Merge V and IV (sustainable use)

subdata <- subdata %>%
  mutate(IUCN_CAT_simple = case_when(IUCN_CAT == "Ia" ~ "Strict protection",
                                     IUCN_CAT == "Ib" ~ "Strict protection",
                                     IUCN_CAT == "II" ~ "Strict protection",
                                     IUCN_CAT == "III" ~ "Strict protection",
                                     IUCN_CAT == "IV" ~ "Allow management",
                                     IUCN_CAT == "V" ~ "Allow management",
                                     IUCN_CAT == "VI" ~ "Sustainable use"))

subdata$IUCN_CAT_simple <- factor(subdata$IUCN_CAT_simple, 
                                  levels = c("Strict protection","Allow management","Sustainable use"))

table(subdata$IUCN_CAT_simple)


Strict protection  Allow management   Sustainable use 
              937              1649               364

10.1 DBI

Code

# country
mod_DBI_c <- glmmTMB(DBI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1|ADM0_A3), 
                     data = subdata)
summary(mod_DBI_c)

 Family: gaussian  ( identity )
Formula:          
DBI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1 | ADM0_A3)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   618.3    672.2   -300.1    600.3     2941 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 ADM0_A3  (Intercept) 0.04781  0.2187  
 Residual             0.06735  0.2595  
Number of obs: 2950, groups:  ADM0_A3, 100

Dispersion estimate for gaussian family (sigma^2): 0.0674 

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)  
(Intercept)                     -0.056140   0.040144  -1.399   0.1620  
IUCN_CAT_simpleAllow management  0.011561   0.012766   0.906   0.3651  
IUCN_CAT_simpleSustainable use   0.005667   0.018992   0.298   0.7654  
MPA_area_log                    -0.001411   0.001420  -0.994   0.3202  
INCOME_GRPUpper middle income    0.074988   0.057630   1.301   0.1932  
INCOME_GRPLower middle income    0.063263   0.073477   0.861   0.3892  
INCOME_GRPLow income             0.261656   0.121683   2.150   0.0315 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# region
mod_DBI_r <- glmmTMB(DBI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1|REGION_WB), 
                     data = subdata)
summary(mod_DBI_r)

 Family: gaussian  ( identity )
Formula:          
DBI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1 | REGION_WB)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  1368.6   1422.5   -675.3   1350.6     2941 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 REGION_WB (Intercept) 0.02200  0.1483  
 Residual              0.09175  0.3029  
Number of obs: 2950, groups:  REGION_WB, 7

Dispersion estimate for gaussian family (sigma^2): 0.0917 

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                      0.009754   0.060329   0.162 0.871552    
IUCN_CAT_simpleAllow management -0.057042   0.012679  -4.499 6.83e-06 ***
IUCN_CAT_simpleSustainable use  -0.041910   0.019907  -2.105 0.035266 *  
MPA_area_log                    -0.001913   0.001449  -1.321 0.186637    
INCOME_GRPUpper middle income    0.041598   0.023306   1.785 0.074278 .  
INCOME_GRPLower middle income    0.240237   0.028443   8.446  < 2e-16 ***
INCOME_GRPLow income             0.336921   0.096011   3.509 0.000449 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# compare
AIC(mod_DBI_c, mod_DBI_r)

          df      AIC
mod_DBI_c  9  618.285
mod_DBI_r  9 1368.573

10.2 SMI

Code

# country
mod_SMI_c <- glmmTMB(SMI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1|ADM0_A3), 
                     data = subdata)
summary(mod_SMI_c)

 Family: gaussian  ( identity )
Formula:          
SMI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1 | ADM0_A3)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  -387.7   -333.8    202.8   -405.7     2941 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 ADM0_A3  (Intercept) 0.03953  0.1988  
 Residual             0.04770  0.2184  
Number of obs: 2950, groups:  ADM0_A3, 100

Dispersion estimate for gaussian family (sigma^2): 0.0477 

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)  
(Intercept)                     -0.071335   0.036013  -1.981   0.0476 *
IUCN_CAT_simpleAllow management  0.020879   0.010756   1.941   0.0522 .
IUCN_CAT_simpleSustainable use   0.008084   0.016004   0.505   0.6135  
MPA_area_log                    -0.002570   0.001197  -2.148   0.0318 *
INCOME_GRPUpper middle income    0.045123   0.051675   0.873   0.3825  
INCOME_GRPLower middle income    0.031330   0.065600   0.478   0.6329  
INCOME_GRPLow income             0.251763   0.107876   2.334   0.0196 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# region
mod_SMI_r <- glmmTMB(SMI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1|REGION_WB), 
                     data = subdata)
summary(mod_SMI_r)

 Family: gaussian  ( identity )
Formula:          
SMI ~ IUCN_CAT_simple + MPA_area_log + INCOME_GRP + (1 | REGION_WB)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   525.2    579.1   -253.6    507.2     2941 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 REGION_WB (Intercept) 0.01630  0.1277  
 Residual              0.06894  0.2626  
Number of obs: 2950, groups:  REGION_WB, 7

Dispersion estimate for gaussian family (sigma^2): 0.0689 

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -0.016648   0.051956  -0.320  0.74864    
IUCN_CAT_simpleAllow management -0.032643   0.010991  -2.970  0.00298 ** 
IUCN_CAT_simpleSustainable use  -0.033267   0.017256  -1.928  0.05388 .  
MPA_area_log                    -0.002589   0.001256  -2.062  0.03924 *  
INCOME_GRPUpper middle income    0.025892   0.020215   1.281  0.20024    
INCOME_GRPLower middle income    0.225445   0.024657   9.143  < 2e-16 ***
INCOME_GRPLow income             0.339839   0.083415   4.074 4.62e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# compare
AIC(mod_SMI_c, mod_SMI_r)

          df       AIC
mod_SMI_c  9 -387.6952
mod_SMI_r  9  525.2297

10.3 Visualize effects

10.3.0.1 Country models

10.3.0.1.0.1 DBI

Code

emmeans_DBI <- emmeans(mod_DBI_c,
                       type = "response",
                       specs = "IUCN_CAT_simple") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = IUCN_CAT_simple, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

emmeans_DBI <- emmeans(mod_DBI_c,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

No significant effects for DBI in models using country as random intercept…

10.3.0.1.0.2 SMI

Code

emmeans_SMI <- emmeans(mod_SMI_c,
                       type = "response",
                       specs = "IUCN_CAT_simple") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = IUCN_CAT_simple, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

emmeans_SMI <- emmeans(mod_SMI_c,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

10.3.0.2 Region models

10.3.0.2.0.1 DBI

Code

emmeans_DBI <- emmeans(mod_DBI_r,
                       type = "response",
                       specs = "IUCN_CAT_simple") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = IUCN_CAT_simple, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

emmeans_DBI <- emmeans(mod_DBI_r,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_DBI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "DBI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Significant effects using region as random intercept…

Protected areas in categories more strict protection models show positive DBIs.

Protected areas in middle to low and low income countries show positive DBIs.

10.3.0.2.0.2 SMI

Code

emmeans_SMI <- emmeans(mod_SMI_r,
                       type = "response",
                       specs = "IUCN_CAT_simple") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = IUCN_CAT_simple, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "IUCN categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Code

emmeans_SMI <- emmeans(mod_SMI_r,
                       type = "response",
                       specs = "INCOME_GRP") %>%
  as.data.frame() %>%
  mutate(Test = ifelse(sign(lower.CL)==sign(upper.CL),"significant","non-significant"))

ggplot(emmeans_SMI, aes(x = INCOME_GRP, y = emmean, fill = Test)) +
  scale_fill_manual(values = c("non-significant"="white","significant"="black")) +
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL),
                 position = position_dodge(.5))+
  geom_point(position = position_dodge(.5), size = 3, shape = 21) +
  labs(x = "Income categoty", y = "SMI")+
  geom_hline(yintercept = 0)+
  theme_classic()

Same pattern for SMI.