---
title: "Data Analysis - European first-year university students accept evolution but lack substantial knowledge about it: a standardized European cross-country assessment"
author: "Dr. Alex Bergmann"
date: "22/12/2020"
output:
  html_document: default
  pdf_document: default
---



```{r, library & data setup}

knitr::opts_chunk$set(echo = FALSE)
knitr::opts_chunk$set(warning = FALSE)

rm(list=ls())

#install.packages(c("openxlsx","tidyverse", "mosaic", "plotly","plyr","naniar","ggrepel", "EnvStats", "VIM", "mice", "psych", "sjPlot", "lavaan", "lavaanPlot", "semPlot", "lme4", "lmerTest", "jtools", "robustlmm", "Hmisc", "reshape2", "car", "effects", "glmmTMB", "ggpubr"))

library(openxlsx) # reading
library(tidyverse) # data manipulation
library(mosaic) # descriptive analysis (inspect/favstats, gf_plots)
library(plotly)
library(plyr)
library(naniar) # NA handling
library(ggrepel)
library(EnvStats)
library(VIM) # MVA
library(mice) # MVA
library(psych) #Item & scale analysis
library(sjPlot) #Item & scale analysis, model diagnostics
library(lavaan) # CFA
library(lavaanPlot) # CFA
library(semPlot) # CFA
library(lme4) # multilevel modeling
library(lmerTest)# multilevel modeling
library(jtools)# multilevel modeling
library(robustlmm)# multilevel modeling
library(Hmisc) # correlation matrices
library(reshape2) # data wrangling (wide/long)
library(car) # multico
library(effects) # model diagnostics
library(glmmTMB) # model diagnostics
library(ggpubr) # arrange plots

setwd("/Users/Alex/Documents/Projekte/W&E zu Evolution")
set.seed(2308)

data_long <- read.xlsx("data_long_PK_AB.xlsx")
str(data_long)


```


```{r, functions setup}

#Source: https://github.com/tidyverse/ggplot2/blob/eecc450f7f13c5144069705ef22feefe0b8f53f7/R/geom-violin.r#L102

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}



# Source: https://stats.stackexchange.com/questions/233800/how-can-i-get-confidence-intervals-for-fixed-effects-using-the-rlmer-function-r

confint.rlmerMod <- function(object,parm,level=0.95) {
     beta <- fixef(object)
     if (missing(parm)) parm <- names(beta)
     se <- sqrt(diag(vcov(object)))
     z <- qnorm((1+level)/2)
     ctab <- cbind(beta-z*se,beta+z*se)
     colnames(ctab) <- stats:::format.perc(c((1-level)/2,(1+level)/2),
                                           digits=3)
     return(ctab[parm,])
 }
```




\

\

# 0. Data Preparation and Cleaning

```{r}

# set NA values correctly
data_long <- data_long %>%
  replace_with_na_all(condition = ~.x == 99)

# set variable types and levels correctly
data_long$sex <- as.factor(data_long$sex)
data_long$bio <- as.factor(data_long$bio)
data_long$bio <- revalue(data_long$bio, c(
                        "1"="Bio",
                        "2"="Non-Bio"))


data_long$country <- as.factor(data_long$country)
data_long$country <- revalue(data_long$country, c(
                        "1"="Austria",
                        "2"="Belgium",
                        "3"="Bosnia-Herzegovina",
                        "4"="Bulgaria",
                        "5"="Croatia",
                        "6"="Cyprus",
                        "7"="Czech Republic",
                        "8"="Finland",
                        "9"="France",
                        "10"="Germany",
                        "11"="Greece",
                        "12"="Hungary",
                        "13"="Italy",
                        "14"="Latvia",
                        "15"="Macedonia",
                        "16"="Netherlands",
                        "17"="Poland",
                        "18"="Portugal",
                        "19"="Romania",
                        "20"="Serbia",
                        "21"="Slovakia",
                        "22"="Slovenia",
                        "23"="Spain",
                        "24" = "Sweden",
                        "25" = "Switzerland",
                        "26" = "Turkey",
                        "27" = "Ukraine"
                        ))

data_long$denom_summarized <- as.factor(data_long$denom_summarized)
data_long$denom_summarized <- revalue(data_long$denom_summarized, c(
   "0"="No answer",
  "1"="Protestant",
  "2"="Christian free churches",
  "3"="Catholic",
  "4"="Orthodox",
  "5"="Muslim",
  "11"="None",
  "12" = "Other"
))



```

\

\


# 1. Missing Value Analysis

