---
title: "R Notebook"
output: html_notebook
---

# Packages

```{r}
# install.packages("fastglm")
library(survey)
library(glmnet)
library(tidyverse)
library(xtable)
library(future)
library(vcd)
library(doFuture)
library(fastglm)

source("../codes/adalasso-chen-phd.R")
```

# Data

```{r}
totals <- readRDS("../data/gus-woj-sek-boot-totals.rds") %>%
  filter(substr(kod2,1,1) != 6) %>%
  mutate(kod2 = ifelse(kod2 == 95, 96, kod2),
         kod2 = ifelse(kod2 == 92, 91, kod2),
         kod2 = ifelse(kod2 == 53, 54, kod2),
         #kod1 = as.character(kod1),
         kod2 = as.character(kod2),
         rok = as.character(rok)) 

head(totals)

final_data <- readRDS("../data/bkl-final.rds") %>%
  filter(zawod1 != 6, rok!=2012, nace %in% unique(totals$sekcja)) %>%
  mutate(zawod2 = ifelse(zawod2 == 95, 96, zawod2),
         zawod2 = ifelse(zawod2 == 92, 91, zawod2),
         zawod2 = ifelse(zawod2 == 99, 96, zawod2),
         zawod2 = ifelse(zawod2 == 53, 54, zawod2),
         zeros = str_extract(zawod6, "0+$"),
         zeros = str_count(zeros, "0"),
         zeros = ifelse(is.na(zeros), 6, 6-zeros),
         rok = as.character(rok)) %>%
  rename(kod1 = zawod1,
         kod2 = zawod2,
         kod6 = zawod6,
         sekcja = nace) %>%
  filter(zeros > 1) 
```


# Correlations


```{r}
final_data %>%
  select(kod2, woj, sekcja, komp_techniczne:komp_biurowe) %>%
  gather(comps, vals, komp_techniczne:komp_biurowe) -> for_corr

for_corr  %>%
  count(kod2, comp=comps, vals) %>%
  group_by(comp) %>%
  do(occupancy = xtabs(n~vals + kod2, data = .) %>% assocstats(.) %>% .$cramer) %>%
  unnest() %>%
  left_join(
    for_corr  %>%
  count(sekcja, comp=comps, vals) %>%
  group_by(comp) %>%
  do(NACE = xtabs(n~vals + sekcja, data = .) %>% assocstats(.) %>% .$cramer) %>%
  unnest()) %>%
  left_join(for_corr  %>%
  count(woj, comp=comps, vals) %>%
  group_by(comp) %>%
  do(Voivodeship = xtabs(n~vals + woj, data = .) %>% assocstats(.) %>% .$cramer) %>%
  unnest()) %>%
  mutate(comp = case_when(comp == "komp_kulturalne" ~ "Artistic",
                            comp == "komp_dyspozycyjne"~ "Availability",
                            comp == "komp_kognitywne" ~ "Cognitive",
                            comp == "komp_komputerowe" ~ "Computer",
                            comp == "komp_interpersonalne" ~ "Interpersonal",
                            comp == "komp_kierownicze" ~ "Managerial",
                            comp == "komp_matematyczne" ~ "Mathematical",
                            comp == "komp_biurowe"~ "Office",
                            comp == "komp_fizyczne" ~ "Physical",
                            comp == "komp_indywidualne"~ "Self-organization",
                            comp =="komp_techniczne"  ~ "Technical")) %>%
  arrange(comp) %>%
  xtable(digits = 2) %>%
  print.xtable(include.rownames = F)
```

```{r}
vcd::assocstats(xtabs(~komp_indywidualne + kod1 + rok, data = final_data)) %>% map("cramer")
vcd::assocstats(xtabs(~komp_indywidualne + kod2 + rok, data = final_data)) %>% map("cramer")
```

# Calibration

- zawod1 
- zawod1 + woj
- zawod1 + sekcja
- zawod1 + woj + sekcja

```{r}
des <- svydesign(ids = ~1, data = final_data)
```

### Calibration based on kod2

