GISD_Generate_2017_revision2020.Rmd

---
title: "GISD - German Index of Socio-Economic Deprivation"
author: "Niels Michalski (Überarbeitung der Syntax von Lars Kroll)"
date: "21 April 2020"
output: 
  bookdown::html_document2:
    keep_md: true
    code_folding: hide
    toc: true
    toc_float: false
    toc_depth: 2
    number_sections: false
    fig_caption: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Intro

Auf diesem Codeblog stelle ich die Generierung des German Index of Socio-Economic Deprivation vor. Dabei handelt es sich um einen Index sozioökonomischer Deprivation auf regionalräumlicher Ebene, der seit 2012 im Fachgebiet Soziale Determinanten der Gesundheit am RKI entwickelt wurde und seither jährlich aktualisiert wird. Für die Generierung werden Indikatoren der INKAR-Datenbank des BBSR verwendet. Es wird die Revision 2020 vorgestellt, die Daten aus den Jahren 1998 bis 2017 verwendet. 
In einer früheren Revision für 2019 tauchten zum Teil starke Abweichungen zur Revision 2018 auf. Der Grund dafür war, dass bei der Addition der Teildimensionen die Bildungsdimension in umgekehrter Richtung in den GISD einging. Ausschlaggebend für diesen Fehler war eine negative Korrelation der Bildungsdimension mit dem Anteil Arbeitsloser, die im Code der vergangenen Revisionen zu einer Umpolung des Teilscores führt. Der Anteil Arbeitsloser wird als Markerindikator verwendet. Die folgende Darstellung zeigt und diskutiert die Problematik und stellt eine Lösung vor.

# Kommentierte Darstellung der Syntax

Das folgende Kapitel stellt die Syntax zur Generierung der Daten mit R vor. Der Code ist optional darstellbar und enthält detaillierte Kommentare.


## 0. Benötigte Pakete

Der Code nutzt im Wesentlichen die Pakete des Tidyverse. 

```{r Libraries und Pfade, message=FALSE, warning=FALSE}

library("tidyverse") # Tidyverse Methods
library("bookdown") 
library("readxl") # Read Excel
library("zoo")
library("imputeTS") # Impute Missing Features
library("haven") # write Stata-dta
library("sf") # write Stata-dta
library(pastecs) # descriptive stats

# Create Output directories in working directory if necessary
dir.create("Outfiles", showWarnings=F)
dir.create("Outfiles/2020", showWarnings=F)
dir.create("Outfiles/2020/Bund", showWarnings=F)
dir.create("Outfiles/2020/Other", showWarnings=F)
dir.create("Outfiles/2020/Stata", showWarnings=F)
```

## I.  Generierung eines ID-Datensatzes

Zunächst muss ein Datensatz generiert werden in dem den kleinsten regionalen Einheiten (Gemeinden) alle übergeordneten regionalen Einheiten und deren Regionalkennziffern zugeordnet werden. Datenquelle ist die Gebietsstandsreferenz von Destatis Stand 31.12.2015.
Es wird insbesondere geprüft, ob die Referenzdaten Missings auf den Regionalkennziffern oder den Namen der Gebietsstände aufweisen.
An diesen ID-Datensatz werden später die Indikatoren angespielt. 

```{r ID-Datensatz aus den Referenzdaten generieren , message=FALSE, warning=FALSE}

print_missings = function(data) {
  df = data[-1,]; 
  if(sum(is.na(df))>0){print("Missing observations: "); print(df[!complete.cases(df),])}; 
  df}
# Eine Funktion die später häufiger verwendet wird. Erklärung der Befehle:
# 1. Die erste Zeile des Datensatzes wird entfernt, da es sich dabei nicht um 
# eine Beobachtung, sondern nur um eine Variablenbeschreibung handelt
# 2. Wenn es im Datensatz fehlende Werte gibt, werden die dazugehörigen 
# Beobachtungen ausgegeben
# 3. Zuletzt wird der Datensatz aufgerufen damit er in Pipes weiterverarbeitet 
# werden kann

load_dataset = function(sheet) {
  suppressMessages(
    read_excel("Data/Referenz/Referenz_1998_2017.xlsx", sheet = sheet, na = "NA")
  )
}
# Eine Funktion die später häufiger verwendet wird. Sie dient dazu die einzelnen 
# Sheets aus Excel einzulesen. Außerdem werden Warnmeldungen beim Laden der 
# Daten unterdrückt. Zum Aufrufen der Funktion muss der Name des gewünschten
# Excel-Blattes als Argument angegeben werden.

Gemeinden_INKAR <- load_dataset("Gemeinden-GVB") %>% 
  print_missings() %>% na.omit() %>%
  mutate(Kennziffer=as.numeric(gem17),"Kennziffer Gemeindeverband"=vbgem17, fl17=as.numeric(fl17))
# Pipes: 
# 1. Wenn es fehlende Werte gibt, wird man hierdurch benachrichtigt
# 2. Gemeinden ohne fehlende Beobachtungen
# 3. Rename von zwei Variablen; " um Leerzeichen zu berücksichtigen 

Gemeindeverbaende_INKAR <- load_dataset("Gemeindeverbände") %>% 
  print_missings() %>% na.omit() %>% 
  select("Kennziffer Gemeindeverband"=gvb17,"Name des Gemeindeverbands"=gvb17name)
# Das ganze nochmal für Gemeindeverbände  
# Pipes: 
# 1. Wenn es fehlende Werte gibt, wird man hierdurch benachrichtigt
# 2. Missings herausfiltern
# 3. Nur die Variablen gvb17 und Name des Gemeindeverbands ausgewählt

Kreise_INKAR <- load_dataset("KRS") %>%
  print_missings() %>% na.omit() %>% 
  mutate(krs17= as.numeric(krs17)/1000, fl17 = as.numeric(fl17))
# ... und für Kreise
# Pipes: 
# 1. Wenn es fehlende Werte gibt, wird man hierdurch benachrichtigt
# 2. Missings herausfiltern
# 3. Neue Variable generieren, die die Kreisvariable auf Fünfsteller reduzieren

# Die drei Datensätze werden nun ausgehend vom Gemeindedatensatz zu einem ID-Datensatz zusammmengefügt
id_dataset <- Gemeinden_INKAR %>% 
              select(Gemeindekennziffer=Kennziffer,"Name der Gemeinde"=gem17name,"Kennziffer Gemeindeverband") %>% 
              mutate(Kreiskennziffer=floor(Gemeindekennziffer/1000)) %>%
              left_join(.,Kreise_INKAR %>% select("Kreiskennziffer"=krs17,
                                                  "Name des Kreises"=krs17name,
                                                  "Raumordnungsregion Nr"=ROR11,
                                                  Raumordnungsregion=ROR11name,
                                                  NUTS2,
                                                  "NUTS2 Name"=NUTS2name,
                                                  "Bundesland"=...28),by="Kreiskennziffer") %>%
              left_join(.,Gemeindeverbaende_INKAR, by="Kennziffer Gemeindeverband")
# Pipes:  1. (select) Variablenauswahl (gkz, Gemeindename, Gemeindeverband)
#         2. die Kreiskennziffer wird aus der Gemeindekennziffer generiert
#         3. leftjoin spielt Kreisdaten über Kreiskennziffer an
#         3.1 select wdhlt, die anzupielenden Variablen aus, darunter auch NUTS und ROR und Bundesland, dessen Variablenname beim Einlesen zu lang war (...24)
#         3.2 die Kreiskennziffer wurde vor dem leftjoin im Using-Datensatz generiert
#         4. als letztes werden die Gemeindeverbandskennziffern angespielt

```