```{r}
# remove known systematic missing values (three Spanish courses, that were not provided with questions concernign their religious faith)

data_long <- data_long%>%
  dplyr::filter(is.na(course) | course != "4" & 
           course != "5" &
           course != "6")

aggr_plot <- aggr(data_long, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(df), cex.axis=.5, gap=2, ylab=c("Histogram of missing data","Pattern"))

# missing values on main variables
kaevo_count <- data_long%>%
  dplyr::select(KAEVO.A1,
                KAEVO.A2,
                KAEVO.A3,
                KAEVO.A4,
                KAEVO.A5,
                KAEVO.A6,
                KAEVO.A7,
                KAEVO.A8,
                KAEVO.A9.1,
                KAEVO.A9.2,
                KAEVO.A10,
                KAEVO.A11)
kaevo_count$na_count <- rowSums(is.na(kaevo_count))
table(kaevo_count$na_count)

perf_count <- data_long%>%
  dplyr::select(PERF.F1:PERF.F10)

perf_count$na_count <- rowSums(is.na(perf_count))
table(perf_count$na_count)

atevo_count <- data_long%>%
  dplyr::select(ATEVO.E1:ATEVO.E8)
atevo_count$na_count <- rowSums(is.na(atevo_count))
table(atevo_count$na_count)

# summary main scales
main_count <- data.frame(kaevo_count,atevo_count,perf_count)
main_count$na_count_all <- rowSums(is.na(main_count))
table(main_count$na_count_all)

# summary all variables
data_long$na_count_all <- rowSums(is.na(data_long))
table(data_long$na_count_all)

one_missing <- data_long%>%
  dplyr::filter(na_count_all == 10)

aggr_plot <- aggr(one_missing, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(df), cex.axis=.5, gap=2, ylab=c("Histogram of missing data","Pattern"))

data_clean <- data_long%>%
  dplyr::select(-na_count_all )

```

\

\

# 2. Item- and scale analysis

## 2.1. EFA: Religious faith & acceptance of evolution

```{r}

perf_fa <- data_clean%>%
  dplyr::select(PERF.F1:PERF.F10)

items.parallel <- fa.parallel(perf_fa, fa="both", quant = .95)
```

\

\
```{r}
cortest.bartlett(perf_fa)
KMO(perf_fa)
```
\

\

```{r}
items.pa <- fa(perf_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 1,
               rotate="varimax")
```

\

\
```{r}
print(items.pa)
```

\

\
```{r}
items.pa <- fa(perf_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 2,
               rotate="varimax")
```

\

\
```{r}
print(items.pa)
sjt.itemanalysis(perf_fa)
```


```{r}
kaevo_fa <- data_clean%>%
  dplyr::select(KAEVO.A1:KAEVO.A11)

items.parallel <- fa.parallel(kaevo_fa, fa="both", quant = .95)
cortest.bartlett(kaevo_fa)
```

```{r}
KMO(kaevo_fa)
```


```{r}
items.pa <- fa(kaevo_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 4,
               rotate="varimax")
```


```{r}
print(items.pa)
```

```{r}
items.pa <- fa(kaevo_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 3,
               rotate="varimax")
```


```{r}
print(items.pa)
```

```{r}
items.pa <- fa(kaevo_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 2,
               rotate="varimax")
```


```{r}
print(items.pa)
```



```{r}
items.pa <- fa(kaevo_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 1,
               rotate="varimax")
```


```{r}
print(items.pa)
```



```{r}
atevo_fa <- data_clean%>%
  dplyr::select(ATEVO.E1:ATEVO.E8)

items.parallel <- fa.parallel(atevo_fa, fa="both", quant = .95)
cortest.bartlett(atevo_fa)
```


```{r}
KMO(atevo_fa)
```


```{r}
items.pa <- fa(atevo_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 1,
               rotate="varimax")
```


```{r}
print(items.pa)
```


```{r, include = F}
items.pa <- fa(atevo_fa,
               SMC=TRUE,
               fm="pa",
               nfactors = 2,
               rotate="varimax")
```
\

\
```{r}
print(items.pa)
sjt.itemanalysis(atevo_fa)
```
\



## 2.2. CFA: Knowledge about evolution

```{r}

model <-
    'ada =~ KAEVO.A1 + KAEVO.A3 + KAEVO.A5 + KAEVO.A6
    her =~ KAEVO.A7 + KAEVO.A8
    tre =~ KAEVO.A9.1 + KAEVO.A9.2
    spe =~ KAEVO.A4 + KAEVO.A10'
  
  fit <- cfa(model, data=kaevo_fa,std.lv=TRUE)
summary(fit,fit.measures=TRUE,standardized=TRUE)
parameterEstimates(fit)

```

