#' Prepare inputs for TMB model.
#'
#' @param naomi_data Naomi data object
#' @param report_likelihood Option to report likelihood in fit object (default true).
#' @param anchor_home_district Option to include random effect home district attractiveness to retain residents on ART within home districts (default true).
#'
#' @return Inputs ready for TMB model
#'
#' @seealso [select_naomi_data]
#' @export
prepare_tmb_inputs <- function(naomi_data,
report_likelihood = 1L) {
stopifnot(methods::is(naomi_data, "naomi_data"))
stopifnot(methods::is(naomi_data, "naomi_mf"))
## ANC observation aggregation matrices
##
## TODO: Refactor code to make the function create_artattend_Amat() more generic.
## Should not refer to 'ART' specific; also useful for ANC attendance,
## fertility, etc.
create_anc_Amat <- function(anc_obs_dat) {
df_attend_anc <- naomi_data$mf_model %>%
dplyr::select(reside_area_id = area_id,
attend_area_id = area_id,
sex,
age_group,
idx)
dat <- dplyr::rename(anc_obs_dat, attend_area_id = area_id)
Amat <- create_artattend_Amat(
dat,
age_groups = naomi_data$age_groups,
sexes = naomi_data$sexes,
area_aggregation = naomi_data$area_aggregation,
df_art_attend = df_attend_anc,
by_residence = FALSE
)
Amat
}
create_survey_Amat <- function(survey_dat) {
df_attend_survey <- naomi_data$mf_model %>%
dplyr::select(reside_area_id = area_id,
attend_area_id = area_id,
sex,
age_group,
idx)
survey_dat$attend_area_id <- survey_dat$area_id
Amat <- create_artattend_Amat(
survey_dat,
age_groups = naomi_data$age_groups,
sexes = naomi_data$sexes,
area_aggregation = naomi_data$area_aggregation,
df_art_attend = df_attend_survey,
by_residence = FALSE,
by_survey = TRUE
)
Amat
}
A_anc_clients_t2 <- create_anc_Amat(naomi_data$anc_clients_t2_dat)
A_anc_prev_t1 <- create_anc_Amat(naomi_data$anc_prev_t1_dat)
A_anc_prev_t2 <- create_anc_Amat(naomi_data$anc_prev_t2_dat)
A_anc_artcov_t1 <- create_anc_Amat(naomi_data$anc_artcov_t1_dat)
A_anc_artcov_t2 <- create_anc_Amat(naomi_data$anc_artcov_t2_dat)
A_prev <- create_survey_Amat(naomi_data$prev_dat)
A_artcov <- create_survey_Amat(naomi_data$artcov_dat)
A_vls <- create_survey_Amat(naomi_data$vls_dat)
A_recent <- create_survey_Amat(naomi_data$recent_dat)
## ART attendance aggregation
# Default model for ART attending: Anchor home district = add random effect for home district
if(naomi_data$model_options$anchor_home_district) {
Xgamma <- sparse_model_matrix(~0 + attend_area_idf, naomi_data$mf_artattend)
} else {
Xgamma <- sparse_model_matrix(~0 + attend_area_idf:as.integer(jstar != 1),
naomi_data$mf_artattend)
}
if(naomi_data$artattend_t2) {
Xgamma_t2 <- Xgamma
} else {
Xgamma_t2 <- sparse_model_matrix(~0, naomi_data$mf_artattend)
}
df_art_attend <- naomi_data$mf_model %>%
dplyr::rename(reside_area_id = area_id) %>%
dplyr::left_join(naomi_data$mf_artattend, by = "reside_area_id",
multiple = "all",
relationship = "many-to-many") %>%
dplyr::mutate(attend_idf = forcats::as_factor(attend_idx),
idf = forcats::as_factor(idx))
Xart_gamma <- sparse_model_matrix(~0 + attend_idf, df_art_attend)
Xart_idx <- sparse_model_matrix(~0 + idf, df_art_attend)
A_artattend_t1 <- create_artattend_Amat(artnum_df = dplyr::rename(naomi_data$artnum_t1_dat, attend_area_id = area_id),
age_groups = naomi_data$age_groups,
sexes = naomi_data$sexes,
area_aggregation = naomi_data$area_aggregation,
df_art_attend = df_art_attend,
by_residence = FALSE)
A_artattend_t2 <- create_artattend_Amat(artnum_df = dplyr::rename(naomi_data$artnum_t2_dat, attend_area_id = area_id),
age_groups = naomi_data$age_groups,
sexes = naomi_data$sexes,
area_aggregation = naomi_data$area_aggregation,
df_art_attend = df_art_attend,
by_residence = FALSE)
A_artattend_mf <- create_artattend_Amat(artnum_df = dplyr::select(naomi_data$mf_model, attend_area_id = area_id, sex, age_group, artnum_idx = idx),
age_groups = naomi_data$age_groups,
sexes = naomi_data$sexes,
area_aggregation = naomi_data$area_aggregation,
df_art_attend = df_art_attend,
by_residence = FALSE)
A_art_reside_attend <- naomi_data$mf_artattend %>%
dplyr::transmute(
reside_area_id,
attend_area_id,
sex = "both",
age_group = "Y000_999"
) %>%
create_artattend_Amat(age_groups = naomi_data$age_groups,
sexes = naomi_data$sexes,
area_aggregation = naomi_data$area_aggregation,
df_art_attend = df_art_attend,
by_residence = TRUE)
## Construct TMB data and initial parameter vectors
df <- naomi_data$mf_model
X_15to49 <- Matrix::t(sparse_model_matrix(~-1 + area_idf:age15to49, naomi_data$mf_model))
## Paediatric prevalence from 15-49 female ratio
X_15to49f <- Matrix::t(Matrix::sparse.model.matrix(~0 + area_idf:age15to49:as.integer(sex == "female"), df))
df$bin_paed_rho_model <- 1 - df$bin_rho_model
X_paed_rho_ratio <- sparse_model_matrix(~-1 + area_idf:paed_rho_ratio:bin_paed_rho_model, df)
paed_rho_ratio_offset <- 0.5 * df$bin_rho_model
X_paed_lambda_ratio_t1 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t1, df)
X_paed_lambda_ratio_t2 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t2, df)
X_paed_lambda_ratio_t3 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t3, df)
X_paed_lambda_ratio_t4 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t4, df)
X_paed_lambda_ratio_t5 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t5, df)
f_rho_a <- if(all(is.na(df$rho_a_fct))) ~0 else ~0 + rho_a_fct
f_alpha_a <- if(all(is.na(df$alpha_a_fct))) ~0 else ~0 + alpha_a_fct
if (naomi_data$rho_paed_x_term) {
f_rho_xa <- ~0 + area_idf
} else {
f_rho_xa <- ~0
}
## Ratio of paediatric incidence rate to 15-49 female prevalence
## If no sex stratified prevalence data, don't estimate spatial variation in
## sex odds ratio
if ( ! all(c("male", "female") %in% naomi_data$prev_dat$sex)) {
f_rho_xs <- ~0
} else {
f_rho_xs <- ~0 + area_idf
}
## If no sex stratified ART coverage data, don't estimate spatial variation in
## sex odds ratio
if ( ! all(c("male", "female") %in% naomi_data$artcov_dat$sex) &&
! all(c("male", "female") %in% naomi_data$artnum_t1_dat$sex) &&
! all(c("male", "female") %in% naomi_data$artnum_t2_dat$sex) ) {
f_alpha_xs <- ~0
} else {
f_alpha_xs <- ~0 + area_idf
}
## If flag **and** has ART by sex data at both times, estimate time x district x
## sex ART odds ratio.
if (naomi_data$alpha_xst_term) {
if (!all(c("male", "female") %in% naomi_data$artnum_t1_dat$sex) &&
!all(c("male", "female") %in% naomi_data$artnum_t2_dat$sex)) {
stop(paste("Sex-stratified ART data are required at both Time 1 and Time 2",
"to estimate district x sex x time interaction for ART coverage"))
}
f_alpha_xst <- ~0 + area_idf
} else {
f_alpha_xst <- ~0
}
## If no ART data at both time points, do not fit a change in ART coverage. Use
## logit difference in ART coverage from Spectrum.
## T1 ART data may be either survey or programme
##
has_t1_art <- nrow(naomi_data$artcov_dat) > 0 | nrow(naomi_data$artnum_t1_dat) > 0
has_t2_art <- nrow(naomi_data$artnum_t2_dat) > 0
if( !has_t1_art | !has_t2_art ) {
f_alpha_t2 <- ~0
f_alpha_xt <- ~0
logit_alpha_t1t2_offset <- naomi_data$mf_model$logit_alpha_t1t2_offset
} else {
f_alpha_t2 <- ~1
f_alpha_xt <- ~0 + area_idf
logit_alpha_t1t2_offset <- numeric(nrow(naomi_data$mf_model))
}
## Paediatric ART coverage random effects
artnum_t1_dat <- naomi_data$artnum_t1_dat %>%
dplyr::left_join(get_age_groups(), by = "age_group") %>%
dplyr::mutate(age_group_end = age_group_start + age_group_span - 1)
artnum_t2_dat <- naomi_data$artnum_t2_dat %>%
dplyr::left_join(get_age_groups(), by = "age_group") %>%
dplyr::mutate(age_group_end = age_group_start + age_group_span - 1)
has_t1_paed_art <- any(artnum_t1_dat$age_group_end < 15)
has_t2_paed_art <- any(artnum_t2_dat$age_group_end < 15)
if(has_t1_paed_art | has_t2_paed_art) {
f_alpha_xa <- ~0 + area_idf
} else {
f_alpha_xa <- ~0
}
if(has_t1_paed_art & has_t2_paed_art) {
f_alpha_t2 <- ~1 + age_below15
f_alpha_xat <- ~0 + area_idf
} else {
f_alpha_xat <- ~0
}
## If no recent infection data, do not estimate incidence sex ratio or
## district random effects
if(nrow(naomi_data$recent_dat) == 0) {
f_lambda <- ~0
f_lambda_x <- ~0
} else {
f_lambda <- ~ 1 + female_15plus
f_lambda_x <- ~0 + area_idf
}
dtmb <- list(
population_t1 = df$population_t1,
population_t2 = df$population_t2,
Lproj_hivpop_t1t2 = naomi_data$Lproj_t1t2$Lproj_hivpop,
Lproj_incid_t1t2 = naomi_data$Lproj_t1t2$Lproj_incid,
Lproj_paed_t1t2 = naomi_data$Lproj_t1t2$Lproj_paed,
X_rho = as.matrix(sparse_model_matrix(~female_15plus, df, "bin_rho_model", TRUE)),
X_alpha = stats::model.matrix(~female_15plus, df),
X_alpha_t2 = stats::model.matrix(f_alpha_t2, df),
X_lambda = stats::model.matrix(f_lambda, df),
X_asfr = stats::model.matrix(~1, df),
X_ancrho = stats::model.matrix(~1, df),
X_ancalpha = stats::model.matrix(~1, df),
Z_x = sparse_model_matrix(~0 + area_idf, df),
Z_rho_x = sparse_model_matrix(~0 + area_idf, df, "bin_rho_model", TRUE),
Z_rho_xs = sparse_model_matrix(f_rho_xs, df, "female_15plus", TRUE),
Z_rho_a = sparse_model_matrix(f_rho_a, df, "bin_rho_model", TRUE),
Z_rho_as = sparse_model_matrix(f_rho_a, df, "female_15plus", TRUE),
Z_rho_xa = sparse_model_matrix(f_rho_xa, df, "age_below15"),
Z_alpha_x = sparse_model_matrix(~0 + area_idf, df),
Z_alpha_xs = sparse_model_matrix(f_alpha_xs, df, "female_15plus", TRUE),
Z_alpha_a = sparse_model_matrix(f_alpha_a, df),
Z_alpha_as = sparse_model_matrix(f_alpha_a, df, "female_15plus", TRUE),
Z_alpha_xt = sparse_model_matrix(f_alpha_xt, df),
Z_alpha_xa = sparse_model_matrix(f_alpha_xa, df, "age_below15"),
Z_alpha_xat = sparse_model_matrix(f_alpha_xat, df, "age_below15"),
Z_alpha_xst = sparse_model_matrix(f_alpha_xst, df, "female_15plus", TRUE),
Z_lambda_x = sparse_model_matrix(f_lambda_x, df),
## Z_xa = Matrix::sparse.model.matrix(~0 + area_idf:age_group_idf, df),
Z_asfr_x = sparse_model_matrix(~0 + area_idf, df),
Z_ancrho_x = sparse_model_matrix(~0 + area_idf, df),
Z_ancalpha_x = sparse_model_matrix(~0 + area_idf, df),
log_asfr_t1_offset = log(df$asfr_t1),
log_asfr_t2_offset = log(df$asfr_t2),
log_asfr_t3_offset = log(df$asfr_t3),
logit_anc_rho_t1_offset = log(df$frr_plhiv_t1),
logit_anc_rho_t2_offset = log(df$frr_plhiv_t2),
logit_anc_rho_t3_offset = log(df$frr_plhiv_t3),
logit_anc_alpha_t1_offset = log(df$frr_already_art_t1),
logit_anc_alpha_t2_offset = log(df$frr_already_art_t2),
logit_anc_alpha_t3_offset = log(df$frr_already_art_t3),
##
logit_rho_offset = naomi_data$mf_model$logit_rho_offset * naomi_data$mf_model$bin_rho_model,
logit_alpha_offset = naomi_data$mf_model$logit_alpha_offset,
logit_alpha_t1t2_offset = logit_alpha_t1t2_offset,
##
unaware_untreated_prop_t1 = df$spec_unaware_untreated_prop_t1,
unaware_untreated_prop_t2 = df$spec_unaware_untreated_prop_t2,
unaware_untreated_prop_t3 = df$spec_unaware_untreated_prop_t3,
##
Q_x = methods::as(naomi_data$Q, "dgCMatrix"),
Q_x_rankdef = ncol(naomi_data$Q) - as.integer(Matrix::rankMatrix(naomi_data$Q)),
n_nb = naomi_data$mf_areas$n_neighbors,
adj_i = naomi_data$mf_artattend$reside_area_idx - 1L,
adj_j = naomi_data$mf_artattend$attend_area_idx - 1L,
Xgamma = Xgamma,
Xgamma_t2 = Xgamma_t2,
log_gamma_offset = naomi_data$mf_artattend$log_gamma_offset,
Xart_idx = Xart_idx,
Xart_gamma = Xart_gamma,
##
omega = naomi_data$omega,
OmegaT0 = naomi_data$rita_param$OmegaT0,
sigma_OmegaT = naomi_data$rita_param$sigma_OmegaT,
betaT0 = naomi_data$rita_param$betaT0,
sigma_betaT = naomi_data$rita_param$sigma_betaT,
ritaT = naomi_data$rita_param$ritaT,
##
logit_nu_mean = naomi_data$logit_nu_mean,
logit_nu_sd = naomi_data$logit_nu_sd,
##
X_15to49 = X_15to49,
log_lambda_t1_offset = naomi_data$mf_model$log_lambda_t1_offset,
log_lambda_t2_offset = naomi_data$mf_model$log_lambda_t2_offset,
##
X_15to49f = X_15to49f,
X_paed_rho_ratio = X_paed_rho_ratio,
paed_rho_ratio_offset = paed_rho_ratio_offset,
##
X_paed_lambda_ratio_t1 = X_paed_lambda_ratio_t1,
X_paed_lambda_ratio_t2 = X_paed_lambda_ratio_t2,
X_paed_lambda_ratio_t3 = X_paed_lambda_ratio_t3,
X_paed_lambda_ratio_t4 = X_paed_lambda_ratio_t4,
X_paed_lambda_ratio_t5 = X_paed_lambda_ratio_t5,
##
## Household survey input data
x_prev = naomi_data$prev_dat$x_eff,
n_prev = naomi_data$prev_dat$n_eff,
A_prev = A_prev,
x_artcov = naomi_data$artcov_dat$x_eff,
n_artcov = naomi_data$artcov_dat$n_eff,
A_artcov = A_artcov,
x_vls = naomi_data$vls_dat$x_eff,
n_vls = naomi_data$vls_dat$n_eff,
A_vls = A_vls,
x_recent = naomi_data$recent_dat$x_eff,
n_recent = naomi_data$recent_dat$n_eff,
A_recent = A_recent,
##
## ANC testing input data
x_anc_clients_t2 = naomi_data$anc_clients_t2_dat$anc_clients_x,
offset_anc_clients_t2 = naomi_data$anc_clients_t2_dat$anc_clients_pys_offset,
A_anc_clients_t2 = A_anc_clients_t2,
x_anc_prev_t1 = naomi_data$anc_prev_t1_dat$anc_prev_x,
n_anc_prev_t1 = naomi_data$anc_prev_t1_dat$anc_prev_n,
A_anc_prev_t1 = A_anc_prev_t1,
x_anc_artcov_t1 = naomi_data$anc_artcov_t1_dat$anc_artcov_x,
n_anc_artcov_t1 = naomi_data$anc_artcov_t1_dat$anc_artcov_n,
A_anc_artcov_t1 = A_anc_artcov_t1,
x_anc_prev_t2 = naomi_data$anc_prev_t2_dat$anc_prev_x,
n_anc_prev_t2 = naomi_data$anc_prev_t2_dat$anc_prev_n,
A_anc_prev_t2 = A_anc_prev_t2,
x_anc_artcov_t2 = naomi_data$anc_artcov_t2_dat$anc_artcov_x,
n_anc_artcov_t2 = naomi_data$anc_artcov_t2_dat$anc_artcov_n,
A_anc_artcov_t2 = A_anc_artcov_t2,
##
## Number on ART input data
A_artattend_t1 = A_artattend_t1,
x_artnum_t1 = naomi_data$artnum_t1_dat$art_current,
A_artattend_t2 = A_artattend_t2,
x_artnum_t2 = naomi_data$artnum_t2_dat$art_current,
A_artattend_mf = A_artattend_mf,
A_art_reside_attend = A_art_reside_attend,
##
## Time 3 projection inputs
population_t3 = df$population_t3,
Lproj_hivpop_t2t3 = naomi_data$Lproj_t2t3$Lproj_hivpop,
Lproj_incid_t2t3 = naomi_data$Lproj_t2t3$Lproj_incid,
Lproj_paed_t2t3 = naomi_data$Lproj_t2t3$Lproj_paed,
logit_alpha_t2t3_offset = df$logit_alpha_t2t3_offset,
log_lambda_t3_offset = df$log_lambda_t3_offset,
##
## Time 4 projection inputs
population_t4 = df$population_t4,
Lproj_hivpop_t3t4 = naomi_data$Lproj_t3t4$Lproj_hivpop,
Lproj_incid_t3t4 = naomi_data$Lproj_t3t4$Lproj_incid,
Lproj_paed_t3t4 = naomi_data$Lproj_t3t4$Lproj_paed,
logit_alpha_t3t4_offset = df$logit_alpha_t3t4_offset,
log_lambda_t4_offset = df$log_lambda_t4_offset,
##
## Time 5 projection inputs
population_t5 = df$population_t5,
Lproj_hivpop_t4t5 = naomi_data$Lproj_t4t5$Lproj_hivpop,
Lproj_incid_t4t5 = naomi_data$Lproj_t4t5$Lproj_incid,
Lproj_paed_t4t5 = naomi_data$Lproj_t4t5$Lproj_paed,
logit_alpha_t4t5_offset = df$logit_alpha_t4t5_offset,
log_lambda_t5_offset = df$log_lambda_t5_offset,
##
A_out = naomi_data$A_out,
A_anc_out = naomi_data$A_anc_out,
calc_outputs = 1L,
report_likelihood = report_likelihood
)
ptmb <- list(
beta_rho = numeric(ncol(dtmb$X_rho)),
beta_alpha = numeric(ncol(dtmb$X_alpha)),
beta_alpha_t2 = numeric(ncol(dtmb$X_alpha_t2)),
beta_lambda = numeric(ncol(dtmb$X_lambda)),
beta_asfr = numeric(1),
beta_anc_rho = numeric(1),
beta_anc_alpha = numeric(1),
beta_anc_rho_t2 = numeric(1),
beta_anc_alpha_t2 = numeric(1),
u_rho_x = numeric(ncol(dtmb$Z_rho_x)),
us_rho_x = numeric(ncol(dtmb$Z_rho_x)),
u_rho_xs = numeric(ncol(dtmb$Z_rho_xs)),
us_rho_xs = numeric(ncol(dtmb$Z_rho_xs)),
u_rho_a = numeric(ncol(dtmb$Z_rho_a)),
u_rho_as = numeric(ncol(dtmb$Z_rho_as)),
u_rho_xa = numeric(ncol(dtmb$Z_rho_xa)),
ui_asfr_x = numeric(ncol(dtmb$Z_asfr_x)),
ui_anc_rho_x = numeric(ncol(dtmb$Z_ancrho_x)),
ui_anc_alpha_x = numeric(ncol(dtmb$Z_ancalpha_x)),
ui_anc_rho_xt = numeric(ncol(dtmb$Z_ancrho_x)),
ui_anc_alpha_xt = numeric(ncol(dtmb$Z_ancalpha_x)),
##
u_alpha_x = numeric(ncol(dtmb$Z_alpha_x)),
us_alpha_x = numeric(ncol(dtmb$Z_alpha_x)),
u_alpha_xs = numeric(ncol(dtmb$Z_alpha_xs)),
us_alpha_xs = numeric(ncol(dtmb$Z_alpha_xs)),
u_alpha_a = numeric(ncol(dtmb$Z_alpha_a)),
u_alpha_as = numeric(ncol(dtmb$Z_alpha_as)),
u_alpha_xt = numeric(ncol(dtmb$Z_alpha_xt)),
u_alpha_xa = numeric(ncol(dtmb$Z_alpha_xa)),
u_alpha_xat = numeric(ncol(dtmb$Z_alpha_xat)),
u_alpha_xst = numeric(ncol(dtmb$Z_alpha_xst)),
##
log_sigma_lambda_x = log(1.0),
ui_lambda_x = numeric(ncol(dtmb$Z_lambda_x)),
##
logit_phi_rho_a = 0,
log_sigma_rho_a = log(2.5),
logit_phi_rho_as = 2.582,
log_sigma_rho_as = log(2.5),
logit_phi_rho_x = 0,
log_sigma_rho_x = log(2.5),
logit_phi_rho_xs = 0,
log_sigma_rho_xs = log(2.5),
log_sigma_rho_xa = log(0.5),
##
logit_phi_alpha_a = 0,
log_sigma_alpha_a = log(2.5),
logit_phi_alpha_as = 2.582,
log_sigma_alpha_as = log(2.5),
logit_phi_alpha_x = 0,
log_sigma_alpha_x = log(2.5),
logit_phi_alpha_xs = 0,
log_sigma_alpha_xs = log(2.5),
log_sigma_alpha_xt = log(2.5),
log_sigma_alpha_xa = log(2.5),
log_sigma_alpha_xat = log(2.5),
log_sigma_alpha_xst = log(2.5),
##
OmegaT_raw = 0,
log_betaT = 0,
logit_nu_raw = 0,
##
log_sigma_asfr_x = log(0.5),
log_sigma_ancrho_x = log(2.5),
log_sigma_ancalpha_x = log(2.5),
log_sigma_ancrho_xt = log(2.5),
log_sigma_ancalpha_xt = log(2.5),
##
log_or_gamma = numeric(ncol(dtmb$Xgamma)),
log_sigma_or_gamma = log(2.5),
log_or_gamma_t1t2 = numeric(ncol(dtmb$Xgamma_t2)),
log_sigma_or_gamma_t1t2 = log(2.5)
)
v <- list(data = dtmb,
par_init = ptmb)
class(v) <- "naomi_tmb_input"
v
}
sparse_model_matrix <- function(formula, data, binary_interaction = 1,
drop_zero_cols = FALSE) {
if(is.character(binary_interaction))
binary_interaction <- data[[binary_interaction]]
stopifnot(length(binary_interaction) %in% c(1, nrow(data)))
mm <- Matrix::sparse.model.matrix(formula, data)
mm <- mm * binary_interaction
mm <- Matrix::drop0(mm)
if(drop_zero_cols)
mm <- mm[ , apply(mm, 2, Matrix::nnzero) > 0]
mm
}
make_tmb_obj <- function(data, par, calc_outputs = 1L, inner_verbose = FALSE,
progress = NULL) {
data$calc_outputs <- as.integer(calc_outputs)
obj <- TMB::MakeADFun(data = data,
parameters = par,
DLL = "naomi",
silent = !inner_verbose,
random = c("beta_rho",
"beta_alpha", "beta_alpha_t2",
"beta_lambda",
"beta_asfr",
"beta_anc_rho", "beta_anc_alpha",
"beta_anc_rho_t2", "beta_anc_alpha_t2",
"u_rho_x", "us_rho_x",
"u_rho_xs", "us_rho_xs",
"u_rho_a", "u_rho_as",
"u_rho_xa",
##
"u_alpha_x", "us_alpha_x",
"u_alpha_xs", "us_alpha_xs",
"u_alpha_a", "u_alpha_as",
"u_alpha_xt",
"u_alpha_xa", "u_alpha_xat", "u_alpha_xst",
##
"ui_lambda_x",
"logit_nu_raw",
##
"ui_asfr_x",
"ui_anc_rho_x", "ui_anc_alpha_x",
"ui_anc_rho_xt", "ui_anc_alpha_xt",
##
"log_or_gamma", "log_or_gamma_t1t2"))
if (!is.null(progress)) {
obj$fn <- report_progress(obj$fn, progress)
}
obj
}
report_progress <- function(fun, progress) {
fun <- match.fun(fun)
function(...) {
progress$iterate_fit()
fun(...)
}
}
#' Fit TMB model
#'
#' @param tmb_input Model input data
#' @param outer_verbose If TRUE print function and parameters every iteration
#' @param inner_verbose If TRUE then disable tracing information from TMB
#' @param max_iter maximum number of iterations
#' @param progress Progress printer, if null no progress printed
#'
#' @return Fit model.
#' @export
fit_tmb <- function(tmb_input,
outer_verbose = TRUE,
inner_verbose = FALSE,
max_iter = 250,
progress = NULL
) {
stopifnot(inherits(tmb_input, "naomi_tmb_input"))
obj <- make_tmb_obj(tmb_input$data, tmb_input$par_init, calc_outputs = 0L,
inner_verbose, progress)
trace <- if (outer_verbose) 1 else 0
f <- withCallingHandlers(
stats::nlminb(obj$par, obj$fn, obj$gr,
control = list(trace = trace,
iter.max = max_iter)),
warning = function(w) {
if (grepl("NA/NaN function evaluation", w$message)) {
invokeRestart("muffleWarning")
}
}
)
if(f$convergence != 0)
warning(paste("convergence error:", f$message))
if(outer_verbose)
message(paste("converged:", f$message))
f$par.fixed <- f$par
f$par.full <- obj$env$last.par
objout <- make_tmb_obj(tmb_input$data, tmb_input$par_init, calc_outputs = 1L, inner_verbose)
f$mode <- objout$report(f$par.full)
val <- c(f, obj = list(objout))
class(val) <- "naomi_fit"
val
}
#' Calculate Posterior Mean and Uncertainty Via TMB `sdreport()`
#'
#' @param naomi_fit Fitted TMB model.
#'
#' @export
report_tmb <- function(naomi_fit) {
stopifnot(methods::is(fit, "naomi_fit"))
naomi_fit$sdreport <- TMB::sdreport(naomi_fit$obj, naomi_fit$par,
getReportCovariance = FALSE,
bias.correct = TRUE)
naomi_fit
}
#' Sample TMB fit
#'
#' @param fit The TMB fit
#' @param nsample Number of samples
#' @param rng_seed seed passed to set.seed.
#' @param random_only Random only
#' @param verbose If TRUE prints additional information.
#'
#' @return Sampled fit.
#' @export
sample_tmb <- function(fit, nsample = 1000, rng_seed = NULL,
random_only = TRUE, verbose = FALSE) {
set.seed(rng_seed)
stopifnot(methods::is(fit, "naomi_fit"))
stopifnot(nsample > 1)
to_tape <- TMB:::isNullPointer(fit$obj$env$ADFun$ptr)
if (to_tape)
fit$obj$retape(FALSE)
if(!random_only) {
if(verbose) print("Calculating joint precision")
hess <- sdreport_joint_precision(fit$obj, fit$par.fixed)
if(verbose) print("Inverting precision for joint covariance")
cov <- solve(hess)
if(verbose) print("Drawing sample")
## TODO: write a version of rmvnorm that uses precision instead of covariance
smp <- mvtnorm::rmvnorm(nsample, fit$par.full, cov)
} else {
r <- fit$obj$env$random
par_f <- fit$par.full[-r]
par_r <- fit$par.full[r]
hess_r <- fit$obj$env$spHess(fit$par.full, random = TRUE)
smp_r <- rmvnorm_sparseprec(nsample, par_r, hess_r)
smp <- matrix(0, nsample, length(fit$par.full))
smp[ , r] <- smp_r
smp[ ,-r] <- matrix(par_f, nsample, length(par_f), byrow = TRUE)
colnames(smp)[r] <- colnames(smp_r)
colnames(smp)[-r] <- names(par_f)
}
if(verbose) print("Simulating outputs")
sim <- apply(smp, 1, fit$obj$report)
r <- fit$obj$report()
if(verbose) print("Returning sample")
fit$sample <- Map(vapply, list(sim), "[[", lapply(lengths(r), numeric), names(r))
is_vector <- vapply(fit$sample, inherits, logical(1), "numeric")
fit$sample[is_vector] <- lapply(fit$sample[is_vector], matrix, nrow = 1)
names(fit$sample) <- names(r)
fit
}
rmvnorm_sparseprec <- function(
n,
mean = rep(0, nrow(prec)),
prec = diag(length(mean))
) {
z = matrix(stats::rnorm(n * length(mean)), ncol = n)
L_inv = Matrix::Cholesky(prec)
v <- mean + Matrix::solve(methods::as(L_inv, "pMatrix"), Matrix::solve(Matrix::t(methods::as(L_inv, "Matrix")), z))
as.matrix(Matrix::t(v))
}
create_artattend_Amat <- function(artnum_df, age_groups, sexes, area_aggregation,
df_art_attend, by_residence = FALSE, by_survey = FALSE) {
## If by_residence = TRUE, merge by reside_area_id, else aggregate over all
## reside_area_id
by_vars <- c("attend_area_id", "sex", "age_group")
if (by_residence) {
by_vars <- c(by_vars, "reside_area_id")
}
id_vars <- by_vars
if (by_survey) {
id_vars <- c(id_vars, "survey_id")
}
if(!("artnum_idx" %in% colnames(artnum_df))) {
artnum_df$artnum_idx <- seq_len(nrow(artnum_df))
}
A_artnum <- artnum_df %>%
dplyr::select(tidyselect::all_of(id_vars), artnum_idx) %>%
dplyr::rename(artdat_age_group = age_group,
artdat_sex = sex) %>%
dplyr::left_join(
get_age_groups() %>%
dplyr::transmute(
artdat_age_group = age_group,
artdat_age_start = age_group_start,
artdat_age_end = age_group_start + age_group_span
),
by = "artdat_age_group",
relationship = "many-to-many"
) %>%
## Note: this would be much faster with tree data structure for age rather than crossing...
tidyr::crossing(
get_age_groups() %>%
dplyr::filter(age_group %in% age_groups)
) %>%
dplyr::filter(
artdat_age_start <= age_group_start,
age_group_start + age_group_span <= artdat_age_end
) %>%
dplyr::left_join(
data.frame(artdat_sex = c("male", "female", "both", "both", "both"),
sex = c("male", "female", "male", "female", "both"),
stringsAsFactors = FALSE) %>%
dplyr::filter(sex %in% sexes),
by = "artdat_sex",
multiple = "all",
relationship = "many-to-many"
)
## Map artattend_area_id to model_area_id
A_artnum <- A_artnum %>%
dplyr::left_join(
area_aggregation,
by = c("attend_area_id" = "area_id"),
multiple = "all",
relationship = "many-to-many"
) %>%
dplyr::mutate(attend_area_id = model_area_id,
model_area_id = NULL)
## Check no areas with duplicated reporting
art_duplicated_check <- A_artnum %>%
dplyr::group_by_at(id_vars) %>%
dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
dplyr::filter(n > 1)
if (nrow(art_duplicated_check)) {
stop(paste("ART or ANC data multiply reported for some age/sex strata in areas:",
paste(unique(art_duplicated_check$attend_area_id), collapse = ", ")))
}
## Merge to ART attendance data frame
df_art_attend <- df_art_attend %>%
dplyr::select(tidyselect::all_of(by_vars)) %>%
dplyr::mutate(
Aidx = dplyr::row_number(),
value = 1
)
A_artnum <- dplyr::left_join(A_artnum, df_art_attend, by = by_vars, multiple = "all")
A_artnum <- A_artnum %>%
{
Matrix::spMatrix(nrow(artnum_df),
nrow(df_art_attend),
.$artnum_idx,
.$Aidx,
.$value)
}
A_artnum
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.