## II. Erzeugen eines Datensatzes mit Kennziffern als ID unabhängig von der Ebene 

In diesem Code-Abschnitt werden die INKAR-Daten zu den Indikatoren in einem Datensatz zusammengeführt. Die Information für die Indikatoren, die für die Berechnung des GISD verwendet werden, liegt auf unterschiedlichen Ebenen vor. Die Faktorenanalysen sollen später auf Gemeindeebene durchgeführt werden, weshalb Information der Kreisebene an alle Gemeinden dieser Kreise angespielt wird. Percentile des Indexes können so später für jede regionale Ebene separat berechnet werden. Datenbasis sind die INKAR-Daten der jeweiligen Indikatoren im Excel-Format, die zu jeder Revision aus der INKAR-Datenbank heruntergeladen werden. Tabelle \@ref(tab:indicators) stellt die Indikatoren dar.

```{r Indikatoren einlesen und mit dem ID-Datensatz Zusammeführen, message=FALSE}

# Basis erzeugen: Ausgangspunkt Kreisdaten
# Es werden Indikatoren allen Ebenen angespielt, als erstes die Kreise.
Basedata <- Kreise_INKAR %>% select(Kennziffer=krs17) %>% mutate(Jahr=2017)
# Datensatz zum Anspielen der Daten generieren
# Ausgangspunkt Kreisdatensatz
# Pipes:  1. nur Kreiskennzifern ausgewählt
#         2. Jahresvariable generiert (2017)

# Liste der Variablen generieren
inputdataset <- list.files("Data/INKAR_1998_2017/") # Variablenliste der Dateinamen im Ordner

# Einlesen der einzelnen Excelfiles zu den Daten (Schleife) 
# for testing file<-inputdataset[1]
for(file in inputdataset){
  suppressMessages(myimport <- read_excel(paste0("Data/INKAR_1998_2017/",file), skip = 1, sheet = "Daten"))
  names(myimport)[1] <- "Kennziffer"
  myimport[2:3] <- NULL
  myimport <- myimport %>% gather(key = "Jahr", value = "Value" , -"Kennziffer", convert=T, na.rm = T) %>%
    mutate(Kennziffer=as.numeric(as.character(Kennziffer)), Value=as.numeric(Value))
  names(myimport)[3] <- strsplit(strsplit(file,"_")[[1]][2],"[.]")[[1]][1]
  Basedata <- full_join(Basedata, myimport, by=c("Kennziffer", "Jahr"))
}
# Schleife für jedes Excel-File
# 1. Einlesen der Exceldatei; jeweils das Sheet "Daten"; erste Zeile wird geskippt, die Daten werden als Text eingelesen
# 2. Die erste Splate wird als Kennziffer benannt
# 3. Die zweite und dritte Zeile werden gelöscht
# 4. Die Daten werde reshaped, um die Jahresinfos im langen Format zu speichern; convert konvertiert die Jahreszahlen zu Integern;
# rm.na entfert Zeilen, bei denen der "Value" fehlt; -"Kennziffer" sorgt dafür, dass jeder Kennziffer mehrere Jahre zugeordnet werden
# 5. von innen nach außen 
# 5.1 das innere strsplit(file, "_") teilt den Filenamen inkl. Dateiendung beim "_"; mit [[1]][2] wird das zweite Element davon ausgewählt
# 5.3 das äußere strsplit teilt dies dann beim ".", sodass nur noch der Dateiname übrig bleibt, der mit [[1]][1] ausgewählt wird
# 5.5 names(import)[3] nimmt diesen Dateinamen als Variablennamen für die dritte Spalte
# 6. Jedes file der Schleife wird an Basedata gejoint über Kennziffer und Jahr; full_join übernimmt dabei jede Zeile und Spalte jeder Seite,
# auch wenn die Werte auf einer Seite missing enthalten

rm(inputdataset) 


# Liste der Indikatoren erstellen
listofdeterminants <- names(Basedata)[3:length(Basedata)]

# Regionale Tiefe der Indikatoren 
ind_level <- c("Gemeindeverband","Gemeindeverband","Kreis", "Kreis", "Kreis", "Kreis", "Kreis", "Gemeindeverband", "Kreis", "Kreis")
level_table <- cbind(listofdeterminants,ind_level)
# Tabelle der Indikatoren mit regionaler Tiefe
ind_col = c("Indikator","Tiefe des Indikators")

# Datensatz für die Gemeindeverbandsebene generieren
Basedata_Gemeindeverbandsebene <- Basedata %>% select(Kennziffer,Jahr,Arbeitslosigkeit,Beschaeftigtenquote,Einkommensteuer) %>%   
  gather(key,value,3:5) %>% filter(!is.na(value)) %>% spread(key,value) %>% filter(Jahr>=1998) %>% rename("Gemeindeverband"=Kennziffer)
# Pipes:  1. Auswahl der Variablen 
#         2. Reshape der Daten von wide nach long      
#         3. Auswahl von Non-Missing 
#         4. Reshape von long nach wide 
#         5. Auswahl der Daten Jahr>=1998
#         6. Umbenennung der Kennziffervariable

# Datensatz für die Kreisebene generieren 
Basedata_Kreisebene <- Basedata %>% select(krs15=Kennziffer,Jahr,listofdeterminants) %>% 
  select(-Arbeitslosigkeit,-Einkommensteuer,-Beschaeftigtenquote) %>% rename(Kreis=krs15)
# Pipes:  1. neben der Kennziffer, die einen anderen Namen bekommt wird das Jahr und die Variablenliste ausgewählt
#         2. drei Variablen werden aus der Auswahl ausgeschlossen
#         3. die Kreisvariable wird in Kreis umbenannt, weil im nächsten Schritt Kreisinfos an die Gemeinden angespielt werden

# Hinweis: für die Gemeindeebene wird kein Basedata-Datensatz erstell, da es keine Infos auf der Gemeindeebene gibt

# Join different levels
# Nun werden die Daten bezogen auf die Ebenen gemergt
# Dazu wird erstmal ein Leerdatensatz im Longformat erstellt, der Fälle für alle Gemeinden für jedes Jahr generiert
Workfile <- as.data.frame(expand.grid("Kennziffer"=Gemeinden_INKAR %>% pull(Kennziffer),"Jahr"=seq(min(Basedata$Jahr):max(Basedata$Jahr)) + min(Basedata$Jahr)-1)) %>% mutate(Kreiskennziffer=floor(as.numeric(Kennziffer)/1000)) %>% as_tibble() %>%
   left_join(. , Gemeinden_INKAR,by=c("Kennziffer"))  %>%
   select(Gemeindekennziffer=Kennziffer,Kreis=Kreiskennziffer,Gemeindeverband="Kennziffer Gemeindeverband",Jahr,Bevoelkerung=bev17) %>% mutate(Gemeindeverband=as.numeric(Gemeindeverband), Bevoelkerung=as.numeric(Bevoelkerung)) %>% 
  arrange(Gemeindekennziffer,Jahr) %>% # Join Metadata
   left_join(. , Basedata_Kreisebene,by=c("Kreis","Jahr")) %>% # Hier wird über Kreis gematched
   left_join(. , Basedata_Gemeindeverbandsebene,by=c("Gemeindeverband","Jahr")) %>%  # Join Indicators for Level: Gemeindeverband 
   filter(Jahr>=1998)
# als erstes wird ein data.frame erzeugt (Workfile); der alle Gemeindewellen (1998-201x) in den Zeilen stehen hat
# 1. expand.grid erzeugt ein tibble mit allen Kombinationen von Kennziffern und Jahren
#     pull erzeugt einen Vektor für die Variablenwerte von Kennziffer aus dem Datensatz
#     + min(...) wird zu der Sequenz von Jahren aus dem Basedata addiert (1 bis X) damit auch Jahreswerte weitergeben werden[ist das nötig?]
#     stringAsFactors sorgt dafür, dass die Kennziffern nicht als Factors sondern als Strings geladen werden und es damit keine Probleme bei der Weiterverarbeitung gibt
# 2. mutate generiert eine Kreiskennziffer
# 3. as_tibble erzeugt einen tibble, damit left_join genutzt werden kann
# 4. erstes left_join spielt die Gemeindedaten über Kennziffer an, das geht so, weil Gemeinden_INKAR als tibble gespeichert ist
# 5. select, wählt die inhaltlichen Variablen aus, und ändert die Variablennamen
# 6. arrange im select sortiert nach Gemeindekennziffer und Jahr
# 7. zweites left_join spielt die Daten der Kreisebene via Kreis und Jahr an
# 8. drittes left_join spielt die Daten der Gemeindeverbandsebene via Gemeindeverband und Jahr an
# Notiz: . in den Befehlen bezieht sich auf den tibble bzw. data.frame der in der Pipe bearbeitet wird

rm(myimport)

# Stata-Datensatz rausschreiben
write_dta(Workfile, paste0("S:/OE/FG28/205 Regionale Unterschiede/GISD/Plausibilitätschecks/workfile.dta"))

# Ende Generierung Basisdatensatz
```

