## 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 )
## Warning: package 'ggplot2' was built under R version 4.5.2
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: Matrix
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
##
## Attaching package: 'ordinal'
## The following object is masked from 'package:dplyr':
##
## slice
# Load additional required packages
library(glmmTMB)
# Function to rescale ordinal data to (0,1) for beta regression
rescale_to_beta <- function(x, min_val = 1, max_val = 5) {
# Rescale from [min, max] to (0, 1) using the transformation:
# (y - min + c) / (max - min + 2c) where c = 0.5
c <- 0.5
(x - min_val + c) / (max_val - min_val + 2 * c)
}
# Function to summarize beta regression models with effect sizes
summarise_beta_model <- function(model, params) {
print(summary(model))
cat("\n")
aic_value <- AIC(model)
cat(sprintf("AIC: %.1f\n\n", aic_value))
# For beta regression, we'll use conditional R-squared from performance package
# or just skip if it causes issues
tryCatch({
r2 <- performance::r2(model)
cat(sprintf("Conditional R²: %.3f\n", r2$R2_conditional))
cat(sprintf("Marginal R²: %.3f\n\n", r2$R2_marginal))
}, error = function(e) {
cat("R² calculation skipped\n\n")
})
# Extract coefficients for specified parameters
summ <- summary(model)
coef_table <- coef(summ)$cond # conditional model coefficients
# For beta regression, coefficients are interpreted differently than odds ratios
# We'll report exp(coef) as multiplicative effects on the odds
effects <- exp(coef_table[params, "Estimate"])
names(effects) <- params
# Calculate confidence intervals
ci_lower <- exp(coef_table[params, "Estimate"] - 1.96 * coef_table[params, "Std. Error"])
ci_upper <- exp(coef_table[params, "Estimate"] + 1.96 * coef_table[params, "Std. Error"])
ci <- cbind(ci_lower, ci_upper)
rownames(ci) <- params
p_values <- coef_table[params, "Pr(>|z|)"]
names(p_values) <- params
cat("Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:\n")
for (param in params) {
cat(sprintf("%s: Exp(B) = %.2f, 95%% CI [%.2f, %.2f], p = %s\n",
param, effects[param], ci[param, 1], ci[param, 2],
format.pval(p_values[param], digits = 3, eps = 0.001)))
}
cat("\n")
result <- list(
effects = effects,
confidence_intervals = ci,
p_values = p_values,
aic = aic_value
)
invisible(result)
}
# Create comparison table for beta models
create_beta_comparison_table <- function(model_results, param) {
n_models <- length(model_results)
comparison_table <- data.frame(
Model = paste("Model", 1:n_models),
Effect_CI = sapply(1:n_models, function(i) {
eff <- model_results[[i]]$effects[param]
ci_lower <- model_results[[i]]$confidence_intervals[param, 1]
ci_upper <- model_results[[i]]$confidence_intervals[param, 2]
sprintf("%.2f [%.2f, %.2f]", eff, ci_lower, ci_upper)
}),
p = sapply(1:n_models, function(i) {
sprintf("%.3f", model_results[[i]]$p_values[param])
}),
AIC = round(sapply(model_results, function(x) x$aic), 1)
)
return(comparison_table)
}
## ESI Beta Regression Models
cat("=== ESI BETA REGRESSION MODELS ===\n\n")
## === ESI BETA REGRESSION MODELS ===
# Initialize list to store results
model_results_ESI_beta <- list()
# DEFINE THE ORIGINAL FORMULAS (these are needed for standardised_complete_cases)
formula_ESI_1 <- ESI_factor ~ ESI_norm
formula_ESI_2 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI
formula_ESI_3 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy)
formula_ESI_4 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy) + ESI_norm:stakeholderCategoryForESI
formula_ESI_5 <- ESI_factor ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy) + (0 + ESI_norm|caseStudy)
# Model 1: Simple model with only ESI_norm
cat("--- ESI Model 1 ---\n")
## --- ESI Model 1 ---
data_ESI_1 <- standardised_complete_cases(dF, formula_ESI_1)
# Convert factor to numeric properly: factor -> character -> numeric
data_ESI_1$ESI_factor_beta <- rescale_to_beta(as.numeric(as.character(data_ESI_1$ESI_factor)))
formula_ESI_1_beta <- ESI_factor_beta ~ ESI_norm
model_ESI_1_beta <- glmmTMB(
formula_ESI_1_beta,
data = data_ESI_1,
family = beta_family()
)
model_results_ESI_beta[[1]] <- summarise_beta_model(model_ESI_1_beta, "ESI_norm")
## Family: beta ( logit )
## Formula: ESI_factor_beta ~ ESI_norm
## Data: data_ESI_1
##
## AIC BIC logLik -2*log(L) df.resid
## -381.4 -370.1 193.7 -387.4 316
##
##
## Dispersion parameter for beta family (): 6.65
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.24059 0.04826 25.707 < 2e-16 ***
## ESI_norm 0.11352 0.04372 2.597 0.00941 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -381.4
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## ESI_norm: Exp(B) = 1.12, 95% CI [1.03, 1.22], p = 0.00941
# Model 2: Add stakeholder category
cat("--- ESI Model 2 ---\n")
## --- ESI Model 2 ---
data_ESI_2 <- standardised_complete_cases(dF, formula_ESI_2)
data_ESI_2$ESI_factor_beta <- rescale_to_beta(as.numeric(as.character(data_ESI_2$ESI_factor)))
formula_ESI_2_beta <- ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI
model_ESI_2_beta <- glmmTMB(
formula_ESI_2_beta,
data = data_ESI_2,
family = beta_family()
)
model_results_ESI_beta[[2]] <- summarise_beta_model(model_ESI_2_beta, "ESI_norm")
## Family: beta ( logit )
## Formula: ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI
## Data: data_ESI_2
##
## AIC BIC logLik -2*log(L) df.resid
## -382.6 -352.5 199.3 -398.6 311
##
##
## Dispersion parameter for beta family (): 6.9
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 1.12641 0.07576 14.867
## ESI_norm 0.11274 0.04481 2.516
## stakeholderCategoryForESICommercial 0.03524 0.18287 0.193
## stakeholderCategoryForESIEnvironment protectors 0.17124 0.15203 1.126
## stakeholderCategoryForESIGovernment 0.04312 0.17292 0.249
## stakeholderCategoryForESITourists 0.37986 0.12250 3.101
## stakeholderCategoryForESIWild resource users 0.04247 0.13333 0.319
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## ESI_norm 0.01188 *
## stakeholderCategoryForESICommercial 0.84717
## stakeholderCategoryForESIEnvironment protectors 0.26002
## stakeholderCategoryForESIGovernment 0.80310
## stakeholderCategoryForESITourists 0.00193 **
## stakeholderCategoryForESIWild resource users 0.75007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -382.6
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## ESI_norm: Exp(B) = 1.12, 95% CI [1.03, 1.22], p = 0.0119
# Model 3: Add random intercept for case study
cat("--- ESI Model 3 ---\n")
## --- ESI Model 3 ---
data_ESI_3 <- standardised_complete_cases(dF, formula_ESI_3)
data_ESI_3$ESI_factor_beta <- rescale_to_beta(as.numeric(as.character(data_ESI_3$ESI_factor)))
formula_ESI_3_beta <- ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy)
model_ESI_3_beta <- glmmTMB(
formula_ESI_3_beta,
data = data_ESI_3,
family = beta_family()
)
model_results_ESI_beta[[3]] <- summarise_beta_model(model_ESI_3_beta, "ESI_norm")
## Family: beta ( logit )
## Formula:
## ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI + (1 | caseStudy)
## Data: data_ESI_3
##
## AIC BIC logLik -2*log(L) df.resid
## -391.2 -357.3 204.6 -409.2 310
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.07947 0.2819
## Number of obs: 319, groups: caseStudy, 7
##
## Dispersion parameter for beta family (): 7.46
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 1.24639 0.15340 8.125
## ESI_norm 0.11759 0.04566 2.575
## stakeholderCategoryForESICommercial -0.15705 0.20949 -0.750
## stakeholderCategoryForESIEnvironment protectors 0.21684 0.18696 1.160
## stakeholderCategoryForESIGovernment -0.30196 0.25564 -1.181
## stakeholderCategoryForESITourists -0.05541 0.28373 -0.195
## stakeholderCategoryForESIWild resource users 0.14178 0.19050 0.744
## Pr(>|z|)
## (Intercept) 4.46e-16 ***
## ESI_norm 0.010 *
## stakeholderCategoryForESICommercial 0.453
## stakeholderCategoryForESIEnvironment protectors 0.246
## stakeholderCategoryForESIGovernment 0.238
## stakeholderCategoryForESITourists 0.845
## stakeholderCategoryForESIWild resource users 0.457
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -391.2
##
## Conditional R²: 0.776
## Marginal R²: 0.252
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## ESI_norm: Exp(B) = 1.12, 95% CI [1.03, 1.23], p = 0.01
# Model 4: Add interaction
cat("--- ESI Model 4 ---\n")
## --- ESI Model 4 ---
data_ESI_4 <- standardised_complete_cases(dF, formula_ESI_4)
data_ESI_4$ESI_factor_beta <- rescale_to_beta(as.numeric(as.character(data_ESI_4$ESI_factor)))
formula_ESI_4_beta <- ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy) +
ESI_norm:stakeholderCategoryForESI
model_ESI_4_beta <- glmmTMB(
formula_ESI_4_beta,
data = data_ESI_4,
family = beta_family()
)
model_results_ESI_beta[[4]] <- summarise_beta_model(model_ESI_4_beta, "ESI_norm")
## Family: beta ( logit )
## Formula:
## ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI + (1 |
## caseStudy) + ESI_norm:stakeholderCategoryForESI
## Data: data_ESI_4
##
## AIC BIC logLik -2*log(L) df.resid
## -384.9 -332.2 206.5 -412.9 305
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.06168 0.2484
## Number of obs: 319, groups: caseStudy, 7
##
## Dispersion parameter for beta family (): 7.52
##
## Conditional model:
## Estimate Std. Error
## (Intercept) 1.22119 0.14482
## ESI_norm 0.12845 0.07179
## stakeholderCategoryForESICommercial -0.18761 0.20207
## stakeholderCategoryForESIEnvironment protectors 0.17158 0.19322
## stakeholderCategoryForESIGovernment -0.19775 0.26798
## stakeholderCategoryForESITourists 0.05177 0.28721
## stakeholderCategoryForESIWild resource users 0.14749 0.19016
## ESI_norm:stakeholderCategoryForESICommercial -0.14678 0.15554
## ESI_norm:stakeholderCategoryForESIEnvironment protectors 0.09846 0.18875
## ESI_norm:stakeholderCategoryForESIGovernment 0.07499 0.16644
## ESI_norm:stakeholderCategoryForESITourists -0.09142 0.11693
## ESI_norm:stakeholderCategoryForESIWild resource users 0.12954 0.15316
## z value Pr(>|z|)
## (Intercept) 8.432 <2e-16 ***
## ESI_norm 1.789 0.0736 .
## stakeholderCategoryForESICommercial -0.928 0.3532
## stakeholderCategoryForESIEnvironment protectors 0.888 0.3745
## stakeholderCategoryForESIGovernment -0.738 0.4606
## stakeholderCategoryForESITourists 0.180 0.8570
## stakeholderCategoryForESIWild resource users 0.776 0.4380
## ESI_norm:stakeholderCategoryForESICommercial -0.944 0.3454
## ESI_norm:stakeholderCategoryForESIEnvironment protectors 0.522 0.6019
## ESI_norm:stakeholderCategoryForESIGovernment 0.451 0.6523
## ESI_norm:stakeholderCategoryForESITourists -0.782 0.4343
## ESI_norm:stakeholderCategoryForESIWild resource users 0.846 0.3977
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -384.9
##
## Conditional R²: 0.750
## Marginal R²: 0.291
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## ESI_norm: Exp(B) = 1.14, 95% CI [0.99, 1.31], p = 0.0736
# Model 5: Add random slope for ESI_norm
cat("--- ESI Model 5 ---\n")
## --- ESI Model 5 ---
data_ESI_5 <- standardised_complete_cases(dF, formula_ESI_5)
data_ESI_5$ESI_factor_beta <- rescale_to_beta(as.numeric(as.character(data_ESI_5$ESI_factor)))
formula_ESI_5_beta <- ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI + (1|caseStudy) +
(0 + ESI_norm|caseStudy)
model_ESI_5_beta <- glmmTMB(
formula_ESI_5_beta,
data = data_ESI_5,
family = beta_family()
)
model_results_ESI_beta[[5]] <- summarise_beta_model(model_ESI_5_beta, "ESI_norm")
## Family: beta ( logit )
## Formula:
## ESI_factor_beta ~ ESI_norm + stakeholderCategoryForESI + (1 |
## caseStudy) + (0 + ESI_norm | caseStudy)
## Data: data_ESI_5
##
## AIC BIC logLik -2*log(L) df.resid
## -390.0 -352.3 205.0 -410.0 309
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.08205 0.2864
## caseStudy.1 ESI_norm 0.01027 0.1014
## Number of obs: 319, groups: caseStudy, 7
##
## Dispersion parameter for beta family (): 7.58
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 1.24296 0.15416 8.063
## ESI_norm 0.12967 0.06296 2.060
## stakeholderCategoryForESICommercial -0.17848 0.20792 -0.858
## stakeholderCategoryForESIEnvironment protectors 0.18348 0.18895 0.971
## stakeholderCategoryForESIGovernment -0.29385 0.25063 -1.172
## stakeholderCategoryForESITourists -0.05208 0.27773 -0.188
## stakeholderCategoryForESIWild resource users 0.17577 0.19543 0.899
## Pr(>|z|)
## (Intercept) 7.47e-16 ***
## ESI_norm 0.0394 *
## stakeholderCategoryForESICommercial 0.3907
## stakeholderCategoryForESIEnvironment protectors 0.3315
## stakeholderCategoryForESIGovernment 0.2410
## stakeholderCategoryForESITourists 0.8512
## stakeholderCategoryForESIWild resource users 0.3684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -390.0
##
## Conditional R²: 0.787
## Marginal R²: 0.264
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## ESI_norm: Exp(B) = 1.14, 95% CI [1.01, 1.29], p = 0.0394
# Create comparison table
comparison_table_ESI_beta <- create_beta_comparison_table(model_results_ESI_beta, "ESI_norm")
## WTS Beta Regression Models
cat("\n=== WTS BETA REGRESSION MODELS ===\n\n")
##
## === WTS BETA REGRESSION MODELS ===
# Initialize list to store results
model_results_WTS_beta <- list()
# DEFINE THE ORIGINAL FORMULAS (these are needed for standardised_complete_cases)
formula_WTS_1 <- WTS_factor ~ WTS_norm
formula_WTS_2 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS
formula_WTS_3 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy)
formula_WTS_4 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy) + WTS_norm:stakeholderCategoryForWTS
formula_WTS_5 <- WTS_factor ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy) + (0 + WTS_norm|caseStudy)
formula_WTS_6 <- WTS_factor ~ WTS_norm + ESI + stakeholderCategoryForWTS + (1|caseStudy)
# Model 1: Simple model with only WTS_norm
cat("--- WTS Model 1 ---\n")
## --- WTS Model 1 ---
data_WTS_1 <- standardised_complete_cases(dF, formula_WTS_1)
data_WTS_1$WTS_factor_beta <- rescale_to_beta(as.numeric(as.character(data_WTS_1$WTS_factor)))
formula_WTS_1_beta <- WTS_factor_beta ~ WTS_norm
model_WTS_1_beta <- glmmTMB(
formula_WTS_1_beta,
data = data_WTS_1,
family = beta_family()
)
model_results_WTS_beta[[1]] <- summarise_beta_model(model_WTS_1_beta, "WTS_norm")
## Family: beta ( logit )
## Formula: WTS_factor_beta ~ WTS_norm
## Data: data_WTS_1
##
## AIC BIC logLik -2*log(L) df.resid
## -68.5 -58.4 37.3 -74.5 215
##
##
## Dispersion parameter for beta family (): 3.51
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.47824 0.06456 7.408 1.28e-13 ***
## WTS_norm 0.27560 0.06328 4.355 1.33e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -68.5
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## WTS_norm: Exp(B) = 1.32, 95% CI [1.16, 1.49], p = <0.001
# Model 2: Add stakeholder category
cat("--- WTS Model 2 ---\n")
## --- WTS Model 2 ---
data_WTS_2 <- standardised_complete_cases(dF, formula_WTS_2)
data_WTS_2$WTS_factor_beta <- rescale_to_beta(as.numeric(as.character(data_WTS_2$WTS_factor)))
formula_WTS_2_beta <- WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS
model_WTS_2_beta <- glmmTMB(
formula_WTS_2_beta,
data = data_WTS_2,
family = beta_family()
)
model_results_WTS_beta[[2]] <- summarise_beta_model(model_WTS_2_beta, "WTS_norm")
## Family: beta ( logit )
## Formula: WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS
## Data: data_WTS_2
##
## AIC BIC logLik -2*log(L) df.resid
## -62.9 -35.8 39.4 -78.9 210
##
##
## Dispersion parameter for beta family (): 3.58
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 0.60622 0.10347 5.859
## WTS_norm 0.28016 0.06656 4.209
## stakeholderCategoryForWTSCommercial -0.41013 0.26419 -1.552
## stakeholderCategoryForWTSEnvironment protectors -0.30562 0.20722 -1.475
## stakeholderCategoryForWTSGovernment -0.05146 0.21595 -0.238
## stakeholderCategoryForWTSLocal community -0.23296 0.25074 -0.929
## stakeholderCategoryForWTSWild resource users -0.15029 0.16635 -0.903
## Pr(>|z|)
## (Intercept) 4.66e-09 ***
## WTS_norm 2.57e-05 ***
## stakeholderCategoryForWTSCommercial 0.121
## stakeholderCategoryForWTSEnvironment protectors 0.140
## stakeholderCategoryForWTSGovernment 0.812
## stakeholderCategoryForWTSLocal community 0.353
## stakeholderCategoryForWTSWild resource users 0.366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -62.9
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## WTS_norm: Exp(B) = 1.32, 95% CI [1.16, 1.51], p = <0.001
# Model 3: Add random intercept for case study
cat("--- WTS Model 3 ---\n")
## --- WTS Model 3 ---
data_WTS_3 <- standardised_complete_cases(dF, formula_WTS_3)
data_WTS_3$WTS_factor_beta <- rescale_to_beta(as.numeric(as.character(data_WTS_3$WTS_factor)))
formula_WTS_3_beta <- WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy)
model_WTS_3_beta <- glmmTMB(
formula_WTS_3_beta,
data = data_WTS_3,
family = beta_family()
)
model_results_WTS_beta[[3]] <- summarise_beta_model(model_WTS_3_beta, "WTS_norm")
## Family: beta ( logit )
## Formula:
## WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS + (1 | caseStudy)
## Data: data_WTS_3
##
## AIC BIC logLik -2*log(L) df.resid
## -62.0 -31.5 40.0 -80.0 209
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.04548 0.2133
## Number of obs: 218, groups: caseStudy, 6
##
## Dispersion parameter for beta family (): 3.7
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 0.40165 0.22228 1.807
## WTS_norm 0.26787 0.06726 3.982
## stakeholderCategoryForWTSCommercial -0.20644 0.31492 -0.656
## stakeholderCategoryForWTSEnvironment protectors -0.01406 0.34066 -0.041
## stakeholderCategoryForWTSGovernment 0.08683 0.24412 0.356
## stakeholderCategoryForWTSLocal community -0.08503 0.35084 -0.242
## stakeholderCategoryForWTSWild resource users 0.09093 0.30517 0.298
## Pr(>|z|)
## (Intercept) 0.0708 .
## WTS_norm 6.82e-05 ***
## stakeholderCategoryForWTSCommercial 0.5121
## stakeholderCategoryForWTSEnvironment protectors 0.9671
## stakeholderCategoryForWTSGovernment 0.7221
## stakeholderCategoryForWTSLocal community 0.8085
## stakeholderCategoryForWTSWild resource users 0.7657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -62.0
##
## Conditional R²: 0.485
## Marginal R²: 0.314
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## WTS_norm: Exp(B) = 1.31, 95% CI [1.15, 1.49], p = <0.001
# Model 4: Add interaction
cat("--- WTS Model 4 ---\n")
## --- WTS Model 4 ---
data_WTS_4 <- standardised_complete_cases(dF, formula_WTS_4)
data_WTS_4$WTS_factor_beta <- rescale_to_beta(as.numeric(as.character(data_WTS_4$WTS_factor)))
formula_WTS_4_beta <- WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy) +
WTS_norm:stakeholderCategoryForWTS
model_WTS_4_beta <- glmmTMB(
formula_WTS_4_beta,
data = data_WTS_4,
family = beta_family()
)
model_results_WTS_beta[[4]] <- summarise_beta_model(model_WTS_4_beta, "WTS_norm")
## Family: beta ( logit )
## Formula:
## WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS + (1 |
## caseStudy) + WTS_norm:stakeholderCategoryForWTS
## Data: data_WTS_4
##
## AIC BIC logLik -2*log(L) df.resid
## -57.5 -10.2 42.8 -85.5 204
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 0.06106 0.2471
## Number of obs: 218, groups: caseStudy, 6
##
## Dispersion parameter for beta family (): 3.82
##
## Conditional model:
## Estimate Std. Error
## (Intercept) 0.371515 0.231895
## WTS_norm 0.241155 0.106784
## stakeholderCategoryForWTSCommercial -0.339635 0.447828
## stakeholderCategoryForWTSEnvironment protectors -0.119542 0.364173
## stakeholderCategoryForWTSGovernment 0.109098 0.241015
## stakeholderCategoryForWTSLocal community -0.007051 0.381881
## stakeholderCategoryForWTSWild resource users 0.160250 0.336558
## WTS_norm:stakeholderCategoryForWTSCommercial -0.177535 0.381936
## WTS_norm:stakeholderCategoryForWTSEnvironment protectors 0.303535 0.196651
## WTS_norm:stakeholderCategoryForWTSGovernment -0.021695 0.198483
## WTS_norm:stakeholderCategoryForWTSLocal community 0.271646 0.311653
## WTS_norm:stakeholderCategoryForWTSWild resource users -0.151551 0.180839
## z value Pr(>|z|)
## (Intercept) 1.602 0.1091
## WTS_norm 2.258 0.0239 *
## stakeholderCategoryForWTSCommercial -0.758 0.4482
## stakeholderCategoryForWTSEnvironment protectors -0.328 0.7427
## stakeholderCategoryForWTSGovernment 0.453 0.6508
## stakeholderCategoryForWTSLocal community -0.018 0.9853
## stakeholderCategoryForWTSWild resource users 0.476 0.6340
## WTS_norm:stakeholderCategoryForWTSCommercial -0.465 0.6421
## WTS_norm:stakeholderCategoryForWTSEnvironment protectors 1.544 0.1227
## WTS_norm:stakeholderCategoryForWTSGovernment -0.109 0.9130
## WTS_norm:stakeholderCategoryForWTSLocal community 0.872 0.3834
## WTS_norm:stakeholderCategoryForWTSWild resource users -0.838 0.4020
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -57.5
##
## Conditional R²: 0.561
## Marginal R²: 0.361
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## WTS_norm: Exp(B) = 1.27, 95% CI [1.03, 1.57], p = 0.0239
# Model 5: Add random slope for WTS_norm
cat("--- WTS Model 5 ---\n")
## --- WTS Model 5 ---
data_WTS_5 <- standardised_complete_cases(dF, formula_WTS_5)
data_WTS_5$WTS_factor_beta <- rescale_to_beta(as.numeric(as.character(data_WTS_5$WTS_factor)))
formula_WTS_5_beta <- WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS + (1|caseStudy) +
(0 + WTS_norm|caseStudy)
model_WTS_5_beta <- glmmTMB(
formula_WTS_5_beta,
data = data_WTS_5,
family = beta_family()
)
model_results_WTS_beta[[5]] <- summarise_beta_model(model_WTS_5_beta, "WTS_norm")
## Family: beta ( logit )
## Formula:
## WTS_factor_beta ~ WTS_norm + stakeholderCategoryForWTS + (1 |
## caseStudy) + (0 + WTS_norm | caseStudy)
## Data: data_WTS_5
##
## AIC BIC logLik -2*log(L) df.resid
## -60.0 -26.2 40.0 -80.0 208
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 4.548e-02 2.133e-01
## caseStudy.1 WTS_norm 7.094e-11 8.423e-06
## Number of obs: 218, groups: caseStudy, 6
##
## Dispersion parameter for beta family (): 3.7
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 0.40165 0.22228 1.807
## WTS_norm 0.26787 0.06726 3.982
## stakeholderCategoryForWTSCommercial -0.20644 0.31492 -0.656
## stakeholderCategoryForWTSEnvironment protectors -0.01406 0.34066 -0.041
## stakeholderCategoryForWTSGovernment 0.08683 0.24412 0.356
## stakeholderCategoryForWTSLocal community -0.08503 0.35084 -0.242
## stakeholderCategoryForWTSWild resource users 0.09093 0.30517 0.298
## Pr(>|z|)
## (Intercept) 0.0708 .
## WTS_norm 6.82e-05 ***
## stakeholderCategoryForWTSCommercial 0.5121
## stakeholderCategoryForWTSEnvironment protectors 0.9671
## stakeholderCategoryForWTSGovernment 0.7221
## stakeholderCategoryForWTSLocal community 0.8085
## stakeholderCategoryForWTSWild resource users 0.7657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -60.0
## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Decrease the `tolerance` level to force the calculation of random effect
## variances, or impose priors on your random effects parameters (using
## packages like `brms` or `glmmTMB`).
## Random effect variances not available. Returned R2 does not account for random effects.
## Conditional R²: NA
## Marginal R²: 0.385
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## WTS_norm: Exp(B) = 1.31, 95% CI [1.15, 1.49], p = <0.001
# Create comparison table
comparison_table_WTS_beta <- create_beta_comparison_table(model_results_WTS_beta, "WTS_norm")
# Model 6: WTS with ESI included (matching your final model)
cat("--- WTS Model 6 (with ESI) ---\n")
## --- WTS Model 6 (with ESI) ---
data_WTS_6 <- standardised_complete_cases(dF, formula_WTS_6)
data_WTS_6$WTS_factor_beta <- rescale_to_beta(as.numeric(as.character(data_WTS_6$WTS_factor)))
formula_WTS_6_beta <- WTS_factor_beta ~ WTS_norm + ESI + stakeholderCategoryForWTS + (1|caseStudy)
model_WTS_6_beta <- glmmTMB(
formula_WTS_6_beta,
data = data_WTS_6,
family = beta_family()
)
summarise_beta_model(model_WTS_6_beta, c("WTS_norm", "ESI"))
## Family: beta ( logit )
## Formula:
## WTS_factor_beta ~ WTS_norm + ESI + stakeholderCategoryForWTS +
## (1 | caseStudy)
## Data: data_WTS_6
##
## AIC BIC logLik -2*log(L) df.resid
## -80.0 -46.3 50.0 -100.0 205
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## caseStudy (Intercept) 1.493e-10 1.222e-05
## Number of obs: 215, groups: caseStudy, 6
##
## Dispersion parameter for beta family (): 3.98
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 0.57750 0.10002 5.774
## WTS_norm 0.24883 0.06540 3.805
## ESI 0.29446 0.06309 4.667
## stakeholderCategoryForWTSCommercial -0.48632 0.26172 -1.858
## stakeholderCategoryForWTSEnvironment protectors -0.25337 0.20164 -1.257
## stakeholderCategoryForWTSGovernment -0.00707 0.20850 -0.034
## stakeholderCategoryForWTSLocal community 0.09920 0.25765 0.385
## stakeholderCategoryForWTSWild resource users -0.06577 0.16257 -0.405
## Pr(>|z|)
## (Intercept) 7.74e-09 ***
## WTS_norm 0.000142 ***
## ESI 3.06e-06 ***
## stakeholderCategoryForWTSCommercial 0.063145 .
## stakeholderCategoryForWTSEnvironment protectors 0.208906
## stakeholderCategoryForWTSGovernment 0.972951
## stakeholderCategoryForWTSLocal community 0.700220
## stakeholderCategoryForWTSWild resource users 0.685808
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## AIC: -80.0
## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Decrease the `tolerance` level to force the calculation of random effect
## variances, or impose priors on your random effects parameters (using
## packages like `brms` or `glmmTMB`).
## Random effect variances not available. Returned R2 does not account for random effects.
## Conditional R²: NA
## Marginal R²: 0.579
##
## Exponentiated Coefficients (multiplicative effects) and 95% Confidence Intervals:
## WTS_norm: Exp(B) = 1.28, 95% CI [1.13, 1.46], p = <0.001
## ESI: Exp(B) = 1.34, 95% CI [1.19, 1.52], p = <0.001
# Print comparison tables
{
cat("\n=== COMPARISON TABLES ===\n\n")
cat("ESI Models:\n")
print(comparison_table_ESI_beta)
cat("\n")
cat("WTS Models:\n")
print(comparison_table_WTS_beta)
}
##
## === COMPARISON TABLES ===
##
## ESI Models:
## Model Effect_CI p AIC
## 1 Model 1 1.12 [1.03, 1.22] 0.009 -381.4
## 2 Model 2 1.12 [1.03, 1.22] 0.012 -382.6
## 3 Model 3 1.12 [1.03, 1.23] 0.010 -391.2
## 4 Model 4 1.14 [0.99, 1.31] 0.074 -384.9
## 5 Model 5 1.14 [1.01, 1.29] 0.039 -390.0
##
## WTS Models:
## Model Effect_CI p AIC
## 1 Model 1 1.32 [1.16, 1.49] 0.000 -68.5
## 2 Model 2 1.32 [1.16, 1.51] 0.000 -62.9
## 3 Model 3 1.31 [1.15, 1.49] 0.000 -62.0
## 4 Model 4 1.27 [1.03, 1.57] 0.024 -57.5
## 5 Model 5 1.31 [1.15, 1.49] 0.000 -60.0