inst/doc/convergence_issues.R

In rcbayes: Estimate Rogers-Castro Migration Age Schedules with Bayesian Models

## ----setup, include = FALSE, eval=FALSE---------------------------------------
# knitr::opts_chunk$set(
#   collapse = TRUE,
#   comment = "#>",
#   message = FALSE,
#   fig.width = 8,
#   fig.height = 4,
#   fig.align = "center"
# )

## ----include=FALSE------------------------------------------------------------
set.seed(123)

## -----------------------------------------------------------------------------
library(rcbayes)
library(tibble)
library(ggplot2)

ages <- 0:80
migrants <- c(11804, 10606, 9845, 9244, 8471, 7940, 7348, 6885, 6431,
              6055, 5454, 4997, 4845, 4596, 4397, 4814, 4592, 4646, 5386,
              7180, 11374, 14713, 17195, 18937, 19223, 19091, 18507,
              17615, 16929, 15693, 15246, 14152, 13365, 12340, 11609,
              10278, 9547, 8992, 8438, 7883, 7315, 6909, 6730, 6272,
              5994, 6087, 5896, 5592, 5487, 5237, 6021, 5933, 5577,
              5674, 5503, 4916, 5008, 4822, 4824, 4696, 4086, 4019,
              4139, 4054, 4134, 3625, 3871, 4238, 4306, 4440, 3118,
              2980, 2885, 2845, 2795, 2085, 2076, 2035, 2030, 1986, 2037)
pop <- c(105505, 105505, 105505, 105505, 105505, 106126, 106126, 106126,
         106126, 106126, 100104, 100104, 100104, 100104, 100104, 114880,
         114880, 114880, 114880, 114880, 136845, 136845, 136845, 136845,
         136845, 136582, 136582, 136582, 136582, 136582, 141935, 141935,
         141935, 141935, 141935, 134097, 134097, 134097, 134097, 134097,
         130769, 130769, 130769, 130769, 130769, 133718, 133718, 133718,
         133718, 133718, 154178, 154178, 154178, 154178, 154178, 145386,
         145386, 145386, 145386, 145386, 126270, 126270, 126270, 126270,
         126270, 108314, 108314, 108314, 108314, 108314, 79827, 79827,
         79827, 79827, 79827, 59556, 59556, 59556, 59556, 59556, 59556)

## ----eval=FALSE---------------------------------------------------------------
# rc_res <- mig_estimate_rc(
#   ages, migrants, pop,
#   pre_working_age = TRUE,
#   working_age = TRUE,
#   retirement = TRUE,
#   post_retirement = TRUE
# )

## ----echo=FALSE---------------------------------------------------------------
rc_res = list(check_converge = c(), pars_df = c(), fit_df = c())

rc_res[['check_converge']] <- matrix(
  c(
    0.113605451436614,0.00779593103661137,2.15774607128397,3.70790733560849,
    0.111764533173589,0.00797819680947733,2.18281974401689,3.4397804041693,
    0.200212843072318,0.0476641666967773,2.42785538791501,2.33364604863079,
    0.0860134303532162,0.0029298454713518,2.4448397222948,2.36661380639638,
    0.206149088628833,0.00434306449617798,2.62766374221068,2.08317823108183,
    0.011611407086162,0.00465354274179327,2.09013982530816,4.7075734702239,
    0.0111607261406471,0.00204934887455636,7.11402640579904,1.21806784550916,
    21.1903479818411,0.110825507910137,2.42600816481533,2.40555966266519,
    65.8824541176426,0.527963547098352,2.99425165755755,1.73668210871192,
    0.391006779161944,0.00755861905096238,2.98501838004935,1.75178240898006,
    0.564868092984617,0.211502782189408,2.35090673278088,2.66219716061898,
    0.00291870586352614,0.00744665931477216,2.96818132803141,1.90946024259034,
    0.0140624198563911,0.00535502339245935,2.76119187029469,1.92225295693418),
  ncol=4,
  dimnames = list(c("alpha1[1]", "alpha2[1]", "alpha3[1]", 
                    "a1[1]", "a2[1]", "a3[1]","a4[1]",
                    "mu2[1]", "mu3[1]", 
                    "lambda2[1]", "lambda3[1]","lambda4[1]", "c"),
                  c("mean", "se_mean", "n_eff", "Rhat")),
  byrow = TRUE)