```{r indicators, echo=FALSE}
knitr::kable(level_table, col.names = ind_col, caption = "Liste der Indikatoren")
```

Es gibt noch einige Probleme bei der Auswahl der Indikatoren, die erst später zum Tragen kommen. Insbesondere der Bildungsindikator Schulabgänger ohne Abschluss macht Probleme, weil der nicht mit dem Anteil Beschäftigter mit akademischem Bildungsabschluss korreliert. 
Eine Betrachtung des alternativen Indikators Schulabgänger mit Hochschulreife zeigt, dass dieser besser mit dem Anteil Beschäftigter mit akademischem Bildungsabschluss korreliert, aber dafür recht stark negativ mit dem Anteil der Beschäftigten ohne Abschluss. Das kann daran liegen, dass dort wo der Anteil von Personen mit akademischem Abschluss hoch ist, auch der Druck für Schulabgänger einen Abschluss zu machen geringer ist.


## III.Imputation fehlender Werte

Das bisherige Imputationsmodell nutzt Arbeitslosigkeit als Prädiktoren. In den Daten für 2016 fehlen für diese Variable 6 Werte.
Der Einfachheit halber werden diese interpoliert. 

```{r Vereinzelte Missings auf den Imputationsvariablen interpolieren}

# Anzahl der Missings für die Indikatoren
missings_table = as.data.frame(expand.grid("Jahr"=1998:max(Basedata$Jahr)))
predictors_list = data.frame(Variable=character(), Missings=double(), stringsAsFactors = FALSE)
for (column in level_table[,1]){
  for (year in 1998:max(Basedata$Jahr)){
    missings_table[year-1997,column] = Workfile %>% filter(Jahr==year, Bevoelkerung>0, is.na(Workfile[,column])) %>% nrow()
  }
  predictors_list[nrow(predictors_list) + 1,] = c(column, Workfile %>% filter(Bevoelkerung>0, is.na(Workfile[,column])) %>% nrow())
}
predictors_list = predictors_list %>% mutate(Missings=as.integer(Missings))
predictors_list = predictors_list[order(predictors_list$Missings),]
predictors_list

Missing_on_Imputationsvars <- Workfile %>%  filter(Jahr>=1998, Bevoelkerung>0, is.na(Arbeitslosigkeit) | is.na(SchulabgaengermitHochschulreife))
Missing_on_Imputationsvars
# das betrifft nur Arbeitslosigkeit in 6 Gemeinden in 2016
# für Schulabgänger mit Hochschulreife 36 Gemeinden in 2016

# Fälle betrachten: Beispiel 5154028 
TimeSeries_for_Missing <- Workfile %>%  filter(Gemeindekennziffer==5154028 ) %>% select(Gemeindekennziffer, Jahr, Arbeitslosigkeit, SchulabgaengermitHochschulreife) %>% arrange(Gemeindekennziffer, Jahr)
TimeSeries_for_Missing

# Interpolation der fehlenden Werte über die Zeitreihe (Mittelwert: Vorjahr, Nachjahr)
Workfile <- Workfile %>%  filter(Jahr>=1998, Bevoelkerung>0) %>%  group_by(Gemeindeverband) %>% mutate(impu_arblos = na.approx(Arbeitslosigkeit), impu_mA = na.approx(SchulabgaengermitHochschulreife)) %>% select(-Arbeitslosigkeit, -SchulabgaengermitHochschulreife) %>% rename(Arbeitslosigkeit=impu_arblos, SchulabgaengermitHochschulreife=impu_mA) %>% ungroup()

# zweites Problem: drei Landkreise (Bamberg et al.) mit 0.0 auf SchulabgaengermitHochschulreife (kein Problem, wenn Variable erstmal nicht benutzt wird)
#   mutate(SchulabgaengermitHochschulreife = na_if(SchulabgaengermitHochschulreife, 0.0)

# Check der Interpolation
TimeSeries <- Workfile %>%  filter(Gemeindekennziffer==5154028) %>% select(Gemeindekennziffer, Jahr, Arbeitslosigkeit, SchulabgaengerohneAbschluss) %>% arrange(Gemeindekennziffer, Jahr)
TimeSeries


```

Die Missings auf der Arbeitslosigkeit finden sich in unterschiedlichen Gemeinden verschiedener Landkreise in NRW.

