#Loading packages
#---------------------------------------------------------
{
require(ggplot2)
require(dplyr)
require(reshape2)
require(reshape)
require(lubridate)
require(plyr)
require(vegan)
require(psy)
require(psych)
require(lavaan)
require(fitdistrplus)
require(ltm)
require(tidySEM)
}

#Loading the data
#---------------------------------------------------------
df <- read.csv(file="~/ResultsSurvey.csv")
head(df)

#Renaming columns
#---------------------------------------------------------
df <- df %>% dplyr::rename(., c(mammals = job.iasmammals, birds = job.iasbirds, 
                              amphibians = job.iasamphibians, reptiles = job.iasreptiles,
                              freshwater.fish = job.iasfreshfish, freshwater.invertebrates = job.iasfreshinv,
                              terrestrial.invertebrates = job.iasterrinv, marine.organisms = job.iasmarine, 
                              plants = job.iasplants))

#Rescaling bipolar scales
#---------------------------------------------------------
idx <- c("ias.natconserv", "ias.costbenefit", "ias.nextgen", "ias.nopriority", "ias.collateral", 
         "ias.interfere", "ias.notenough", "nature.humwellbeing", "nature.plantanimals", "nature.priorhum",
         "nature.conssacrifice")
df[,which(names(df)%in%idx)] <- as.data.frame(lapply(df[,which(names(df)%in%idx)],
                                                     function(x)as.numeric(factor(x, levels=c("Fortemente in disaccordo", "In disaccordo", "Parzialmente in disaccordo",             
                                                                                              "Incerto/a", "Parzialmente d'accordo", "D'accordo", "Completamente d'accordo")))-4))

idx <- c("mammals", "birds", "amphibians", "reptiles", "freshwater.fish", 
         "freshwater.invertebrates", "terrestrial.invertebrates", "marine.organisms", "plants")
df[,which(names(df)%in%idx)] <- as.data.frame(lapply(df[,which(names(df)%in%idx)],
                                                     function(x)as.numeric(factor(x, levels=c("No", "Si")))-1))

df$practitioner <- apply(df[, c("practitioner.as", "practitioner.em")], 1, sum)
table(df$practitioner)

df <- dplyr::mutate(df, practitioner.dummy = ifelse(practitioner == 0, 0, 1))
df$practitioner.dummy <- factor(df$practitioner.dummy, levels = c(0,1), labels = c("Researcher", "Practitioner"))

#Characterizing background
#---------------------------------------------------------
df[which(df$mammals==1|df$birds==1|df$amphibians==1|df$reptiles==1),"terrestrial.vertebrates"] <- 1
df[which(is.na(df$terrestrial.vertebrates)),"terrestrial.vertebrates"] <- 0
table(df$terrestrial.vertebrates)

idx <- c("terrestrial.vertebrates", "freshwater.fish", "freshwater.invertebrates", 
         "terrestrial.invertebrates", "marine.organisms", "plants")
df$taxa <- apply(df[,which(names(df)%in%idx)], 1, sum)
d <- df[-which(df$taxa==0),]

head(d)
apply(d %>% dplyr::select(c(mammals, birds, amphibians, reptiles, freshwater.fish, 
                      freshwater.invertebrates, terrestrial.invertebrates, 
                      marine.organisms, plants)), 2, sum)/nrow(d)

#Data exploration
#---------------------------------------------------------
table(d$taxa)
d[which(d$taxa ==1),]

d %>% dplyr::filter(taxa == 1) %>% 
  dplyr::select(terrestrial.vertebrates, freshwater.fish, 
                freshwater.invertebrates, terrestrial.invertebrates, marine.organisms, plants)

d[which(d$taxa ==1), "taxa.type"] <- idx[apply(df %>% dplyr::filter(taxa == 1) %>% 
                                                 dplyr::select(terrestrial.vertebrates, freshwater.fish, 
                                                               freshwater.invertebrates, terrestrial.invertebrates, marine.organisms, plants)
                                               , 1, \(x) which(x > 0)[1])]
d[which(d$taxa !=1), "taxa.type"] <- "multiple"
table(d$taxa.type) #Frew freshwater experts
d$taxa.type <- plyr::revalue(d$taxa.type, c("freshwater.fish" = "freshwater + marine", 
                                            "freshwater.invertebrates" = "freshwater + marine",
                                            "marine.organisms" = "freshwater + marine"))