```{r}
registerDoFuture()
plan(multiprocess)

tic <- Sys.time()
res_calib <- foreach(i = 1:500, .export = c("wynik")) %dopar% {
  set.seed(i)
  
  totals_kod2 <- subset(totals, subset = b == i)

  final_data_kod2 <- final_data %>%  group_by(rok) %>%  
      sample_frac(1, replace = T) %>% ungroup() %>%
      add_count(rok, name = "m") %>%
      left_join( totals_kod2 %>% count(rok, wt = hat_wolne, name = "total")) %>%
      mutate(weight = total / m )
  
  des_kod2 <- svydesign(ids = ~1, weight = ~weight, data = final_data_kod2)
  
  tot_kod2 <- xtabs(hat_wolne~rok + kod2, data = totals_kod2)
  tot_sek <- xtabs(hat_wolne~rok + sekcja, data = totals_kod2)
  cal_kod2 <- calibrate(design = des_kod2, formula = list(~rok + kod2), 
                      population = list(tot_kod2))
  cal_kod2_sek <- rake(design = des_kod2, 
                     sample.margins = list(~rok + kod2, ~rok + sekcja), 
                     population.margins = list(tot_kod2, tot_sek), control = list(maxit = 50))
  
  svyby(formula = ~komp_techniczne + komp_matematyczne + komp_kulturalne + komp_komputerowe + 
        komp_kognitywne + komp_kierownicze + komp_interpersonalne + 
        komp_indywidualne + komp_fizyczne + komp_dyspozycyjne + komp_biurowe, 
      by = ~ rok, 
      FUN = svymean, 
      design = cal_kod2) %>%
  select(rok, komp_techniczne:komp_biurowe) -> wyn1 
  
  svyby(formula = ~komp_techniczne + komp_matematyczne + komp_kulturalne + komp_komputerowe + 
        komp_kognitywne + komp_kierownicze + komp_interpersonalne + 
        komp_indywidualne + komp_fizyczne + komp_dyspozycyjne + komp_biurowe, 
      by = ~ rok, 
      FUN = svymean, 
      design = cal_kod2_sek) %>%
  select(rok, komp_techniczne:komp_biurowe) -> wyn2
  
  wynik <- list(calib_kod2=wyn1, 
                calib_kod2_w = weights(cal_kod2),
                calib_kod2_sek=wyn2,
                calib_kod2_sek_w = weights(cal_kod2_sek))
  
}

Sys.time() - tic
```

```{r}
saveRDS(res_calib, file = "../results/res_calib.rds")
```


### Model-calibration - GLM

```{r}
komps <- c("komp_techniczne", "komp_matematyczne", "komp_kulturalne", "komp_komputerowe", 
                             "komp_kognitywne", "komp_kierownicze", "komp_interpersonalne", 
                             "komp_indywidualne", "komp_fizyczne", "komp_dyspozycyjne", "komp_biurowe")

tic <- Sys.time()
res_glm_kod2 <- foreach(i = 1:500, .export = c("wynik")) %dopar% {

  set.seed(i)
  
  totals_boot <- totals %>% filter(b == i)
  
  final_data_kod2 <- final_data %>%  group_by(rok) %>%  
      sample_frac(1, replace = T) %>% ungroup() %>%
      add_count(rok, name = "m") %>%
      left_join( totals_boot %>% count(rok, wt = hat_wolne, name = "total")) %>%
      mutate(weight = total / m ) %>%
      as.data.frame()

    X <- Matrix::fac2sparse(final_data_kod2$kod2) %>% t() 

    wynik_est <- list()
    wynik_wt <- list()
    wynik_model <- list()
    for (k in 1:length(komps)) {
      
        m1 <- fastglm(x = as.matrix(cbind(1,X[,-1])), y = as.matrix(final_data_kod2[,komps[k]]), 
                      family = binomial(), method = 2)
  
        X_tot <- Matrix::fac2sparse(totals_boot$kod2) %>% t() 
        lin_pred  <- as.numeric(cbind(1, X_tot[,-1]) %*% m1$coefficients)
        t1 <- totals_boot
        t2 <- totals_boot
        t1$pred <- exp(lin_pred) / (1 + exp(lin_pred))
        t2$pred <- 0
        t1$flag <- "1"
        t2$flag <- "0"
  
        totals_cal <- as.data.frame(rbind(t1,t2))
        totals_cal$wt  <- ifelse(totals_cal$flag  == 1, totals_cal$pred*totals_cal$hat_wolne, totals_cal$hat_wolne)
        tab_totals <- xtabs(wt~ rok + flag, data = totals_cal)
        tab_totals[,1] <- tab_totals[,1] - tab_totals[,2]
        final_data_kod2$flag <- as.character(final_data_kod2[,komps[k]])
        glm_kod2 <- svydesign(ids = ~1, weight = ~ weight, data = final_data_kod2)
        cal_glm_kod2 <- calibrate(design = glm_kod2, 
                                    formula = list(~rok + flag), 
                                    population = list(tab_totals))
  
      svyby(formula = as.formula(paste("~", komps[k])), 
            by = ~ rok, 
            FUN = svymean, 
            design = cal_glm_kod2) %>%
        dplyr::select(-se) -> wyn

      wynik_est[[k]] <- wyn
      wynik_wt[[k]] <- weights(cal_glm_kod2)
      wynik_model[[k]] <- m1
}


wynik <- list(wynik_est = bind_cols(wynik_est) %>% select(rok, starts_with("komp")), 
              wynik_wt = bind_cols(wynik_wt),
              wynik_model = wynik_model)


}

toc <- Sys.time() - tic
toc
```

