## Testing an intervention co-designed with stakeholders for altering wildlife consumption: Health messaging matters 

## data sets
health <-read.csv("D:/work-sysu/papers/6 (Revision) My_Amy_BiM_intervention/survey/analysis/health.csv",head=T,sep=",",fill=TRUE, row.names= 1)
identi <-read.csv("D:/work-sysu/papers/6 (Revision) My_Amy_BiM_intervention/survey/analysis/identi.csv",head=T,sep=",",fill=TRUE, row.names= 1)
tip_identi <-read.csv("D:/work-sysu/papers/6 (Revision) My_Amy_BiM_intervention/survey/analysis/tip_identi.csv",head=T,sep=",",fill=TRUE, row.names= 1)

library(lme4)
library(Matrix)
library(rstatix)
library(ggpubr)

## set the variables
health$gender<-as.factor(health$gender)
health$education<-as.factor(health$education)
health$age<-as.factor(health$age)
health$consumer<-as.factor(health$consumer)

health$age <- relevel(health$age, "26~35")
health$education <- relevel(health$education, "undergraduate")
health$consumer <- relevel(health$consumer, "control")

############## Mann-Whitney U-Test
### health
wilcox.test(risk_health2 ~ consumer,data=health)
wilcox.test(risk_health4 ~ consumer,data=health)
wilcox.test(risk_health5 ~ consumer,data=health)
wilcox.test(risk_health7_rev ~ consumer,data=health)
wilcox.test(intention.use.wildanimal ~ consumer,data=health)

### Legality
wilcox.test(efficiency_post ~ consumer,data=identi)
wilcox.test(capacity_post ~ consumer,data=identi)
wilcox.test(risk_law_illegal ~ consumer,data=identi)
wilcox.test(risk_law_legal ~ consumer,data=identi)
wilcox.test(risk_punishment ~ consumer,data=identi)
wilcox.test(intention.use.wildanimal ~ consumer,data=identi)

#### tip_legality
wilcox.test(efficiency_post ~ public,data=tip_identi)
wilcox.test(capacity_post ~ public,data=tip_identi)
wilcox.test(risk_law_illegal ~ public,data=tip_identi)
wilcox.test(risk_law_legal ~ public,data=tip_identi)
wilcox.test(risk_punishment ~ public,data=tip_identi)
wilcox.test(intention.tipoff ~ public,data=tip_identi)

############ effect size 
## r (rank biserial correlation)
## 0.30 (small effect), 0.30 - < 0.50 (moderate effect) and >= 0.50 (large effect).
library(rcompanion)
## health
wilcoxonR(x = health$risk_health2, g = health$consumer)
wilcoxonR(x = health$risk_health4, g = health$consumer)
wilcoxonR(x = health$risk_health5, g = health$consumer)
wilcoxonR(x = health$risk_health7_rev, g = health$consumer)
wilcoxonR(x = health$intention.use.wildanimal, g = health$consumer)

## legality
wilcoxonR(x = identi$efficiency_post, g = identi$consumer)
wilcoxonR(x = identi$capacity_post, g = identi$consumer)
wilcoxonR(x = identi$risk_law_illegal, g = identi$consumer)
wilcoxonR(x = identi$risk_law_legal, g = identi$consumer)
wilcoxonR(x = identi$risk_punishment, g = identi$consumer)
wilcoxonR(x = identi$intention.use.wildanimal, g = identi$consumer)

## tip_legality
wilcoxonR(x = tip_identi$efficiency_post, g = tip_identi$public)
wilcoxonR(x = tip_identi$capacity_post, g = tip_identi$public)
wilcoxonR(x = tip_identi$risk_law_illegal, g = tip_identi$public)
wilcoxonR(x = tip_identi$risk_law_legal, g = tip_identi$public)
wilcoxonR(x = tip_identi$risk_punishment, g = tip_identi$public)
wilcoxonR(x = tip_identi$intention.tipoff, g = tip_identi$public)

############# ordered logistic regression 
########## health
health$intention.use.wildanimal <- as.factor(health$intention.use.wildanimal)

install.packages("ordinal")
install.packages("rcompanion")
library(ordinal)
library(rcompanion)

### y = intention.use.wildanimal
### model 1
lm_health1 <- clm(as.factor(intention.use.wildanimal) ~ consumer,
                  data = health)
summary(lm_health1)
nagelkerke(lm_health1)
### model 2
lm_health2 <- clm(as.factor(intention.use.wildanimal) ~ consumer +
                    gender + age + education + 
                    field_ecoen + scale(v_use) + scale(v_protect)+
                    used.bearbile,
                  data = health)
summary(lm_health2)
nagelkerke(lm_health2)
head(health)
## risk_health2=risk_health2_wild; risk_health4=risk_health4_source;
## risk_health5=risk_health5_doctor; risk_improper use=risk_health7_rev
### model 3
lm_health3 <- clm(as.factor(intention.use.wildanimal) ~ consumer +
                    gender + age + education + 
                    field_ecoen + scale(v_use) + scale(v_protect)+
                    used.bearbile +
                    scale(risk_health2) + scale(risk_health4) + 
                    scale(risk_health5) + scale(risk_health7_rev),
                  data = health)
summary(lm_health3)
nagelkerke(lm_health3)

### y = risk_health2_wild
### model 1
lm_hrw1 <- clm(as.factor(risk_health2) ~ consumer,
               data = health)
summary(lm_hrw1)
nagelkerke(lm_hrw1)
### model 2
lm_hrw2 <- clm(as.factor(risk_health2) ~ consumer +
                 gender + age +
                 #education +
                 field_ecoen + scale(v_use) + scale(v_protect)+
                 used.bearbile,
               data = health)
