## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, fig.width = 6, fig.height = 4, fig.align = "center" ) ## ----------------------------------------------------------------------------- library(rcbayes) library(tibble) library(ggplot2) pars <- c(a1= 0.09, alpha1= 0.1, a2= 0.2, alpha2= 0.1, mu2= 21, lambda2= 0.4, a3= 0.02, alpha3= 0.25, mu3= 67, lambda3= 0.6, a4= 0.01, lambda4= 0.01, c= 0.01) ages <- 0:100 mx <- mig_calculate_rc(ages = ages, pars = pars) # plot to see what the schedule looks like df <- tibble(age = ages, mx = mx) df %>% ggplot(aes(age, mx)) + geom_line() + ggtitle("Rogers-Castro age-specific migration schedule (13-parameter)") ## ----------------------------------------------------------------------------- pars <- c(a1= 0.09, alpha1= 0.1, a2= 0.2, alpha2= 0.05, mu2= 25, lambda2= 0.4, c= 0.01) ages <- 0:100 mx <- mig_calculate_rc(ages = ages, pars = pars) # plot to see what the schedule looks like df <- tibble(age = ages, mx = mx) df %>% ggplot(aes(age, mx)) + geom_line() + ggtitle("Rogers-Castro age-specific migration schedule (9-parameter)") ## ---- include=FALSE----------------------------------------------------------- set.seed(123) ## ----------------------------------------------------------------------------- fl_ages <- 0:80 fl_migrants <- c(49, 48, 48, 52, 50, 45, 42, 46, 45, 44, 47, 55, 57, 59, 67, 69, 71, 78, 93, 88, 116, 106, 102, 104, 102, 123, 112, 102, 112, 105, 100, 83, 81, 77, 78, 77, 66, 64, 65, 64, 68, 52, 59, 51, 54, 55, 52, 58, 64, 53, 68, 53, 57, 67, 71, 78, 75, 77, 77, 83, 88, 80, 84, 79, 77, 83, 71, 59, 65, 67, 64, 63, 56, 50, 43, 46, 46, 38, 32, 28, 29) fl_pop <- c(2028, 2193, 2271, 2370, 2403, 2160, 2109, 2206, 2456, 2334, 2392, 2534, 2542, 2601, 2526, 2416, 2420, 2344, 2606, 2355, 2867, 2589, 2426, 2390, 2377, 2909, 2753, 2633, 2847, 2819, 2979, 2608, 2708, 2602, 2745, 2883, 2624, 2607, 2677, 2637, 2964, 2414, 2481, 2464, 2510, 2695, 2552, 2711, 2794, 2683, 2888, 2439, 2631, 2814, 2854, 2999, 2959, 2852, 2957, 2985, 2970, 2882, 2839, 2737, 2782, 2799, 2710, 2527, 2512, 2530, 2505, 2521, 2551, 2125, 1838, 2057, 2037, 1804, 1542, 1470, 1452) df <- tibble(age = fl_ages, mx = fl_migrants / fl_pop) df %>% ggplot(aes(age, mx)) + geom_point() + ggtitle("Observed migration rates") ## ---- eval=FALSE-------------------------------------------------------------- # rc_res <- mig_estimate_rc( # ages=fl_ages, migrants=fl_migrants, pop=fl_pop, # pre_working_age = TRUE, # working_age = TRUE, # retirement = TRUE, # post_retirement = FALSE, # # (optional) arguments for Stan # chains = 4, # iter = 2000, # control = list(adapt_delta = 0.8, max_treedepth = 10) # ) ## ---- echo=FALSE-------------------------------------------------------------- rc_res = list(check_converge = c(), pars_df = c(), fit_df = c()) rc_res[['check_converge']] <- matrix( c(0.870226738894234,0.0106440162045495,2978.76440313297,0.99947658219495, 0.190530488368833,0.00054434025065072,1786.98524371334,1.000853893877, 0.205526356090562,0.00149377245438701,1706.25886324742,1.00109410733695, 0.00472626420715214,6.11293188266945e-05,2180.68698683192,1.00187646979496, 0.0583736582490384,0.000117918917152783,1821.75108753928,1.00032592498145, 0.019992049679246,0.00010495874412877,1671.26255704968,1.00289212248805, 24.9853085965459,0.02063815493269,2198.85598448598,1.00027683051002, 64.8492122289137,0.0192261420687394,2799.85127600686,1.00071123654953, 0.140144946218822,0.000320010302131112,1860.81851649046,1.00209202390757, 0.12811347166421,0.000663115574100773,1639.28943120517,1.00037088664524, 0.0200589900501721,3.5975555236367e-05,1002.95302272389,1.0015353838189), ncol=4, dimnames = list(c("alpha1[1]", "alpha2[1]", "alpha3[1]", "a1[1]", "a2[1]", "a3[1]", "mu2[1]", "mu3[1]", "lambda2[1]", "lambda3[1]", "c"), c("mean", "se_mean", "n_eff", "Rhat")), byrow=TRUE) rc_res[['pars_df']] <- tibble(variable = c("a1","a2","a3","alpha1","alpha2","alpha3","c","lambda2","lambda3","mu2","mu3"), median = c(0.004388728,0.058483556,0.019955846,0.751815037,0.189520468,0.196011541,0.020194227, 0.138963085,0.124622576,24.984743365,64.834417626), lower = c(4.158378e-04,4.851475e-02, 1.158713e-02, 6.763524e-02, 1.478663e-01, 1.130389e-01, 1.718623e-02, 1.157654e-01, 8.611379e-02, 2.308518e+01,6.291426e+01), upper = c(0.01106290, 0.06797684, 0.02883773, 2.25796406, 0.23777804, 0.35320899, 0.02185859, 0.16919224, 0.19132286, 26.91274982, 66.84958398)) rc_res[['fit_df']] <- tibble(ages = 0:80, data = fl_migrants/fl_pop, median = c(0.02445917, 0.02212957, 0.02123731, 0.02083668, 0.02063243, 0.02050311, 0.02042182, 0.02038452, 0.02039545, 0.02048247, 0.02069491, 0.02113354, 0.02188968, 0.02309565, 0.02488725, 0.02723605, 0.03001204, 0.03295246, 0.03584502, 0.03844924, 0.04053613, 0.04205060, 0.04286610, 0.04304004, 0.04260980, 0.04172190, 0.04044603, 0.03891856, 0.03727069, 0.03555740, 0.03387381, 0.03229206, 0.03078721, 0.02940696, 0.02814990, 0.02702949, 0.02604895, 0.02518098, 0.02442590, 0.02378689, 0.02322727, 0.02275823, 0.02236749, 0.02203849, 0.02176177, 0.02154706, 0.02139283, 0.02131295, 0.02129494, 0.02137630, 0.02155494, 0.02186323, 0.02233507, 0.02302582, 0.02390995, 0.02484880, 0.02578982, 0.02669928, 0.02746516, 0.02808168, 0.02848546, 0.02864050, 0.02857913, 0.02831552, 0.02791748, 0.02741569, 0.02683604, 0.02623238, 0.02562902, 0.02503111, 0.02447808, 0.02395072, 0.02346067, 0.02302026, 0.02263513, 0.02229433, 0.02200634, 0.02176027, 0.02153727, 0.02135720, 0.02120329), lower = c(0.02063855, 0.02002963, 0.01939085, 0.01908285, 0.01891132, 0.01875816, 0.01866277, 0.01863718, 0.01869357, 0.01888990, 0.01927156, 0.01980051, 0.02055228, 0.02144848, 0.02274630, 0.02459381, 0.02706351, 0.03002805, 0.03311143, 0.03588292, 0.03798246, 0.03936151, 0.04008736, 0.04026842, 0.03995621, 0.03921550, 0.03812486, 0.03676263, 0.03525482, 0.03361536, 0.03193636, 0.03032246, 0.02890042, 0.02756915, 0.02644255, 0.02541254, 0.02453467, 0.02377473, 0.02312531, 0.02255027, 0.02203500, 0.02157904, 0.02119271, 0.02084455, 0.02056494, 0.02035391, 0.02020902, 0.02013073, 0.02013357, 0.02015945, 0.02028732, 0.02046103, 0.02075511, 0.02113461, 0.02162741, 0.02227860, 0.02317744, 0.02426672, 0.02524938, 0.02593666, 0.02630505, 0.02648453, 0.02645882, 0.02626886, 0.02594994, 0.02549666, 0.02493988, 0.02434173, 0.02367642, 0.02310194, 0.02260649, 0.02216076, 0.02179043, 0.02142710, 0.02110892, 0.02084890, 0.02060545, 0.02039068, 0.02015061, 0.01993089, 0.01975857), upper = c(0.03112025, 0.02513024, 0.02345860, 0.02269542, 0.02235882, 0.02215090, 0.02204610, 0.02199635, 0.02197468, 0.02199267, 0.02209632, 0.02249465, 0.02355002, 0.02522402, 0.02740956, 0.02995684, 0.03269014, 0.03554723, 0.03836364, 0.04098622, 0.04329547, 0.04494095, 0.04580337, 0.04591680, 0.04536656, 0.04421209, 0.04277925, 0.04113415, 0.03937585, 0.03759591, 0.03588578, 0.03424492, 0.03272378, 0.03129877, 0.02997162, 0.02875385, 0.02768118, 0.02671604, 0.02585786, 0.02511999, 0.02447777, 0.02395213, 0.02352872, 0.02319783, 0.02291037, 0.02271109, 0.02258587, 0.02251641, 0.02257225, 0.02280738, 0.02326907, 0.02389398, 0.02456863, 0.02534185, 0.02615691, 0.02704329, 0.02793647, 0.02883292, 0.02971780, 0.03038933, 0.03089967, 0.03112968, 0.03105258, 0.03067914, 0.03009405, 0.02951902, 0.02879010, 0.02812087, 0.02744805, 0.02686322, 0.02629736, 0.02577599, 0.02529982, 0.02485982, 0.02441999, 0.02402877, 0.02364263, 0.02330680, 0.02303527, 0.02278987, 0.02260402), diff_sq = c(8.846661e-08, 5.843961e-08, 1.025078e-08, 1.219354e-06, 3.058745e-08, 1.090453e-07, 2.572242e-07, 2.187438e-07, 4.297242e-06, 2.659236e-06, 1.094279e-06, 3.263517e-07, 2.847359e-07, 1.697953e-07, 2.679424e-06, 1.751785e-06, 4.531883e-07, 1.049723e-07, 2.500966e-08, 1.170591e-06, 5.733003e-09, 1.228009e-06, 6.749902e-07, 2.252466e-07, 9.086339e-08, 3.143570e-07, 5.610471e-08, 3.221258e-08, 4.280630e-06, 2.855579e-06, 9.333020e-08, 2.179993e-07, 7.670842e-07, 3.446848e-08, 7.043755e-08, 1.031716e-07, 8.037304e-07, 3.990348e-07, 2.102037e-08, 2.333985e-07, 8.139379e-08, 1.481629e-06, 1.997268e-06, 1.796781e-06, 6.141552e-08, 1.297080e-06, 1.033592e-06, 6.620450e-09, 2.596243e-06, 2.631823e-06, 3.963146e-06, 1.769188e-08, 4.493109e-07, 6.141955e-07, 9.358831e-07, 1.345288e-06, 1.966191e-07, 8.958878e-08, 2.031352e-06, 7.616767e-08, 1.309121e-06, 7.779241e-07, 1.017574e-06, 3.005237e-07, 5.738324e-08, 5.007556e-06, 4.054851e-07, 8.320548e-06, 6.090081e-08, 2.105700e-06, 1.146669e-06, 1.080280e-06, 2.275552e-06, 2.592394e-07, 5.774014e-07, 4.670081e-09, 3.316458e-07, 4.843763e-07, 6.162222e-07, 5.334144e-06, 1.514953e-06)) ## ----------------------------------------------------------------------------- rc_res[['check_converge']] ## ----------------------------------------------------------------------------- rc_res[['pars_df']] rc_res[['fit_df']] ## ----------------------------------------------------------------------------- rc_res[["fit_df"]] %>% ggplot(aes(ages, data)) + geom_point(aes(color = "data")) + geom_line(aes(x = ages, y = median, color = "fit")) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) + scale_color_manual(name = "", values = c(data = "red", fit = "black")) + ylab("migration rate")