```{r}
# Anzahl der Missings über die Indikatoren
# summary(Workfile %>% select(all_of(listofdeterminants)))
# sapply(Workfile  %>% select(listofdeterminants) , function(x) sum(is.na(x)))

# Imputation
imputationsliste <- subset(listofdeterminants , !(listofdeterminants %in%
                                                    c('Arbeitslosigkeit','SchulabgaengermitHochschulreife','SchulabgaengerohneAbschluss')))
# Variablenliste für die Regressionsimputation wird erstellt
# das betrifft alle Variablen, außer die im angebenen Vektor
# letztere sind frei von Missings und können im Imputationsmodell genutzt werden 


Impdata <-  Workfile %>%  dplyr::filter(Jahr>=1998, Bevoelkerung>0) %>% 
  gather(key,value,6:15) %>% mutate(value=ifelse(value<0.00001,NA,value)) %>% spread(key,value)
# Imputationsdatensatz generieren: Jahr>=1998, Bevoelkerung>0 
# gather und spread identifiziern key-Variablen automatisch 
# es geht hier vor allem darum Werten<0 ein NA zuzuordnen
# vor allem auf der Variablen Schulabgänger mit Hochschulreife gibt es ebenfalls viele 0-Werte,
# Es ist möglich, dass es in manchen Kreisen keine Schulen mit gymnasialer Oberstufe gibt, 
# das sollte allerdings kein Grund sein, diese im GISD abzustufen. Deshalb werden auch diese Wert auf NA kodiert und 
# im Folgenden imputiert
 

# sapply(Impdata  %>% select(listofdeterminants) , function(x) sum(is.na(x)))
# Einige Missings basierten auf Gebietsständen ohne Bevölkerung, diese sind entfernt 


# Als erstes wird die Imputationsfunktion erstellt (hier werden noch keine Daten generiert)
# Impute_function (NOT FOR GROUPED DATA!)
my_ts_imputer <- function(data,outcome_name){
  mydata   <- data %>% group_by(Gemeindekennziffer) %>%
    select(Gemeindekennziffer,Jahr,Arbeitslosigkeit,SchulabgaengerohneAbschluss,"Outcome"=paste(outcome_name)) %>% 
    mutate(MEAN=mean(Outcome , na.rm=T)) %>% ungroup()
  mymodell <- lm(Outcome ~
                  I(Jahr*Jahr*MEAN) + I(Jahr*MEAN) + Arbeitslosigkeit + 
                   SchulabgaengerohneAbschluss,
                   data = mydata  , na.action="na.exclude")
  mydata %>% select(Outcome) %>% mutate(Imputed = predict(mymodell, newdata =mydata )) %>%
    mutate(Outcome=ifelse(is.na(Outcome),Imputed,Outcome)) %>% 
    mutate(Outcome=ifelse(Outcome<0,0,Outcome)) %>% pull(Outcome)
  }
# Hier wird eine Funktion generiert, die im Datensatz (data) fehlende Daten für ausgewählte Variablen (outcome_name) imputiert
# 1. zunächst werden Mittelwerte für das Outcome (siehe select) jeweils für die Gemeinde generiert, d.h. über alle Wellen aggregiert
# 2. mymodell definiert das Modell (lm); "I()" sichert ab, dass der Operator * erkannt wird und dass ein Spaltenvektor in die Formel eingeht
# 3. zweites mydata: es wird eine Variable Imputed generiert, die sich aus der prediction aus mymodell ergibt
#    während der vorherige Befehl (mymodell) die Koeffizienten generiert, werden nun auf Basis dieses Modells predictions generiert, 
#    und zwar auch für Fälle mit Missing auf den Outcomes
# 4. fehlende Werte in den Outcomes werden durch Werte auf der Variable Imputed ersetzt
# 5. Für einige Fälle erzeugt die prediction unplausible Werte (negative Outcomes), diese werden auf 0 gesetzt
# 6. pull kreiert einen Vektor (hier Variable Outcome), die im nächsten Befehl verwendet wird

# Test Function if necessary
# Impdata %>% mutate(Test=my_ts_imputer(.,"Bruttoverdienst")) %>% select(Gemeindekennziffer,Jahr,Bruttoverdienst,Test) %>% head()

Impdata.imputed <- Impdata %>% mutate(
  Beschaeftigtenquote=my_ts_imputer(.,"Beschaeftigtenquote"),
  Bruttoverdienst=my_ts_imputer(.,"Bruttoverdienst"),
  BeschaeftigtemitakadAbschluss=my_ts_imputer(.,"BeschaeftigtemitakadAbschluss"),
  BeschaeftigteohneAbschluss=my_ts_imputer(.,"BeschaeftigteohneAbschluss"),
  Einkommensteuer=my_ts_imputer(.,"Einkommensteuer"),
  Haushaltseinkommen=my_ts_imputer(.,"Haushaltseinkommen"),
  Schuldnerquote=my_ts_imputer(.,"Schuldnerquote"),
  SchulabgaengermitHochschulreife=my_ts_imputer(.,"SchulabgaengermitHochschulreife")    
  )
# hier wird der Datensatz mit den imputierten Werten generiert. Die Funktion my_ts_imputer wird auf jeden Indikator mit Missings angewendet

# Result of Imputation
summary(as.data.frame(Impdata.imputed) %>% ungroup()  %>% select(listofdeterminants))



# Percentile für ausgewählte Variablen bilden

Impdata.imputed <- Impdata.imputed %>% 
          mutate(Bruttoverdienst2 = findInterval(Bruttoverdienst, quantile(Bruttoverdienst, probs=0:100/100 , type=9)),
          Haushaltseinkommen2 = findInterval(Haushaltseinkommen, quantile(Haushaltseinkommen, probs=0:100/100 , type=9)),
          Einkommensteuer2 = findInterval(Einkommensteuer, quantile(Einkommensteuer, probs=0:100/100 , type=9)),
          SchulabgaengermitHochschulreife2 = findInterval(SchulabgaengermitHochschulreife, quantile(SchulabgaengermitHochschulreife, probs=0:100/100 , type=9)),
          SchulabgaengerohneAbschluss2 = findInterval(SchulabgaengerohneAbschluss, quantile(SchulabgaengerohneAbschluss, probs=0:100/100 , type=9)),
          BeschaeftigtemitakadAbschluss2 = findInterval(BeschaeftigtemitakadAbschluss, quantile(BeschaeftigtemitakadAbschluss, probs=0:100/100 , type=9))) 



# Stata-Datensatz rausschreiben
# write_dta(Impdata.imputed, paste0("Outfiles/2020/Stata/impdata.dta"))
```



```{r Histograms for List of Determinants, eval=FALSE, message=TRUE, include=FALSE, results="HIDE"}

library(ggplot2)
for(i in listofdeterminants){
hist_over_year <- ggplot(data = Impdata.imputed) + 
  geom_histogram(mapping =  aes_string(x = i)) + 
  facet_wrap(~Jahr)
print(hist_over_year)
}

hist_over_all <- ggplot(data = Impdata.imputed) + geom_histogram(mapping = aes_string(x = "BeschaeftigteohneAbschluss")) 
print(hist_over_all)


hist_over_westost <- ggplot() +
  geom_histogram(data = Impdata.imputed[Impdata.imputed[,"Kreis"]<11000,], aes_string(x = "BeschaeftigteohneAbschluss"), fill ='darkblue')  + 
  geom_histogram(data = Impdata.imputed[Impdata.imputed[,"Kreis"]>=11000,], aes_string(x = "BeschaeftigteohneAbschluss"), fill ='darkred')  
print(hist_over_westost)

hist_over_year_all <- ggplot(data = Impdata.imputed]) + geom_histogram(mapping = aes_string(x = "BeschaeftigteohneAbschluss")) + facet_wrap(~Jahr)
print(hist_over_year_all)

# Es gibt eine bimodale Verteilung bei den Beschäftigten ohne Abschluss, die in Ost- und Westdeutschland jedoch unimodal ist
```

## IV. Faktorenanalyse (Hauptkomponentenanalyse) inklusive Generierung der Faktorscores
```{r Tibbles für die Teilscores generieren, echo=FALSE}

# Variablenliste für die Faktorenanalyse 
print(listofdeterminants)
TS_Arbeitswelt <- Impdata.imputed  %>% ungroup() %>% dplyr::select(Beschaeftigtenquote,Arbeitslosigkeit,Bruttoverdienst) 
TS_Einkommen   <- Impdata.imputed %>% dplyr::select(Einkommensteuer,Haushaltseinkommen,Schuldnerquote) 
# für den Vergleich der Ergebnisse wird zunächst ein Datensatz für die Variablenauswahl der Revision 2019 generiert
TS_Bildung <- Impdata.imputed %>% dplyr::select(BeschaeftigtemitakadAbschluss,BeschaeftigteohneAbschluss,SchulabgaengerohneAbschluss) 
# Check dieser Lösung für das 2014er Sample 
TS_Bildung_r2014 <- Impdata.imputed %>% filter(Jahr<2015) %>%  dplyr::select(BeschaeftigtemitakadAbschluss,BeschaeftigteohneAbschluss,SchulabgaengerohneAbschluss) 

TS_Bildung_4items <- Impdata.imputed %>% dplyr::select(BeschaeftigtemitakadAbschluss,BeschaeftigteohneAbschluss,SchulabgaengerohneAbschluss, SchulabgaengermitHochschulreife) 

TS_Bildung_4items_without_BoA <- Impdata.imputed %>% dplyr::select(BeschaeftigtemitakadAbschluss,SchulabgaengerohneAbschluss, SchulabgaengermitHochschulreife) 

```