table(d$taxa.type)

table(d$taxa.type)
prop.table(table(d$practitioner.dummy))
d[which(d$practitioner.dummy == "Practitioner"), "name"]

################################################################################
#Modelling attitudes towards IAS
################################################################################

#Exploring the distribution of answers (attitudes towards IAS)
#---------------------------------------------------------
idx <- c("ias.natconserv", "ias.costbenefit", "ias.nextgen", "ias.nopriority", "ias.collateral", "ias.interfere",
         "ias.notenough")

question.idx <- c("Le specie aliene invasive sono una minaccia rilevante per la conservazione della natura",
                  "I benefici derivanti dal controllo delle specie aliene invasive superano i danni che questo può comportare per le specie autoctone",
                  "Un investimento serio nel controllo delle specie aliene invasive andrebbe a vantaggio delle prossime generazioni",
                  "Il controllo delle specie aliene invasive è meno urgente di altri problemi legati alla conservazione della natura (rev)",
                  "Il controllo delle specie aliene invasive può avere effetti collaterali imprevedibili  (rev)",
                  "Il controllo delle specie aliene invasive interferisce con la natura (rev)",
                  "Fino ad oggi,  non è stato fatto abbastanza per controllare le specie aliene invasive")
for(i in seq_along(idx)){
  print(barplot(table(d[,c(idx[i])]), main = question.idx[i]))
}

#Reversing negative questions (only for Cronbach's alpha, otherwise all loadings are positive)
#---------------------------------------------------------
#d$ias.nopriority <- d$ias.nopriority*-1
#d$ias.collateral <- d$ias.collateral*-1
#d$ias.interfere <- d$ias.interfere*-1

#Cronbach's alpha (> 0.7, not bad)
#---------------------------------------------------------
for(i in seq_along(idx)){
  print(barplot(table(d[,c(idx[i])]), main = question.idx[i]))
}
cronbach.alpha(d[, c("ias.natconserv", "ias.costbenefit", "ias.nextgen", "ias.nopriority", "ias.collateral", "ias.interfere",
                     "ias.notenough")])

#A basic confirmatory factor analysis with a robust estimator (Satorra/Bentler correction)
#---------------------------------------------------------
mod.a <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
               ias.collateral + ias.interfere + ias.notenough'
fit.a <- cfa(mod.a, data=d, std.ov= TRUE, std.lv=TRUE, estimator="MLM")

summary(fit.a, fit.measures=TRUE)
fitmeasures(fit.a, fit.measures = c("fmin", "chisq", "df", "pvalue", "cfi", "tli", "srmr", "rmsea"))

fitmeasures(fit.a, "cfi") #Should be higher than 0.95
fitmeasures(fit.a, "tli") #Should be higher than 0.95
fitmeasures(fit.a, "srmr") #Should be smaller than 0.08
fitmeasures(fit.a, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.a, "rmsea.ci.lower") #Its CI should be within 0.6-0.8
fitmeasures(fit.a, "rmsea.ci.upper") #Its CI should be within 0.6-0.8
fit.a.cov <- fitted(fit.a)$cov #model-implied covariances
fit.a.cor <- cov2cor(fit.a.cov) # transform model-implied covariances to correlations
fit.a.cor # original correlations
residuals(fit.a, type = "cor") #ias.interfere and ias.collateral covary

#Allowing ias.collateral and ias.interfere to be correlated
#---------------------------------------------------------
mod.b <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
          #############################################
          ias.collateral ~~ ias.interfere'
fit.b <- cfa(mod.b, data=d, std.ov= TRUE, std.lv=TRUE, estimator="MLM") #Satorra/Bentler correction

fitmeasures(fit.b, fit.measures = c("fmin", "chisq", "df", "pvalue"))
fitmeasures(fit.b, "cfi") #Should be higher than 0.95
fitmeasures(fit.b, "tli") #Should be higher than 0.95
fitmeasures(fit.b, "srmr") #Should be smaller than 0.08
fitmeasures(fit.b, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.b, "rmsea.ci.lower") #Its CI should be within 0.6-0.8
fitmeasures(fit.b, "rmsea.ci.upper") #Its CI should be within 0.6-0.8
fit.b.cov <- fitted(fit.b)$cov #model-implied covariances
fit.b.cor <- cov2cor(fit.b.cov) # transform model-implied covariances to correlations
fit.b.cor # original correlations
residuals(fit.b, type="cor") #ias.interfere and ias.collateral covary
str(fitted(fit.b))