rc_res[['fit_df']] <- tibble(age = 0:80,
       data = migrants/pop,
       median = c(0.11107413, 0.10207562, 0.09392648, 0.08667344, 0.08025095, 0.07451119, 
                  0.06935217, 0.06473068, 0.06059616, 0.05690437, 0.05363063, 0.05073012, 
                  0.04815361, 0.04582636, 0.04373611, 0.04189403, 0.04047357, 0.04090839, 
                  0.04697327, 0.06189863, 0.08392274, 0.10695466, 0.12542619, 0.13688149, 
                  0.14149304, 0.14075261, 0.13643752, 0.13001158, 0.12257184, 0.11481948, 
                  0.10718039, 0.09990469, 0.09312533, 0.08687808, 0.08115890, 0.07595000, 
                  0.07122528, 0.06694617, 0.06308948, 0.05961437, 0.05648621, 0.05366093, 
                  0.05111955, 0.04883263, 0.04678603, 0.04495830, 0.04332511, 0.04186645, 
                  0.04056364, 0.03940147, 0.03836733, 0.03744961, 0.03664648, 0.03594472, 
                  0.03532742, 0.03476267, 0.03425199, 0.03379800, 0.03340861, 0.03307914, 
                  0.03279761, 0.03256650, 0.03236997, 0.03221412, 0.03225893, 0.03356208, 
                  0.03610573, 0.03875565, 0.03972603, 0.03947603, 0.03873864, 0.03777993, 
                  0.03683512, 0.03600479, 0.03562471, 0.03506771, 0.03456024, 0.03413901, 
                  0.03382805, 0.03361884, 0.03349211),
       lower = c(0.10955889, 0.10097274, 0.09318158, 0.08588609, 0.07923866, 0.07338551, 
                 0.06821759, 0.06368101, 0.05970448, 0.05620466, 0.05309501, 0.05028441, 
                 0.04768464, 0.04534359, 0.04322938, 0.04133397, 0.03981681, 0.04017004, 
                 0.04617863, 0.06108230, 0.08293964, 0.10596553, 0.12447152, 0.13592858, 
                 0.14058096, 0.13989090, 0.13564572, 0.12925827, 0.12181410, 0.11405511, 
                 0.10645688, 0.09922248, 0.09247883, 0.08626839, 0.08058662, 0.07542530, 
                 0.07076293, 0.06652748, 0.06270095, 0.05924734, 0.05613392, 0.05330736, 
                 0.05075355, 0.04842025, 0.04629537, 0.04440086, 0.04273021, 0.04124867, 
                 0.03995819, 0.03883248, 0.03787124, 0.03703954, 0.03633343, 0.03569461, 
                 0.03507823, 0.03450711, 0.03399487, 0.03351492, 0.03305895, 0.03264036, 
                 0.03225994, 0.03191384, 0.03162075, 0.03139179, 0.03142643, 0.03237460, 
                 0.03424834, 0.03437457, 0.03452557, 0.03467760, 0.03484352, 0.03502117, 
                 0.03520709, 0.03535974, 0.03483698, 0.03431782, 0.03391187, 0.03354558, 
                 0.03316710, 0.03275000, 0.03232411),
       upper = c(0.11341196, 0.10324869, 0.09465255, 0.08733431, 0.08086402, 0.07510810, 
                 0.06994563, 0.06533161, 0.06118399, 0.05747201, 0.05414250, 0.05118015,
                 0.04861464, 0.04649246, 0.04464422, 0.04301912, 0.04171672, 0.04194436, 
                 0.04780617, 0.06280798, 0.08490775, 0.10794241, 0.12647352, 0.13794525, 
                 0.14251906, 0.14171389, 0.13726329, 0.13077593, 0.12327584, 0.11550739, 
                 0.10784189, 0.10052317, 0.09368927, 0.08740651, 0.08168648, 0.07651546, 
                 0.07181934, 0.06755270, 0.06365780, 0.06011183, 0.05690655, 0.05400963, 
                 0.05145158, 0.04918898, 0.04716633, 0.04535589, 0.04372723, 0.04227591, 
                 0.04096728, 0.03979244, 0.03873892, 0.03778779, 0.03694441, 0.03620099, 
                 0.03560906, 0.03520374, 0.03488964, 0.03464470, 0.03446929, 0.03434677, 
                 0.03427521, 0.03423904, 0.03424497, 0.03428457, 0.03435298, 0.03469481, 
                 0.03713969, 0.03975826, 0.04080284, 0.04042104, 0.03946999, 0.03845057, 
                 0.03761703, 0.03686507, 0.03617822, 0.03629931, 0.03654233, 0.03679272, 
                 0.03707064, 0.03734264, 0.03761274))