# Faktorenanalyse basierend auf Hauptkomponentenanalyse für jede der drei Subscalen
```{r PCA für die Teilscores, echo=FALSE}



# PCA für die Arbeitsweltdimension
TS_Arbeitswelt.pca <- prcomp(TS_Arbeitswelt, center = TRUE, scale. = TRUE, retx=TRUE)
TS_Arbeitswelt.pca
	# Option retx erzeugt rotierte Lösung

head(TS_Arbeitswelt.pca$sdev)
# nur die erste Komponente mit Eigenwert über 1
	# (prcomp gibt standardmäßig Sdev statt Varianz aus)
plot(TS_Arbeitswelt.pca)
	# screeplot - bei nur drei Variablen wird ein Balkendiagramm angezeigt
TS_Arbeitswelt.pca
# die Faktorladungen der drei Hauptkomponenten für Arbeitswelt 
# die Ladungen der ersten Komponente enstprechen der Erwartung

TS_Arbeitswelt.pca <- prcomp(TS_Arbeitswelt, center = TRUE, scale. = TRUE, retx=TRUE, rank. = 1)
# die Option rank erlaubt die Beschränkung der Anzahl an Komponenten (Faktoren)
TS_Arbeitswelt.pca

# PCA für die Einkommensdimension
TS_Einkommen.pca <- prcomp(TS_Einkommen, center = TRUE, scale. = TRUE, retx=TRUE) 
plot(TS_Einkommen.pca)
TS_Einkommen.pca <- prcomp(TS_Einkommen, center = TRUE, scale. = TRUE, retx=TRUE, rank. = 1) 
TS_Einkommen.pca

# PCA für die Bildungsdimension
TS_Bildung.pca <- prcomp(TS_Bildung, center = TRUE, scale. = TRUE, retx=TRUE) 
TS_Bildung.pca

# Alternativ Bildungskomponente mit BeschaeftigtemitakadAbschluss,SchulabgaengermitHochschulreife,SchulabgaengerohneAbschluss
TS_Bildung_new.pca <- prcomp(TS_Bildung, center = TRUE, scale. = TRUE, retx=TRUE, rank. = 1) 
# für die Bildung deutet die Analyse eher auf zwei Komponenten hin
# die Faktorladung für SchulabgaengerohneAbschluss ist auf dem zweiten Faktor schwach
# es wird die Komponente ausgewählt, bei der Beschaeftigte mit akad Abschluss positiv korreliert und 
# BeschaeftigteohneAbschluss und SchulabgaengerohneAbschluss negativ
# regionale Deprivation als Merkmal geringer Anteile von Akademikern bei gleichzeitigen hohen Anteilen 
# von Beschaeftigten ohne Abschluss und Schulabgaengern ohne Abschluss


# Check der Bildungskomponente in Revision 2018 (Daten für 2014)
TS_Bildung_r2014.pca <- prcomp(TS_Bildung_r2014, center = TRUE, scale. = TRUE, retx=TRUE) 
TS_Bildung_r2014.pca

TS_Bildung_4items.pca <- prcomp(TS_Bildung_4items, center = TRUE, scale. = TRUE, retx=TRUE) 
plot(TS_Bildung_4items.pca)
TS_Bildung_4items.pca


# für den JoHM Beitrag wird die Bildungsdimension mit den zwei Indikatoren generiert

# TS_Bildung_r2012.pca <- prcomp(TS_Bildung_r2012, center = TRUE, scale. = TRUE, retx=TRUE)
# TS_Bildung_r2012.pca

# Die Bildungsdimension wird als Index generiert: der Datensatz w|rde so aussehen

TS_Bildung <- TS_Bildung %>% mutate(z_BaA = scale(BeschaeftigtemitakadAbschluss),
                                           z_SoA = scale(SchulabgaengerohneAbschluss),
                                           Bildung_Index = (2*z_BaA - z_SoA))   %>% 
                                    cor(use="pairwise.complete.obs") 


# Es wurde außerdem eine Komponentenanalyse mit allen vier Bildungsindikatoren durchgeführt. 
# Aber auch hier bestand das Problem, der inkonsistenten Korrelationen zwischen den Teildimensionen.
```

# Nun wird die Generierung der Faktorscores vorbereitet.
```{r Generierung der Faktorscores, echo=FALSE}
# Componentoverview
GISD_Komponents <- cbind("Teildimension"="Arbeitswelt","Anteil"=TS_Arbeitswelt.pca$rotation^2,"Score"=TS_Arbeitswelt.pca$rotation) 
# cbind erstellt Spaltenvektoren mit den Infos aus Teildimension, den (rotierten) Faktorladungen und den Components; 
GISD_Komponents <- rbind(GISD_Komponents,cbind("Teildimension"="Einkommen","Anteil"=TS_Einkommen.pca$rotation^2,"Score"=TS_Einkommen.pca$rotation))

# Matrix der "Faktorladungen" für die Bildungsdimension erstellen
x <- c("Bildung", 0.66^2, 0.66)
y <- c("Bildung", (-0.33)^2, -0.33) 
row.TS_Bildung <- rbind(x, y) 
rownames(row.TS_Bildung) <- c("BeschaeftigtemitakadAbschluss", "SchulabgaengerohneAbschluss")
colnames(row.TS_Bildung) <- c("Teildimension", "PC1", "PC1")

# rbind erstellt Zeilenvektoren, diese werden hier in die bereits vorhandenen Spaltenvektoren eingebunden
GISD_Komponents <- rbind(GISD_Komponents,row.TS_Bildung) 
# auch für die Teildimension Bildung werden Zeilenvektoren eingebunden
GISD_Komponents <- cbind("Variables"=as.data.frame(rownames(GISD_Komponents)),as.data.frame(GISD_Komponents))
# als letztes wird die Matrix in einen Dataframe übersetzt

rownames(GISD_Komponents) <- NULL
# die überflüssigen Zeilennamen werden gestrichen
colnames(GISD_Komponents) <- c("Variable","Dimension","Anteil","Score")
# aussagekräftige Spaltennamen vergeben
GISD_Komponents$GISD <- "GISD"
# eine weitere Spalte wird eingef|gt mit dem String "GISD" in jeder Zeile
GISD_Komponents$Proportion <- round(as.numeric(as.character(GISD_Komponents$Anteil))*100,digits=1)
# eine weitere Spalte Proportion wird eingef|gt mit prozentualen Anteilswerten (eine Nachkommastelle)

# Hier findet die Prediction der Scores statt
Resultdataset <- Impdata.imputed
Resultdataset$TS_Arbeitswelt <- as.numeric(predict(TS_Arbeitswelt.pca, newdata = Impdata.imputed))
Resultdataset$TS_Einkommen <- as.numeric(predict(TS_Einkommen.pca , newdata = Impdata.imputed))
Resultdataset <- Resultdataset %>% mutate(TS_Bildung = as.numeric((2*scale(BeschaeftigtemitakadAbschluss) - scale(SchulabgaengerohneAbschluss))))


summary(Resultdataset %>% dplyr::select(TS_Arbeitswelt, TS_Einkommen, TS_Bildung))
descs <- stat.desc(Resultdataset[, -5])

# Korrelationen überprüfen
Resultdataset %>% dplyr::select(Arbeitslosigkeit,TS_Arbeitswelt,TS_Einkommen,TS_Bildung)  %>% cor( use="pairwise.complete.obs")  
# die Richtung der Skala der Scores ist nach der Generierung willkürlich 
# sie werden nun anhand der Variable Arbeitslosigkeit ausgerichtet,
# d.h. sie werden so gepolt, dass sie positiv mit Arbeitslosigkeit korrelieren, um Deprivation abzubilden:

if (cor(Resultdataset$Arbeitslosigkeit, Resultdataset$TS_Bildung,use="pairwise.complete.obs")<0) {
   Resultdataset$TS_Bildung <- Resultdataset$TS_Bildung*-1
   }
if (cor(Resultdataset$Arbeitslosigkeit, Resultdataset$TS_Arbeitswelt,use="pairwise.complete.obs")<0) {
  Resultdataset$TS_Arbeitswelt <- Resultdataset$TS_Arbeitswelt*-1
  }
if (cor(Resultdataset$Arbeitslosigkeit, Resultdataset$TS_Einkommen,use="pairwise.complete.obs")<0) {
  Resultdataset$TS_Einkommen <- Resultdataset$TS_Einkommen*-1
}

# Korrelationen erneut überprüfen
Resultdataset %>% dplyr::select(Arbeitslosigkeit,TS_Arbeitswelt,TS_Einkommen,TS_Bildung) %>% cor( use="pairwise.complete.obs")
# nun sind alle Korrelationen positiv
# wenngleich die Korrelation der Bildungsdimension mit Arbeitslosigkeit sehr gering ist
# inhaltlich ist das nicht unplausibel (hvhere Abiturquoten in strukturschwachen Regionen)

GISD_Komponents

# Tabelle der Komponenten mit den Anteilen ausgeben und gespeichert
save(GISD_Komponents, file="Outfiles/2019/GISD_Komponents.RData")
# [Warum das hier plötzlich?]

# Normalization
Resultdataset$TS_Arbeitswelt <- (Resultdataset$TS_Arbeitswelt -min(Resultdataset$TS_Arbeitswelt ))/(max(Resultdataset$TS_Arbeitswelt )-min(Resultdataset$TS_Arbeitswelt ))
Resultdataset$TS_Einkommen <- (Resultdataset$TS_Einkommen -min(Resultdataset$TS_Einkommen ))/(max(Resultdataset$TS_Einkommen )-min(Resultdataset$TS_Einkommen ))
Resultdataset$TS_Bildung <- (Resultdataset$TS_Bildung -min(Resultdataset$TS_Bildung ))/(max(Resultdataset$TS_Bildung )-min(Resultdataset$TS_Bildung ))


# GISD
Resultdataset$GISD_Score <- Resultdataset$TS_Arbeitswelt+Resultdataset$TS_Einkommen+Resultdataset$TS_Bildung
Resultdataset$GISD_Score <- (Resultdataset$GISD_Score -min(Resultdataset$GISD_Score ))/(max(Resultdataset$GISD_Score )-min(Resultdataset$GISD_Score ))

# Result
summary(Resultdataset %>% select(TS_Arbeitswelt,TS_Einkommen,TS_Bildung,GISD_Score))
str(Resultdataset %>% select(TS_Arbeitswelt,TS_Einkommen,TS_Bildung,GISD_Score))

# Teilscores und GISD-Score in Datensatz speichern
Resultdataset <- Resultdataset %>% select(Gemeindekennziffer,Jahr,Bevoelkerung,contains("TS_"),contains("GISD_Score"))


```