#Comparing the difference in the T statistics’ values from the two models, which
#follows a chi-squared distribution (lavaan accounts for correction of robust estimators)
#---------------------------------------------------------
lavTestLRT(fit.a, fit.b) #Significant differences

################################################################################
#Modelling attitudes towards IAS: n.taxa as a group
################################################################################

#Measurement invariance: the multi-group factor model with constraints exhibits
#an adequate fit to the data according to the AFIs; and (b) the change in the CFI
#and MFI values from one model to a more restrictive (nested) model is negligible.
#---------------------------------------------------------
table(d$taxa.type)
d$taxa.type <- as.character(d$taxa.type)
d$taxa.type <- plyr::revalue(d$taxa.type,
                              c("freshwater + marine" = "single", "multiple" = "multiple", 
                                "plants" = "single", "terrestrial.invertebrates" = "single",
                                terrestrial.vertebrates = "single"))

#Configural invariance (LVMs have the same structure in the group: same LVs,
#same indicator variables, same pattern of constrained and estimated parameters.
#Estimated coefficients might differ between groups, and thus LVs might not measure
#the same construct in each group)
#---------------------------------------------------------
mod <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
                      #---------------------------------------------------------
                 ias.collateral ~~ ias.interfere'

fit.conf.inv <- cfa(mod, data = d, std.ov = TRUE, std.lv=TRUE, estimator="MLM",
                    group = "taxa.type")

#Weak invariance (loadings are the same across groups. Latent variances are allowed
#to vary among groups: weak argument that LVs are equivalent between groups, but
#only latent variances and covariances could be compared)
#---------------------------------------------------------
fit.weak.inv <- cfa(mod, data = d, std.ov = TRUE, std.lv=TRUE, estimator="MLM",
                    group = "taxa.type", group.equal = "loadings")

#Strong invariance (For a given indicator variable, intercepts are the same among groups
#LV means are the same among groups. Individuals at the same LV score have the same
#expected values of indicator variables. LV means, variances, and covariances can be 
#compared. Group differences on strongly-invariant LVs, reflect real differences in LVs, 
#as opposed to the difference caused by differente measures)
#---------------------------------------------------------
fit.strong.inv <- cfa(mod, data = d, std.ov = TRUE, std.lv=TRUE, estimator="MLM",
                    group = "taxa.type", group.equal = c("loadings", "intercepts"))

#Testing changes (models behave similarly, strong invariance achieved)
#---------------------------------------------------------
lavTestLRT(fit.conf.inv, fit.weak.inv, fit.strong.inv)
fitmeasures(fit.conf.inv, c("cfi", "mfi", "rmsea", "srmr"))
fitmeasures(fit.weak.inv, c("cfi", "mfi", "rmsea", "srmr"))
fitmeasures(fit.strong.inv, c("cfi", "mfi", "rmsea", "srmr"))

#Comparing means and loadings between groups
#---------------------------------------------------------
bestmod <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
          #############################################
          ias.collateral ~~ ias.interfere'