\

## 2.3. Indexing

```{r}

perf_data <- data_clean%>%
  dplyr::select(PERF.F1:PERF.F10)
perf_data$perf_na <- rowSums(is.na(perf_data))

perf_data <- perf_data%>%
  mutate(perf = ifelse(perf_na == 0, (PERF.F1 + PERF.F2 + PERF.F3 + PERF.F4+ PERF.F5+ PERF.F6+ PERF.F7+ PERF.F8 + PERF.F9 + PERF.F10), "NA"))%>%
  select(perf)

atevo_data <- data_clean%>%
  dplyr::select(ATEVO.E1:ATEVO.E8)
atevo_data$atevo_na <- rowSums(is.na(atevo_data))

atevo_data <- atevo_data%>%
  mutate(atevo = ifelse(atevo_na == 0, (ATEVO.E1 + ATEVO.E2 + ATEVO.E3 + ATEVO.E4 + ATEVO.E5 + ATEVO.E6 + ATEVO.E7 + ATEVO.E8), "NA"))%>%
  select(atevo)

kaevo_data <- data_clean%>%
  dplyr::select(KAEVO.A1,
                KAEVO.A2,
                KAEVO.A3,
                KAEVO.A4,
                KAEVO.A5,
                KAEVO.A6,
                KAEVO.A7,
                KAEVO.A8,
                KAEVO.A9.1,
                KAEVO.A9.2,
                KAEVO.A10,
                KAEVO.A11)
kaevo_data$kaevo_na <- rowSums(is.na(kaevo_data))

kaevo_data <- kaevo_data%>%
  mutate(kaevo = ifelse(kaevo_na == 0, (KAEVO.A1 +
                KAEVO.A2 + 
                KAEVO.A3 + 
                KAEVO.A4 + 
                KAEVO.A5 + 
                KAEVO.A6 + 
                KAEVO.A7 + 
                KAEVO.A8 + 
                KAEVO.A9.1 + 
                KAEVO.A9.2 + 
                KAEVO.A10 + 
                KAEVO.A11), "NA"))%>%
  select(kaevo)

data <- data.frame(data_clean, perf_data, atevo_data, kaevo_data)

data <- data%>%
  dplyr::select(ID, country, bio, course, age, sex, bio_classes, interest_bio, meaning_evo, learn_evo, denomination, denom_summarized, perf, kaevo, atevo)

data$perf <- as.numeric(data$perf)
data$atevo <- as.numeric(data$atevo)
data$kaevo <- as.numeric(data$kaevo)
data$sex <- as.factor(data$sex)

```


\

# 3. Descriptive statistics 


**Data preparation**
```{r}

# preparing cat data set & relevel factors

cat_data_atevo <- data%>%
  dplyr::select(ID, atevo, bio, sex)%>%
  drop_na(atevo)%>%
  mutate(atevo_cat = ifelse(atevo %in% 8:13, "reject",
                            ifelse(atevo %in% 14:19, "rather reject",
                                   ifelse(atevo %in% 20:28, "neutral",
                                            ifelse(atevo %in% 29:34, "rather accept","accept")))))
cat_data_atevo$atevo_cat <- factor(cat_data_atevo$atevo_cat, levels = c("reject", "rather reject", "neutral", "rather accept", "accept"))




cat_data_kaevo <- data%>%
    dplyr::select(ID, kaevo, bio, sex)%>%
  drop_na(kaevo)%>%
  mutate(kaevo_cat = ifelse(kaevo %in% 0:5, "very low",
                                            ifelse(kaevo %in% 6:7, "low",
                                                   ifelse(kaevo %in% 8:9, "moderate",
                                                          ifelse(kaevo %in% 10:11, "rather high","high")))))

cat_data_kaevo$kaevo_cat <- factor(cat_data_kaevo$kaevo_cat, levels = c("very low", "low", "moderate", "rather high", "high"))






cat_data_perf <- data%>%
      dplyr::select(ID, perf, bio, sex)%>%
  drop_na(perf)%>%
  mutate(perf_cat = ifelse(perf %in% 10:17, "not religious at all",
                                            ifelse(perf %in% 8:25, "not religious",
                                                   ifelse(perf %in% 26:34, "neutral",
                                                          ifelse(perf %in% 35:42, "religious","very religious")))))


cat_data_perf$perf_cat <- factor(cat_data_perf$perf_cat, levels = c("not religious at all", "not religious", "neutral", "religious", "very religious"))



```