## V.  Datenexport - Erstellung der Datensätze 
```{r echo=FALSE}
# Merge IDs to Resultdataset
RawResult <- left_join(Resultdataset,id_dataset,by="Gemeindekennziffer")



exportlist<- NULL
exportlist$Kennziffern <- c("Gemeindekennziffer","Kreiskennziffer","Kennziffer Gemeindeverband","Raumordnungsregion Nr","NUTS2")
exportlist$Namen <- c("Name der Gemeinde","Name des Kreises","Name des Gemeindeverbands","Raumordnungsregion","NUTS2 Name")
exportlist$Label <- c("Gemeinde","Kreis","Gemeindeverband","Raumordnungsregion","NUTS2")
# exportlist$Kennziffern <- c("Gemeindekennziffer") # for testing
  

# Es folgt eine sehr lange Schleife
# für alle Regionalkennziffern (siehe Vektor) werden Datensätze generiert und in Ordnern abgelegt
for(mykennziffer in exportlist$Kennziffern) {
  myname <-  exportlist$Namen[exportlist$Kennziffern==mykennziffer]
  mylabel<-  exportlist$Label[exportlist$Kennziffern==mykennziffer]
  print(paste("Level:",myname,"Label:",mylabel))
  
  # Datensatzerstellung
  outputdata <- RawResult 
  
  outputdata$Group <- outputdata[[mykennziffer]]
  mergedataset  <- outputdata %>% dplyr::select(ID=mykennziffer,myname,Bundesland) %>% 
    group_by(ID) %>% filter(row_number()==1) %>% ungroup() 
  names(mergedataset)[1]=mykennziffer
  
  # Aggregation
  outputdata.agg <- outputdata %>% 
    group_by(Group,Jahr) %>% 
    dplyr::select(Group,Jahr,"Bevoelkerung",GISD_Score, TS_Bildung, TS_Einkommen, TS_Arbeitswelt) %>% 
    summarise(GISD_Score = weighted.mean(GISD_Score, Bevoelkerung), 
              TS_Bildung = weighted.mean(TS_Bildung, Bevoelkerung), 
              TS_Einkommen = weighted.mean(TS_Einkommen, Bevoelkerung),
              TS_Arbeitswelt = weighted.mean(TS_Arbeitswelt, Bevoelkerung),
              Bevoelkerung = sum(Bevoelkerung))
  
  
  
  # hier werden die bevoelkerungsgewichteten Mittelwerte über die regionalen Einheiten gebildet
  # Achtung: Referenzrahmen für den Bevölkerungsstand ist das Referenzjahr. Die Varianz der Bevölkerung über die Jahre wird nicht berücksichtigt.
  
  # Daten bereinigen
  names(outputdata.agg)[1] <- mykennziffer
  outputdata.agg <- merge(outputdata.agg,mergedataset,by=mykennziffer) %>%  
    dplyr::select(mykennziffer,myname,Jahr,Bundesland,"Bevoelkerung",GISD_Score, TS_Bildung, TS_Einkommen, TS_Arbeitswelt) %>%
    group_by(Jahr) %>% as_tibble()
  
  # Rekodierung
  # hier wird der GISD-Score neu normalisiert und die Quintile gebildet
  outputdata.agg <- outputdata.agg %>%  mutate(GISD_Score = round((GISD_Score -min(GISD_Score ))/(max(GISD_Score )-min(GISD_Score )), digits=6)) %>% 
                                        group_by(Jahr) %>% mutate(GISD_5 = findInterval(GISD_Score, quantile(GISD_Score,   probs=0:5/5 , type=9)),
                                        GISD_5 = findInterval(GISD_5, c(1:5)),
                                        GISD_10 = findInterval(GISD_Score, quantile(GISD_Score, probs=0:10/10 , type=9)),
                                        GISD_10 = findInterval(GISD_10, c(1:10)),
                                        GISD_k = findInterval(GISD_5, c(1,2,5))) %>% ungroup()
                                       
  summary(outputdata.agg %>% select(contains("GISD")))  
                                       
  
  # Aktuelles Referenzmodell 

  
  # Ausgabe Bund
  dir.create("Outfiles/", showWarnings=F)
  dir.create("Outfiles/2020/", showWarnings=F)
  dir.create("Outfiles/2020/Bund/", showWarnings=F)  
  dir.create(paste0("Outfiles/2020/Bund/",mylabel), showWarnings=F)  
  mydata <- outputdata.agg %>% ungroup() %>% dplyr::select(-Bundesland)
  write.csv(mydata, paste0("Outfiles/2020/Bund/",mylabel,"/",mylabel,".csv"))
  
  names(mydata) <- gsub("\\.","_",make.names(names(mydata)))
  names(mydata) <- gsub("\\?","oe",names(mydata))
  names(mydata) <- gsub("\\?","ae",names(mydata))
  names(mydata) <- gsub("\\?","ue",names(mydata))
  names(mydata) <- gsub("\\?","ss",names(mydata))
  write_dta(mydata, paste0("Outfiles/2020/Bund/",mylabel,"/",mylabel,"_long.dta"))
  
  # Ausgabe Bundeslandspezifisch ohne Stadtstaaten und nur für Ebenen Kreis und Gemeindeverband
  if (mylabel %in% c("Gemeindeverband","Kreis")) {
  outputdata.agg <- outputdata.agg %>% ungroup() %>% filter(!(Bundesland %in% c("Bremen","Hamburg","Berlin"))) %>% dplyr::select(-GISD_k,-GISD_5,-GISD_10)   %>% group_by(Jahr,Bundesland) 
  
  # Rekodierung Bundesland
  outputdata.agg <- outputdata.agg %>%  mutate(GISD_Score = round((GISD_Score -min(GISD_Score ))/(max(GISD_Score )-min(GISD_Score )), digits=6)) %>% 
                                        group_by(Jahr) %>% mutate(GISD_5 = findInterval(GISD_Score, quantile(GISD_Score,   probs=0:5/5 , type=9)),
                                        GISD_5 = findInterval(GISD_5, c(1:5)),
                                        GISD_10 = findInterval(GISD_Score, quantile(GISD_Score, probs=0:10/10 , type=9)),
                                        GISD_10 = findInterval(GISD_10, c(1:10)),
                                        GISD_k = findInterval(GISD_5, c(1,2,5))) %>% ungroup()
                         
  summary(outputdata)


  # Ausgabe Bundesländer
  ListeBula <- unique(outputdata$Bundesland)
  dir.create("Outfiles/2020/Bundesland", showWarnings=F)  
  for(myland in ListeBula) {
  dir.create( paste0("Outfiles/2020/Bundesland/",myland), showWarnings=F)  
    dir.create( paste0("Outfiles/2020/Bundesland/",myland,"/",mylabel), showWarnings=F)  
    mydata <- outputdata %>% filter(Bundesland==myland) %>% ungroup() %>% dplyr::select(-Bundesland)
    write.csv(mydata, paste0("Outfiles/2020/Bundesland/",myland,"/",mylabel,"/",mylabel,".csv"))
    
    mydata <- outputdata %>% filter(Bundesland==myland)
    names(mydata) <- gsub("\\.","_",make.names(names(mydata)))
    names(mydata) <- gsub("\\?","oe",names(mydata))
    names(mydata) <- gsub("\\?","ae",names(mydata))
    names(mydata) <- gsub("\\?","ue",names(mydata))
    names(mydata) <- gsub("\\?","ss",names(mydata))
    write_dta(mydata, paste0("Outfiles/2020/Bundesland/",myland,"/",mylabel,"/",mylabel,".dta"))
  }
  }  
}

```