bestfit <- cfa(mod.b, data=d, std.ov= TRUE, std.lv=TRUE, estimator="MLM", group = "taxa.type") #Satorra/Bentler correction
fitmeasures(bestfit, fit.measures = c("fmin", "chisq", "df", "pvalue"))
fitmeasures(bestfit, "cfi") #Should be higher than 0.95
fitmeasures(bestfit, "tli") #Should be higher than 0.95
fitmeasures(bestfit, "srmr") #Should be smaller than 0.08
fitmeasures(bestfit, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(bestfit, "rmsea.ci.lower") #Its CI should be within 0.6-0.8
fitmeasures(bestfit, "rmsea.ci.upper") #Its CI should be within 0.6-0.8
residuals(bestfit, type="cor") #ias.interfere and ias.collateral covary
str(predict(bestfit))

#Plot scores and loadings
#---------------------------------------------------------
d[which(d$taxa.type == "single"), "score"] <- as.numeric(predict(bestfit)$single)
d[which(d$taxa.type == "multiple"), "score"] <- as.numeric(predict(bestfit)$multiple)
d$score <- as.numeric(d$score)
foo <- d
foo$taxa.type <- factor(foo$taxa.type, levels = c("single", "multiple"), labels = c("Homogeneous", "Heterogeneous"))
ggplot() + 
  geom_violin(data = foo, aes(x = taxa.type, y = score, fill = taxa.type)) +
  geom_boxplot(data = foo, aes(x = taxa.type, y = score), width = 0.1, color="black", alpha=0.6) +
  scale_fill_manual(values = c("steelblue", "gold")) +
  theme_bw() + labs(y = "PMA scores (Z-scores)", x = "IAS taxa", fill = "Professional background") +
  guides(fill = FALSE) +
  theme(
    axis.text.x = element_text(colour="black",size=11,angle=0,hjust=.5,vjust=.5,face="plain"),
    axis.text.y = element_text(colour="black",size=11,angle=0,hjust=1,vjust=.5,face="plain"),
    panel.background = element_rect(colour = "black", size = 0.4, linetype = "solid"))

inspect(bestfit, what = "std")$single$lambda
inspect(bestfit, what = "std")$multiple$lambda
summary(bestfit)

item <- c("ias.natconserv", "ias.costbenefit", "ias.nextgen", 
          "ias.nopriority", "ias.collateral", "ias.interfere", 
          "ias.notenough")
mean <- c(0.688, 0.497, 0.890, -0.490, -0.372, -0.437, 0.426,
           0.719, 0.408, 0.759, -0.586, -0.204, -0.401, 0.415)
se <- c(0.128, 0.074, 0.109, 0.085, 0.078, 0.089, 0.095,
        0.115, 0.076, 0.104, 0.070, 0.088, 0.068, 0.107)
cv <- se/mean
foo <- data.frame(item, mean, se, cv)
foo$item <- factor(foo$item, levels = c("ias.natconserv", "ias.costbenefit", "ias.nextgen", 
                                        "ias.nopriority", "ias.collateral", "ias.interfere", 
                                        "ias.notenough"))
foo$group <- c(rep("single", 7), rep("multiple", 7))
foo$group <- factor(foo$group, levels = c("single", "multiple"), labels = c("Single group", "Multiple groups"))

ggplot() + 
  geom_errorbar(data = foo, aes(x = item, ymin = mean - se, ymax = mean + se, color = group), 
                position = position_dodge(width = 0.5), width = 0.1) +
  geom_point(data = foo, aes(x = item, y = mean, fill = group), 
             shape = 21, size = 5, color = "black",  position = position_dodge(width = 0.5)) +
  geom_hline(aes(yintercept = 0), linetype = "dashed", size = 0.3) +
  scale_color_manual(values = c("black", "black")) + 
  #scale_fill_manual(values = c("#FF9966", "#66CC33")) +
  scale_fill_manual(values = c("steelblue", "gold")) +
  scale_x_discrete(labels = c("Nature \nconservation", "Costs and\nbenefits", "Next \ngenerations",
                            "Not a \npriority", "Collateral \ndamages", "Interference \nwith nature", "Not enough \ncontrol")) +
  scale_y_continuous(lim = c(-1, 1)) +
  theme_bw() + labs(y = "Factor loadings (mean pm S.E.)", fill = "IAS Taxa") +
  guides(color = FALSE) +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_text(colour="black",size=11,angle=0,hjust=.5,vjust=.5,face="plain"),
    axis.text.y = element_text(colour="black",size=11,angle=0,hjust=1,vjust=.5,face="plain"),
    panel.background = element_rect(colour = "black", size = 0.4, linetype = "solid"),
    legend.position = c(0.15, 0.15),
    legend.background = element_rect(fill = "white", color = "black", size = 0.3),
    legend.text = element_text(colour="black",size=11,angle=0, face="plain"),
    legend.title = element_text(colour="black",size=11,angle=0, face="bold"))