\

\

**Cats religios faith**
```{r}
cat_data_perf%>%
   group_by(perf_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))

cat_data_perf%>%
   group_by(bio, perf_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))


cat_data_perf%>%
  filter(bio == "Bio")%>%
   group_by(perf_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))


cat_data_perf%>%
  filter(bio == "Non-Bio")%>%
   group_by(perf_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))
```

\

\


**Cats knowledge**
```{r}
cat_data_kaevo%>%
   group_by(kaevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))

cat_data_kaevo%>%
   group_by(bio, kaevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))

cat_data_kaevo%>%
  filter(bio == "Bio")%>%
   group_by(kaevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))


cat_data_kaevo%>%
  filter(bio == "Non-Bio")%>%
   group_by(kaevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))

```

\

\

**Cats acceptance**
```{r}
cat_data_atevo%>%
   group_by(atevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))

cat_data_atevo%>%
   group_by(bio, atevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))

cat_data_atevo%>%
  filter(bio == "Bio")%>%
   group_by(atevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))


cat_data_atevo%>%
  filter(bio == "Non-Bio")%>%
   group_by(atevo_cat) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))

favstats(atevo ~ bio, data = cat_data_atevo)
favstats(kaevo ~ bio, data = cat_data_kaevo)
favstats(perf ~ bio, data = cat_data_perf)


favstats(~ atevo, data = cat_data_atevo)
favstats(~ kaevo, data = cat_data_kaevo)
favstats(~ perf, data = cat_data_perf)


cat_data_atevo%>%
  filter(atevo_cat == "reject")%>%
   group_by(bio) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))


cat_data_atevo%>%
  filter(atevo_cat == "rather reject")%>%
   group_by(bio) %>%
   dplyr::summarise(count = n()) %>%
   mutate(prop = count/sum(count))


```

\

\

**Favstats & Frequencies**
```{r}

tally(~ bio, data = data)
tally(~ bio, data = data, format = "proportion")
tally(~ sex, data = data)
tally(~ country, data = data)

tally(~ bio | country, data = data)
tally(~ bio | country, data = data, format = "proportion")
tally(~ sex| country, data = data)

favstats( ~ atevo, data = data)
favstats(atevo ~ sex, data = data)
favstats(atevo ~ bio, data = data)
favstats(atevo ~ country, data = data)

favstats( ~ kaevo, data = data)
favstats(kaevo ~ sex, data = data)
favstats(kaevo ~ bio, data = data)
favstats(kaevo ~ country, data = data)

favstats( ~ perf, data = data)
favstats(perf ~ sex, data = data)
favstats(perf ~ bio, data = data)
favstats(perf ~ country, data = data)



```

\

\

```{r}
cor_data <- data%>%
  dplyr::select(perf, atevo, kaevo)
```

*Correlations*
```{r}

rcorr(as.matrix(cor_data), type = "pearson")

```


**t-tests (multilevel)**
```{r}

model_0_a <-lmer(atevo ~ 1 + (1 | country), REML = F, data = data)
summary(model_0_a)
summ(model_0_a)

model_atevo <- lmer(atevo ~ bio + (1 | country), REML = F, data = data)
summary(model_atevo)
summ(model_atevo)

#f2 = 0.01 (very small)

model_0_k <-lmer(kaevo ~ 1 + (1 | country), REML = F, data = data)
summ(model_0_k)
model_kaevo <- lmer(kaevo ~ bio + (1 | country), REML = F, data = data)
summary(model_kaevo)
summ(model_kaevo)
#f2 = 0.04 (medium)

model_0_p <-lmer(perf ~ 1 + (1 | country), REML = F, data = data)
summ(model_0_p)
model_perf <- lmer(perf ~ bio + (1 | country), REML = F, data = data)
summary(model_perf)
summ(model_perf)

#f2 = 0.00 (no effect)

```

\

\

# 4. Descriptive Plotting (main scales)

