## Setup ###################################################################################################
# To run at the command line, replace this with a function that calls setwd()
# to the local directory where you have the data and R files.
setwd_PROCOAST_pluralistic_ignorance_1()
# Results in packages and own-written analysis functions loaded into memory
source( "analysisFunctions.R", local=T )
# Results in a data frame dF loaded into memory
# dF only has data points that are complete cases
# (both personal and perceived attitude) for at least one of ESI and WTS
# dFWithIncompleteCases is also loaded, which has all data
source( "dataLoadAndPrep.R", local=T )
## Sample description
# Number of survey responses
nrow(dFWithIncompleteCases)
## [1] 453
# Number complete cases for at least one of ESI and WTS
nrow(dF)
## [1] 337
# Proportion with university degrees
sum( table(dF$education)[c("6","8")] ) / sum( table(dF$education) )
## [1] 0.6497696
# Proportions in age bands
dF %>%
mutate(age_group = case_when(
age <= 30 ~ "18-30",
age <= 50 ~ "31-50",
age <= 70 ~ "51-70",
age > 70 ~ "71+"
)) %>%
count(age_group) %>%
mutate(proportion = n / sum(n))
## # A tibble: 5 × 3
## age_group n proportion
## <chr> <int> <dbl>
## 1 18-30 93 0.276
## 2 31-50 95 0.282
## 3 51-70 71 0.211
## 4 71+ 36 0.107
## 5 <NA> 42 0.125
# Proportions of genders
table(dF$gender)/sum(table(dF$gender))
##
## man nonBinary woman
## 0.48986486 0.01013514 0.50000000
# Mean economic situation
mean(dF$econSituation,na.rm=T)
## [1] 3.75
# Describe stakeholder group sample sizes
# Note, stakeholder groups are specific to each case study,
# unlike stakeholder categories, which are generally applicable.
summarise_stakeholders(dF)
## Stakeholder groups - Total (ESI complete, WTS complete):
##
## Estonia, Environment protectors, National Nature Conservation Institutions: 2 (2, 2)
## Estonia, Government, Local Government: 2 (2, 0)
## Estonia, Local community, Local (Permanent) Community Representatives: 10 (10, 0)
## Estonia subtotal: 14 (14, 2)
##
## Ireland, Commercial, Farmer/ Landowner: 6 (6, 6)
## Ireland, Environment protectors, NGO: 1 (1, 1)
## Ireland, Environment protectors, NPWS (National Park & Wildlife Services): 1 (1, 1)
## Ireland, Government, SCC (Sligo County Council): 6 (6, 6)
## Ireland subtotal: 14 (14, 14)
##
## Italy, Environment protectors, NGO: 7 (7, 0)
## Italy, Government, Local Government: 1 (1, 0)
## Italy, Local community, Higher education: 16 (16, 0)
## Italy, Local community, Local citizens: 10 (10, 0)
## Italy subtotal: 34 (34, 0)
##
## Malta, Commercial, Tour operator, Boat operator: 6 (6, 6)
## Malta, Government, Environment Ministry & Relevant Managing Authorities: 11 (11, 10)
## Malta, Government, Malta Tourism Authority Ministry: 3 (3, 3)
## Malta, Government, Ministry of Gozo & Local Council: 1 (1, 1)
## Malta, Tourists, Local visitors / Foreign tourists: 73 (73, 72)
## Malta subtotal: 94 (94, 92)
##
## Norway, Environment protectors, Divers and others who clean up underwater: 15 (15, 15)
## Norway, Environment protectors, Environmental organisations: 5 (5, 5)
## Norway, Environment protectors, Possible financiers for cleanup operations: 1 (1, 1)
## Norway, Government, Directorate of Fisheries: 1 (1, 1)
## Norway, Wild resource users, Fishers: 31 (30, 31)
## Norway subtotal: 53 (52, 53)
##
## Romania, Commercial, Small Business Owners: 1 (0, 1)
## Romania, Government, Local Authorities: 3 (0, 3)
## Romania, Tourists, Tourists: 12 (0, 12)
## Romania subtotal: 16 (0, 16)
##
## Slovenia, Commercial, Farmers: 9 (9, 0)
## Slovenia, Environment protectors, Protected area practitioners: 3 (3, 0)
## Slovenia, Local community, Local community: 59 (59, 0)
## Slovenia subtotal: 71 (71, 0)
##
## UK, Commercial, Farmer/Forester: 2 (2, 2)
## UK, Environment protectors, Nature conservation professional: 3 (3, 3)
## UK, Local community, Community member: 16 (15, 16)
## UK, Wild resource users, Hunter/Fisher: 3 (3, 3)
## UK, Wild resource users, Recreational wildlife user: 17 (17, 17)
## UK subtotal: 41 (40, 41)
##
## Stakeholder category totals across all case studies:
## Commercial: 24 (23, 15)
## Environment protectors: 38 (38, 28)
## Government: 28 (25, 24)
## Local community: 111 (110, 16)
## Tourists: 85 (73, 84)
## Wild resource users: 51 (50, 51)
##
## GRAND TOTAL: 337 (319, 218)
##
## Number of stakeholder groups: 32
## Median size: 5.5
## Mean size: 10.5
## Reliability
#ESI
cor.test(dF$envCit1,dF$envCit2)
##
## Pearson's product-moment correlation
##
## data: dF$envCit1 and dF$envCit2
## t = 14.954, df = 331, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5661261 0.6950130
## sample estimates:
## cor
## 0.6349674
sum(complete.cases(dF$envCit1, dF$envCit2))
## [1] 333
#WTS
cor.test(dF$willToSac1,dF$willToSac2)
##
## Pearson's product-moment correlation
##
## data: dF$willToSac1 and dF$willToSac2
## t = 13.246, df = 228, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5795478 0.7268216
## sample estimates:
## cor
## 0.6594665
sum(complete.cases(dF$willToSac1, dF$willToSac2))
## [1] 230
## Pluralistic analysis #############################################################
## Starting with agreement (H1)
## Variable descriptions
# Grand mean of percent agree with ESI item
grand_mean(dF$ESI >= 4, dF$complete_cases_ESI, dF$caseStudy, function(x) x * 100)
## [1] 86.6744
# Grand mean of estimated percent agree with ESI item
grand_mean(dF$ESI_norm, dF$complete_cases_ESI, dF$caseStudy)
## [1] 59.08489
# Grand mean of percent agree with WTS item
grand_mean(dF$WTS >= 4, dF$complete_cases_WTS, dF$caseStudy, function(x) x * 100)
## [1] 52.54095
# Grand mean of estimated percent agree with WTS item
grand_mean(dF$WTS_norm, dF$complete_cases_WTS, dF$caseStudy)
## [1] 30.33032
## Pre-registered tests
# H1 for ESI
H1_results_ESI <- permutation_sign_test(dF$ESI_EPI_score)
print(H1_results_ESI)
## Permutation Sign Test
## =====================
##
## Observed mean: 25.201
## Alternative hypothesis: two.sided
## Number of permutations: 99999
## Number of observations: 321
## P-value: 0
## Significance: ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# H1 for WTS
H1_results_WTS <- permutation_sign_test(dF$WTS_EPI_score)
print(H1_results_WTS)
## Permutation Sign Test
## =====================
##
## Observed mean: 25.6169
## Alternative hypothesis: two.sided
## Number of permutations: 99999
## Number of observations: 220
## P-value: 0
## Significance: ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Now disagreement (H2)
# Grand mean of percent disagree with ESI item
grand_mean(dF$ESI <= 2, dF$complete_cases_ESI, dF$caseStudy, function(x) x * 100)
## [1] 3.819626
# Grand mean of estimated percent disagree with ESI item
grand_mean(dF$ESI_norm_dis, dF$complete_cases_ESI, dF$caseStudy)
## [1] 36.86891
# Grand mean of percent disagree with WTS item
grand_mean(dF$WTS <= 2, dF$complete_cases_WTS, dF$caseStudy, function(x) x * 100)
## [1] 20.67506
# Grand mean of estimated percent disagree with WTS item
grand_mean(dF$WTS_norm_dis, dF$complete_cases_WTS, dF$caseStudy)
## [1] 42.80778
## Pre-registered tests
# H1 for ESI
H2_results_ESI <- permutation_sign_test(dF$ESI_EPI_score_dis)
print(H2_results_ESI)
## Permutation Sign Test
## =====================
##
## Observed mean: -35.2366
## Alternative hypothesis: two.sided
## Number of permutations: 99999
## Number of observations: 215
## P-value: 0
## Significance: ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# H1 for WTS
H2_results_WTS <- permutation_sign_test(dF$WTS_EPI_score_dis)
print(H2_results_WTS)
## Permutation Sign Test
## =====================
##
## Observed mean: -28.0421
## Alternative hypothesis: two.sided
## Number of permutations: 99999
## Number of observations: 219
## P-value: 0
## Significance: ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Plotting pluralistic ignorance
table_for_plot_ESI <- create_pluralistic_ignorance_group_results_table(dF,"ESI")
table_for_plot_ESI
## # A tibble: 16 × 5
## caseStudy stakeholderCategory propAgree propPercNormAgree n
## <chr> <chr> <dbl> <dbl> <int>
## 1 Malta Commercial 1 0.333 6
## 2 Malta Tourists 0.959 0.616 73
## 3 Malta Government 0.867 0.633 15
## 4 Norway Wild resource users 0.8 0.59 30
## 5 Norway Environment protectors 0.952 0.757 21
## 6 Estonia Local community 0.8 0.53 10
## 7 UK Wild resource users 1 0.67 20
## 8 UK Local community 0.733 0.563 15
## 9 UK Environment protectors 1 0.883 3
## 10 Slovenia Local community 0.847 0.622 59
## 11 Slovenia Commercial 0.444 0.678 9
## 12 Slovenia Environment protectors 1 0.767 3
## 13 Ireland Government 0.833 0.383 6
## 14 Ireland Commercial 0.833 0.217 6
## 15 Italy Local community 1 0.727 26
## 16 Italy Environment protectors 1 0.771 7
table_for_plot_WTS <- create_pluralistic_ignorance_group_results_table(dF,"WTS")
table_for_plot_WTS
## # A tibble: 12 × 5
## caseStudy stakeholderCategory propAgree propPercNormAgree n
## <chr> <chr> <dbl> <dbl> <int>
## 1 Malta Commercial 0.5 0.15 6
## 2 Malta Tourists 0.708 0.344 72
## 3 Malta Government 0.786 0.443 14
## 4 Norway Wild resource users 0.645 0.326 31
## 5 Norway Environment protectors 0.762 0.543 21
## 6 Romania Tourists 0.333 0.35 12
## 7 Romania Government 0.667 0.333 3
## 8 UK Wild resource users 0.75 0.41 20
## 9 UK Local community 0.562 0.325 16
## 10 UK Environment protectors 1 0.75 3
## 11 Ireland Government 0.167 0.15 6
## 12 Ireland Commercial 0.333 0.183 6
combined_plot <- create_pluralistic_ignorance_facet_plot(
data_list = list(table_for_plot_ESI, table_for_plot_WTS),
attitude_vars = c("ESI", "WTS")
)
print(combined_plot)
ggsave("pluralistic_ignorance_plot.png",
plot = combined_plot,
width = 6.7,
height = 6.2,
units = "in",
dpi = 300,
bg = "white")
## Generality analysis
set.seed(123)
mod_EPI_ESI <- lmer( ESI_EPI_score ~ stakeholderCategory + ( 1 | caseStudy), data=dF )
mod_EPI_WTS <- lmer( WTS_EPI_score ~ stakeholderCategory + ( 1 | caseStudy), data=dF )
# Run the bootstrap
boot_results_ESI <- bootMer(mod_EPI_ESI, calc_proportion_positive_groups,
nsim = 99999, use.u = FALSE, type = "parametric")
boot_results_WTS <- bootMer(mod_EPI_WTS, calc_proportion_positive_groups,
nsim = 99999, use.u = FALSE, type = "parametric")
ci_ESI <- quantile(boot_results_ESI$t, c(0.025, 0.975))
cat(" ESI Point estimate:", calc_proportion_positive_groups(mod_EPI_ESI), "\n",
"ESI Bootstrap mean:", mean(boot_results_ESI$t), "\n",
"ESI 95% CI: [", ci_ESI[1], ",", ci_ESI[2], "]\n")
## ESI Point estimate: 0.802537
## ESI Bootstrap mean: 0.7983417
## ESI 95% CI: [ 0.5818582 , 0.9671808 ]
ci_WTS <- quantile(boot_results_WTS$t, c(0.025, 0.975))
cat(" WTS Point estimate:", calc_proportion_positive_groups(mod_EPI_WTS), "\n",
"WTS Bootstrap mean:", mean(boot_results_WTS$t), "\n",
"WTS 95% CI: [", ci_WTS[1], ",", ci_WTS[2], "]\n")
## WTS Point estimate: 0.8794341
## WTS Bootstrap mean: 0.8646631
## WTS 95% CI: [ 0.6107596 , 1 ]
## Modelling of personal attitude ~ perceived norm ######################################################
# Linear models are just not OK with this heavily skewed ordinal data, as revealed by the residual plot.
# So as per prereg, we are going to do ordinal models instead.
formula_ESI_1 <- ESI ~ ESI_norm
model_ESI_1 <- lm(formula_ESI_1, data = standardised_complete_cases( dF, formula_ESI_1 ) )
plot(jitter(fitted(model_ESI_1)), jitter(residuals(model_ESI_1)), xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs Fitted")
## ESI models
# Initialize list to store results of all models
model_results_ESI <- list()
formula_ESI_1 <- ESI_factor ~ ESI_norm
model_ESI_1_ordinal <- clm(formula_ESI_1, data = standardised_complete_cases(dF, formula_ESI_1))
model_results_ESI[[1]] <- summarise_model_with_OR(model_ESI_1_ordinal, "ESI_norm")
## formula: ESI_factor ~ ESI_norm
## data: standardised_complete_cases(dF, formula_ESI_1)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 319 -318.77 647.54 5(0) 5.19e-08 3.1e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## ESI_norm 0.3056 0.1155 2.646 0.00814 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -4.4074 0.5039 -8.746
## 2|3 -3.1922 0.2843 -11.230
## 3|4 -2.0500 0.1759 -11.651
## 4|5 -0.4244 0.1157 -3.667
##
## AIC: 647.5
##
## Odds Ratios and 95% Confidence Intervals:
## ESI_norm: OR = 1.36, 95% CI [1.08, 1.71], p = 0.00814
#Testing the proportional odds assumption for ESI
nominal_test(model_ESI_1_ordinal)
## Tests of nominal effects
##
## formula: ESI_factor ~ ESI_norm
## Df logLik AIC LRT Pr(>Chi)
## <none> -318.77 647.54
## ESI_norm 3 -315.18 646.35 7.191 0.06605 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
formula_ESI_2 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI
model_ESI_2_ordinal <- clm(formula_ESI_2, data = standardised_complete_cases(dF, formula_ESI_2))
model_results_ESI[[2]] <- summarise_model_with_OR(model_ESI_2_ordinal, "ESI_norm")
## formula: ESI_factor ~ ESI_norm + stakeholderCategoryForESI
## data: standardised_complete_cases(dF, formula_ESI_2)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 319 -308.92 637.84 5(0) 5.72e-08 4.5e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## ESI_norm 0.33348 0.12258 2.720 0.00652 **
## stakeholderCategoryForESICommercial 0.04715 0.46188 0.102 0.91870
## stakeholderCategoryForESIEnvironment protectors 0.42467 0.37694 1.127 0.25990
## stakeholderCategoryForESIGovernment 0.57373 0.47844 1.199 0.23046
## stakeholderCategoryForESITourists 1.36436 0.34620 3.941 8.12e-05 ***
## stakeholderCategoryForESIWild resource users 0.07728 0.32408 0.238 0.81153
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -4.11895 0.52125 -7.902
## 2|3 -2.89838 0.31451 -9.216
## 3|4 -1.73504 0.22248 -7.798
## 4|5 -0.04159 0.18618 -0.223
##
## AIC: 637.8
##
## Odds Ratios and 95% Confidence Intervals:
## ESI_norm: OR = 1.40, 95% CI [1.10, 1.78], p = 0.00652
formula_ESI_3 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy)
model_ESI_3_ordinal <- clmm(formula_ESI_3, data = standardised_complete_cases(dF, formula_ESI_3))
model_results_ESI[[3]] <- summarise_model_with_OR(model_ESI_3_ordinal, "ESI_norm")
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1 | caseStudy)
## data: standardised_complete_cases(dF, formula_ESI_3)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 319 -296.89 615.79 639(2534) 3.71e-06 1.3e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 1.263 1.124
## Number of groups: caseStudy 7
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## ESI_norm 0.3522 0.1340 2.630 0.00855 **
## stakeholderCategoryForESICommercial -0.5897 0.5334 -1.105 0.26895
## stakeholderCategoryForESIEnvironment protectors 0.7124 0.5371 1.326 0.18474
## stakeholderCategoryForESIGovernment -0.9145 0.7331 -1.247 0.21222
## stakeholderCategoryForESITourists -0.7876 0.8770 -0.898 0.36917
## stakeholderCategoryForESIWild resource users 0.4081 0.5292 0.771 0.44059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -4.8889 0.7333 -6.667
## 2|3 -3.6527 0.6018 -6.069
## 3|4 -2.4220 0.5496 -4.407
## 4|5 -0.5405 0.5223 -1.035
##
## AIC: 615.8
##
## Odds Ratios and 95% Confidence Intervals:
## ESI_norm: OR = 1.42, 95% CI [1.09, 1.85], p = 0.00855
formula_ESI_4 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy) + ESI_norm:stakeholderCategoryForESI
model_ESI_4_ordinal <- clmm(formula_ESI_4, data = standardised_complete_cases(dF, formula_ESI_4))
model_results_ESI[[4]] <- summarise_model_with_OR(model_ESI_4_ordinal, "ESI_norm")
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1 | caseStudy) + ESI_norm:stakeholderCategoryForESI
## data: standardised_complete_cases(dF, formula_ESI_4)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 319 -295.81 623.61 1259(4897) 1.86e-04 1.6e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 1.137 1.066
## Number of groups: caseStudy 7
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## ESI_norm 0.3263 0.2062 1.582 0.114
## stakeholderCategoryForESICommercial -0.6389 0.5328 -1.199 0.230
## stakeholderCategoryForESIEnvironment protectors 0.5598 0.5509 1.016 0.310
## stakeholderCategoryForESIGovernment -0.6988 0.8751 -0.799 0.425
## stakeholderCategoryForESITourists -0.6106 0.9785 -0.624 0.533
## stakeholderCategoryForESIWild resource users 0.4341 0.5320 0.816 0.415
## ESI_norm:stakeholderCategoryForESICommercial -0.1905 0.4692 -0.406 0.685
## ESI_norm:stakeholderCategoryForESIEnvironment protectors 0.4817 0.5224 0.922 0.356
## ESI_norm:stakeholderCategoryForESIGovernment 0.2101 0.5571 0.377 0.706
## ESI_norm:stakeholderCategoryForESITourists -0.1676 0.3547 -0.473 0.637
## ESI_norm:stakeholderCategoryForESIWild resource users 0.2182 0.4027 0.542 0.588
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -4.8507 0.7319 -6.628
## 2|3 -3.6149 0.6002 -6.023
## 3|4 -2.3816 0.5467 -4.356
## 4|5 -0.4866 0.5171 -0.941
##
## AIC: 623.6
##
## Odds Ratios and 95% Confidence Intervals:
## ESI_norm: OR = 1.39, 95% CI [0.93, 2.08], p = 0.114
formula_ESI_5 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy) + (0 + ESI_norm|caseStudy)
model_ESI_5_ordinal <- clmm(formula_ESI_5, data = standardised_complete_cases(dF, formula_ESI_5))
model_results_ESI[[5]] <- summarise_model_with_OR(model_ESI_5_ordinal, "ESI_norm")
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1 | caseStudy) + (0 + ESI_norm | caseStudy)
## data: standardised_complete_cases(dF, formula_ESI_5)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 319 -296.59 617.17 864(3434) 3.15e-05 1.2e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 1.26000 1.122
## caseStudy ESI_norm 0.06004 0.245
## Number of groups: caseStudy 7
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## ESI_norm 0.3516 0.1710 2.056 0.0398 *
## stakeholderCategoryForESICommercial -0.6055 0.5362 -1.129 0.2588
## stakeholderCategoryForESIEnvironment protectors 0.6750 0.5375 1.256 0.2092
## stakeholderCategoryForESIGovernment -0.8720 0.7351 -1.186 0.2355
## stakeholderCategoryForESITourists -0.7469 0.8764 -0.852 0.3940
## stakeholderCategoryForESIWild resource users 0.4903 0.5413 0.906 0.3651
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -4.8936 0.7340 -6.667
## 2|3 -3.6531 0.6026 -6.062
## 3|4 -2.4114 0.5505 -4.381
## 4|5 -0.5051 0.5250 -0.962
##
## AIC: 617.2
##
## Odds Ratios and 95% Confidence Intervals:
## ESI_norm: OR = 1.42, 95% CI [1.02, 1.99], p = 0.0398
comparison_table_ESI <- create_comparison_table(model_results_ESI, "ESI_norm")
## WTS models
# Initialize list to store results of all models
model_results_WTS <- list()
formula_WTS_1 <- WTS_factor ~ WTS_norm
model_WTS_1_ordinal <- clm(formula_WTS_1, data = standardised_complete_cases(dF, formula_WTS_1))
model_results_WTS[[1]] <- summarise_model_with_OR(model_WTS_1_ordinal, "WTS_norm")
## formula: WTS_factor ~ WTS_norm
## data: standardised_complete_cases(dF, formula_WTS_1)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 218 -311.31 632.62 5(0) 3.46e-13 1.3e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## WTS_norm 0.5872 0.1292 4.545 5.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -2.6406 0.2632 -10.031
## 2|3 -1.3851 0.1702 -8.140
## 3|4 -0.5841 0.1449 -4.031
## 4|5 0.9712 0.1557 6.236
##
## AIC: 632.6
##
## Odds Ratios and 95% Confidence Intervals:
## WTS_norm: OR = 1.80, 95% CI [1.40, 2.33], p = <0.001
#Testing the proportional odds assumption for WTS
nominal_test(model_WTS_1_ordinal)
## Tests of nominal effects
##
## formula: WTS_factor ~ WTS_norm
## Df logLik AIC LRT Pr(>Chi)
## <none> -311.31 632.62
## WTS_norm 3 -305.92 627.85 10.779 0.01299 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
formula_WTS_2 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS
model_WTS_2_ordinal <- clm(formula_WTS_2, data = standardised_complete_cases(dF, formula_WTS_2))
model_results_WTS[[2]] <- summarise_model_with_OR(model_WTS_2_ordinal, "WTS_norm")
## formula: WTS_factor ~ WTS_norm + stakeholderCategoryForWTS
## data: standardised_complete_cases(dF, formula_WTS_2)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 218 -309.07 638.13 5(0) 8.68e-13 6.9e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## WTS_norm 0.60140 0.13780 4.364 1.28e-05 ***
## stakeholderCategoryForWTSCommercial -0.75796 0.49294 -1.538 0.124
## stakeholderCategoryForWTSEnvironment protectors -0.57526 0.40480 -1.421 0.155
## stakeholderCategoryForWTSGovernment -0.03572 0.42435 -0.084 0.933
## stakeholderCategoryForWTSLocal community -0.46171 0.51879 -0.890 0.373
## stakeholderCategoryForWTSWild resource users -0.36321 0.32474 -1.118 0.263
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -2.9229 0.3169 -9.224
## 2|3 -1.6602 0.2428 -6.838
## 3|4 -0.8530 0.2241 -3.807
## 4|5 0.7243 0.2222 3.259
##
## AIC: 638.1
##
## Odds Ratios and 95% Confidence Intervals:
## WTS_norm: OR = 1.82, 95% CI [1.40, 2.40], p = <0.001
formula_WTS_3 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy)
model_WTS_3_ordinal <- clmm(formula_WTS_3, data = standardised_complete_cases(dF, formula_WTS_3))
model_results_WTS[[3]] <- summarise_model_with_OR(model_WTS_3_ordinal, "WTS_norm")
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1 | caseStudy)
## data: standardised_complete_cases(dF, formula_WTS_3)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 218 -308.35 638.69 741(1510) 5.59e-05 2.3e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.1461 0.3823
## Number of groups: caseStudy 6
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## WTS_norm 0.58417 0.14027 4.165 3.12e-05 ***
## stakeholderCategoryForWTSCommercial -0.46197 0.55777 -0.828 0.408
## stakeholderCategoryForWTSEnvironment protectors -0.03871 0.63601 -0.061 0.951
## stakeholderCategoryForWTSGovernment 0.18022 0.46770 0.385 0.700
## stakeholderCategoryForWTSLocal community -0.17985 0.68658 -0.262 0.793
## stakeholderCategoryForWTSWild resource users 0.09603 0.56818 0.169 0.866
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -2.5819 0.4568 -5.652
## 2|3 -1.3055 0.4163 -3.136
## 3|4 -0.4839 0.4137 -1.170
## 4|5 1.1291 0.4304 2.624
##
## AIC: 638.7
##
## Odds Ratios and 95% Confidence Intervals:
## WTS_norm: OR = 1.79, 95% CI [1.36, 2.36], p = <0.001
formula_WTS_4 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy) + WTS_norm:stakeholderCategoryForWTS
model_WTS_4_ordinal <- clmm(formula_WTS_4, data = standardised_complete_cases(dF, formula_WTS_4))
model_results_WTS[[4]] <- summarise_model_with_OR(model_WTS_4_ordinal, "WTS_norm")
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1 | caseStudy) + WTS_norm:stakeholderCategoryForWTS
## data: standardised_complete_cases(dF, formula_WTS_4)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 218 -305.02 642.04 1218(2474) 1.25e-04 2.4e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.1845 0.4296
## Number of groups: caseStudy 6
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## WTS_norm 0.56129 0.22973 2.443 0.0146 *
## stakeholderCategoryForWTSCommercial -0.84458 0.75498 -1.119 0.2633
## stakeholderCategoryForWTSEnvironment protectors -0.28531 0.68177 -0.418 0.6756
## stakeholderCategoryForWTSGovernment 0.20739 0.47069 0.441 0.6595
## stakeholderCategoryForWTSLocal community -0.13604 0.72429 -0.188 0.8510
## stakeholderCategoryForWTSWild resource users 0.23926 0.62011 0.386 0.6996
## WTS_norm:stakeholderCategoryForWTSCommercial -0.47191 0.68013 -0.694 0.4878
## WTS_norm:stakeholderCategoryForWTSEnvironment protectors 0.56993 0.39895 1.429 0.1531
## WTS_norm:stakeholderCategoryForWTSGovernment 0.05693 0.41783 0.136 0.8916
## WTS_norm:stakeholderCategoryForWTSLocal community 0.59305 0.69114 0.858 0.3909
## WTS_norm:stakeholderCategoryForWTSWild resource users -0.44275 0.37714 -1.174 0.2404
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -2.6025 0.4682 -5.559
## 2|3 -1.2982 0.4286 -3.029
## 3|4 -0.4598 0.4267 -1.078
## 4|5 1.1878 0.4433 2.680
##
## AIC: 642.0
##
## Odds Ratios and 95% Confidence Intervals:
## WTS_norm: OR = 1.75, 95% CI [1.12, 2.75], p = 0.0146
formula_WTS_5 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy) + (0 + WTS_norm|caseStudy)
model_WTS_5_ordinal <- clmm(formula_WTS_5, data = standardised_complete_cases(dF, formula_WTS_5))
model_results_WTS[[5]] <- summarise_model_with_OR(model_WTS_5_ordinal, "WTS_norm")
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1 | caseStudy) + (0 + WTS_norm | caseStudy)
## data: standardised_complete_cases(dF, formula_WTS_5)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 218 -308.35 640.69 775(1578) 2.23e-04 2.1e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.1461 0.3823
## caseStudy WTS_norm 0.0000 0.0000
## Number of groups: caseStudy 6
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## WTS_norm 0.58418 0.14027 4.165 3.12e-05 ***
## stakeholderCategoryForWTSCommercial -0.46198 0.55776 -0.828 0.408
## stakeholderCategoryForWTSEnvironment protectors -0.03871 0.63593 -0.061 0.951
## stakeholderCategoryForWTSGovernment 0.18023 0.46769 0.385 0.700
## stakeholderCategoryForWTSLocal community -0.17987 0.68654 -0.262 0.793
## stakeholderCategoryForWTSWild resource users 0.09602 0.56811 0.169 0.866
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -2.5820 0.4567 -5.653
## 2|3 -1.3055 0.4162 -3.137
## 3|4 -0.4839 0.4135 -1.170
## 4|5 1.1290 0.4302 2.624
##
## AIC: 640.7
##
## Odds Ratios and 95% Confidence Intervals:
## WTS_norm: OR = 1.79, 95% CI [1.36, 2.36], p = <0.001
comparison_table_WTS <- create_comparison_table(model_results_WTS, "WTS_norm")
print(comparison_table_ESI)
## Model OR_CI p AIC
## 1 Model 1 1.36 [1.08, 1.71] 0.008 647.5
## 2 Model 2 1.40 [1.10, 1.78] 0.007 637.8
## 3 Model 3 1.42 [1.09, 1.85] 0.009 615.8
## 4 Model 4 1.39 [0.93, 2.08] 0.114 623.6
## 5 Model 5 1.42 [1.02, 1.99] 0.040 617.2
print(comparison_table_WTS)
## Model OR_CI p AIC
## 1 Model 1 1.80 [1.40, 2.33] 0.000 632.6
## 2 Model 2 1.82 [1.40, 2.40] 0.000 638.1
## 3 Model 3 1.79 [1.36, 2.36] 0.000 638.7
## 4 Model 4 1.75 [1.12, 2.75] 0.015 642.0
## 5 Model 5 1.79 [1.36, 2.36] 0.000 640.7
formula_WTS_6 <- WTS_factor ~ WTS_norm + ESI + stakeholderCategoryForWTS + (1|caseStudy)
model_WTS_6_ordinal <- clmm(formula_WTS_6, data = standardised_complete_cases(dF, formula_WTS_6))
summarise_model_with_OR(model_WTS_6_ordinal, c( "WTS_norm", "ESI" ) )
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: WTS_factor ~ WTS_norm + ESI + stakeholderCategoryForWTS + (1 | caseStudy)
## data: standardised_complete_cases(dF, formula_WTS_6)
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 215 -293.49 610.99 1214(1573) 9.19e-06 2.4e+06
##
## Random effects:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 6.271e-09 7.919e-05
## Number of groups: caseStudy 6
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## WTS_norm 0.56142 0.14091 3.984 6.77e-05 ***
## ESI 0.59175 0.12906 4.585 4.54e-06 ***
## stakeholderCategoryForWTSCommercial -1.00361 0.52093 -1.927 0.054 .
## stakeholderCategoryForWTSEnvironment protectors -0.49176 0.41323 -1.190 0.234
## stakeholderCategoryForWTSGovernment 0.07166 0.43127 0.166 0.868
## stakeholderCategoryForWTSLocal community 0.18235 0.52737 0.346 0.730
## stakeholderCategoryForWTSWild resource users -0.17076 0.33006 -0.517 0.605
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -3.0474 0.3311 -9.205
## 2|3 -1.6675 0.2461 -6.774
## 3|4 -0.8189 0.2261 -3.622
## 4|5 0.8691 0.2262 3.842
##
## AIC: 611.0
##
## Odds Ratios and 95% Confidence Intervals:
## WTS_norm: OR = 1.75, 95% CI [1.33, 2.31], p = <0.001
## ESI: OR = 1.81, 95% CI [1.40, 2.33], p = <0.001
exp(coef(model_WTS_6_ordinal))[c("WTS_norm", "ESI")]
## WTS_norm ESI
## 1.753162 1.807150
exp(confint(model_WTS_6_ordinal))[c("WTS_norm", "ESI"),]
## 2.5 % 97.5 %
## WTS_norm 1.330089 2.310806
## ESI 1.403253 2.327300
# To render:
# rmarkdown::render("YOUR_ABSOLUTE_PATH/analysis.R", output_format = rmarkdown::html_document( pandoc_args = c("--metadata", "author=Anon For Peer Review")))