#Plotting beliefs about collateral impacts between the two groups
#---------------------------------------------------------
foo <- d
foo$ias.collateral <- foo$ias.collateral*-1
fee <- setNames(as.data.frame(prop.table(table(foo[which(foo$taxa.type == "single"), "ias.collateral"]))), c("score", "n"))
faa <- setNames(as.data.frame(prop.table(table(foo[which(foo$taxa.type == "multiple"), "ias.collateral"]))), c("score", "n"))
fee$score <- factor(fee$score)
faa$score <- factor(faa$score)
fee$group <- "single"
faa$group <- "multiple"
dat <- rbind(faa, fee)
dat$group <- factor(dat$group, levels = c("single", "multiple"), labels = c("Single group", "Multiple groups"))

dat %>%mutate(
  n = ifelse(group =="Single group", n*(-1), n*1))%>%
  ggplot(aes(x = score, y = n, fill = group)) +
  geom_bar(stat = "identity", color = "black") +
  geom_hline(aes(yintercept = 0)) +
  labs(fill = "Professional background", x = "Item score", y = "% of answers") +
  guides(fill = FALSE) +
  theme_bw() +
  coord_flip() + scale_fill_manual(values = c("steelblue", "gold")) +
  scale_y_continuous(breaks = c(-0.2, -0.1, 0, 0.1, 0.2), labels = c("20%", "10%", "", "10%", "20%")) +
  scale_x_discrete(position = "top", labels = c("Strongly \ndisagree", "Disagree", "Somehow \ndisagree",
                                                "Uncertain", "Somehow \nagree", "Agree", "Strongly \nagree")) +
  theme(
    axis.title.y = element_blank(),
    axis.text.x = element_text(colour="black",size=11,angle=0,hjust=.5,vjust=.5,face="plain"),
    axis.text.y = element_text(colour="black",size=11,angle=0,hjust=1,vjust=.5,face="plain"),
    panel.background = element_rect(colour = "black", size = 0.4, linetype = "solid"))

################################################################################
#Modelling attitudes towards IAS: practitioners vs researcher
################################################################################

#Configural invariance (LVMs have the same structure in the group: same LVs,
#same indicator variables, same pattern of constrained and estimated parameters.
#Estimated coefficients might differ between groups, and thus LVs might not measure
#the same construct in each group)
#---------------------------------------------------------
dat <- d[which(!is.na(d$practitioner.dummy)),]
mod <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
                      #---------------------------------------------------------
                 ias.collateral ~~ ias.interfere'

fit.conf.inv <- cfa(mod, data = dat, std.ov = TRUE, std.lv=TRUE, estimator="MLM",
                    group = "practitioner.dummy")

#Weak invariance (loadings are the same across groups. Latent variances are allowed
#to vary among groups: weak argument that LVs are equivalent between groups, but
#only latent variances and covariances could be compared)
#---------------------------------------------------------
fit.weak.inv <- cfa(mod, data = dat, std.ov = TRUE, std.lv=TRUE, estimator="MLM",
                    group = "practitioner.dummy", group.equal = "loadings")

#Strong invariance (For a given indicator variable, intercepts are the same among groups
#LV means are the same among groups. Individuals at the same LV score have the same
#expected values of indicator variables. LV means, variances, and covariances can be 
#compared. Group differences on strongly-invariant LVs, reflect real differences in LVs, 
#as opposed to the difference caused by differente measures)
#---------------------------------------------------------
fit.strong.inv <- cfa(mod, data = dat, std.ov = TRUE, std.lv=TRUE, estimator="MLM",
                      group = "practitioner.dummy", group.equal = c("loadings", "intercepts"))

#Testing changes (models behave similarly, strong invariance achieved)
#---------------------------------------------------------
lavTestLRT(fit.conf.inv, fit.weak.inv, fit.strong.inv)
fitmeasures(fit.conf.inv, c("cfi", "mfi", "rmsea", "srmr"))
fitmeasures(fit.weak.inv, c("cfi", "mfi", "rmsea", "srmr"))
fitmeasures(fit.strong.inv, c("cfi", "mfi", "rmsea", "srmr"))

bestmod <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
          #############################################
          ias.collateral ~~ ias.interfere'