```{r}
# Bubbleplot

size_df <- data%>%
  group_by(country)%>%
  dplyr::summarise(n = n())

size_df <- size_df%>%
  select(-country)

data_b <- data%>%
  group_by(country)%>%
  dplyr::summarize(
            perf_mean = mean(perf, na.rm = T),
            perf_sd = sd(perf, na.rm = T),
            perf_n =  n(),
            perf_se = perf_sd / sqrt(perf_n),
            atevo_mean = mean(atevo, na.rm = T),
            atevo_sd = sd(atevo, na.rm = T),
            atevo_n =  n(),
            atevo_se = atevo_sd / sqrt(atevo_n),
            kaevo_mean = mean(kaevo, na.rm = T),
            kaevo_sd = sd(kaevo, na.rm = T),
            kaevo_n =  n(),
            kaevo_se = kaevo_sd / sqrt(kaevo_n))
            

data_bubble <- cbind(data_b, size_df)


y_tilt <- 33.5 - 32.63514
x_tilt <- 4.551887 - 4
mygrey = "black"
mysize = 0.15

p1 <- ggplot(data_bubble, aes(x = kaevo_mean, y = atevo_mean))+
  scale_size_continuous(name = "Sample size",
                        guide = guide_legend(override.aes = list(color = "#7a0177")))+
  scale_colour_gradient(name = "Religious faith score",
                        low = "#66c2a4", high = "#7a0177")+
  xlab("Knowledge score")+
  ylab("Acceptance score")+
  geom_point(aes(size = kaevo_n, colour = perf_mean), alpha = .8)+ 
  geom_errorbarh(aes(
                    xmin= kaevo_mean - kaevo_se, 
                    xmax= kaevo_mean + kaevo_se), 
                width= .05,
                size = .2)+
  geom_errorbar(aes(ymin= atevo_mean - atevo_se,ymax= atevo_mean + atevo_se), 
    width= .05,
    size = .2)+
    theme_bw()+
  scale_x_continuous(limits=c(0, 12), 
                     breaks=c(0,2, 4, 6, 8, 10, 12))+
  scale_y_continuous(limits=c(8, 40))+
  geom_rect(mapping=aes(xmin=2, xmax=8.5, ymin=28, ymax=37), color="black", size = 0.25, alpha=0)+
annotate(geom = "text", x = 2.5, y = 36, label = "b", color = mygrey, fontface = "bold", size = 3, hjust = 1, vjust = .5)+
  annotate(geom = "text", x = 0.5, y = 40, label = "a", color = mygrey, fontface = "bold", size = 7, hjust = 1, vjust = .5)







p2 <- ggplot(data_bubble, aes(x = kaevo_mean, y = atevo_mean))+
  scale_size_continuous(name = "Sample size",
                       guide = guide_legend(override.aes = list(color = "#7a0177")))+
  scale_colour_gradient(name = "Religious faith score",
                        low = "#66c2a4", high = "#7a0177")+
  xlab("Knowledge score")+
  ylab("Acceptance score")+
  geom_point(aes(size = kaevo_n, colour = perf_mean), alpha = .8)+ 
  geom_errorbarh(aes(
                    xmin= kaevo_mean - kaevo_se, 
                    xmax= kaevo_mean + kaevo_se), 
                width= .05,
                size = .25)+
  geom_errorbar(aes(ymin= atevo_mean - atevo_se,ymax= atevo_mean + atevo_se), 
    width= .05,
    size = .25)+
    theme_bw()+
  scale_x_continuous(limits=c(2, 8), 
                     breaks=c(2, 4, 6, 8))+
  scale_y_continuous(limits=c(28, 37), 
                     breaks=c(28, 30, 32, 34, 36))+
  #theme(axis.title.x = element_blank(),
     #   axis.title.y= element_blank())+
  annotate(geom = "text", x = 7.8 + 0.2, y = 36.085, label = "FI", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 7.513447 + 0.2, y = 35.01856, label = "NL", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 7.126316 + 0.2, y = 34.70588 - 0.2, label = "ES", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 6.127418 + 0.2, y = 32.95301 + 0.0, label = "DE", color = mygrey, size = 3, fontface = "bold", hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 6.127321 + 0.2 , y = 33.46875, label = "BE", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 6.830508 + 0.3, y = 32.70312, label = "CH", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 6.375806 + 0.15, y = 31.67107, label = "FR", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
      annotate(geom = "text", x = 5.620427 + 0.15, y = 32.26619, label = "IT", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 5.152975 + 0.2, y = 31.91957, label = "HR", color = mygrey, fontface = "bold",size = 3, hjust = 0, vjust = .50)+
    annotate(geom = "text", x = 4.120623 - 0.25, y = 31.16443 + 0.3, label = "SI", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+

    annotate(geom = "text", x = 3.944444 - 0.35, y = 30.99648 - 0.3, label = "RS", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+

    annotate(geom = "text", x = 4.2 + 0.3, y = 31.12, label = "CY", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 4.5 + 0.2, y = 30.30028, label = "CZ", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 2.88587 + 0.2, y = 29.94144, label = "BA", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 2.801418 - 0.2, y = 30.16867, label = "LV", color = mygrey, fontface = "bold", size = 3, hjust = 1, vjust = .50)+
  annotate(geom = "text", x = 3.032258 + 0.35, y = 30.68085 - 0.2, label = "SK", color = mygrey, fontface = "bold", size = 3, hjust = 1, vjust = .50)+
  annotate(geom = "text", x = 3.034926 - 0.2, y = 30.8281 + 0.1, label = "RO", color = mygrey, fontface = "bold", size = 3, hjust = 1, vjust = .50)+
  annotate(geom = "text", x = 3.125 + 0.2, y = 31.21053 + 0.1, label = "UA", color = mygrey, size = 3, fontface = "bold", hjust = 0, vjust = .50)+
  annotate(geom = "text", x = 2.865672 - 0.2, y = 31.40244, label = "EL", color = mygrey, fontface = "bold", size = 3, hjust = 1, vjust = .50)+
  annotate(geom = "text", x = 4.588235 - 0.2, y = 32.44286 - 0.15, label = "PT", color = mygrey, size = 3, fontface = "bold", hjust = 1, vjust = .50)+
  annotate(geom = "text", x = 3.916667 - 0.2, y = 32.63768 - 0.15, label = "TR", color = mygrey, fontface = "bold", size = 3, hjust = 1, vjust = .50)+
  annotate(geom = "text", x = 3.608108 - 0.2, y = 32.95109, label = "BG", color = mygrey, fontface = "bold", size = 3, hjust = 1, vjust = .50)+
  annotate(geom = "text", x = 4.551887 + 0.2, y = 32.63514 + 0.1, label = "HU", color = mygrey, fontface = "bold", size = 3, hjust = 0, vjust = .50)+
    annotate(geom = "text", x = 4.967254 - 0.2, y = 33.28571, label = "PL", color = mygrey, fontface = "bold",size = 3, hjust = 1, vjust = .50)+
    annotate(geom = "text", x = 5.97931 - 0.2, y = 33.36538 + 0.1, label = "AT", color = mygrey, fontface = "bold", size = 3, hjust = .5, vjust = 0)+
    annotate(geom = "text", x = 5.451613 - 0.15, y = 33.1875 + 0.1, label = "SE",  color = mygrey,fontface = "bold", size = 3, hjust = .5, vjust = 0)+
  annotate(geom = "text", x = 2.3, y = 37, label = "b", color = mygrey, fontface = "bold", size = 7, hjust = 1, vjust = .5)


ggarrange(p1, p2, ncol=2, nrow=1, common.legend = TRUE, legend="right")

ggsave(
  "bubbleplot_final",
  device = "jpeg",
  scale = 1,
  width = 26,
  height = 12,
  units = "cm",
  dpi = 300,
  limitsize = TRUE
)

```