## -----------------------------------------------------------------------------
rc_res[['check_converge']]

## ----fig.width=6--------------------------------------------------------------
rc_res[["fit_df"]] %>%
  ggplot(aes(age, data)) +
  geom_point(aes(color = "data")) +
  geom_line(aes(x = age, 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")

## ----eval=FALSE---------------------------------------------------------------
# res <- mig_estimate_rc(ages, migrants, pop,
#                        pre_working_age = TRUE,
#                        working_age = TRUE,
#                        retirement = TRUE,
#                        post_retirement = TRUE,
#                        #optional inputs into stan
#                        iter = 3000,
#                        control = list(adapt_delta = 0.95, max_treedepth = 10)
#                        )

## ----eval=FALSE---------------------------------------------------------------
# init_vals <- init_rc(ages,
#                      migrants,
#                      pop,
#                      pre_working_age = TRUE,
#                      working_age = TRUE,
#                      retirement = TRUE,
#                      post_retirement = TRUE,
#                      nchains = 4)
# 
# res <- mig_estimate_rc(ages, migrants, pop,
#                        pre_working_age = TRUE,
#                        working_age = TRUE,
#                        retirement = TRUE,
#                        post_retirement = TRUE,
#                        #optional inputs into stan
#                        init = init_vals
#                        )

## ----eval=FALSE---------------------------------------------------------------
# init_vals <- init_rc(ages,
#                      migrants,
#                      pop,
#                      pre_working_age = TRUE,
#                      working_age = TRUE,
#                      retirement = TRUE,
#                      post_retirement = TRUE,
#                      nchains = 4)
# 
# rc_res_fixed <- mig_estimate_rc(ages, migrants, pop,
#                        pre_working_age = TRUE,
#                        working_age = TRUE,
#                        retirement = TRUE,
#                        post_retirement = TRUE,
#                        #optional inputs into stan
#                        control = list(adapt_delta = 0.95, max_treedepth = 10),
#                        init = init_vals
# )

## ----echo=FALSE---------------------------------------------------------------
rc_res_fixed = list(check_converge = c(), pars_df = c(), fit_df = c())

rc_res_fixed[['check_converge']] <- matrix(
  c(0.107162666095038,0.000113905396688904,914.0993478393,1.00234765459143,
    0.103416351522329,5.64917230559274e-05,853.927068031998,1.00281875605411,
    0.210669831871199,0.00140706757932151,988.257923785188,1.00150360579644,
    0.0878408694500156,4.26344484838096e-05,1025.19360646823,1.00125328173498,
    0.202957753238002,5.51135303627361e-05,1096.26976516497,1.00141331160793,
    0.016091838651906,6.14366866169559e-05,1022.81686620825,1.00149544086794,
    0.0107940563428325,0.000202777673594207,770.930151592325,1.00187319942316,
    21.0673530496958,0.00160357875882577,914.9772707573,1.00197537852488,
    66.5184716633019,0.0108069832903751,982.132724052394,1.00157929557214,
    0.395991098815759,0.00017884423246015,1158.66065253182,1.00080606722439,
    0.744339527000941,0.00373523428808142,1503.88004401573,1.00163128898661,
    0.0091660363074148,0.000171349903624075,623.526448086348,1.00399958735813,
    0.0120475931079795,0.000236317198704928,716.808119784298,1.0020118971368),
  ncol=4,
  dimnames = list(c("alpha1[1]", "alpha2[1]", "alpha3[1]", 
                    "a1[1]", "a2[1]", "a3[1]","a4[1]",
                    "mu2[1]", "mu3[1]", 
                    "lambda2[1]", "lambda3[1]","lambda4[1]", "c"),
                  c("mean", "se_mean", "n_eff", "Rhat")),
  byrow = TRUE)

rc_res_fixed[['fit_df']] <- tibble(age = 0:80,
       data = migrants/pop,
       median = c(0.11066577, 0.10183467, 0.09389724, 0.08677859, 0.08039617, 0.07466633,
                  0.06952911, 0.06491991, 0.06079010, 0.05708701, 0.05377364, 0.05079940, 
                  0.04814050, 0.04576222, 0.04363252, 0.04173797, 0.04023607, 0.04059125, 
                  0.04677098, 0.06199728, 0.08416159, 0.10711504, 0.12541505, 0.13670053, 
                  0.14124400, 0.14053808, 0.13631138, 0.12999417, 0.12263667, 0.11494149, 
                  0.10733717, 0.10006698, 0.09326341, 0.08697700, 0.08122151, 0.07598656, 
                  0.07123679, 0.06694529, 0.06307296, 0.05958576, 0.05644594, 0.05362702, 
                  0.05109293, 0.04882184, 0.04678431, 0.04495985, 0.04332688, 0.04186791, 
                  0.04056457, 0.03940128, 0.03836677, 0.03744697, 0.03663026, 0.03590659, 
                  0.03526751, 0.03470475, 0.03421152, 0.03378090, 0.03340640, 0.03308209, 
                  0.03280512, 0.03256866, 0.03237309, 0.03221376, 0.03217486, 0.03312048, 
                  0.03615995, 0.03904294, 0.04014689, 0.03983118, 0.03889165, 0.03782808, 
                  0.03683831, 0.03599620, 0.03530893, 0.03475818, 0.03433355, 0.03401210, 
                  0.03377150, 0.03360034, 0.03348575),
       lower = c(0.10938343, 0.10085839, 0.09315967, 0.08614682, 0.07983181, 0.07413997, 
                 0.06901742, 0.06442625, 0.06031352, 0.05663560, 0.05335241, 0.05039448, 
                 0.04773767, 0.04535350, 0.04321303, 0.04128219, 0.03975127, 0.04009192, 
                 0.04605109, 0.06116725, 0.08333080, 0.10622850, 0.12446871, 0.13582829, 
                 0.14044398, 0.13978288, 0.13557999, 0.12927923, 0.12196578, 0.11434229, 
                 0.10678598, 0.09956663, 0.09280558, 0.08654702, 0.08081728, 0.07559846, 
                 0.07086616, 0.06659222, 0.06273188, 0.05925050, 0.05611922, 0.05331224, 
                 0.05079336, 0.04853664, 0.04651009, 0.04469189, 0.04306340, 0.04161878, 
                 0.04031990, 0.03916481, 0.03813495, 0.03721814, 0.03640352, 0.03567912, 
                 0.03503863, 0.03447041, 0.03397472, 0.03353608, 0.03314746, 0.03280966, 
                 0.03250885, 0.03224883, 0.03202669, 0.03183978, 0.03180072, 0.03236501, 
                 0.03520049, 0.03831441, 0.03940962, 0.03915478, 0.03828029, 0.03723175, 
                 0.03624018, 0.03539759, 0.03471905, 0.03420878, 0.03381553, 0.03348453, 
                 0.03321844, 0.03300026, 0.03281042),
       upper = c(0.11202162, 0.10283789, 0.09466222, 0.08740818, 0.08093686, 0.07517373, 
                 0.07001886, 0.06539917, 0.06126140, 0.05753726, 0.05420355, 0.05121903, 
                 0.04855045, 0.04618349, 0.04406297, 0.04219851, 0.04071295, 0.04111381, 
                 0.04748327, 0.06285125, 0.08503755, 0.10802562, 0.12635079, 0.13758453, 
                 0.14204044, 0.14127723, 0.13701374, 0.13068959, 0.12331207, 0.11556133, 
                 0.10790406, 0.10057925, 0.09373162, 0.08740192, 0.08162284, 0.07636317, 
                 0.07160537, 0.06729989, 0.06341591, 0.05992296, 0.05677216, 0.05393834, 
                 0.05139425, 0.04911275, 0.04706269, 0.04522956, 0.04358717, 0.04211525, 
                 0.04080241, 0.03963269, 0.03859004, 0.03766736, 0.03684888, 0.03612591, 
                 0.03548708, 0.03492513, 0.03443622, 0.03401193, 0.03364986, 0.03333821, 
                 0.03307808, 0.03285901, 0.03268092, 0.03254230, 0.03254832, 0.03405741, 
                 0.03701209, 0.03986310, 0.04093148, 0.04053379, 0.03951791, 0.03839881, 
                 0.03739697, 0.03653334, 0.03584743, 0.03529488, 0.03485950, 0.03452575, 
                 0.03431082, 0.03420497, 0.03418890))

## -----------------------------------------------------------------------------
rc_res_fixed[['check_converge']]

## ----fig.width=6--------------------------------------------------------------
rc_res_fixed[["fit_df"]] %>%
  ggplot(aes(age, data)) +
  geom_point(aes(color = "data")) +
  geom_line(aes(x = age, 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")