bestfit <- cfa(mod.b, data = dat, std.ov = TRUE, std.lv = TRUE, estimator = "MLM", group = "practitioner.dummy")
fitmeasures(bestfit, fit.measures = c("fmin", "chisq", "df", "pvalue"))
fitmeasures(bestfit, "cfi") #Should be higher than 0.95
fitmeasures(bestfit, "tli") #Should be higher than 0.95
fitmeasures(bestfit, "srmr") #Should be smaller than 0.08
fitmeasures(bestfit, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(bestfit, "rmsea.ci.lower") #Its CI should be within 0.6-0.8
fitmeasures(bestfit, "rmsea.ci.upper") #Its CI should be within 0.6-0.8
residuals(bestfit, type="cor") #ias.interfere and ias.collateral covary
str(predict(bestfit))

#Plot scores and loadings
#---------------------------------------------------------
table(dat$practitioner.dummy)
dat[which(dat$practitioner.dummy== "Practitioner"), "score"] <- as.numeric(predict(bestfit)$Practitioner)
dat[which(dat$practitioner.dummy== "Researcher"), "score"] <- as.numeric(predict(bestfit)$Researcher)
dat$score <- as.numeric(dat$score)
boxplot(score ~ practitioner.dummy, data = dat)

ggplot() + 
  geom_violin(data = dat, aes(x = practitioner.dummy, y = score, fill = practitioner.dummy)) +
  geom_boxplot(data = dat, aes(x = practitioner.dummy, y = score), width = 0.1, color="black", alpha=0.6) +
  scale_fill_manual(values = c("steelblue", "gold")) +
  theme_bw() + labs(y = "PMA scores (Z-scores)", x = "Background") +
  scale_y_continuous(position = "right") +
  guides(fill = FALSE) +
  theme(
    axis.text.x = element_text(colour="black",size=11,angle=0,hjust=.5,vjust=.5,face="plain"),
    axis.text.y = element_text(colour="black",size=11,angle=0,hjust=1,vjust=.5,face="plain"),
    panel.background = element_rect(colour = "black", size = 0.4, linetype = "solid"))

inspect(bestfit, what = "std")$Practioner$lambda
inspect(bestfit, what = "std")$Researcher$lambda
summary(bestfit)

item <- c("ias.natconserv", "ias.costbenefit", "ias.nextgen", 
          "ias.nopriority", "ias.collateral", "ias.interfere", 
          "ias.notenough")
mean <- c(0.687, 0.414, 0.864, -0.519, -0.276, -0.368, 0.392,
          0.730, 0.503, 0.761, -0.592, -0.306, -0.511, 0.541)
se <- c(0.094, 0.058, 0.092, 0.066, 0.070, 0.069, 0.090,
        0.152, 0.104, 0.127, 0.109, 0.115, 0.097, 0.124)
cv <- se/mean
foo <- data.frame(item, mean, se, cv)
foo$item <- factor(foo$item, levels = c("ias.natconserv", "ias.costbenefit", "ias.nextgen", 
                                        "ias.nopriority", "ias.collateral", "ias.interfere", 
                                        "ias.notenough"))
foo$group <- c(rep("Practitioner", 7), rep("Researcher", 7))
foo$group <- factor(foo$group, levels = c("Practitioner", "Researcher"))

ggplot() + 
  geom_errorbar(data = foo, aes(x = item, ymin = mean - se, ymax = mean + se, color = group), 
                position = position_dodge(width = 0.5), width = 0.1) +
  geom_point(data = foo, aes(x = item, y = mean, fill = group), 
             shape = 21, size = 5, color = "black",  position = position_dodge(width = 0.5)) +
  geom_hline(aes(yintercept = 0), linetype = "dashed", size = 0.3) +
  scale_color_manual(values = c("black", "black")) + 
  #scale_fill_manual(values = c("#FF9966", "#66CC33")) +
  scale_fill_manual(values = c("steelblue", "gold")) +
  scale_x_discrete(labels = c("Nature \nconservation", "Costs and\nbenefits", "Next \ngenerations",
                              "Not a \npriority", "Collateral \ndamages", "Interference \nwith nature", "Not enough \ncontrol")) +
  scale_y_continuous(lim = c(-1, 1)) +
  theme_bw() + labs(y = "Factor loadings (mean pm S.E.)", fill = "Professional background") +
  guides(color = FALSE) +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_text(colour="black",size=11,angle=0,hjust=.5,vjust=.5,face="plain"),
    axis.text.y = element_text(colour="black",size=11,angle=0,hjust=1,vjust=.5,face="plain"),
    panel.background = element_rect(colour = "black", size = 0.4, linetype = "solid"),
    legend.position = c(0.15, 0.15),
    legend.background = element_rect(fill = "white", color = "black", size = 0.3),
    legend.text = element_text(colour="black",size=11,angle=0, face="plain"),
    legend.title = element_text(colour="black",size=11,angle=0, face="bold"))

