# install.packages("fastglm")
library(survey)
Loading required package: grid
Loading required package: Matrix
Loading required package: survival
Attaching package: ‘survey’
The following object is masked from ‘package:graphics’:
dotchart
library(glmnet)
Loading required package: foreach
Loaded glmnet 2.0-18
library(tidyverse)
[30m── [1mAttaching packages[22m ────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──[39m
[30m[32m✔[30m [34mggplot2[30m 3.2.0 [32m✔[30m [34mpurrr [30m 0.3.2
[32m✔[30m [34mtibble [30m 2.1.3 [32m✔[30m [34mdplyr [30m 0.8.1
[32m✔[30m [34mtidyr [30m 0.8.3 [32m✔[30m [34mstringr[30m 1.4.0
[32m✔[30m [34mreadr [30m 1.3.1 [32m✔[30m [34mforcats[30m 0.4.0[39m
[30m── [1mConflicts[22m ───────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31m✖[30m [34mpurrr[30m::[32maccumulate()[30m masks [34mforeach[30m::accumulate()
[31m✖[30m [34mtidyr[30m::[32mexpand()[30m masks [34mMatrix[30m::expand()
[31m✖[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31m✖[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()
[31m✖[30m [34mpurrr[30m::[32mwhen()[30m masks [34mforeach[30m::when()[39m
library(xtable)
library(future)
Attaching package: ‘future’
The following object is masked from ‘package:survival’:
cluster
library(vcd)
library(doFuture)
Loading required package: globals
Loading required package: iterators
Loading required package: parallel
library(fastglm)
source("../codes/adalasso-chen-phd.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)
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)
Joining, by = "comp"
Joining, by = "comp"
% latex table generated in R 3.5.1 by xtable 1.8-4 package
% Wed Jul 24 14:01:48 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrr}
\hline
comp & occupancy & NACE & Voivodeship \\
\hline
Artistic & 0.22 & 0.11 & 0.05 \\
Availability & 0.15 & 0.14 & 0.05 \\
Cognitive & 0.21 & 0.06 & 0.06 \\
Computer & 0.45 & 0.23 & 0.10 \\
Interpersonal & 0.42 & 0.23 & 0.06 \\
Managerial & 0.34 & 0.15 & 0.04 \\
Mathematical & 0.05 & 0.02 & 0.03 \\
Office & 0.11 & 0.06 & 0.03 \\
Physical & 0.17 & 0.09 & 0.04 \\
Self-organization & 0.34 & 0.19 & 0.04 \\
Technical & 0.31 & 0.11 & 0.07 \\
\hline
\end{tabular}
\end{table}
vcd::assocstats(xtabs(~komp_indywidualne + kod1 + rok, data = final_data)) %>% map("cramer")
$`rok:2011`
[1] 0.2501016
$`rok:2013`
[1] 0.231543
$`rok:2014`
[1] 0.2218449
vcd::assocstats(xtabs(~komp_indywidualne + kod2 + rok, data = final_data)) %>% map("cramer")
$`rok:2011`
[1] 0.3659692
$`rok:2013`
[1] 0.3381546
$`rok:2014`
[1] 0.3601699
des <- svydesign(ids = ~1, data = final_data)
No weights or probabilities supplied, assuming equal probability
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
saveRDS(res_calib, file = "../results/res_calib.rds")
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)
}
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
Joining, by = "rok"
toc <- Sys.time() - tic
toc
Time difference of 2.531316 hours
Kalibracja wg kod2 + nace
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
saveRDS(res_lasso_kod2, file = "../results/res_lasso_kod2.rds")
Kalibracja wg kod2 + nace
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
saveRDS(res_lasso_kod2_nace, file = "../results/res_lasso_kod2_nace.rds")
registerDoFuture()
plan(multiprocess)
[ONE-TIME WARNING] Forked processing ('multicore') is disabled in future (>= 1.13.0) when running R from RStudio, because it is considered unstable. Because of this, plan("multicore") will fall back to plan("sequential"), and plan("multiprocess") will fall back to plan("multisession") - not plan("multicore") as in the past. For more details, how to control forked processing or not, and how to silence this warning in future R sessions, see ?future::supportsMulticore
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:500, .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)
}
numValues: 500, numResults: 0, stopped: TRUE
Error: Failed to retrieve the value of MultisessionFuture (<none>) from cluster SOCKnode #8 (PID 17737 on localhost ‘localhost’). The reason reported was ‘vector memory exhausted (limit reached?)’. Post-mortem diagnostic: A process with this PID exists, which suggests that the localhost worker is still alive.
saveRDS(res_alasso_kod2, file = "../results/res_alasso_kod2.rds")
Kalibracja wg kod2 + nace
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
saveRDS(res_alasso_kod2_nace, file = "../results/res_alasso_kod2_nace.rds")