summary(lm_hrw2)
nagelkerke(lm_hrw2)

#### y = risk_health4_source (unknown source)
#### model 1
lm_hri1 <- clm(as.factor(risk_health4) ~ consumer,
               data = health)
summary(lm_hri1)
nagelkerke(lm_hri1)
### model 2
lm_hri2 <- clm(as.factor(risk_health4) ~ consumer +
                 gender + age + education + 
                 field_ecoen + scale(v_use) + scale(v_protect)+
                 used.bearbile,
               data = health)
summary(lm_hri2)
nagelkerke(lm_hri2)

#### y = risk_health5_doctor (safety)
#### model 1
lm_hrd1 <- clm(as.factor(risk_health5) ~ consumer,
               data = health)
summary(lm_hrd1)
nagelkerke(lm_hrd1)
### model 2
lm_hrd2 <- clm(as.factor(risk_health5) ~ consumer +
                 gender + age + education + 
                 field_ecoen + scale(v_use) + scale(v_protect)+
                 used.bearbile,
               data = health)
summary(lm_hrd2)
nagelkerke(lm_hrd2)

#### y = risk_health7_rev  (improper use)
#### model 1
lm_hru1 <- clm(as.factor(risk_health7_rev)  ~ consumer,
               data = health)
summary(lm_hru1)
nagelkerke(lm_hru1)
### model 2
lm_hru2 <- clm(as.factor(risk_health7_rev) ~ consumer +
                 gender + age + education + 
                 field_ecoen + scale(v_use) + scale(v_protect)+
                 used.bearbile,
               data = health)
summary(lm_hru2)
nagelkerke(lm_hru2)

############ legality 
####  y = intention.use.wildanimal
####### model 1
lm_identi1 <- clm(as.factor(intention.use.wildanimal)  ~ consumer,
                  data = identi)
summary(lm_identi1)
nagelkerke(lm_identi1)
####### model 2
lm_identi2 <- clm(as.factor(intention.use.wildanimal) ~ consumer +
                    gender + age + education + 
                    field_ecoen + scale(v_use) + scale(v_protect)+
                    used.bearbile,
                  data = identi)
summary(lm_identi2)
nagelkerke(lm_identi2)
####### model 3
lm_identi3 <- clm(as.factor(intention.use.wildanimal) ~ consumer +
                    gender + age + education + 
                    field_ecoen + scale(v_use) + scale(v_protect)+
                    used.bearbile +
                    scale(efficiency_post) + scale(capacity_post) + 
                    scale(risk_law_illegal) + scale(risk_law_legal) +scale(risk_punishment),
                  data = identi)
summary(lm_identi3)
nagelkerke(lm_identi3)

####  y = efficiency_post
####### model 1
lm_legeff1 <- clm(as.factor(efficiency_post)  ~ consumer,
                  data = identi)
summary(lm_legeff1)
nagelkerke(lm_legeff1)
####### model 2
lm_legeff2 <- clm(as.factor(efficiency_post) ~ consumer +
                    gender + age + education + 
                    field_ecoen + scale(v_use) + scale(v_protect)+
                    used.bearbile,
                  data = identi)
summary(lm_legeff2)
nagelkerke(lm_legeff2)

####  y = capacity_post
####### model 1
lm_legcap1 <- clm(as.factor(capacity_post)  ~ consumer,
                  data = identi)
summary(lm_legcap1)
nagelkerke(lm_legcap1)
####### model 2
lm_legcap2 <- clm(as.factor(capacity_post) ~ consumer +
                    gender + age + education + 
                    field_ecoen + scale(v_use) + scale(v_protect)+
                    used.bearbile,
                  data = identi)
summary(lm_legcap2)
nagelkerke(lm_legcap2)

####  y = risk_law_illegal
####### model 1
lm_legill1 <- clm(as.factor(risk_law_illegal)  ~ consumer,
                  data = identi)
summary(lm_legill1)
nagelkerke(lm_legill1)
####### model 2
lm_legill2 <- clm(as.factor(risk_law_illegal) ~ consumer +
                    gender + age + education + 
                    field_ecoen + scale(v_use) + scale(v_protect)+
                    used.bearbile,
                  data = identi)
summary(lm_legill2)
nagelkerke(lm_legill2)

#####  y = risk_law_legal
####### model 1
lm_legle1 <- clm(as.factor(risk_law_legal)  ~ consumer,
                 data = identi)
summary(lm_legle1)
nagelkerke(lm_legle1)
####### model 2
lm_legle2 <- clm(as.factor(risk_law_legal) ~ consumer +
                   gender + age + education + 
                   field_ecoen + scale(v_use) + scale(v_protect)+
                   used.bearbile,
                 data = identi)
summary(lm_legle2)
nagelkerke(lm_legle2)

####  y = risk_punishment
####### model 1
lm_legpu1 <- clm(as.factor(risk_punishment)  ~ consumer,
                 data = identi)
summary(lm_legpu1)
nagelkerke(lm_legpu1)
####### model 2
lm_legpu2 <- clm(as.factor(risk_punishment) ~ consumer +
                   gender + age + education + 
                   field_ecoen + scale(v_use) + scale(v_protect)+
                   used.bearbile,
                 data = identi)
summary(lm_legpu2)
nagelkerke(lm_legpu2)


############## mediation analyses
## using the PROCESS command package (version 4.1) in SPSS 25.0 (Hayes 2022).