################################################################################
#Comparing a CFA without groups with a CFA with taxa and background
################################################################################
mod.a <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
          #############################################
          ias.collateral ~~ ias.interfere'

fit.a <- cfa(mod.a, data = dat, std.ov= TRUE, std.lv=TRUE, estimator="MLM") #Satorra/Bentler correction
fit.b <- cfa(mod.a, data = dat, std.ov= TRUE, std.lv=TRUE, estimator="MLM", group = "taxa.type") #Satorra/Bentler correction
fit.c <- cfa(mod.a, data = dat, std.ov= TRUE, std.lv=TRUE, estimator="MLM", group = "practitioner.dummy") #Satorra/Bentler correction

fitmeasures(fit.a, fit.measures = c("fmin", "chisq", "df", "pvalue"))

fitmeasures(fit.a, "cfi", "tli", "srmr", "rmsea") #Should be higher than 0.95
fitmeasures(fit.b, "cfi", "tli", "srmr", "rmsea") #Should be higher than 0.95
fitmeasures(fit.c, "cfi", "tli", "srmr", "rmsea") #Should be higher than 0.95

fitmeasures(fit.a, "tli") #Should be higher than 0.95
fitmeasures(fit.b, "tli") #Should be higher than 0.95
fitmeasures(fit.c, "tli") #Should be higher than 0.95

fitmeasures(fit.a, "srmr") #Should be smaller than 0.08
fitmeasures(fit.b, "srmr") #Should be smaller than 0.08
fitmeasures(fit.c, "srmr") #Should be smaller than 0.08

fitmeasures(fit.a, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.b, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.c, "rmsea") #Its CI should be within 0.6-0.8

################################################################################
#SEM with Age
################################################################################

age <- read.csv(file="~/Authors_Age_20Jun2023.csv")
names(d)
names(age) <- c("n", "name", "year.birth", "year.degree")
age$age <- 2020 - age$year.birth
nrow(age[which(is.na(age$age)),]) #69 missing
nrow(age[which(!is.na(age$age)),]) #255 complete
foo <- dplyr::left_join(d, age, by = c("name"))
nrow(foo[which(!is.na(foo$age)),]) #248 researchers with complete age
d <- foo
dat <- foo[which(!is.na(foo$age)),]  #248 researchers with complete age

mean(dat$age)
sd(dat$age)

mod.sem <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
             attitudes ~ age
          #############################################
          ias.collateral ~~ ias.interfere'
fit.sem <- cfa(mod.sem, data=d, std.ov= TRUE, std.lv=TRUE, estimator="MLM") #Satorra/Bentler correction
fitmeasures(fit.sem, fit.measures = c("fmin", "chisq", "df", "pvalue"))

bestfit <- fit.a

fitmeasures(bestfit, "cfi") #Should be higher than 0.95
fitmeasures(fit.sem, "cfi") #Should be higher than 0.95

fitmeasures(bestfit, "tli") #Should be higher than 0.95
fitmeasures(fit.sem, "tli") #Should be higher than 0.95

fitmeasures(bestfit, "srmr") #Should be smaller than 0.08
fitmeasures(fit.sem, "srmr") #Should be smaller than 0.08

fitmeasures(bestfit, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.sem, "rmsea") #Its CI should be within 0.6-0.8

summary(fit.sem)


################################################################################
#SEM with N.taxa (issues)
################################################################################
mod.sem <- 'attitudes =~ NA*ias.natconserv + ias.costbenefit + ias.nextgen + ias.nopriority +
                      ias.collateral + ias.interfere + ias.notenough
             attitudes ~ taxa
          #############################################
          ias.collateral ~~ ias.interfere'