```{r}
saveRDS(res_glm_kod2, file = "../results/res_glm_kod2.rds")
```


### Lasso


Kalibracja wg kod2 + nace

```{r}
registerDoFuture()
plan(multiprocess)
tic <- Sys.time()
komps <- c("komp_techniczne", "komp_matematyczne", "komp_kulturalne", "komp_komputerowe", 
                             "komp_kognitywne", "komp_kierownicze", "komp_interpersonalne", 
                             "komp_indywidualne", "komp_fizyczne", "komp_dyspozycyjne", "komp_biurowe")

res_lasso_kod2 <- foreach(i = 1:500, .export = c("wynik")) %dopar% {
  
  set.seed(i)
  
  totals_boot <- totals %>% filter(b == i)
  
  final_data_kod2 <- final_data %>%  group_by(rok) %>%  
      sample_frac(1, replace = T) %>% ungroup() %>%
      add_count(rok, name = "m") %>%
      left_join( totals_boot %>% count(rok, wt = hat_wolne, name = "total")) %>%
      mutate(weight = total / m ) %>%
      as.data.frame()

    X <- Matrix::fac2sparse(final_data_kod2$kod2) %>% t() 

    wynik_est <- list()
    wynik_wt <- list()
    wynik_model <- list()
    for (k in 1:length(komps)) {
      
        m1 <- mycv.glmnet(x = X, y = as.matrix(final_data_kod2[,komps[k]]), intercept = T)
        X_tot <- Matrix::fac2sparse(totals_boot$kod2) %>% t() 
        lin_pred  <- as.numeric(cbind(1, X_tot) %*% m1$coef)
        t1 <- totals_boot
        t2 <- totals_boot
        t1$pred <- exp(lin_pred) / (1 + exp(lin_pred))
        t2$pred <- 0
        t1$flag <- "1"
        t2$flag <- "0"
  
        totals_cal <- as.data.frame(rbind(t1,t2))
        totals_cal$wt  <- ifelse(totals_cal$flag  == 1, totals_cal$pred*totals_cal$hat_wolne, totals_cal$hat_wolne)
        tab_totals <- xtabs(wt~ rok + flag, data = totals_cal)
        tab_totals[,1] <- tab_totals[,1] - tab_totals[,2]
        final_data_kod2$flag <- as.character(final_data_kod2[,komps[k]])
        lasso_kod2 <- svydesign(ids = ~1, weight = ~ weight, data = final_data_kod2)
        cal_lasso_kod2 <- calibrate(design = lasso_kod2, 
                                    formula = list(~rok + flag), 
                                    population = list(tab_totals))
  
      svyby(formula = as.formula(paste("~", komps[k])), 
            by = ~ rok, 
            FUN = svymean, 
            design = cal_lasso_kod2) %>%
        select(-se) -> wyn

      wynik_est[[k]] <- wyn
      wynik_wt[[k]] <- weights(cal_lasso_kod2)
      wynik_model[[k]] <- m1
}


wynik <- list(wynik_est = bind_cols(wynik_est) %>% select(rok, starts_with("komp")), 
              wynik_wt = bind_cols(wynik_wt),
              wynik_model = wynik_model)


}

toc <- Sys.time() - tic
toc
```

```{r}
saveRDS(res_lasso_kod2, file = "../results/res_lasso_kod2.rds")
```


Kalibracja wg kod2 + nace