## VI.  Datensätze für PLZ generieren

```{r eval=FALSE, include=FALSE}
#### PLZ Daten noch nicht gecheckt (nm) ### 

# Output Postcode Data
load("Data/SHP/GEM_Zipcode_Intersections_2015.RData") # AGS/Postcode-Intersections-Dataset in sf format


for (mykennziffer in c("PLZ2","PLZ3","PLZ4","PLZ5")) {
  myname <-  paste0(mykennziffer)
  mylabel<-  paste0(mykennziffer)
  print(paste("Level:",myname,"Label:",mylabel))
  
  # Datensatzerstellung # weighted.mean fehlt wg. Fehler Evaluation error: 'x' and 'w' must have the same length
  outputdata <- Resultdataset 
  outputdata <- outputdata %>% dplyr::select(AGS=Gemeindekennziffer,Jahr,GISD_Score)
  outputdata <- left_join(as.data.frame(PLZ.df) %>% ungroup() %>% mutate(AGS=as.numeric(as.character(AGS))),
                          outputdata,by=c("AGS"), all.x = TRUE)
  outputdata <- outputdata %>% filter(!is.na(mykennziffer) & !is.na(EW_Area) & !is.na(Jahr) & EW_Area>0)
  mycol <- which(mykennziffer %in% names(outputdata))
  outputdata <- outputdata %>% group_by(Jahr,AGS) 
  outputdata <- outputdata %>% mutate(GISD_Score = weighted.mean(GISD_Score,EW_Area))
  names(outputdata)[names(outputdata)=="Jahr"]<- "JAHR" # Seltsames Problem Name "Jahr"
  outputdata <- outputdata %>% group_by_at(vars("JAHR",mykennziffer)) %>% 
    summarise(GISD_Score = weighted.mean(GISD_Score,EW_Area), Bevölkerung = sum(EW_Area)) %>%
    group_by(JAHR)
  
  outputdata <- outputdata %>%  mutate(GISD_Score = round((GISD_Score -min(GISD_Score ))/(max(GISD_Score )-min(GISD_Score )), digits=6),
                                       GISD_5 = findInterval(GISD_Score, quantile(GISD_Score,   probs=0:5/5 , type=9)),
                                       GISD_5 = findInterval(GISD_5, c(1:5)),
                                       GISD_10 = findInterval(GISD_Score, quantile(GISD_Score, probs=0:10/10 , type=9)),
                                       GISD_10 = findInterval(GISD_10, c(1:10)),
                                       GISD_k = findInterval(GISD_5, c(1,2,5))) 
  summary(outputdata)            
  head(outputdata)
  ListeJahre <- unique(outputdata$JAHR)
  dir.create( paste0("Revisions/"), showWarnings=F)
  dir.create( paste0("Revisions/2020/"), showWarnings=F) 
  dir.create( paste0("Revisions/2020/Bund/"), showWarnings=F) 
  dir.create( paste0("Revisions/2020/Bund/",mylabel), showWarnings=F) 
  mydata <- outputdata %>% ungroup() 
  write.csv2(mydata, paste0("Revisions/2020/Bund/",mylabel,"/",mylabel,".csv"))
  mydata <- outputdata %>% ungroup() 
  names(mydata) <- gsub("\\.","_",make.names(names(mydata)))
  names(mydata) <- gsub("\\?","oe",names(mydata))
  names(mydata) <- gsub("\\?","ae",names(mydata))
  names(mydata) <- gsub("\\?","ue",names(mydata))
  names(mydata) <- gsub("\\?","ss",names(mydata))
  write_dta(mydata, paste0("Revisions/2020/Bund/",mylabel,"/",mylabel,"_long.dta"))
  }


```

# Allgemeine SOP für die Revision (nach Lars Kroll)
1. Neue Daten und Gebietsstände aus der INKAR-Datenbank herunterladen. Variablennamen und Formate überprüfen.
2. Postleitzahlen in GISD_generate_postcodes.R anhand der Gebietsstandsdatei überprüfen.
3. GISD_Generate.R ausführen

# Anknüpfungungspunkte für eine grundsätzliche Überarbeitung der GISD-Generierung 

Es gibt gute Gründe dafür am Konzept Bildung, Einkommen und Arbeitsweltindikatoren im GISD zu vereinen, auch wenn die Korrelation der Teildimensionen mit Einzelindikatoren der anderen Teildimensionen nur gering korrelieren. 
Es gibt andererseits Möglichkeiten den GISD weiter zu verbessern. Einzelne Schwachstellen sollen hier kurz erwdähnt werden.

