vignettes/ffv_opt_solin_relow.Rmd
ffv_opt_solin_relow.RmdGiven allocation space parameters (A, alpha), test the lower-bounded continuous optimal allocation line by line without function. This demonstrates relative optimality conditions, inverse optimal allocation, and optimal allocation.
Solves the linear optimal allocation problem. Each individual receives input \(N_i\), and we are interested in individual expected outcomes \(H_i\). Production functions or reduced causal form relationships are linear with individual specific slopes and intercepts. Solve for different planners.
library(dplyr)
library(tidyr)
library(REconTools)
library(PrjOptiAlloc)
bl_save_rda = FALSELoad data that is generated by regression ffy_opt_dtgch_cbem4 (vignette).
Alternatively, provide hand input data.
bl_save_rda = FALSE
it_data_which = 2
if (it_data_which == 1) {
# the data here has very little heterogeneity in alpha
# Load Library
ls_opti_alpha_A <- ffy_opt_dtgch_cbem4()
df_raw <- ls_opti_alpha_A$df_raw
df_hw_cebu_m24 <- df_raw
df_esti <- ls_opti_alpha_A$df_esti
# Review dataframes
# raw file
head(df_raw, 10)
head(df_esti, 10)
ar_prot_data = df_hw_cebu_m24$prot
fl_N_agg = sum(ar_prot_data)
# Attach
attach(df_raw)
}
if (it_data_which == 2) {
# generate 50 individuals with arbitrary A, alpha and beta
set.seed(123)
it_N <- 50
ar_alpha <- runif(it_N, min=0.1, max=1)
ar_A <- runif(it_N, min=1, max=2)
ar_beta <- runif(it_N, min=0, max=1)
ar_beta <- ar_beta/sum(ar_beta)
mt_cts <- cbind(seq(1, it_N), ar_A, ar_alpha, ar_beta)
colnames(mt_cts) <- c('indi.id', 'A_lin', 'alpha_lin', 'beta')
df_esti <- as_tibble(mt_cts)
fl_N_agg <- 10*mean(ar_alpha)
}I would like to be able to directly call the dataset generated here in various functions. Save the datafile we just created in the project folder.
df <- df_esti
svr_A_i <- 'A_lin'
svr_alpha_i <- 'alpha_lin'
svr_beta_i <- 'beta'
fl_N_agg <- fl_N_agg
fl_rho <- -10
# a. select only relevant data
df_opti <- df %>% rename(A = !!sym(svr_A_i), alpha = !!sym(svr_alpha_i), beta = !!sym(svr_beta_i))
# b. Generate V4, Rank Index Value, rho specific
# df_opti <- df_opti %>% mutate(!!paste0('rv_', it_rho_ctr) := A/((alpha*beta))^(1/(1-fl_rho)))
df_opti <- df_opti %>% mutate(rank_val = A/((alpha*beta))^(1/(1-fl_rho)))
# c. Generate Rank Index based on y-intercept
df_opti <- df_opti %>% arrange(rank_val) %>% mutate(rank_y_inter_idx = row_number())
# d. Populate lowest index alpha, beta, and A to all rows
df_opti <- df_opti %>% mutate(lowest_rank_A = A[rank_y_inter_idx==1]) %>%
mutate(lowest_rank_alpha = alpha[rank_y_inter_idx==1]) %>%
mutate(lowest_rank_beta = beta[rank_y_inter_idx==1])
# Print
head(df_opti, 10)## # A tibble: 10 x 9
## indi.id A alpha beta rank_val rank_y_inter_idx lowest_rank_A
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl>
## 1 4 1.12 0.895 0.0370 1.53 1 1.12
## 2 24 1.00 0.995 0.00853 1.54 2 1.12
## 3 7 1.13 0.575 0.0355 1.61 3 1.12
## 4 1 1.05 0.359 0.0233 1.62 4 1.12
## 5 26 1.22 0.738 0.0382 1.69 5 1.12
## 6 12 1.09 0.508 0.0117 1.74 6 1.12
## 7 31 1.24 0.967 0.0240 1.75 7 1.12
## 8 30 1.11 0.232 0.0268 1.76 8 1.12
## 9 14 1.27 0.615 0.0368 1.80 9 1.12
## 10 40 1.18 0.308 0.0171 1.89 10 1.12
## # ... with 2 more variables: lowest_rank_alpha <dbl>, lowest_rank_beta <dbl>
summary(df_opti)## indi.id A alpha beta
## Min. : 1.00 Min. :1.001 Min. :0.1222 Min. :0.0004062
## 1st Qu.:13.25 1st Qu.:1.226 1st Qu.:0.3442 1st Qu.:0.0100739
## Median :25.50 Median :1.445 Median :0.5527 Median :0.0195919
## Mean :25.50 Mean :1.477 Mean :0.5681 Mean :0.0200000
## 3rd Qu.:37.75 3rd Qu.:1.743 3rd Qu.:0.8028 3rd Qu.:0.0285548
## Max. :50.00 Max. :1.985 Max. :0.9948 Max. :0.0381932
## rank_val rank_y_inter_idx lowest_rank_A lowest_rank_alpha
## Min. :1.529 Min. : 1.00 Min. :1.122 Min. :0.8947
## 1st Qu.:1.961 1st Qu.:13.25 1st Qu.:1.122 1st Qu.:0.8947
## Median :2.295 Median :25.50 Median :1.122 Median :0.8947
## Mean :2.302 Mean :25.50 Mean :1.122 Mean :0.8947
## 3rd Qu.:2.728 3rd Qu.:37.75 3rd Qu.:1.122 3rd Qu.:0.8947
## Max. :3.034 Max. :50.00 Max. :1.122 Max. :0.8947
## lowest_rank_beta
## Min. :0.03704
## 1st Qu.:0.03704
## Median :0.03704
## Mean :0.03704
## 3rd Qu.:0.03704
## Max. :0.03704## lowest_rank_A lowest_rank_alpha lowest_rank_beta
## Min. :1.122 Min. :0.8947 Min. :0.03704
## 1st Qu.:1.122 1st Qu.:0.8947 1st Qu.:0.03704
## Median :1.122 Median :0.8947 Median :0.03704
## Mean :1.122 Mean :0.8947 Mean :0.03704
## 3rd Qu.:1.122 3rd Qu.:0.8947 3rd Qu.:0.03704
## Max. :1.122 Max. :0.8947 Max. :0.03704
# e. relative slope and relative intercept with respect to lowest index
df_opti <- df_opti %>%
mutate(rela_slope_to_lowest =
((
(lowest_rank_alpha*lowest_rank_beta)/(alpha*beta)
)^(1/(fl_rho-1))
*(lowest_rank_alpha/alpha))
) %>%
mutate(rela_intercept_to_lowest =
(((
(lowest_rank_alpha*lowest_rank_beta)/(alpha*beta)
)^(1/(fl_rho-1))
*(lowest_rank_A/alpha)) - (A/alpha))
)
# h. Relative x-intercept points and ranking based on x-intercepts
df_opti <- df_opti %>%
mutate(rela_x_intercept =
(-1)*(rela_intercept_to_lowest/rela_slope_to_lowest)
) %>% arrange(rela_x_intercept) %>%
mutate(rank_x_inter_idx = row_number())
# In this example, is x-intercept rank and y-intercept rank the same?
sum(df_opti$rank_y_inter_idx-df_opti$rank_x_inter_idx)## [1] 0## rela_slope_to_lowest rela_intercept_to_lowest
## Min. :0.7946 Min. :-6.2816
## 1st Qu.:1.0155 1st Qu.:-1.4794
## Median :1.3717 Median :-0.8238
## Mean :1.8444 Mean :-1.1871
## 3rd Qu.:2.1876 3rd Qu.:-0.5305
## Max. :5.7846 Max. : 0.0000
df_opti %>%
select(A, alpha, beta,
rank_val, rela_intercept_to_lowest, rela_x_intercept,
rank_y_inter_idx, rank_x_inter_idx)## # A tibble: 50 x 8
## A alpha beta rank_val rela_intercept_to_lowest rela_x_intercept
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.12 0.895 0.0370 1.53 0 0
## 2 1.00 0.995 0.00853 1.54 -0.00948 0.0119
## 3 1.13 0.575 0.0355 1.61 -0.0938 0.0630
## 4 1.05 0.359 0.0233 1.62 -0.156 0.0710
## 5 1.22 0.738 0.0382 1.69 -0.155 0.130
## 6 1.09 0.508 0.0117 1.74 -0.266 0.177
## 7 1.24 0.967 0.0240 1.75 -0.163 0.182
## 8 1.11 0.232 0.0268 1.76 -0.635 0.192
## 9 1.27 0.615 0.0368 1.80 -0.310 0.221
## 10 1.18 0.308 0.0171 1.89 -0.733 0.299
## # ... with 40 more rows, and 2 more variables: rank_y_inter_idx <int>,
## # rank_x_inter_idx <int>
# f. cumulative sums
df_opti <- df_opti %>%
mutate(rela_slope_to_lowest_cumsum =
cumsum(rela_slope_to_lowest)
) %>%
mutate(rela_intercept_to_lowest_cumsum =
cumsum(rela_intercept_to_lowest)
)
# Print
summary(df_opti %>% select(rela_slope_to_lowest_cumsum, rela_intercept_to_lowest_cumsum))## rela_slope_to_lowest_cumsum rela_intercept_to_lowest_cumsum
## Min. : 1.00 Min. :-59.356
## 1st Qu.:30.81 1st Qu.:-28.787
## Median :48.36 Median :-15.944
## Mean :46.42 Mean :-19.426
## 3rd Qu.:64.62 3rd Qu.: -7.318
## Max. :92.22 Max. : 0.000
# g. inverting cumulative slopes and intercepts
df_opti <- df_opti %>%
mutate(rela_slope_to_lowest_cumsum_invert =
(1/rela_slope_to_lowest_cumsum)
) %>%
mutate(rela_intercept_to_lowest_cumsum_invert =
((-1)*(rela_intercept_to_lowest_cumsum)/(rela_slope_to_lowest_cumsum))
)
# Print
summary(df_opti %>% select(rela_slope_to_lowest_cumsum_invert, rela_intercept_to_lowest_cumsum_invert))## rela_slope_to_lowest_cumsum_invert rela_intercept_to_lowest_cumsum_invert
## Min. :0.01084 Min. :0.0000
## 1st Qu.:0.01548 1st Qu.:0.2375
## Median :0.02068 Median :0.3296
## Mean :0.06883 Mean :0.3284
## 3rd Qu.:0.03247 3rd Qu.:0.4453
## Max. :1.00000 Max. :0.6436
# # h. Relative x-intercept points
# df_opti <- df_opti %>%
# mutate(rela_x_intercept =
# (-1)*(rela_intercept_to_lowest/rela_slope_to_lowest)
# )
# i. Inverted relative x-intercepts
df_opti <- df_opti %>%
mutate(opti_lowest_spline_knots =
(rela_intercept_to_lowest_cumsum + rela_slope_to_lowest_cumsum*rela_x_intercept)
)
# Print
summary(df_opti %>% select(rela_x_intercept, opti_lowest_spline_knots))## rela_x_intercept opti_lowest_spline_knots
## Min. :0.0000 Min. : 0.000
## 1st Qu.:0.3538 1st Qu.: 3.586
## Median :0.6277 Median :14.415
## Mean :0.6334 Mean :18.657
## 3rd Qu.:0.9831 3rd Qu.:34.778
## Max. :1.2340 Max. :54.438
# j. Sort by order of receiving transfers/subsidies
df_opti <- df_opti %>% arrange(rela_x_intercept)
# k. Find position of subsidy
df_opti <- df_opti %>% arrange(opti_lowest_spline_knots) %>%
mutate(tot_devi = opti_lowest_spline_knots - fl_N_agg) %>%
arrange((-1)*case_when(tot_devi < 0 ~ tot_devi)) %>%
mutate(allocate_lowest =
case_when(row_number() == 1 ~
rela_intercept_to_lowest_cumsum_invert +
rela_slope_to_lowest_cumsum_invert*fl_N_agg)) %>%
mutate(allocate_lowest = allocate_lowest[row_number() == 1])
# Print
summary(df_opti %>% select(tot_devi, allocate_lowest))## tot_devi allocate_lowest
## Min. :-5.681 Min. :0.4198
## 1st Qu.:-2.095 1st Qu.:0.4198
## Median : 8.735 Median :0.4198
## Mean :12.976 Mean :0.4198
## 3rd Qu.:29.097 3rd Qu.:0.4198
## Max. :48.758 Max. :0.4198
# l. Find position of subsidy
df_opti <- df_opti %>%
mutate(opti_allocate =
rela_intercept_to_lowest +
rela_slope_to_lowest*allocate_lowest) %>%
mutate(opti_allocate =
case_when(opti_allocate >= 0 ~ opti_allocate, TRUE ~ 0)) %>%
mutate(allocate_total = sum(opti_allocate, na.rm=TRUE))
# if want 0 to be 0
# case_when(opti_allocate >= 0 ~ opti_allocate, TRUE ~ 0)) %>%
# if want 0 to be NA, makes it easier to see average change for non zeros
# case_when(opti_allocate >= 0 ~ opti_allocate)) %>%
# Print
summary(df_opti %>% select(opti_allocate, allocate_total))## opti_allocate allocate_total
## Min. :0.0000 Min. :5.681
## 1st Qu.:0.0000 1st Qu.:5.681
## Median :0.0000 Median :5.681
## Mean :0.1136 Mean :5.681
## 3rd Qu.:0.1792 3rd Qu.:5.681
## Max. :0.7674 Max. :5.681
# l. Predictes Expected choice: A + alpha*opti_allocate
df_opti <- df_opti %>% mutate(opti_exp_outcome = A + alpha*opti_allocate)
# m. Drop some variables that I do not want to keep even in full df to export
# inpalc = input allocation optimal
# expout = expected outcome given input allocation
df_opti <- df_opti %>% select(-one_of(c('lowest_rank_alpha', 'lowest_rank_beta')))
ar_opti_inpalc <- df_opti %>% pull(opti_allocate)
ar_opti_expout <- df_opti %>% pull(opti_exp_outcome)
# Print
head(df_opti, 10)## # A tibble: 10 x 21
## indi.id A alpha beta rank_val rank_y_inter_idx lowest_rank_A
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl>
## 1 33 1.42 0.722 0.0261 2.03 15 1.12
## 2 48 1.09 0.519 0.00298 1.97 14 1.12
## 3 6 1.21 0.141 0.0346 1.96 13 1.12
## 4 41 1.13 0.229 0.0121 1.93 12 1.12
## 5 35 1.10 0.122 0.0202 1.90 11 1.12
## 6 40 1.18 0.308 0.0171 1.89 10 1.12
## 7 14 1.27 0.615 0.0368 1.80 9 1.12
## 8 30 1.11 0.232 0.0268 1.76 8 1.12
## 9 31 1.24 0.967 0.0240 1.75 7 1.12
## 10 12 1.09 0.508 0.0117 1.74 6 1.12
## # ... with 14 more variables: rela_slope_to_lowest <dbl>,
## # rela_intercept_to_lowest <dbl>, rela_x_intercept <dbl>,
## # rank_x_inter_idx <int>, rela_slope_to_lowest_cumsum <dbl>,
## # rela_intercept_to_lowest_cumsum <dbl>,
## # rela_slope_to_lowest_cumsum_invert <dbl>,
## # rela_intercept_to_lowest_cumsum_invert <dbl>,
## # opti_lowest_spline_knots <dbl>, tot_devi <dbl>, allocate_lowest <dbl>, ...
summary(df_opti)## indi.id A alpha beta
## Min. : 1.00 Min. :1.001 Min. :0.1222 Min. :0.0004062
## 1st Qu.:13.25 1st Qu.:1.226 1st Qu.:0.3442 1st Qu.:0.0100739
## Median :25.50 Median :1.445 Median :0.5527 Median :0.0195919
## Mean :25.50 Mean :1.477 Mean :0.5681 Mean :0.0200000
## 3rd Qu.:37.75 3rd Qu.:1.743 3rd Qu.:0.8028 3rd Qu.:0.0285548
## Max. :50.00 Max. :1.985 Max. :0.9948 Max. :0.0381932
## rank_val rank_y_inter_idx lowest_rank_A rela_slope_to_lowest
## Min. :1.529 Min. : 1.00 Min. :1.122 Min. :0.7946
## 1st Qu.:1.961 1st Qu.:13.25 1st Qu.:1.122 1st Qu.:1.0155
## Median :2.295 Median :25.50 Median :1.122 Median :1.3717
## Mean :2.302 Mean :25.50 Mean :1.122 Mean :1.8444
## 3rd Qu.:2.728 3rd Qu.:37.75 3rd Qu.:1.122 3rd Qu.:2.1876
## Max. :3.034 Max. :50.00 Max. :1.122 Max. :5.7846
## rela_intercept_to_lowest rela_x_intercept rank_x_inter_idx
## Min. :-6.2816 Min. :0.0000 Min. : 1.00
## 1st Qu.:-1.4794 1st Qu.:0.3538 1st Qu.:13.25
## Median :-0.8238 Median :0.6277 Median :25.50
## Mean :-1.1871 Mean :0.6334 Mean :25.50
## 3rd Qu.:-0.5305 3rd Qu.:0.9831 3rd Qu.:37.75
## Max. : 0.0000 Max. :1.2340 Max. :50.00
## rela_slope_to_lowest_cumsum rela_intercept_to_lowest_cumsum
## Min. : 1.00 Min. :-59.356
## 1st Qu.:30.81 1st Qu.:-28.787
## Median :48.36 Median :-15.944
## Mean :46.42 Mean :-19.426
## 3rd Qu.:64.62 3rd Qu.: -7.318
## Max. :92.22 Max. : 0.000
## rela_slope_to_lowest_cumsum_invert rela_intercept_to_lowest_cumsum_invert
## Min. :0.01084 Min. :0.0000
## 1st Qu.:0.01548 1st Qu.:0.2375
## Median :0.02068 Median :0.3296
## Mean :0.06883 Mean :0.3284
## 3rd Qu.:0.03247 3rd Qu.:0.4453
## Max. :1.00000 Max. :0.6436
## opti_lowest_spline_knots tot_devi allocate_lowest opti_allocate
## Min. : 0.000 Min. :-5.681 Min. :0.4198 Min. :0.0000
## 1st Qu.: 3.586 1st Qu.:-2.095 1st Qu.:0.4198 1st Qu.:0.0000
## Median :14.415 Median : 8.735 Median :0.4198 Median :0.0000
## Mean :18.657 Mean :12.976 Mean :0.4198 Mean :0.1136
## 3rd Qu.:34.778 3rd Qu.:29.097 3rd Qu.:0.4198 3rd Qu.:0.1792
## Max. :54.438 Max. :48.758 Max. :0.4198 Max. :0.7674
## allocate_total opti_exp_outcome
## Min. :5.681 Min. :1.133
## 1st Qu.:5.681 1st Qu.:1.357
## Median :5.681 Median :1.471
## Mean :5.681 Mean :1.529
## 3rd Qu.:5.681 3rd Qu.:1.743
## Max. :5.681 Max. :1.985## ar_opti_expout ar_opti_inpalc
## [1,] 1.422612 0.006880531
## [2,] 1.133311 0.076469304
## [3,] 1.257978 0.364864308
## [4,] 1.194809 0.280558877
## [5,] 1.182663 0.653270709
## [6,] 1.266780 0.297369431
## [7,] 1.446481 0.279664218
## [8,] 1.286286 0.753652841
## [9,] 1.449922 0.213404330
## [10,] 1.280823 0.366106406
## [11,] 1.475576 0.346299292
## [12,] 1.321196 0.767417565
## [13,] 1.432987 0.530954043
## [14,] 1.323050 0.324096274
## [15,] 1.497510 0.419810625
## [16,] 1.320373 0.000000000
## [17,] 1.187691 0.000000000
## [18,] 1.439832 0.000000000
## [19,] 1.511505 0.000000000
## [20,] 1.434893 0.000000000
## [21,] 1.442200 0.000000000
## [22,] 1.475317 0.000000000
## [23,] 1.560948 0.000000000
## [24,] 1.665115 0.000000000
## [25,] 1.668056 0.000000000
## [26,] 1.379817 0.000000000
## [27,] 1.374463 0.000000000
## [28,] 1.448516 0.000000000
## [29,] 1.351798 0.000000000
## [30,] 1.466779 0.000000000
## [31,] 1.383970 0.000000000
## [32,] 1.754475 0.000000000
## [33,] 1.753308 0.000000000
## [34,] 1.629221 0.000000000
## [35,] 1.788196 0.000000000
## [36,] 1.653102 0.000000000
## [37,] 1.710182 0.000000000
## [38,] 1.794342 0.000000000
## [39,] 1.798925 0.000000000
## [40,] 1.886469 0.000000000
## [41,] 1.984957 0.000000000
## [42,] 1.656758 0.000000000
## [43,] 1.612771 0.000000000
## [44,] 1.810064 0.000000000
## [45,] 1.895045 0.000000000
## [46,] 1.893051 0.000000000
## [47,] 1.814640 0.000000000
## [48,] 1.343516 0.000000000
## [49,] 1.812390 0.000000000
## [50,] 1.782294 0.000000000
if (bl_save_rda) {
# df_opt_dtgch_ropti: dataframe, opt project, data guat cebu height, cebu edu mother, results relative linear optimal
df_opt_dtgch_cbem4_rrlop <- df_opti
usethis::use_data(df_opt_dtgch_cbem4_rrlop, df_opt_dtgch_cbem4_rrlop, overwrite = TRUE)
}