---
title: "EvoCost_Countryalpha"
author: "Alex Bergmann"
date: "12/11/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)


```

\

\

# 0. Data Preparation and Cleaning

```{r, include = F}

# 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, include = F}
# remove known systematic missings
data_long <- data_long%>%
  dplyr::filter(is.na(course) | course != "4" & 
           course != "5" &
           course != "6")

str(data_long)

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"))

#percentmiss <- function(x) {sum(is.na(x)) / length(x)*100}
#apply(data_long, 2, percentmiss)


# 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))
#str(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))
#str(perf_count)
table(perf_count$na_count)
559/9200*100

atevo_count <- data_long%>%
  dplyr::select(ATEVO.E1:ATEVO.E8)


atevo_count$na_count <- rowSums(is.na(atevo_count))
#str(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)

6579/9200*100
1235/9200*100

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 )

```

\

\

# 1. Item- and scale analysis:  PERF

```{r}
vec <- levels(data_clean$country)


### PERF
for(x in vec){
  
  co <- data_clean%>%
    filter(country == x)%>%
    distinct(country)
  
  nr <- data_clean%>%
    filter(country == x)%>%
    nrow()
  
  data_c <- data_clean%>%
    filter(country == x)%>%
    dplyr::select(PERF.F1:PERF.F10)
  
  a <- (sjt.itemanalysis(data_c))
  print(co)
  print(nr)
  print(a$cronbach.values)
}


```

# 2. Item- and scale analysis: ATEVO

```{r}
### ATEVO
for(x in vec){
  
  co <- data_clean%>%
    filter(country == x)%>%
    distinct(country)
  nr <- data_clean%>%
    filter(country == x)%>%
    nrow()
  
  data_c <- data_clean%>%
    filter(country == x)%>%
    dplyr::select(ATEVO.E1:ATEVO.E8)
  
  
  a <- (sjt.itemanalysis(data_c))
  print(co)
  print(nr)
  print(a$cronbach.values)
}
```