\

\

```{r}
ggplot(data, aes(x = kaevo, y = atevo, colour = perf))+
  geom_jitter(alpha = .5, size = .2, stroke = 0.8)+
  xlab("Knowledge score")+
  ylab("Acceptance score")+
  scale_colour_gradient(name = "Religious faith score",
                        low = "#66c2a4", high = "#7a0177",
                        na.value=NA)+
  theme_bw()

ggsave(
  "scatterplot_final",
  device = "jpeg",
  scale = 1,
  width = 16,
  height = 12,
  units = "cm",
  dpi = 300,
  limitsize = TRUE
)

```

\

\


# 5. Scale Comparison

```{r}

comparison_data <- data%>%
  dplyr::mutate(perf_t = ((perf - 10)/40)*100,
         atevo_t = ((atevo - 8)/32)*100,
         kaevo_t = ((kaevo - 0)/12)*100)%>%
  dplyr::select(ID, perf_t, atevo_t, kaevo_t, bio)

comp_tab <- comparison_data%>%
  dplyr::summarize(
            perf_mean = mean(perf_t, na.rm = T),
            perf_sd = sd(perf_t, na.rm = T),
            perf_n =  n(),
            perf_se = perf_sd / sqrt(perf_n),
            atevo_mean = mean(atevo_t, na.rm = T),
            atevo_sd = sd(atevo_t, na.rm = T),
            atevo_n =  n(),
            atevo_se = atevo_sd / sqrt(atevo_n),
            kaevo_mean = mean(kaevo_t, na.rm = T),
            kaevo_sd = sd(kaevo_t, na.rm = T),
            kaevo_n =  n(),
            kaevo_se = kaevo_sd / sqrt(kaevo_n))

melt_comparison <- melt(comparison_data, id=c("ID", "bio"))

melt_comparison$variable <- factor(melt_comparison$variable , levels=c("atevo_t", "kaevo_t", "perf_t"))


#y = ((x - min(x))/(max(x) - min(x)))*100
mylabels <- c("Acceptance score", "Knowledge score", "Religious faith score")

ggplot(melt_comparison, aes(x = variable, y = value))+
  geom_violin()+
  geom_boxplot(width=0.1)+
  xlab("")+
  stat_n_text(size = 3) +
  ylab("Scaled values (0-100)")+
  scale_x_discrete(labels= mylabels)+
  theme_bw()

favstats(value ~ variable, data = melt_comparison)
favstats(value ~ variable | bio, data = melt_comparison)


melt_comparison$bio <- revalue(melt_comparison$bio, c(
  "Bio" = "biology-related",
  "Non-Bio" = "non-biology"
))


melt_comparison <- transform(melt_comparison,
          bio=revalue(bio,c("biology-related"="biology-related")))

melt_comparison$bio <- ordered(melt_comparison$bio)


#melt_comparison%>%
  #filter(bio =="biology-related")

#melt_comparison%>%
 # filter(bio =="non-biology")

melt_comparison_na <- melt_comparison%>%
  drop_na()




melt_comparison_na%>%
  group_by(variable)%>%
  dplyr::summarise(n = n())

a <- expression("Acceptance score"~"("~italic("n")~"= 8,526"~")")
b <- expression("Knowledge score"~"("~italic("n")~"= 7,805"~")")
c <- expression("Religious faith score"~"("~italic("n")~"= 8,353"~")")

mylabels <- c(a,b,c)


```