1. Missing Data 
* Hoher Anteil an Missing Data in den frühen Wellen
* Umgang mit Missing Data kann verbessert werden
2. Faktorenanalyse
* Bisher wird die Faktorenanalyse per pcf-Verfahren durchgeführt. Hier wäre zu prüfen, ob Common Factor-Verfahren oder konfirmatorische Faktorenanalyse zu einer Verbesserung führen könnten.
3.Indikatorenauswahl 
* Die Struktur der Faktorladungen der Bildungsindikatoren ist nicht robust gegenüber Datenschwankungen. Der erste Faktor bildet die intendierte Kompomente ab. Es gibt einen zweiten Faktor mit Eigenwert über 1. Die Gewichte der Faktorladungen der Indikatoren BeschaeftigteohneAbschluss und Schulgaengerohneabschluss variieren sehr stark zwischen den Revisionen 2018 und 2019. Hier könnte man über eine andere Auswahl von Indikatoren nachdenken. Die bisherigen Indikatoren BeschaeftigteohneAbschluss und BeschaeftigtemitHochschulabschluss dieser Teildimension weisen die höchsten Anteile an MissingData auf (75%). Zudem wurde bisher noch nicht berücksichtigt, dass die  zwischenzeitliche Verkürzung der Schulzeit für das Abitur (G8 Reform) und die spätere Rücknahme dieser Reform in einigen Bundesldndern im Untersuchungszeitraum zu statistischen Artefakten in den Schulabgdngerquoten führt.
4. Methodologische Grundlagen
* Diskussion der dem Messmodell zugrunde liegende Kausalmechanismen 
Dimensionen sozioökonomischer Deprivation auf räumlicher Ebene: Einkommen, Arbeitswelt, Bildung
- Einkommen: (HH-Einkommen, Steueraufkommen, Schuldnerquote)
- Kaufkraft berücksichtigen, Vermögen berücksichtigen
- betrifft Handlungsspielräume der Kommunen, Proxy für Wirtschaftskraft der Kommunen


```{r eval=FALSE, include=FALSE}

# Hier wird eine Faktoranalyse für die Bildungskomponente durchgeführt, ohne die
# Beschäftigung ohne Abschluss zu verwenden, welche sich systematisch zwischen
# West und Ost unterscheidet
TS_Bildung_4items_without_BoA <- Impdata.imputed %>% dplyr::select(BeschaeftigtemitakadAbschluss,SchulabgaengerohneAbschluss, SchulabgaengermitHochschulreife) 
TS_Bildung_4items_without_BoA.pca <- prcomp(TS_Bildung_4items_without_BoA, center = TRUE, scale. = TRUE, retx=TRUE) 
TS_Bildung_4items_without_BoA.pca

GISD_Komponents <- cbind("Teildimension"="Arbeitswelt","Anteil"=TS_Arbeitswelt.pca$rotation^2,"Score"=TS_Arbeitswelt.pca$rotation) 
GISD_Komponents <- rbind(GISD_Komponents,cbind("Teildimension"="Einkommen","Anteil"=TS_Einkommen.pca$rotation^2,"Score"=TS_Einkommen.pca$rotation))
# Code wie oben
GISD_Komponents <- rbind(GISD_Komponents,cbind("Teildimension"="Bildung","Anteil"=TS_Bildung_4items_without_BoA.pca$rotation^2,
                                               "Score"=TS_Bildung_4items_without_BoA.pca$rotation))
# Neuer Code

GISD_Komponents <- cbind("Variables"=as.data.frame(rownames(GISD_Komponents)),as.data.frame(GISD_Komponents))
# als letztes wird die Matrix in einen Dataframe übersetzt

rownames(GISD_Komponents) <- NULL
colnames(GISD_Komponents) <- c("Variable","Dimension","Anteil","Score")
GISD_Komponents$GISD <- "GISD"
GISD_Komponents$Proportion <- round(as.numeric(as.character(GISD_Komponents$Anteil))*100,digits=1)
Resultdataset <- Impdata.imputed
Resultdataset$TS_Arbeitswelt <- as.numeric(predict(TS_Arbeitswelt.pca, newdata = Impdata.imputed))
Resultdataset$TS_Einkommen <- as.numeric(predict(TS_Einkommen.pca , newdata = Impdata.imputed))
# Code wie oben
Resultdataset$TS_Bildung <- as.numeric(predict(TS_Bildung_4items_without_BoA.pca , newdata = Impdata.imputed))
# Neuer Code
summary(Resultdataset %>% dplyr::select(TS_Arbeitswelt, TS_Einkommen, TS_Bildung))

if (cor(Resultdataset$Arbeitslosigkeit, Resultdataset$TS_Bildung,use="pairwise.complete.obs")<0) {
   Resultdataset$TS_Bildung <- Resultdataset$TS_Bildung*-1
   }
if (cor(Resultdataset$Arbeitslosigkeit, Resultdataset$TS_Arbeitswelt,use="pairwise.complete.obs")<0) {
  Resultdataset$TS_Arbeitswelt <- Resultdataset$TS_Arbeitswelt*-1
  }
if (cor(Resultdataset$Arbeitslosigkeit, Resultdataset$TS_Einkommen,use="pairwise.complete.obs")<0) {
  Resultdataset$TS_Einkommen <- Resultdataset$TS_Einkommen*-1
}
Resultdataset$TS_Arbeitswelt <- (Resultdataset$TS_Arbeitswelt -min(Resultdataset$TS_Arbeitswelt ))/(max(Resultdataset$TS_Arbeitswelt )-min(Resultdataset$TS_Arbeitswelt ))
Resultdataset$TS_Einkommen <- (Resultdataset$TS_Einkommen -min(Resultdataset$TS_Einkommen ))/(max(Resultdataset$TS_Einkommen )-min(Resultdataset$TS_Einkommen ))
Resultdataset$TS_Bildung <- (Resultdataset$TS_Bildung -min(Resultdataset$TS_Bildung ))/(max(Resultdataset$TS_Bildung )-min(Resultdataset$TS_Bildung ))
Resultdataset$GISD_Score <- Resultdataset$TS_Arbeitswelt+Resultdataset$TS_Einkommen+Resultdataset$TS_Bildung
Resultdataset$GISD_Score <- (Resultdataset$GISD_Score -min(Resultdataset$GISD_Score ))/(max(Resultdataset$GISD_Score )-min(Resultdataset$GISD_Score ))
Resultdataset <- Resultdataset %>% select(Gemeindekennziffer,Jahr,Bevoelkerung,contains("TS_"),contains("GISD_Score"))
summary(Resultdataset %>% select(TS_Arbeitswelt,TS_Einkommen,TS_Bildung,GISD_Score))
# Code wie oben

# Danach muss der Code oben beginnend beim Datenexport angewandt werden (Achtung! Oben erstellte Exporte werden dann überschrieben)
# Dafür wird das Working Directory temporär geändert
dir.create( paste0("Meik_3Bildungsindikatoren/"), showWarnings=F)
setwd("Meik_3Bildungsindikatoren")
# An dieser Stelle dann die Datenexport und PLZ Chunks laufen lassen
```

```{r eval=FALSE, include=FALSE}
# Das working directory wieder zurücksetzen
setwd("..") 
getwd()

```

```{r eval=FALSE, include=FALSE}

# Und hier noch eine Faktoranalyse mit allen 9 Komponenten ohne BeschaeftigteohneAbschluss
TS <- Impdata.imputed  %>% ungroup() %>% dplyr::select(all_of(listofdeterminants), -BeschaeftigteohneAbschluss) 
TS.pca <- prcomp(TS, center = TRUE, scale. = TRUE, retx=TRUE, rank. = 2) 


```