#############################################################################################################################################################################

###### Repliationscode zur GESIS Survey Guideline
###### Sand, Matthias und Kunz, Tanja (2020). Gewichtung in der Praxis. Mannheim, GESIS – Leibniz-Institut für Sozialwissenschaften (GESIS Survey Guidelines).
###### DOI: 10.15465/gesis-sg_030
###### Teil 2: Kalibrierungsgewichte für ALLBUS (Kapitel 2.2)

#############################################################################################################################################################################

rm(list = ls())
library(foreign)
library(xlsx)
library(doParallel)
library(missForest)
library(dplyr)
library(reshape2)
library(readstata13)

registerDoParallel(cores = 4)

Pfad <- paste0(getwd(),"/")
ALLBUS <- read.dta13(paste0(Pfad,"ZA5250_v2-1-0.dta"),convert.factors = T,encoding = "UTF8",missing.type = T)

## Aus den Daten der Bevölkerungsfortschreibung gemeinsame Verteilung und Margins bilden
Zens <- read.xlsx("Kopie von 12411-0012.xlsx",header = T,sheetIndex = 2,encoding = "UTF8")
row.names(Zens) <- Zens$Alter
Zens2 <- melt(Zens,id.vars = "Alter")
Zens2$rel <- Zens2$value/sum(Zens2$value)
Zens2_Alter <- tapply(Zens2$value,Zens2$Alter,function(x) sum(x)/sum(Zens2$value))
Zens2_BL <- tapply(Zens2$value,Zens2$variable,function(x) sum(x)/sum(Zens2$value))

### Daten vorbereiten ###

Allbus2 <- data.frame(respid = ALLBUS$respid,agec = ALLBUS$agec, land = ALLBUS$land, wghtpew = ALLBUS$wghtpew)
Allbus2$land[Allbus2$land == "EHEM. BERLIN-OST"] <- "EHEM. BERLIN-WEST"
levels(Allbus2$land)[levels(Allbus2$land) == "EHEM. BERLIN-WEST"] <- "Berlin"
s <- which(Allbus2$agec == "NICHT GENERIERBAR")
Allbus2 <- Allbus2[-s,] # Da nur zu Demonstrationszwecken verwendet, werden NAs einfach gelöscht anstelle einer ordentlichen Imputation
Allbus2 <- droplevels(Allbus2)
Allbus2$land <- factor(Allbus2$land,levels(Allbus2$land)[c(8,9,11,12,4,2,6,13,3,5,7,10,14,15,1,16)])


### Daten Aggregieren ###

Allb <- aggregate(Allbus2$wghtpew, by = list(agec = Allbus2$agec,land = Allbus2$land), sum)


##### Kalibrierungsmethoden im Vergleich

### 1. Poststratifikation ###

## Da einige Zellen unbesetzt sind, muss weiter aggregiert werden. Dies betrifft einige BLs gekreuzt mit 89+
## Daher werden 75 - 89 und  älter 89 zusammengefasst

Allbus3 <- Allbus2
Allbus3$agec[Allbus3$agec == "UEBER 89 JAHRE"] <- "75-89 JAHRE"
Allbus3 <- droplevels(Allbus3)
Allb2 <- aggregate(Allbus3$wghtpew, by = list(agec = Allbus3$agec,land = Allbus3$land), sum)
Allb2$rel <- Allb2$x/sum(Allb2$x)
## Auch für Fortschreibung ###
Zens3 <- Zens
Zens3[5,] <- Zens3[5,]+Zens3[6,]
Zens3 <- Zens3[-6,]
Zens3$Alter[5] <- factor("75 +")
Zens4 <- melt(Zens3,id.vars = "Alter")
Zens4$rel <- Zens4$value/sum(Zens4$value)


### Gewichtung -> Überprüfen, ob Anordnung gleich ist ###


Allb2$weight <- Zens4$rel/Allb2$rel


### 2. Raking ###

## Hier gehen die ursprünglichen Kategorien

w_0 <- rep(1,times = nrow(Allb))
times <- 0