```{r}
tic <- Sys.time()
res_lasso_kod2_nace <- foreach(i = 1:500, .export = c("wynik")) %dopar% {
  
  set.seed(i)
  
  totals_boot <- totals %>% filter(b == i)
  
  final_data_kod2 <- final_data %>%  group_by(rok) %>%  
      sample_frac(1, replace = T) %>% ungroup() %>%
      add_count(rok, name = "m") %>%
      left_join( totals_boot %>% count(rok, wt = hat_wolne, name = "total")) %>%
      mutate(weight = total / m ) %>%
      as.data.frame()

    
    Xkod <- Matrix::fac2sparse(final_data_kod2$kod2) %>% t() 
    Xsek <- Matrix::fac2sparse(final_data_kod2$sekcja) %>% t() 
    X <- cbind(Xkod, Xsek) 

    wynik_est <- list()
    wynik_wt <- list()
    wynik_model <- list()
    for (k in 1:length(komps)) {
      
        m1 <- mycv.glmnet(x = X, y = as.matrix(final_data_kod2[,komps[k]]), intercept = T)
        Xkod_tot <- Matrix::fac2sparse(totals_boot$kod2) %>% t() 
        Xsek_tot <- Matrix::fac2sparse(totals_boot$sekcja) %>% t() 
        X_tot <- cbind(Xkod_tot, Xsek_tot)
        lin_pred  <- as.numeric(cbind(1, X_tot) %*% m1$coef)
        t1 <- totals_boot
        t2 <- totals_boot
        t1$pred <- exp(lin_pred) / (1 + exp(lin_pred))
        t2$pred <- 0
        t1$flag <- "1"
        t2$flag <- "0"
  
        totals_cal <- as.data.frame(rbind(t1,t2))
        totals_cal$wt  <- ifelse(totals_cal$flag  == 1, totals_cal$pred*totals_cal$hat_wolne, totals_cal$hat_wolne)
        tab_totals <- xtabs(wt~ rok + flag, data = totals_cal)
        tab_totals[,1] <- tab_totals[,1] - tab_totals[,2]
        final_data_kod2$flag <- as.character(final_data_kod2[,komps[k]])
        lasso_kod2 <- svydesign(ids = ~1, weight = ~ weight, data = final_data_kod2)
        cal_lasso_kod2 <- calibrate(design = lasso_kod2, 
                                    formula = list(~rok + flag), 
                                    population = list(tab_totals))
  
      svyby(formula = as.formula(paste("~", komps[k])), 
            by = ~ rok, 
            FUN = svymean, 
            design = cal_lasso_kod2) %>%
        select(-se) -> wyn

      wynik_est[[k]] <- wyn
      wynik_wt[[k]] <- weights(cal_lasso_kod2)
      wynik_model[[k]] <- m1
}


wynik <- list(wynik_est = bind_cols(wynik_est) %>% select(rok, starts_with("komp")), 
              wynik_wt = bind_cols(wynik_wt),
              wynik_model = wynik_model)
}
toc <- Sys.time() - tic
toc
```

```{r}
saveRDS(res_lasso_kod2_nace, file = "../results/res_lasso_kod2_nace.rds")
```

### Adaptive lasso


```{r}
registerDoFuture()
plan(multiprocess)
```

```{r}

tic <- Sys.time()
komps <- c("komp_techniczne", "komp_matematyczne", "komp_kulturalne", "komp_komputerowe", 
                             "komp_kognitywne", "komp_kierownicze", "komp_interpersonalne", 
                             "komp_indywidualne", "komp_fizyczne", "komp_dyspozycyjne", "komp_biurowe")

res_alasso_kod2 <- foreach(i = 1:10, .export = c("wynik"), .verbose = TRUE) %dopar% {
  
  set.seed(i)
  
  totals_boot <- totals %>% filter(b == i)
  
  final_data_kod2 <- final_data %>%  group_by(rok) %>%  
      sample_frac(1, replace = T) %>% ungroup() %>%
      add_count(rok, name = "m") %>%
      left_join( totals_boot %>% count(rok, wt = hat_wolne, name = "total")) %>%
      mutate(weight = total / m ) %>%
      as.data.frame()

    X <- Matrix::fac2sparse(final_data_kod2$kod2) %>% t() 

    wynik_est <- list()
    wynik_wt <- list()
    wynik_model <- list()
    for (k in 1:length(komps)) {
      
        m1 <- mycv.glmnet(x = X, y = as.matrix(final_data_kod2[,komps[k]]), intercept = T, alpha = 0)
        m1 <- mycv.glmnet(x = X, y = as.matrix(final_data_kod2[,komps[k]]), intercept = T, penalty.factor = 1/abs(m1$coef[-1]))
        X_tot <- Matrix::fac2sparse(totals_boot$kod2) %>% t() 
        lin_pred  <- as.numeric(cbind(1, X_tot) %*% m1$coef)
        t1 <- totals_boot
        t2 <- totals_boot
        t1$pred <- exp(lin_pred) / (1 + exp(lin_pred))
        t2$pred <- 0
        t1$flag <- "1"
        t2$flag <- "0"
  
        totals_cal <- as.data.frame(rbind(t1,t2))
        totals_cal$wt  <- ifelse(totals_cal$flag  == 1, totals_cal$pred*totals_cal$hat_wolne, totals_cal$hat_wolne)
        tab_totals <- xtabs(wt~ rok + flag, data = totals_cal)
        tab_totals[,1] <- tab_totals[,1] - tab_totals[,2]
        final_data_kod2$flag <- as.character(final_data_kod2[,komps[k]])
        lasso_kod2 <- svydesign(ids = ~1, weight = ~ weight, data = final_data_kod2)
        cal_lasso_kod2 <- calibrate(design = lasso_kod2, 
                                    formula = list(~rok + flag), 
                                    population = list(tab_totals))
  
      svyby(formula = as.formula(paste("~", komps[k])), 
            by = ~ rok, 
            FUN = svymean, 
            design = cal_lasso_kod2) %>%
        select(-se) -> wyn

      wynik_est[[k]] <- wyn
      wynik_wt[[k]] <- weights(cal_lasso_kod2)
      wynik_model[[k]] <- m1
}


wynik <- list(wynik_est = bind_cols(wynik_est), 
              wynik_wt = bind_cols(wynik_wt),
              wynik_model = wynik_model)


}

toc <- Sys.time() - tic
toc
```