fit.sem <- cfa(mod.sem, data = dat, std.ov= TRUE, std.lv=TRUE, estimator="MLM") #Satorra/Bentler correction

bestfit <- fit.a

fitmeasures(bestfit, "cfi") #Should be higher than 0.95
fitmeasures(fit.sem, "cfi") #Should be higher than 0.95

fitmeasures(bestfit, "tli") #Should be higher than 0.95
fitmeasures(fit.sem, "tli") #Should be higher than 0.95

fitmeasures(bestfit, "srmr") #Should be smaller than 0.08
fitmeasures(fit.sem, "srmr") #Should be smaller than 0.08

fitmeasures(bestfit, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.sem, "rmsea") #Its CI should be within 0.6-0.8

summary(fit.sem)


################################################################################
#Exploring the distribution of answers (attitudes towards conservation, CA = 0.58)
################################################################################
{
#Descriptive statistics
#---------------------------------------------------------
idx <- c("nature.humwellbeing", "nature.plantanimals", "nature.priorhum", "nature.conssacrifice")
question.idx <- c("La natura è importante per il contributo che può fornire al piacere e al benessere degli esseri umani",
                  "Dovremmo proteggere la natura per il benessere di piante ed animali invece che per il benessere degli esseri umani",
                  "La felicità e la sopravvivenza del genere umano sono meno importanti di un pianeta sano",
                  "La conservazione della natura è importante anche quando riduce lo standard di vita delle persone")
for(i in seq_along(idx)){
  print(barplot(table(d[,c(idx[i])]), main = question.idx[i]))
}

d$nature.humwellbeing <- d$nature.humwellbeing*-1
cronbach.alpha(d[,c("nature.humwellbeing", "nature.plantanimals", "nature.priorhum", "nature.conssacrifice")])
cronbach.alpha(d[,c("nature.plantanimals", "nature.priorhum", "nature.conssacrifice")]) #0.66
cronbach.alpha(d[,c("nature.humwellbeing", "nature.priorhum", "nature.conssacrifice")])

#############################################
#ias.collateral ~~ ias.interfere'

mod.a <- 'attitudes =~ NA*nature.plantanimals + nature.priorhum + nature.conssacrifice'
fit.a <- cfa(mod.a, data=d, std.ov= TRUE, std.lv=TRUE, estimator="MLM") #Satorra/Bentler correction

fitmeasures(fit.a, fit.measures = c("fmin", "chisq", "df", "pvalue"))
fitmeasures(fit.a, "cfi") #Should be higher than 0.95
fitmeasures(fit.a, "tli") #Should be higher than 0.95
fitmeasures(fit.a, "srmr") #Should be smaller than 0.08
fitmeasures(fit.a, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.a, "rmsea.ci.lower") #Its CI should be within 0.6-0.8
fitmeasures(fit.a, "rmsea.ci.upper") #Its CI should be within 0.6-0.8
fit.a.cov <- fitted(fit.a)$cov #model-implied covariances
fit.a.cor <- cov2cor(fit.a.cov) # transform model-implied covariances to correlations
fit.a.cor # original correlations

mod.b <- 'attitudes =~ NA*nature.plantanimals + nature.priorhum + nature.conssacrifice
            nature.plantanimals ~~ nature.priorhum'
fit.b <- cfa(mod.b, data=d, std.ov= TRUE, std.lv=TRUE, estimator="MLM") #Satorra/Bentler correction

fitmeasures(fit.a, fit.measures = c("fmin", "chisq", "df", "pvalue"))
fitmeasures(fit.a, "cfi") #Should be higher than 0.95
fitmeasures(fit.a, "tli") #Should be higher than 0.95
fitmeasures(fit.a, "srmr") #Should be smaller than 0.08
fitmeasures(fit.a, "rmsea") #Its CI should be within 0.6-0.8
fitmeasures(fit.a, "rmsea.ci.lower") #Its CI should be within 0.6-0.8
fitmeasures(fit.a, "rmsea.ci.upper") #Its CI should be within 0.6-0.8
fit.a.cov <- fitted(fit.a)$cov #model-implied covariances
fit.a.cor <- cov2cor(fit.a.cov) # transform model-implied covariances to correlations
fit.a.cor # original correlations
}