while(times <= 10000){
  Allb$d <- with(Allb, x*w_0)
  #### Anpassung nach Alter #######
  Allb$d <- with(Allb, x*w_0)
  agec <- tapply(Allb$d,Allb$agec,function(x) sum(x)/sum(Allb$d))
  w_1 <- Zens2_Alter/agec
  w_1 <- w_1[as.numeric(Allb$agec)]
  ### Anpassen nach Bundesland ####
  Allb$d2 <- with(Allb, x*w_0*w_1)
  land <- tapply(Allb$d2,Allb$land,function(x) sum(x)/sum(Allb$d2))
  w_2 <- Zens2_BL/land
  w_2 <- w_2[as.numeric(Allb$land)]
  ### final ###
  w_3 <- w_0*w_1*w_2
  w_0 <- w_3
  ###Break
  if(!all(agec==Zens2_Alter)&
     !all(land==Zens2_BL)
  )
  {times<-times+1}
  else {break} 
  cat("iteration",times-1,"\n")
} #iteration 8

Allb$w0 <- w_0

### 3. GREG  ###

Allbus4 <- ALLBUS
Allbus4$land[Allbus4$land == "EHEM. BERLIN-OST"] <- "EHEM. BERLIN-WEST"
levels(Allbus4$land)[levels(Allbus4$land) == "EHEM. BERLIN-WEST"] <- "Berlin"
s <- which(Allbus4$agec == "NICHT GENERIERBAR")
Allbus4 <- Allbus4[-s,] # Da nur zu Demonstrationszwecken verwendet, werden NAs einfach gelöscht anstelle einer ordentlichen Imputation
Allbus4 <- droplevels(Allbus4)
Allbus4$land <- factor(Allbus4$land,levels(Allbus4$land)[c(8,9,11,12,4,2,6,13,3,5,7,10,14,15,1,16)])


Allbus4$pii <- nrow(Allbus4)/67626000 * 1/Allbus4$wghtpew

pii <- Allbus4$pii
XS <- model.matrix(~Allbus4$agec+Allbus4$land)
XS2 <- XS[,-1]
Zens2_Alter * sum(Zens2$value)
#   18 - 29  30 - 44  45 - 59  60 - 74  75 - 89    90 +  
#   11691962 15101099 19467345 13499518  8541870   749597 
Zens2_BL * sum(Zens2$value)
# Baden.WÃ.rttemberg                 Bayern                 Berlin            Brandenburg 
#             9096555               10808692                2997278                2112080 
# Bremen                Hamburg                 Hessen Mecklenburg.Vorpommern 
# 570376                1513250                5177556                1370491 
# Niedersachsen    Nordrhein.Westfalen        Rheinland.Pfalz               Saarland 
#       6615277               14906990                3410015                 851242 
# Sachsen         Sachsen.Anhalt     Schleswig.Holstein             ThÃ.ringen 
# 3456717                1917048                2409422                1838402 

TS <- c(69051391,15101099, 19467345, 13499518,8541870, 749597,10808692,2997278,2112080,570376,1513250,5177556,1370491,
        6615277,14906990,3410015,851242,3456717,1917048,2409422,1838402) # Model Matrix nimmt erste Ausprägung als Ref-Value +
# intercept = 1 ->  erste Kat weg und am Anfang die Populationsgöße
greg.unit <- function(XS,TX,piiS){
  yS       <- rep(1,length(piiS))
  XS.tilde <- colSums(XS/piiS)
  Bhat     <- solve(t(XS)%*%(XS/piiS))%*%t(XS)
  greg     <- 1/piiS*(1+t(Bhat)%*%(TX-XS.tilde))
  return(greg)
}
gr<-greg.unit(XS,TS,pii)

Allbus4$gr <- gr[,1]
Allbus4$agec2 <- Allbus4$agec
Allbus4$agec2[Allbus4$agec2 == "UEBER 89 JAHRE"] <- "75-89 JAHRE"
Allbus4$ID2 <- paste0(Allbus4$agec2,Allbus4$land)
Allbus4$ID3 <- paste0(Allbus4$agec,Allbus4$land)
Allb2$ID2 <- paste0(Allb2$agec,Allb2$land)
Allb$ID3 <- paste0(Allb$agec,Allb$land)
Allbus5 <- merge(Allbus4,Allb2[,c("ID2","weight")],by = "ID2",all.x = T)
Allbus5 <- merge(Allbus5,Allb[,c("ID3","w0")],by = "ID3",all.x = T)