```{r}
saveRDS(res_alasso_kod2, file = "../results/res_alasso_kod2.rds")
```



Kalibracja wg kod2 + nace


```{r}
tic <- Sys.time()
res_alasso_kod2_nace <- foreach(i = 1:500, .export = c("wynik")) %dopar% {
  
  set.seed(i)
  
  totals_boot <- totals %>% filter(b == i)
  
  final_data_kod2 <- final_data %>%  group_by(rok) %>%  
      sample_frac(1, replace = T) %>% ungroup() %>%
      add_count(rok, name = "m") %>%
      left_join( totals_boot %>% count(rok, wt = hat_wolne, name = "total")) %>%
      mutate(weight = total / m ) %>%
      as.data.frame()

    
    Xkod <- Matrix::fac2sparse(final_data_kod2$kod2) %>% t() 
    Xsek <- Matrix::fac2sparse(final_data_kod2$sekcja) %>% t() 
    X <- cbind(Xkod, Xsek) 

    wynik_est <- list()
    wynik_wt <- list()
    wynik_model <- list()
    for (k in 1:length(komps)) {
      
        m1 <- mycv.glmnet(x = X, y = as.matrix(final_data_kod2[,komps[k]]), intercept = T, alpha = 0)
        m1 <- mycv.glmnet(x = X, y = as.matrix(final_data_kod2[,komps[k]]), intercept = T, penalty.factor = 1/abs(m1$coef[-1]))
        Xkod_tot <- Matrix::fac2sparse(totals_boot$kod2) %>% t() 
        Xsek_tot <- Matrix::fac2sparse(totals_boot$sekcja) %>% t() 
        X_tot <- cbind(Xkod_tot, Xsek_tot)
        lin_pred  <- as.numeric(cbind(1, X_tot) %*% m1$coef)
        t1 <- totals_boot
        t2 <- totals_boot
        t1$pred <- exp(lin_pred) / (1 + exp(lin_pred))
        t2$pred <- 0
        t1$flag <- "1"
        t2$flag <- "0"
  
        totals_cal <- as.data.frame(rbind(t1,t2))
        totals_cal$wt  <- ifelse(totals_cal$flag  == 1, totals_cal$pred*totals_cal$hat_wolne, totals_cal$hat_wolne)
        tab_totals <- xtabs(wt~ rok + flag, data = totals_cal)
        tab_totals[,1] <- tab_totals[,1] - tab_totals[,2]
        final_data_kod2$flag <- as.character(final_data_kod2[,komps[k]])
        lasso_kod2 <- svydesign(ids = ~1, weight = ~ weight, data = final_data_kod2)
        cal_lasso_kod2 <- calibrate(design = lasso_kod2, 
                                    formula = list(~rok + flag), 
                                    population = list(tab_totals))
  
      svyby(formula = as.formula(paste("~", komps[k])), 
            by = ~ rok, 
            FUN = svymean, 
            design = cal_lasso_kod2) %>%
        select(-se) -> wyn

      wynik_est[[k]] <- wyn
      wynik_wt[[k]] <- weights(cal_lasso_kod2)
      wynik_model[[k]] <- m1
}


wynik <- list(wynik_est = bind_cols(wynik_est), 
              wynik_wt = bind_cols(wynik_wt),
              wynik_model = wynik_model)
}
toc <- Sys.time() - tic
toc
```

```{r}
saveRDS(res_alasso_kod2_nace, file = "../results/res_alasso_kod2_nace.rds")
```