\

\

```{r}
ggplot(melt_comparison_na)+
  geom_split_violin(aes(x = variable, y = value, fill = bio), alpha = .85)+
  geom_boxplot(aes(x = variable, y = value, fill = bio), width=0.2, outlier.alpha = .05, outlier.fill = "white", outlier.size = 0.8, alpha = 1)+
  scale_fill_manual(name ="",
                    labels = c("biology-related", "non-biology"),
                      values = c("#0072B2", "#56b4df"))+
  xlab("")+
  ylab("Scaled values (0-100)")+
  scale_x_discrete(labels= mylabels)+
  theme_bw()

ggsave(
  "comparison_final",
  device = "jpeg",
  scale = 1,
  width = 20,
  height = 15,
  units = "cm",
  dpi = 300,
  limitsize = TRUE
)

```

\

\

# 6. Descriptive analysis for modelling procedure + data cleaning 2.0.

```{r}


# remove countries with less than 150 obervations
data <- data %>% 
  dplyr::group_by(country) %>% 
  dplyr::mutate(freq = n()) %>% 
  ungroup() %>% 
  filter(freq >= 150)%>%
  dplyr::select(-freq)


# Descriptives & further Cleaning
gf_jitter(atevo ~ perf, data = data)
gf_jitter(atevo ~ kaevo, data = data)
gf_jitter(atevo ~ bio, data = data)
gf_jitter(atevo ~ interest_bio, data = data)
gf_jitter(atevo ~ age, data = data)
gf_jitter(atevo ~ sex, data = data)


data <- data%>%
    filter(sex != 3)%>%
    filter(bio != "NA")%>%
  drop_na() 




data %>% 
  ggplot(aes(x = atevo,
             y = perf)) +
  geom_smooth(aes(group = country),
              method = "lm",
              se     = FALSE,      
              size   = 0.3) +   
  geom_smooth(method = "lm",
              se     = FALSE,
              color = "red",       
              size   = 2) +        
  theme_bw()

```

\


```{r}
ggplot(data, aes(x = atevo, y = perf))+
  geom_point()+
  geom_smooth(method = "lm")+
  facet_grid(~country)
```

\


```{r}
tally(~country, data = data)

greece_data <- data%>%
  filter(country == "Greece")%>%
  filter(atevo > 20)

ggplot(greece_data, aes(x = atevo, y = perf))+
  geom_point()+
  geom_smooth(method = "lm")
```

\


```{r}
greece_data <- data%>%
  filter(country == "Greece")%>%
  filter(atevo > 20)

ggplot(greece_data, aes(x = atevo, y = perf))+
  geom_point()+
  geom_smooth(method = "lm")
```


\


```{r}
data %>% 
  ggplot(aes(x = atevo,
             y = kaevo)) +
  geom_smooth(aes(group = country),
              method = "lm",
              se     = FALSE,      
              size   = 0.3) +   
  geom_smooth(method = "lm",
              se     = FALSE,
              color = "red",       
              size   = 2) +        
  theme_bw() 
```

\

```{r}
gf_jitter(atevo ~ perf, data = data)
gf_jitter(atevo ~ kaevo, data = data)
gf_jitter(atevo ~ bio, data = data)
gf_jitter(atevo ~ interest_bio, data = data)
gf_jitter(atevo ~ age, data = data)
gf_jitter(atevo ~ sex, data = data)

```


\

\

# 7. Multilevel Modeling


```{r}
names(data)

model_0 <-lmer(atevo ~ 1 + (1 | country), REML = F, data = data)
summary(model_0)
summ(model_0)
ranova(model_0)
```

\

\


```{r}
model_1 <-lmer(atevo ~ 1 + age + sex + (1 | country), REML = F, data = data)
summary(model_1)
summ(model_1)
ranova(model_1)


confint(model_1, method = "boot", nsim = 10000,
        parallel = "multicore", ncpus = 4)

model_2 <-lmer(atevo ~ 1 + age + sex + bio + interest_bio + (1 | country), REML = F, data = data)
summary(model_2)
summ(model_2)
ranova(model_2)

confint(model_2, method = "boot", nsim = 10000,
        parallel = "multicore", ncpus = 4)


model_3 <-lmer(atevo ~ 1 + age + sex + bio + interest_bio + kaevo + (1 | country), REML = F, data = data)
summary(model_3)
summ(model_3)
ranova(model_3)

confint(model_3, method = "boot", nsim = 10000,
        parallel = "multicore", ncpus = 4)


model_4 <-lmer(atevo ~ 1 + age + sex + bio + interest_bio +  kaevo + perf + (1  | country), REML = F,  data = data)

summary(model_4)
summ(model_4)
ranova(model_4) 

confint(model_4, method = "boot", nsim = 10000,
        parallel = "multicore", ncpus = 4)

data$denom_summarized <- factor(data$denom_summarized, levels = c("None", "Protestant", "Catholic", "Orthodox", "Christian free churches", "Muslim", "Other"))

model_4_denom <-lmer(atevo ~ 1 + age + sex + bio + interest_bio +  kaevo + perf + denom_summarized + (1  | country), REML = F,  data = data)

summary(model_4_denom)
summ(model_4_denom)


confint(model_4_denom, method = "boot", nsim = 10000,
        parallel = "multicore", ncpus = 4)



model_5 <-lmer(atevo ~ 1 + age + sex + bio + interest_bio +  kaevo*perf + (1  | country), REML = F,  data = data)
summary(model_5) 
summ(model_5)
ranova(model_5)





```


\


*Modelselection via LRT*
```{r}

anova(model_0,
      model_1,
      model_2,
      model_3,
      model_4,
      model_4_denom) #--> Model 4 it is, da Model 5 VIF-Anaylse MulticolinearitÃ¤t nahelegt
```

\


*Multicollinearity and model diagnostics*

```{r}
vif(model_4) 
plot_model(model_4, type='diag')
```


\

*Robust Alternative*
```{r}

r_model_4 <-rlmer(atevo ~ 1 + age + sex + bio + interest_bio + perf + kaevo + (1  | country), REML = F,  data = data) 
#saveRDS(r_model_4, "r_model_4.rds")
#r_model_4 <- readRDS("r_model_4.rds")
summary(r_model_4)
confint(r_model_4)
```


\

*Additional tests: Random Intercept + Random Slope*
```{r}

data$age_s <- scale(data$age, center = T, scale = T)
data$interest_bio_s <- scale(data$interest_bio, center = T, scale = T)
data$perf_s <- scale(data$perf, center = T, scale = T)
data$kaevo_s <- scale(data$kaevo, center = T, scale = T)


model_4_rs <-lmer(atevo ~ 1 + age_s + sex + bio + interest_bio_s + perf_s + kaevo_s + (perf_s | country), REML = F,  data = data)

summary(model_4_rs)
summ(model_4_rs)
plot_model(model_4_rs, type='diag')
anova(model_0,
      model_1,
      model_2,
      model_3,
      model_4,
      #model_5,
      model_4_rs)

```

```{r}
sessionInfo()
```

