diff --git a/DESCRIPTION b/DESCRIPTION index d3f3bb5..4766acb 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: serosv Type: Package Title: Model Infectious Disease Parameters from Serosurveys -Version: 1.2.0 +Version: 1.2.0.9000 Authors@R: c( person("Anh", "Phan Truong Quynh", email = "anhptq@oucru.org", role = c("aut", "cre"), comment = c(ORCID = "0009-0000-2129-435X")), person("Nguyen", "Pham Nguyen The", email = "nguyenpnt@oucru.org", role = c("aut"), comment = c(ORCID = "0000-0002-0356-2776")), @@ -40,6 +40,7 @@ Imports: rstan (>= 2.18.1), rstantools (>= 2.4.0), boot, + pROC, stats4 Suggests: covr, @@ -49,16 +50,13 @@ Suggests: testthat (>= 3.0.0) Collate: 'data.R' - 'mseir_model.R' - 'sir_basic_model.R' - 'sir_static_model.R' - 'sir_subpops_model.R' 'fractional_polynomial_models.R' 'polynomial_models.R' 'utils.R' 'compare_models.R' 'correct_prevalence.R' 'weibull_model.R' + 'farrington_model.R' 'nonparametric.R' 'semiparametric_models.R' 'mixture_model.R' @@ -69,6 +67,8 @@ Collate: 'plots.R' 'compute_ci.R' 'age_time_model.R' + 'predict.R' + 'print.R' Config/testthat/edition: 3 URL: https://oucru-modelling.github.io/serosv/, https://github.com/OUCRU-Modelling/serosv VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 80b1540..ff45bcb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,13 @@ # Generated by roxygen2: do not edit by hand +S3method(compute_ci,age_time_model) +S3method(compute_ci,default) +S3method(compute_ci,fp_model) +S3method(compute_ci,hierarchical_bayesian_model) +S3method(compute_ci,lp_model) +S3method(compute_ci,mixture_model) +S3method(compute_ci,penalized_spline_model) +S3method(compute_ci,weibull_model) S3method(plot,age_time_model) S3method(plot,estimate_from_mixture) S3method(plot,farrington_model) @@ -7,34 +15,38 @@ S3method(plot,fp_model) S3method(plot,hierarchical_bayesian_model) S3method(plot,lp_model) S3method(plot,mixture_model) -S3method(plot,mseir_model) S3method(plot,penalized_spline_model) S3method(plot,polynomial_model) -S3method(plot,sir_basic_model) -S3method(plot,sir_static_model) -S3method(plot,sir_subpops_model) S3method(plot,weibull_model) +S3method(predict,age_time_model) +S3method(predict,farrington_model) +S3method(predict,fp_model) +S3method(predict,hierarchical_bayesian_model) +S3method(predict,lp_model) +S3method(predict,penalized_spline_model) +S3method(predict,polynomial_model) +S3method(predict,weibull_model) +S3method(print,age_time_model) +S3method(print,estimate_from_mixture) +S3method(print,farrington_model) +S3method(print,fp_model) +S3method(print,hierarchical_bayesian_model) +S3method(print,lp_model) +S3method(print,mixture_model) +S3method(print,penalized_spline_model) +S3method(print,polynomial_model) +S3method(print,weibull_model) export(add_thresholds) export(age_time_model) export(compare_models) -export(compute_ci) -export(compute_ci.age_time_model) -export(compute_ci.fp_model) -export(compute_ci.hierarchical_bayesian_model) -export(compute_ci.lp_model) -export(compute_ci.mixture_model) -export(compute_ci.penalized_spline_model) -export(compute_ci.weibull_model) export(correct_prevalence) export(est_foi) export(estimate_from_mixture) export(farrington_model) -export(find_best_fp_powers) export(fp_model) export(hierarchical_bayesian_model) export(lp_model) export(mixture_model) -export(mseir_model) export(pava) export(penalized_spline_model) export(plot_corrected_prev) @@ -43,9 +55,6 @@ export(plot_standard_curve) export(plot_titer_qc) export(polynomial_model) export(set_plot_style) -export(sir_basic_model) -export(sir_static_model) -export(sir_subpops_model) export(standardize_data) export(to_titer) export(transform_data) @@ -55,15 +64,15 @@ import(assertthat) import(dplyr) import(ggplot2) import(graphics) +import(locfit) import(methods) +import(pROC) import(patchwork) -import(purrr) import(scam) import(tidyr) importFrom(RcppParallel,RcppParallelLibs) importFrom(assertthat,assert_that) importFrom(boot,inv.logit) -importFrom(deSolve,ode) importFrom(dplyr,group_by) importFrom(dplyr,mutate) importFrom(dplyr,n) @@ -87,7 +96,17 @@ importFrom(mixdist,mix) importFrom(mixdist,mixgroup) importFrom(mixdist,mixparam) importFrom(mvtnorm,rmvnorm) +importFrom(purrr,as_mapper) +importFrom(purrr,compact) importFrom(purrr,imap_dfr) +importFrom(purrr,map) +importFrom(purrr,map2) +importFrom(purrr,map_dbl) +importFrom(purrr,map_dfc) +importFrom(purrr,map_dfr) +importFrom(purrr,pmap) +importFrom(purrr,pmap_dfr) +importFrom(purrr,walk) importFrom(rstan,sampling) importFrom(rstan,summary) importFrom(rstantools,rstan_config) @@ -103,6 +122,9 @@ importFrom(stats,predict.glm) importFrom(stats,prop.test) importFrom(stats,qnorm) importFrom(stats,qt) +importFrom(stats4,AIC) +importFrom(stats4,BIC) +importFrom(stats4,logLik) importFrom(stats4,mle) importFrom(stringr,str_detect) useDynLib(serosv, .registration=TRUE) diff --git a/NEWS.md b/NEWS.md index 76a8a94..f8c5a13 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,10 @@ +# serosv (development version) +* clean up documentation + # serosv 1.2.0 +* add to_titer() function to convert assay reading to titer +* update the docs +* update article for to_titer() on the website # serosv 1.1.0 * add correct_prevalence() function to estimate real prevalence from imperfect test diff --git a/R/age_time_model.R b/R/age_time_model.R index 9e75d98..a2c9a64 100644 --- a/R/age_time_model.R +++ b/R/age_time_model.R @@ -3,15 +3,21 @@ #' #' @description Fit age-stratified seroprevalence across multiple time points. Also try to monotonize age (or birth cohort) - specific seroprevalence. #' -#' @param data - input data, must have`age`, `status`, time, group columns, where group column determines how data is aggregated -#' @param time_col - name of the column for time (default to `date`) -#' @param grouping_col - name of the column for time (default to `group`) -#' @param age_correct - a boolean, if `TRUE`, monotonize age-specific prevalence. Monotonize birth cohort-specific seroprevalence otherwise. -#' @param le - number of bins to generate age grid, used when monotonizing data -#' @param ci - confidence interval for smoothing -#' @param monotonize_method - either "pava" or "scam" +#' @param data input data, must have age, status, time, group columns, where group column determines how data is aggregated +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). +#' @param time_col name of the column for time (default to "date") +#' @param grouping_col name of the column for time (default to "group") +#' @param age_correct a boolean, if `TRUE`, monotonize age-specific prevalence. Monotonize birth cohort-specific seroprevalence otherwise. +#' @param le number of bins to generate age grid, used when monotonizing data +#' @param ci confidence interval for smoothing +#' @param monotonize_method either "pava" or "scam" +#' #' @import scam assertthat #' @importFrom mgcv gam predict.gam betar +#' @importFrom purrr map map_dbl map2 #' #' @return a list of class time_age_model with 4 items #' \item{out}{a data.frame with dimension n_group x 9, where columns `info`, `sp`, `foi` store output for non-monotonized @@ -20,7 +26,10 @@ #' \item{age_correct}{a boolean indicating whether the data is monotonized across age or cohort} #' \item{datatype}{whether the input data is aggregated or line-listing data} #' @export -age_time_model <- function(data, time_col="date", grouping_col="group", age_correct=F, le=512, ci = 0.95, monotonize_method = "pava"){ +age_time_model <- function(data, + age_col="age", status_col="status", pos_col="pos", tot_col="tot", + time_col="date", grouping_col="group", + age_correct=F, le=512, ci = 0.95, monotonize_method = "pava"){ # work around to resolve no visible binding note NOTE during check() x <- label <- family <- fit <- se.fit <- ymin <- ymax <- y <- mean_time <- prevalence <- sim_data <- NULL age <- ys <- shift_no <- cohort <- col_time <- monotonized_mod <- df <- info <- sp <- monotonized_info <- monotonized_sp <- NULL @@ -82,7 +91,8 @@ age_time_model <- function(data, time_col="date", grouping_col="group", age_corr model <- list() # --- preprocess data ------ - check_input <- check_input(data) + check_input <- check_input(data, + stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) age_range <- range(data$age) age_grid <- seq(age_range[1], age_range[2], length.out = le) @@ -230,8 +240,10 @@ age_time_model <- function(data, time_col="date", grouping_col="group", age_corr ) model$out <- out + model$monotonize_method <- monotonize_method model$grouping_col <- grouping_col model$age_correct <- age_correct + model$ci <- ci class(model) <- "age_time_model" diff --git a/R/compare_models.R b/R/compare_models.R index 681b521..fb6100c 100644 --- a/R/compare_models.R +++ b/R/compare_models.R @@ -1,37 +1,151 @@ #' Compare models #' +#' @param data input data to fit into the models +#' @param method method to compare models. Can be one of the built-in methods or a function to compute the returned metrics (see Details). #' @param ... models to be compared. Must be models created by serosv. If models' names are not provided, indices will be used instead for the `model` column in the returned data.frame. #' #' #' @return #' a data.frame of 4 columns -#' \item{model}{name or index of the model} +#' \item{label}{name or index of the model} #' \item{type}{model type of the given model (a serosv model name)} #' \item{AIC}{AIC value for the model (lower value indicates better fit)} #' \item{BIC}{BIC value for the model (lower value indicates better fit)} #' +#' @details +#' Built-in comparison methods include: +#' - computing AIC and BIC, which returns AIC, BIC values of the model if available +#' - cross validation, which returns MSE and logloss (negative ) +#' #' @importFrom magrittr %>% -#' @importFrom purrr imap_dfr +#' @importFrom purrr imap_dfr as_mapper #' @importFrom stringr str_detect +#' @importFrom assertthat assert_that #' #' @export -compare_models <- function(...){ +compare_models <- function(data, method="AIC/BIC",...){ list(...) %>% imap_dfr(~ { # return error if input contains non-serosv models - if(!all(str_detect(class(.x), "_model"))) { - stop("Inputs must be serosv models") + # if(!all(str_detect(class(.x), "_model"))) { + # stop("Inputs must be serosv models") + # } + + # get function to compute comparison metrics + metric_func <- if(is.character(method)){ + switch( + method, + "AIC/BIC" = aic_bic, + "CV" = cv, + method + ) + }else{ + method } - # TODO: apply cross validation to get matrix instead + assert_that(is.function(metric_func), + msg = "Function to compute the metrics must be provided") + + out <- metric_func(data, as_mapper(.x)) - data.frame( - model = .y, - type = class(.x), - AIC = AIC(.x$info), - BIC = BIC(.x$info) + assert_that("data.frame" %in% class(out), + msg = "Function to compute the metrics must return a data.frame") + + out %>% mutate( + label = .y, + .before = 1 ) + }) } +# function returning goodness-of-fit metrics such as AIC/BIC, likelihood (with degree-of-freedom) +#' @importFrom stats4 logLik AIC BIC +aic_bic <- function(dat, mod_func){ + out <- mod_func(dat) + + aic <- tryCatch(stats4::AIC(out$info), + error = \(e){NULL}) + bic <- tryCatch(stats4::BIC(out$info), + error = \(e){NULL}) + ll <- tryCatch(stats4::logLik(out$info), + error = \(e){NULL}) + + tibble( + type = class(out), + AIC = aic, + BIC = bic, + logLik = if (!is.null(ll)) as.numeric(ll) else NA, + df = if (!is.null(ll) && !is.null(attr(ll, "df"))) attr(ll, "df") else NA, + mod_out = list(out), + plots = list(plot(out)+ggtitle(paste("Fitted model using", class(out)))) + ) +} + + +# function to compute metrics from cross validation +# assess the generalization/prediction of the model +#' @importFrom stats4 logLik AIC BIC +#' @importFrom stats predict.glm +#' @import tidyr dplyr pROC +cv <- function(dat, mod_func, k=4){ + # assign each row of data to each fold + idx_fold <- sort(rep(1:k, length.out=nrow(dat))) + + metrics <- lapply(1:k, \(fold){ + curr_metric <- list() + + # split data + fit_dat <- dat[idx_fold != fold, ] + test_dat <- dat[idx_fold == fold, ] + + # get model info + out <- mod_func(fit_dat) + curr_metric$type <- class(out) + # generate prediction + pred <- predict(out, data.frame(age=test_dat[,1]), type="response") + + if(out$datatype == "aggregated"){ + # if data is aggregated + seroprev_obs <- test_dat$pos/test_dat$tot + + # MSE + curr_metric$mse <- sum((pred - seroprev_obs)**2)/nrow(test_dat) + # compute logloss (negative binomial loglikelihood) + curr_metric$logloss <- -sum( + dbinom(test_dat$pos, test_dat$tot, prob=pred, log=TRUE), + na.rm = TRUE) + + }else{ + # if data is linelisting + # make sure pred is slightly higher than 0 and lower than 1 to avoid log(0) + eps <- .Machine$double.eps # smallest positive floating point number ~ 2.2e-16 + pred <- pmax(eps, pmin(1 - eps, pred)) + + # compute logloss (negative bernoulli loglikelihood) + curr_metric$logloss <- -sum(test_dat$status*log(pred) + (1 - test_dat$status)*log(1-pred), + na.rm = TRUE) + # and estimate auc + curr_metric$auc <- as.numeric(pROC::auc(test_dat$status, pred)) + } + + curr_metric + }) %>% + bind_rows() %>% + summarise( + # compute average of the metrics + across(where(is.numeric), mean), + # for type, simply get the first one + type = first(type)) + + # also return the model when it is fitted using the whole data + out <- mod_func(dat) + + metrics %>% + mutate( + mod_out = list(out), + plots = list(plot(out)+ggtitle(paste("Fitted model using", class(out)))) + ) +} + diff --git a/R/compute_ci.R b/R/compute_ci.R index 54228e5..21bc310 100644 --- a/R/compute_ci.R +++ b/R/compute_ci.R @@ -2,12 +2,12 @@ compute_ci <- function(x, ci = 0.95, le = 100, ...){ UseMethod("compute_ci") } -#' Compute confidence interval +#' Compute confidence interval for a model of serosv #' -#' @param x - serosv models -#' @param ci - confidence interval -#' @param le - number of data for computing confidence interval -#' @param ... - arbitrary argument +#' @param x serosv models +#' @param ci confidence interval +#' @param le number of data for computing confidence interval +#' @param ... arbitrary argument #' #' @importFrom stats qt predict.glm #' @import dplyr @@ -15,7 +15,7 @@ compute_ci <- function(x, ci = 0.95, le = 100, ...){ #' @return confidence interval dataframe with 4 variables, x and y for the fitted values and ymin and ymax for the confidence interval #' #' @export -compute_ci <- function(x, ci = 0.95, le = 100, ...){ +compute_ci.default <- function(x, ci = 0.95, le = 100, ...){ # resolve no visible binding issue with CRAN check fit <- se.fit <- NULL @@ -23,10 +23,10 @@ compute_ci <- function(x, ci = 0.95, le = 100, ...){ link_inv <- x$info$family$linkinv dataset <- x$info$data n <- nrow(dataset) - length(x$info$coefficients) - age_range <- range(dataset$Age) + age_range <- range(dataset$age) ages <- seq(age_range[1], age_range[2], le = le) - mod1 <- predict.glm(x$info,data.frame(Age = ages), se.fit = TRUE) + mod1 <- predict.glm(x$info,data.frame(age = ages), se.fit = TRUE) n1 <- mod1 %>% as_tibble() %>% select(fit, se.fit) %>% mutate(age = ages ) %>% mutate(lwr = link_inv(fit + qt( p, n) * se.fit), @@ -40,10 +40,10 @@ compute_ci <- function(x, ci = 0.95, le = 100, ...){ #' Compute confidence interval for fractional polynomial model #' -#' @param x - serosv models -#' @param ci - confidence interval -#' @param le - number of data for computing confidence interval -#' @param ... - arbitrary argument +#' @param x serosv models +#' @param ci confidence interval +#' @param le number of data for computing confidence interval +#' @param ... arbitrary argument #' #' @import dplyr #' @return confidence interval dataframe with 4 variables, x and y for the fitted values and ymin and ymax for the confidence interval @@ -74,9 +74,9 @@ compute_ci.fp_model <- function(x, ci = 0.95, le = 100, ...){ #' Compute confidence interval for Weibull model #' -#' @param x - serosv models -#' @param ci - confidence interval -#' @param ... - arbitrary argument +#' @param x serosv models +#' @param ci confidence interval +#' @param ... arbitrary argument #' #' @import dplyr #' @return confidence interval dataframe with 4 variables, x and y for the fitted values and ymin and ymax for the confidence interval @@ -109,15 +109,15 @@ compute_ci.weibull_model <- function(x, ci = 0.95, ...){ #' Compute confidence interval for local polynomial model #' -#' @param x - serosv models -#' @param ci - confidence interval -#' @param ... - arbitrary arguments +#' @param x serosv models +#' @param ci confidence interval +#' @param ... arbitrary arguments #' @return confidence interval dataframe with 4 variables, x and y for the fitted values and ymin and ymax for the confidence interval #' @export compute_ci.lp_model <- function(x,ci = 0.95, ...){ ages <- x$df$age - crit<- crit(x$pi,cov = ci)$crit.val - mod1 <- predict(x$pi, data.frame(a = ages),se.fit = TRUE) + crit<- crit(x$info,cov = ci)$crit.val + mod1 <- predict(x$info, data.frame(a = ages),se.fit = TRUE) out.DF <- data.frame(x = ages, y = mod1$fit,ymin= mod1$fit-crit*(mod1$se.fit/100), ymax= mod1$fit+crit*(mod1$se.fit/100)) out.DF @@ -126,9 +126,9 @@ compute_ci.lp_model <- function(x,ci = 0.95, ...){ #' Compute confidence interval for penalized_spline_model #' -#' @param x - serosv models -#' @param ci - confidence interval -#' @param ... - arbitrary arguments +#' @param x serosv models +#' @param ci confidence interval +#' @param ... arbitrary arguments #' @importFrom mgcv predict.gam #' @import dplyr #' @@ -181,8 +181,8 @@ compute_ci.penalized_spline_model <- function(x,ci = 0.95, ...){ #' Compute 95\% credible interval for hierarchical Bayesian model #' -#' @param x - serosv models -#' @param ... - arbitrary arguments +#' @param x serosv models +#' @param ... arbitrary arguments #' @importFrom mgcv predict.gam #' @import dplyr #' @@ -233,9 +233,9 @@ compute_ci.hierarchical_bayesian_model <- function(x, ...){ #' Compute confidence interval for mixture model #' -#' @param x - serosv mixture_model object -#' @param ci - confidence interval -#' @param ... - arbitrary arguments +#' @param x serosv mixture_model object +#' @param ci confidence interval +#' @param ... arbitrary arguments #' @importFrom stats qnorm #' #' @return list of confidence interval for susceptible and infected. Each confidence interval is a list with 2 items for lower and upper bound of the interval. @@ -258,10 +258,10 @@ compute_ci.mixture_model <- function(x,ci = 0.95, ...){ #' Compute confidence interval for time age model #' -#' @param x - serosv models -#' @param ci - confidence interval -#' @param le - number of data for computing confidence interval -#' @param ... - arbitrary argument +#' @param x serosv models +#' @param ci confidence interval +#' @param le number of data for computing confidence interval +#' @param ... arbitrary argument #' #' @importFrom mgcv predict.gam #' @import dplyr @@ -340,6 +340,8 @@ compute_ci.age_time_model <- function(x, ci=0.95, le = 100, ...){ }) ) %>% select(!!sym(x$grouping_col), sp_df, foi_df) + + out } diff --git a/R/correct_prevalence.R b/R/correct_prevalence.R index 2bc626d..9659749 100644 --- a/R/correct_prevalence.R +++ b/R/correct_prevalence.R @@ -1,6 +1,6 @@ #' Estimate the true sero prevalence using Frequentist/Bayesian estimation #' -#' @param data the input data frame, must either have `age`, `pos`, `tot` columns (for aggregated data) OR `age`, `status` for (linelisting data) +#' @param data the input data frame, must either have columns for `age`, `pos`, `tot` (for aggregated data) OR `age`, `status` (for linelisting data) #' @param bayesian whether to adjust sero-prevalence using the Bayesian or frequentist approach. If set to `TRUE`, true sero-prevalence is estimated using MCMC. #' @param init_se sensitivity of the serological test #' @param init_sp specificity of the serological test @@ -9,6 +9,10 @@ #' @param chains (applicable when `bayesian=TRUE`) number of Markov chains #' @param warmup (applicable when `bayesian=TRUE`) number of warm up runs #' @param iter (applicable when `bayesian=TRUE`) number of iterations +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). #' #' @importFrom rstan sampling summary #' @importFrom dplyr mutate @@ -17,7 +21,7 @@ #' #' @return a list of 3 items #' \item{info}{estimated parameters (when `bayesian = TRUE`) or formula to compute corrected prevalence (when `bayesian = FALSE`)} -#' \item{df}{data.frame of input data (in aggregated form)} +#' \item{df}{data.frame of input data (in aggregated form) with the 95\% confidence interval for apparent (i.e. observed) seroprevalence} #' \item{corrected_sero}{data.frame containing age, the corresponding estimated seroprevalance with 95\% confidence/credible interval, and adjusted tot and pos} #' @export #' @@ -25,6 +29,7 @@ #' data <- rubella_uk_1986_1987 #' correct_prevalence(data) correct_prevalence <- function(data, bayesian=TRUE, + age_col="age",pos_col="pos", tot_col="tot", status_col="status", init_se = 0.95, init_sp = 0.8, study_size_se = 1000, study_size_sp = 1000, chains = 1, warmup = 1000, iter = 2000){ # resolve no visible binding note @@ -32,18 +37,26 @@ correct_prevalence <- function(data, bayesian=TRUE, output <- list() - data <- check_input(data) + data <- check_input(data, stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) age <- data$age pos <- data$pos tot <- data$tot if (data$type == "linelisting"){ - transform_df <- transform_data(age, pos) + transform_df <- transform_data(data, stratum_col="age", status_col="pos") age <- transform_df$t pos <- transform_df$pos tot <- transform_df$tot } + # quantify CI for apparent seroprev + apparent_thetas <- data.frame( + # use prop.test to get confidence interval + lower = mapply(\(pos, tot){prop.test(pos, tot, conf.level = 0.95)$conf.int[1]}, pos, tot), + fit = pos/tot, + upper = mapply(\(pos, tot){prop.test(pos, tot, conf.level = 0.95)$conf.int[2]}, pos, tot) + ) + if (bayesian){ # format data @@ -68,12 +81,7 @@ correct_prevalence <- function(data, bayesian=TRUE, upper = summary(fit)$summary[3:(length(age) + 2), "97.5%"] ) }else{ - thetas <- data.frame( - # use prop.test to get confidence interval - lower = mapply(\(pos, tot){prop.test(pos, tot, conf.level = 0.95)$conf.int[1]}, pos, tot), - fit = pos/tot, - upper = mapply(\(pos, tot){prop.test(pos, tot, conf.level = 0.95)$conf.int[2]}, pos, tot) - ) %>% + thetas <- apparent_thetas %>% mutate( # estimate true prevalence and lower, upper bound lower = pmax(0, (lower - 1 + init_sp)/(init_se + init_sp - 1)), @@ -88,14 +96,17 @@ correct_prevalence <- function(data, bayesian=TRUE, output$df <- data.frame( age = age, pos = pos, - tot = tot + tot = tot, + sero = apparent_thetas$fit, + sero_lwr = apparent_thetas$lower, + sero_upr = apparent_thetas$upper ) output$corrected_se <- data.frame( age = age, - sero = thetas$fit, pos = thetas$fit*tot, #adjusted pos with estimated sero tot = tot, + sero = thetas$fit, sero_lwr = thetas$lower, sero_upr = thetas$upper ) diff --git a/R/data.R b/R/data.R index 50931cb..ea9248d 100644 --- a/R/data.R +++ b/R/data.R @@ -7,6 +7,9 @@ # data are normalized and cleaned for the purpose of this package +# ===== Hepatitis A ======= + + #' Hepatitis A serological data from Belgium in 1993 and 1994 (aggregated) #' #' A study of the prevalence of HAV antibodies conducted in the Flemish @@ -20,14 +23,12 @@ #' } #' #' @examples -#' # Reproduce Fig 4.1 (upper left panel), p. 63 -#' age <- hav_be_1993_1994$age -#' pos <- hav_be_1993_1994$pos -#' tot <- hav_be_1993_1994$tot -#' plot( -#' age, pos / tot, -#' pty = "s", cex = 0.06 * tot, pch = 16, xlab = "age", -#' ylab = "seroprevalence", xlim = c(0, 86), ylim = c(0, 1) +#' with(hav_be_1993_1994, +#' plot( +#' age, pos / tot, +#' pty = "s", cex = 0.34 * sqrt(tot), pch = 16, xlab = "age", +#' ylab = "seroprevalence", xlim = c(0, 86), ylim = c(0, 1) +#' ) #' ) #' #' @source Beutels, M., Van Damme, P., Aelvoet, W. et al. Prevalence of @@ -48,15 +49,18 @@ #' } #' #' @examples -#' # Reproduce Fig 4.1 (upper right panel), p. 63 #' library(dplyr) #' df <- hav_be_2002 %>% #' group_by(age) %>% #' summarise(pos = sum(seropositive), tot = n()) -#' plot( -#' df$age, df$pos / df$tot, -#' pty = "s", cex = 0.06 * df$tot, pch = 16, xlab = "age", -#' ylab = "seroprevalence", xlim = c(0, 86), ylim = c(0, 1) +#' +#' with( +#' df, +#' plot( +#' age, pos / tot, +#' pty = "s", cex = 0.34 * sqrt(tot), pch = 16, xlab = "age", +#' ylab = "seroprevalence", xlim = c(0, 86), ylim = c(0, 1) +#' ) #' ) #' #' @source Thiry, N., Beutels, P., Shkedy, Z. et al. The seroepidemiology @@ -78,14 +82,13 @@ #' } #' #' @examples -#' # Reproduce Fig 4.1 (lower panel), p. 63 -#' age <- hav_bg_1964$age -#' pos <- hav_bg_1964$pos -#' tot <- hav_bg_1964$tot -#' plot( +#' with( +#' hav_bg_1964, +#' plot( #' age, pos / tot, -#' pty = "s", cex = 0.08 * tot, pch = 16, xlab = "age", +#' pty = "s", cex = 0.6 * sqrt(tot), pch = 16, xlab = "age", #' ylab = "seroprevalence", xlim = c(0, 86), ylim = c(0, 1) +#' ) #' ) #' #' @source Keiding, Niels. "Age-Specific Incidence and Prevalence: A @@ -94,6 +97,8 @@ #' \doi{doi:10.2307/2983150} "hav_bg_1964" +# ===== Hepatitis B ======= + #' Hepatitis B serological data from Russia in 1999 (aggregated) #' #' A seroprevalence study conducted in St. Petersburg @@ -108,7 +113,6 @@ #' } #' #' @examples -#' # Reproduce Fig 4.2, p. 65 #' library(dplyr) #' hbv_ru_1999$age <- trunc(hbv_ru_1999$age / 1) * 1 #' hbv_ru_1999$age[hbv_ru_1999$age > 40] <- trunc( @@ -117,9 +121,10 @@ #' df <- hbv_ru_1999 %>% #' group_by(age) %>% #' summarise(pos = sum(pos), tot = sum(tot)) +#' #' plot( #' df$age, df$pos / df$tot, -#' cex = 0.05 * df$tot, pch = 16, xlab = "age", +#' cex = 0.32 * sqrt(df$tot), pch = 16, xlab = "age", #' ylab = "seroprevalence", xlim = c(0, 72) #' ) #' @@ -129,6 +134,7 @@ #' St Petersburg, 2000. "hbv_ru_1999" +# ===== Hepatitis C ======= #' Hepatitis C serological data from Belgium in 2006 (line listing) #' #' A study of HCV infection among injecting drug users. All injecting drug users @@ -143,7 +149,6 @@ #' } #' #' @examples -#' # Reproduce Fig 4.3, p. 66 #' library(dplyr) #' # snapping age to aggregated age group #' # (credit: https://stackoverflow.com/a/12861810) @@ -162,11 +167,13 @@ #' group_by(dur) %>% #' summarise(tot = n(), pos = sum(seropositive)) #' -#' plot( -#' hcv_be_2006$dur, hcv_be_2006$pos / hcv_be_2006$tot, -#' cex = 0.1 * hcv_be_2006$tot, pch = 16, -#' xlab = "duration of injection (years)", -#' ylab = "seroprevalence", xlim = c(0, 25), ylim = c(0, 1) +#' with(hcv_be_2006, +#' plot( +#' dur, pos / tot, +#' cex = 0.42 * sqrt(tot), pch = 16, +#' xlab = "duration of injection (years)", +#' ylab = "seroprevalence", xlim = c(0, 25), ylim = c(0, 1) +#' ) #' ) #' #' @source Mathei, C., Shkedy, Z., Denis, B., Kabali, C., Aerts, @@ -177,6 +184,8 @@ #' \doi{doi:10.1111/j.1365-2893.2006.00725.x} "hcv_be_2006" +# ===== Mumps ======= + #' Mumps serological data from the UK in 1986 and 1987 (aggregated) #' #' a large survey of prevalence of antibodies to mumps and rubella viruses in @@ -185,19 +194,17 @@ #' #' @format A data frame with 3 variables: #' \describe{ -#' \item{age}{Age group} +#' \item{age}{Midpoint of the age group (e.g. 1.5 = 1-2 years old, 2.5 = 2-3 years old)} #' \item{pos}{Number of seropositive individuals} #' \item{tot}{Total number of individuals surveyed} #' } #' #' @examples -#' # Reproduce Fig 4.4 (left panel), p. 67 -#' age <- mumps_uk_1986_1987$age -#' pos <- mumps_uk_1986_1987$pos -#' tot <- mumps_uk_1986_1987$tot -#' plot(age, pos / tot, -#' cex = 0.008 * tot, pch = 16, xlab = "age", ylab = "seroprevalence", -#' xlim = c(0, 45), ylim = c(0, 1) +#' with(mumps_uk_1986_1987, +#' plot(age, pos / tot, +#' cex = 0.1 * sqrt(tot), pch = 16, xlab = "age", ylab = "seroprevalence", +#' xlim = c(0, 45), ylim = c(0, 1) +#' ) #' ) #' #' @source Morgan-Capner P, Wright J, Miller C L, Miller E. Surveillance of @@ -205,6 +212,8 @@ #' 1988; 297 :770 \doi{doi:10.1136/bmj.297.6651.770} "mumps_uk_1986_1987" +# ===== Parvo 19 ======= + #' Parvo B19 serological data from Belgium from 2001-2003 (line listing) #' #' A seroprevalence survey testing for parvovirus B19 IgG antibody, performed on @@ -222,13 +231,12 @@ #' } #' #' @examples -#' # Reproduce Fig 4.5 (left upper panel), p. 68 #' library(dplyr) #' df <- parvob19_be_2001_2003 %>% #' group_by(age) %>% #' summarise(pos = sum(seropositive), tot = n()) #' plot(df$age, df$pos / df$tot, -#' cex = 0.02 * df$tot, pch = 16, xlab = "age", ylab = "seroprevalence", +#' cex = 0.3 * sqrt(df$tot), pch = 16, xlab = "age", ylab = "seroprevalence", #' xlim = c(0, 82), ylim = c(0, 1) #' ) #' @@ -256,17 +264,17 @@ #' } #' #' @examples -#' # Reproduce Fig 4.5 (right upper panel), p. 68 -#' # NB: This figure will look different to that of in the book, since we +#' # Note: This figure will look different to that of in the book, since we #' # believe that the original authors has made some errors in specifying #' # the sample size of the dots. +#' #' library(dplyr) #' df <- parvob19_ew_1996 %>% #' mutate(age = round(age)) %>% #' group_by(age) %>% #' summarise(pos = sum(seropositive), tot = n()) #' plot(df$age, df$pos / df$tot, -#' cex = 0.02 * df$tot, pch = 16, xlab = "age", ylab = "seroprevalence", +#' cex = 0.3 * sqrt(df$tot), pch = 16, xlab = "age", ylab = "seroprevalence", #' xlim = c(0, 82), ylim = c(0, 1) #' ) #' @@ -294,17 +302,17 @@ #' } #' #' @examples -#' # Reproduce Fig 4.5 (left bottom panel), p. 68 -#' # NB: This figure will look different to that of in the book, since we +#' # Note: This figure will look different to that of in the book, since we #' # believe that the original authors has made some errors in specifying #' # the sample size of the dots. +#' #' library(dplyr) #' df <- parvob19_fi_1997_1998 %>% #' mutate(age = round(age)) %>% #' group_by(age) %>% #' summarise(pos = sum(seropositive), tot = n()) #' plot(df$age, df$pos / df$tot, -#' cex = 0.07 * df$tot, pch = 16, xlab = "age", ylab = "seroprevalence", +#' cex = 0.4 * sqrt(df$tot), pch = 16, xlab = "age", ylab = "seroprevalence", #' xlim = c(0, 82), ylim = c(0, 1) #' ) #' @@ -332,16 +340,12 @@ #' } #' #' @examples -#' # Reproduce Fig 4.5 (middle bottom panel), p. 68 -#' # NB: This figure will look different to that of in the book, since we -#' # believe that the original authors has made some errors in specifying -#' # the sample size of the dots. #' library(dplyr) #' df <- parvob19_it_2003_2004 %>% #' group_by(age) %>% #' summarise(pos = sum(seropositive), tot = n()) #' plot(df$age, df$pos / df$tot, -#' cex = 0.07 * df$tot, pch = 16, xlab = "age", ylab = "seroprevalence", +#' cex = 0.32 * sqrt(df$tot), pch = 16, xlab = "age", ylab = "seroprevalence", #' xlim = c(0, 82), ylim = c(0, 1) #' ) #' @@ -369,17 +373,17 @@ #' } #' #' @examples -#' # Reproduce Fig 4.5 (right bottom panel), p. 68 -#' # NB: This figure will look different to that of in the book, since we +#' # Note: This figure will look different to that of in the book, since we #' # believe that the original authors has made some errors in specifying #' # the sample size of the dots. +#' #' library(dplyr) #' df <- parvob19_pl_1995_2004 %>% #' mutate(age = round(age)) %>% #' group_by(age) %>% #' summarise(pos = sum(seropositive), tot = n()) #' plot(df$age, df$pos / df$tot, -#' cex = 0.07 * df$tot, pch = 16, xlab = "age", ylab = "seroprevalence", +#' cex = 0.32 * sqrt(df$tot), pch = 16, xlab = "age", ylab = "seroprevalence", #' xlim = c(0, 82), ylim = c(0, 1) #' ) #' @@ -390,6 +394,9 @@ #' \doi{doi:10.1017/S0950268807009661} "parvob19_pl_1995_2004" + +# ===== Rubella ======= + #' Rubella serological data from the UK in 1986 and 1987 (aggregated) #' #' Prevalence of rubella in the UK, obtained from a large survey of prevalence @@ -397,19 +404,17 @@ #' #' @format A data frame with 3 variables: #' \describe{ -#' \item{age}{Age group} +#' \item{age}{Midpoint of the age group (e.g. 1.5 = 1-2 years old, 2.5 = 2-3 years old)} #' \item{pos}{Number of seropositive individuals} #' \item{tot}{Total number of individuals surveyed} #' } #' #' @examples -#' # Reproduce Fig 4.4 (middle panel), p. 67 -#' age <- rubella_uk_1986_1987$age -#' pos <- rubella_uk_1986_1987$pos -#' tot <- rubella_uk_1986_1987$tot -#' plot(age, pos / tot, -#' cex = 0.008 * tot, pch = 16, xlab = "age", ylab = "seroprevalence", -#' xlim = c(0, 45), ylim = c(0, 1) +#' with(rubella_uk_1986_1987, +#' plot(age, pos / tot, +#' cex = 0.2 * sqrt(tot), pch = 16, xlab = "age", ylab = "seroprevalence", +#' xlim = c(0, 45), ylim = c(0, 1) +#' ) #' ) #' #' @source Morgan-Capner P, Wright J, Miller C L, Miller E. Surveillance of @@ -417,6 +422,8 @@ #' 1988; 297 :770 \doi{doi:10.1136/bmj.297.6651.770} "rubella_uk_1986_1987" +# ===== Tuberculosis ======= + #' Tuberculosis serological data from the Netherlands 1966-1973 (aggregated) #' #' A study of tuberculosis conducted in the Netherlands. Schoolchildren, aged @@ -432,19 +439,17 @@ #' } #' #' @examples -#' # Reproduce Fig 4.6, p.70 -#' age <- tb_nl_1966_1973$age -#' birthyr <- tb_nl_1966_1973$birthyr -#' pos <- tb_nl_1966_1973$pos -#' tot <- tb_nl_1966_1973$tot -#' # left panel -#' plot(age, pos / tot, -#' pch = 16, cex = 0.00005 * tot, xlab = "age", -#' ylab = "prevalence", xlim = c(6, 18) +#' with(tb_nl_1966_1973, +#' plot(age, pos / tot, +#' pch = 16, cex = 0.01 * sqrt(tot), xlab = "age", +#' ylab = "prevalence", xlim = c(6, 18) +#' ) #' ) -#' # right panel -#' plot(birthyr, pos / tot, -#' pch = 16, cex = 0.00005 * tot, xlab = "year", ylab = "prevalence" +#' +#' with(tb_nl_1966_1973, +#' plot(birthyr, pos / tot, +#' pch = 16, cex = 0.01 * sqrt(tot), xlab = "year", ylab = "prevalence" +#' ) #' ) #' #' @source Nagelkerke, N., Heisterkamp, S., Borgdorff, M., Broekmans, J. and @@ -453,6 +458,9 @@ #' \doi{doi:10.1002/(SICI)1097-0258(19990215)18:3<307::AID-SIM15>3.0.CO;2-Z} "tb_nl_1966_1973" + +# ===== VZV ======= + #' VZV serological data from Belgium (Flanders) from 1999-2000 (aggregated) #' #' Age-specific seroprevalence of VZV antibodies, assessed in Flanders (Belgium) @@ -467,15 +475,12 @@ #' } #' #' @examples -#' # Reproduce Fig 4.7 (left panel), p.71 -#' age <- vzv_be_1999_2000$age -#' pos <- vzv_be_1999_2000$pos -#' tot <- vzv_be_1999_2000$tot -#' plot(age, pos / tot, -#' cex = 0.036 * tot, pch = 19, xlab = "age", ylab = "seroprevalence", -#' xlim = c(0, 45), ylim = c(0, 1) +#' with(vzv_be_1999_2000, +#' plot(age, pos / tot, +#' cex = 0.1 * sqrt(tot), pch = 19, xlab = "age", ylab = "seroprevalence", +#' xlim = c(0, 45), ylim = c(0, 1) +#' ) #' ) -#' #' @source Thiry, N., Beutels, P., Shkedy, Z. et al. The seroepidemiology #' of primary varicella-zoster virus infection in Flanders (Belgium). #' Eur J Pediatr 161, 588-593 (2002). @@ -496,14 +501,13 @@ #' } #' #' @examples -#' # Reproduce Fig 4.7 (right panel), p.71 #' library(dplyr) #' df <- vzv_be_2001_2003 %>% #' mutate(age = round(age)) %>% #' group_by(age) %>% #' summarise(pos = sum(seropositive), tot = n()) #' plot(df$age, df$pos / df$tot, -#' cex = 0.036 * df$tot, pch = 19, xlab = "age", ylab = "seroprevalence", +#' cex = 0.1 * sqrt(df$tot), pch = 19, xlab = "age", ylab = "seroprevalence", #' xlim = c(0, 45), ylim = c(0, 1) #' ) #' @@ -514,6 +518,9 @@ #' \doi{doi:10.1017/S0950268807009661} "vzv_be_2001_2003" + +# ===== Rubella - Mumps ======= + #' Rubella - Mumps data from the UK (aggregated) #' #' @@ -532,6 +539,7 @@ #' 1988; 297 :770 \doi{doi:10.1136/bmj.297.6651.770} "rubella_mumps_uk" +# ===== VZV - Parvo B19 ======= #' VZV and Parvovirus B19 serological data in Belgium (line listing) #' @format A data frame with 7 variables: diff --git a/R/farrington_model.R b/R/farrington_model.R new file mode 100644 index 0000000..a73d94d --- /dev/null +++ b/R/farrington_model.R @@ -0,0 +1,93 @@ +#' The Farrington (1990) model. +#' +#' @description Fit age-stratified seroprevalence data using the Farrington (1990) model, which assumes the +#' force of infection increases linearly with age and subsequently decreases exponentially. +#' +#' @details +#' The force of infection is defined as followed +#' +#' \deqn{ +#' \lambda(a) = (\alpha a - \gamma)e^{-\beta a} + \gamma +#' } +#' Where \eqn{\gamma} is called the long term residual for FOI, +#' as \eqn{a \rightarrow \infty} , \eqn{\lambda (a) \rightarrow \gamma} +#' +#' The seroprevalence can thus be estimated using the non-linear model +#' \deqn{ +#' \pi(a) = 1 - exp\{ \frac{\alpha}{\beta}ae^{-\beta a} + +#' \frac{1}{\beta}(\frac{\alpha}{\beta} - +#' \gamma)(e^{-\beta a} - 1) -\gamma a \} +#' } +#' +#' Refer to section 6.1.2. of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. +#' +#' +#' @param data the input data frame, must either have columns for `age`, `pos`, `tot` (for aggregated data) OR `age`, `status` (for linelisting data) +#' @param start Named list of vectors or single vector. +#' Initial values for optimizer. +#' @param fixed Named list of vectors or single vector. +#' Parameter values to keep fixed during optimization. +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). +#' +#' @return a list of class farrington_model with 5 items +#' \item{datatype}{type of datatype used for model fitting (aggregated or linelisting)} +#' \item{df}{the dataframe used for fitting the model} +#' \item{info}{fitted "mle" object} +#' \item{sp}{seroprevalence} +#' \item{foi}{force of infection} +#' @seealso [stats4::mle()] for more information on the fitted mle object +#' +#' @examples +#' df <- rubella_uk_1986_1987 +#' model <- farrington_model( +#' df, +#' start=list(alpha=0.07,beta=0.1,gamma=0.03) +#' ) +#' plot(model) +#' +#' @importFrom stats4 mle +#' +#' @export +farrington_model <- function(data, start, fixed=list(), + age_col="age",pos_col="pos", tot_col="tot", status_col="status") +{ + model <- list() + + # check input whether it is line-listing or aggregated data + data <- check_input(data, stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) + age <- data$age + pos <- data$pos + tot <- data$tot + model$datatype <- data$type + + farrington <- function(alpha,beta,gamma) { + p=1-exp((alpha/beta)*age*exp(-beta*age) + +(1/beta)*((alpha/beta)-gamma)*(exp(-beta*age)-1)-gamma*age) + ll=pos*log(p)+(tot-pos)*log(1-p) + return(-sum(ll)) + } + + model$info <- mle(farrington, fixed=fixed, start=start) + alpha <- model$info@coef[1] + beta <- model$info@coef[2] + gamma <- model$info@coef[3] + model$sp <- 1-exp( + (alpha/beta)*age*exp(-beta*age) + +(1/beta)*((alpha/beta)-gamma)*(exp(-beta*age)-1) + -gamma*age) + model$foi <- (alpha*age-gamma)*exp(-beta*age)+gamma + model$df <- list(age=age, pos=pos, tot=tot) + + class(model) <- "farrington_model" + model +} diff --git a/R/fractional_polynomial_models.R b/R/fractional_polynomial_models.R index 202fe76..759c758 100644 --- a/R/fractional_polynomial_models.R +++ b/R/fractional_polynomial_models.R @@ -20,14 +20,15 @@ formulate <- function(p) { equation } -#' Returns the powers of the GLM fitted model which has the lowest deviance score. +#' Returns the powers of the fractional polynomial model which has the lowest deviance score. #' -#' Refers to section 6.2. +#' Return the best powers for a given degree #' -#' @param data the input data frame, must either have `age`, `pos`, `tot` columns (for aggregated data) OR `age`, `status` for (linelisting data) -#' @param p a powers sequence. +#' @param data the input data frame, must either have columns for `age`, `pos`, `tot` (for aggregated data) OR +#' `age`, `status` (for linelisting data) +#' @param p a powers sequence to be tested. #' @param mc indicates if the returned model should be monotonic. -#' @param degree the degree of the model. Recommended to be <= 2. +#' @param degree the maximum degree (i.e. number of power terms) to search for the best model. Recommended to be <= 2. #' @param link the link function. Defaulted to "logit". #' #' @return list of 3 elements: @@ -35,88 +36,126 @@ formulate <- function(p) { #' \item{deviance}{Deviance of the best fitted model.} #' \item{model}{The best model fitted} #' -#' @examples -#' df <- hav_be_1993_1994 -#' best_p <- find_best_fp_powers( -#' df, -#' p=seq(-2,3,0.1), mc=FALSE, degree=2, link="cloglog" -#' ) -#' best_p -#' #' @importFrom stats glm binomial as.formula -#' -#' @export -find_best_fp_powers <- function(data, p, mc, degree, link="logit"){ - data <- check_input(data) +#' @import dplyr tidyr +#' @importFrom purrr pmap_dfr +find_best_fp_powers <- function(data, + p, mc, degree, link="logit"){ age <- data$age pos <- data$pos tot <- data$tot - glm_best <- NULL - d_best <- NULL - p_best <- NULL - #---- - min_p <- 1 - max_p <- length(p) - state <- rep(min_p, degree) - i <- degree - #---- - - get_cur_p <- function(cur_state) { - cur_p <- c() - for (i in 1:degree) { - cur_p <- c(cur_p, p[cur_state[i]]) - } - cur_p - } + best_mod <- NULL # best model + best_p <- NULL # best powers (p vector) for the given degree - repeat { - if ( - (i < degree && state[i] == max_p) - || (i == degree && state[i] == max_p+1) - ) { - if (i-1 == 0) break - if (state[i-1] < max_p) { - state[i-1] <- state[i-1]+1 - for (j in i:degree) state[j] <- state[i-1] - i <- degree - } else { - i <- i-1 - next - } - } - #------ iteration implementation ------- - p_cur <- get_cur_p(state) + # Starting from the lowest degree + # Get the best combinations of powers p + # Try increment m and only accept when there is a statistically significant improvement (determined by LRT) + for (curr_deg in 1:degree){ + # generate combinations of powers + p_combis <- expand.grid(rep(list(p), curr_deg)) + # filter the rows that are increasing in order + p_combis <- p_combis[apply(p_combis, 1, \(r){all(diff(r)>=0)}), ,drop=FALSE] + # arrange by increasing order of power in first->second degree and so on + p_combis <- p_combis[do.call(order,p_combis),,drop=FALSE] - glm_cur <- glm( - as.formula(formulate(p_cur)), - family=binomial(link=link) - ) - if (glm_cur$converged == TRUE) { - # d_cur <- deviance(glm_cur) - d_cur <- glm_cur$deviance - if (is.null(glm_best) || d_cur < d_best) { - if ((mc && is_monotone(glm_cur)) | !mc) { - glm_best <- glm_cur - d_best <- d_cur - p_best <- p_cur - } + # fit model with all the combinations of p and degree + mods <- p_combis %>% + pmap_dfr(\(...){ + curr_p <- as.numeric(c(...)) + + curr_mod <- glm( + as.formula(formulate(curr_p)), + family=binomial(link=link) + ) + + # only accept the parameters if the model converged + if(curr_mod$converged==FALSE) return(NULL) + + # only accept the parameters if the model is monotonic (if enforced) + if(mc && !is_monotone(curr_mod)) return(NULL) + + tibble( + p = list(curr_p), + mod = list(curr_mod), + deviance = curr_mod$deviance + ) + }) + + # handle scenario where no models converged at current degree + if(nrow(mods) == 0) next + # get the best model and powers for the current degree + best_idx <- which.min(mods$deviance) + curr_deg_mod <- mods[best_idx, ][["mod"]][[1]] + curr_deg_p <- mods[best_idx, ][["p"]][[1]] + + # check if the best model with current degree is better than the last degree + if(!is.null(best_mod)){ + # perform LRT + lrt_out <- anova(best_mod, curr_deg_mod, test="LRT") + lrt_out <- as.data.frame(lrt_out) + + # Royston&Altman (1994) suggests significance level of 0.1 + if(!is.na(lrt_out$`Pr(>Chi)`[2]) && lrt_out$`Pr(>Chi)`[2] < 0.1){ + best_mod <- curr_deg_mod + best_p <- curr_deg_p } + }else{ + best_mod <- curr_deg_mod + best_p <- curr_deg_p } - #--------------------------------------- - if (sum(state != max_p) == 0) break - state[i] <- state[i]+1 } - return(list(p=p_best, deviance=d_best, model=glm_best)) + + if(is.null(best_mod)) stop("Cannot find a converged model with the given degree and power") + + list(p=best_p, model=best_mod) } #' A fractional polynomial model. #' -#' Refers to section 6.2. +#' @description Fractional polynomial model is a generalization of polynomial models +#' where the power of the terms can be fractions, allowing more flexibility and better +#' fit for data where asymptotic behavior is expected. #' +#' @details +#' Instead of a polynomial, the linear predictor is now defined as +#' \deqn{ +#' \eta_m(a, \beta, p_1, p_2, ...,p_m) = \Sigma^m_{i=0} \beta_i H_i(a) +#' } +#' Where \eqn{m} is an integer, \eqn{p_1 \le p_2 \le... \le p_m} is a sequence of powers, +#' and \eqn{H_i(a)} is a transformation given by +#' +#' \deqn{ +#' H_i = \begin{cases} +#' a^{p_i} & \text{ if } p_i \neq p_{i-1}, +#' \\ H_{i-1}(a) \times log(a) & \text{ if } p_i = p_{i-1}, +#' \end{cases} +#' } +#' +#' Refers to section 6.2. of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. #' @param data the input data frame, must either have `age`, `pos`, `tot` columns (for aggregated data) OR `age`, `status` for (linelisting data) -#' @param p the powers of the predictor. +#' @param p is either: +#' \itemize{ +#' \item{a numeric vector specifying the powers to apply to the predictors} +#' \item{ +#' a named list with two elements, \code{"p_range"} and \code{"degree"}. \code{"p_range"} +#' is a sequence of powers and \code{"degree"} is the maximum degree. +#' In which case the package will search for the best degree and power combinations +#' } +#' } #' @param link the link function for model. Defaulted to "logit". +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). +#' @param monotonic whether the returned model should be monotonic (if a search is specified) #' #' @importFrom stats predict as.formula #' @@ -126,7 +165,10 @@ find_best_fp_powers <- function(data, p, mc, degree, link="logit"){ #' \item{info}{a fitted glm model} #' \item{sp}{seroprevalence} #' \item{foi}{force of infection} -#' @seealso [stats::glm()] for more information on glm object +#' @seealso +#' [stats::glm()] for more information on glm object +#' +#' [polynomial_models()] #' #' @examples #' df <- hav_be_1993_1994 @@ -136,20 +178,36 @@ find_best_fp_powers <- function(data, p, mc, degree, link="logit"){ #' plot(model) #' #' @export -fp_model <- function(data,p, link="logit") { +fp_model <- function(data,p,monotonic=FALSE,link="logit", + age_col="age",pos_col="pos", tot_col="tot", status_col="status") { model <- list() - data <- check_input(data) + data <- check_input(data, stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) age <- data$age pos <- data$pos tot <- data$tot - model$datatype <- data$type - model$info <- glm( - as.formula(formulate(p)), - family=binomial(link=link) - ) + # handle powers input here + if(is.numeric(p)){ + model$info <- glm( + as.formula(formulate(p)), + family=binomial(link=link) + ) + model$p <- p + }else if(is.list(p) && all(c("p_range", "degree") %in% names(p))){ + out <- find_best_fp_powers( + data = data.frame(age=age, pos=pos, tot=tot), + p = p$p_range, degree = p$degree, mc = monotonic, link = link + ) + model$p <- out$p + model$info <- out$model + }else{ + stop("Invalid value for `p`: either a numeric vector or a named list with + 2 elements `p_range` and `degree`") + } + + model$sp <- model$info$fitted.values model$foi <- est_foi( t=age, diff --git a/R/hierarchical_bayesian_model.R b/R/hierarchical_bayesian_model.R index c62d032..de6b02c 100644 --- a/R/hierarchical_bayesian_model.R +++ b/R/hierarchical_bayesian_model.R @@ -1,12 +1,69 @@ #' Hierarchical Bayesian Model #' -#' Refers to section 10.3 #' -#' @param data the input data frame, must either have `age`, `pos`, `tot` columns (for aggregated data) OR `age`, `status` for (linelisting data) +#' @description Fit age-stratified seroprevalence to parametric hierarchical Bayesian models. +#' Supported models including Farrington model (2 and 3 parameters variants) +#' and Log Logistic model +#' +#' @details Consider a model for prevalence that has a parametric form +#' \eqn{\pi(a_i, \alpha)} where \eqn{\alpha} is a parameter vector +#' +#' Under a Bayesian framework, we can constraint the parameter space of the prior distribution \eqn{P(\alpha)} +#' to achieve monotonicity of the posterior distribution \eqn{P(\pi_1, \pi_2, ..., \pi_m|y,n)} +#' +#' Where: +#' +#' - \eqn{n = (n_1, n_2, ..., n_m)} and \eqn{n_i} is the sample size at age \eqn{a_i} +#' +#' - \eqn{y = (y_1, y_2, ..., y_m)} and \eqn{y_i} is the number of infected individual from the \eqn{n_i} sampled subjects +#' +#' For \bold{Farrington} model with 3 parameters, prevalence is formulated as follow +#' +#' \deqn{ +#' \pi (a) = 1 - exp\{ \frac{\alpha_1}{\alpha_2}ae^{-\alpha_2 a} + +#' \frac{1}{\alpha_2}(\frac{\alpha_1}{\alpha_2} - \alpha_3)(e^{-\alpha_2 a} - 1) -\alpha_3 a \} +#' } +#' +#' The likelihood model is defined as \eqn{y_i \sim Bin(n_i, \pi_i), \text{ for } i = 1,2,3,...m} +#' +#' The constraint on the parameter space can be incorporated by assuming +#' truncated normal distribution for the components of \eqn{\alpha}, +#' \eqn{\alpha = (\alpha_1, \alpha_2, \alpha_3)} in \eqn{\pi_i = \pi(a_i,\alpha)} +#' +#' The flat hyperpriors are defined as \eqn{\mu_j \sim \mathcal{N}(0, 10000)} and +#' \eqn{\tau^{-2}_j \sim \Gamma(100,100)} +#' +#' For \bold{Farrington} model with 2 parameters, it is equivalent to the previous model with \eqn{\alpha_3 = 0} +#' +#' For \bold{Log logistic model}, seroprevalence is instead defined as +#' +#' \deqn{\pi(a) = \frac{\beta a^\alpha}{1 + \beta a^\alpha}, \text{ } \alpha, \beta > 0} +#' +#' The likelihood is similarly defined as \eqn{y_i \sim Bin(n_i, \pi_i))} +#' +#' The prior model of \eqn{\alpha_1} is specified as \eqn{\alpha_1 \sim \text{truncated } \mathcal{N}(\mu_1, \tau_1)} +#' with flat hyperpriors as in Farrington model +#' +#' \eqn{\beta} is constrained to be positive by specifying \eqn{\alpha_2 \sim \mathcal{N}(\mu_2, \tau_2)} +#' +#' Refer to section Chapter 10.3 of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. +#' +#' @param data the input data frame, must either have columns for `age`, `pos`, `tot` (for aggregated data) OR `age`, `status` (for linelisting data) #' @param type type of model ("far2", "far3" or "log_logistic") #' @param chains number of Markov chains #' @param warmup number of warmup runs #' @param iter number of iterations +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). #' #' @importFrom rstan sampling summary #' @importFrom boot inv.logit @@ -30,11 +87,12 @@ #' plot(model) #' } hierarchical_bayesian_model <- function(data, + age_col="age",pos_col="pos", tot_col="tot", status_col="status", type="far3",chains = 1,warmup = 1500,iter = 5000){ model <- list() # check input whether it is line-listing or aggregated data - data <- check_input(data) + data <- check_input(data, stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) model$datatype <- data$type age <- data$age pos <- data$pos diff --git a/R/mixture_model.R b/R/mixture_model.R index ce3f7ee..1604c9d 100644 --- a/R/mixture_model.R +++ b/R/mixture_model.R @@ -1,12 +1,47 @@ #' Fit a mixture model to classify serostatus #' -#' Refers to section 11.1 - 11.4 +#' @description Fit the antibody level data to a 2-component Gaussian mixture model #' -#' @param antibody_level - vector of the corresponding raw antibody level -#' @param breaks - number of intervals which the antibody_level are grouped into -#' @param pi - proportion of susceptible, infected -#' @param mu - a vector of means of component distributions (vector of 2 numbers in ascending order) -#' @param sigma - a vector of standard deviations of component distributions (vector of 2 number) +#' @details +#' Antibody level (denoted \eqn{Z}) is modeled using a 2-component Gaussian +#' mixture model. Each component \eqn{Z_j} (\eqn{j \in \{I, S\}}) represents the +#' antibody level of the latent Infected and Susceptible sub-populations, following density +#' \eqn{f_j(z_j|\theta_j)} +#' +#' Let \eqn{\pi_{\text{TRUE}}(a)} denotes the age-dependent mixing probability +#' (i.e., the true prevalence), the density of the mixture is formulated as +#' +#' \deqn{f(z|z_I, z_S,a) = (1-\pi_{\text{TRUE}}(a))f_S(z_S|\theta_S)+\pi_{\text{TRUE}}(a)f_I(z_I|\theta_I)} +#' +#' The mean \eqn{E(Z|a)} thus equals +#' \deqn{\mu(a) = (1-\pi_{\text{TRUE}}(a))\mu_S+\pi_{\text{TRUE}}(a)\mu_I} +#' +#' From which true prevalence can be computed as +#' \deqn{\pi_{\text{TRUE}}(a) = \frac{\mu(a) - \mu_S}{\mu_I - \mu_S}} +#' +#' And FOI can then be inferred as +#' \deqn{\lambda_{TRUE} = \frac{\mu'(a)}{\mu_I - \mu(a)}} +#' +#' Function [serosv::mixture_model()] fits antibody level data to \eqn{f_S(z_S|\theta_S)} and +#' \eqn{f_I(z_I|\theta_I)} +#' +#' Function [serosv::estimate_mixture()] will then estimate age-specific antibody level \eqn{\mu(a)} +#' and infer the estimation for \eqn{\pi_{\text{TRUE}}(a)} and \eqn{\lambda_{TRUE}} +#' +#' Refer to section 11.3. of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. +#' +#' @param antibody_level vector of the corresponding raw antibody level +#' @param breaks number of intervals which the antibody_level are grouped into +#' @param pi proportion of susceptible, infected +#' @param mu a vector of means of component distributions (vector of 2 numbers in ascending order) +#' @param sigma a vector of standard deviations of component distributions (vector of 2 number) #' #' @importFrom mixdist mix mixgroup mixparam #' @importFrom stats fitted @@ -41,17 +76,53 @@ mixture_model <- function (antibody_level, breaks=40, pi=c(0.2, 0.8), mu=c(2,6), model } -#' Estimate seroprevalence and foi by combining mixture model and regression +#' Estimate seroprevalence and FOI from a fixed mixture model +#' +#' @description Estimate age-specific seroprevalence and FOI given +#' a fitted mixture model (generated by [serosv::mixture_model()]) +#' +#' @details +#' Antibody level (denoted \eqn{Z}) is modeled using a 2-component Gaussian +#' mixture model. Each component \eqn{Z_j} (\eqn{j \in \{I, S\}}) represents the +#' antibody level of the latent Infected and Susceptible sub-populations, following density +#' \eqn{f_j(z_j|\theta_j)} +#' +#' Let \eqn{\pi_{\text{TRUE}}(a)} denotes the age-dependent mixing probability +#' (i.e., the true prevalence), the density of the mixture is formulated as +#' +#' \deqn{f(z|z_I, z_S,a) = (1-\pi_{\text{TRUE}}(a))f_S(z_S|\theta_S)+\pi_{\text{TRUE}}(a)f_I(z_I|\theta_I)} +#' +#' The mean \eqn{E(Z|a)} thus equals +#' \deqn{\mu(a) = (1-\pi_{\text{TRUE}}(a))\mu_S+\pi_{\text{TRUE}}(a)\mu_I} +#' +#' From which true prevalence can be computed as +#' \deqn{\pi_{\text{TRUE}}(a) = \frac{\mu(a) - \mu_S}{\mu_I - \mu_S}} +#' +#' And FOI can then be inferred as +#' \deqn{\lambda_{TRUE} = \frac{\mu'(a)}{\mu_I - \mu(a)}} +#' +#' Function [serosv::mixture_model()] fits antibody level data to \eqn{f_S(z_S|\theta_S)} and +#' \eqn{f_I(z_I|\theta_I)} +#' +#' Function [serosv::estimate_mixture()] will then estimate age-specific antibody level \eqn{\mu(a)} +#' and infer the estimation for \eqn{\pi_{\text{TRUE}}(a)} and \eqn{\lambda_{TRUE}} +#' +#' Refer to section 11.3. of the the book by Hens et al. (2012) for further details. #' -#' Refers to section 11.2 - 11.4 +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. #' -#' @param age - vector of age -#' @param antibody_level - vector of the corresponding raw antibody level -#' @param mixture_model - mixture_model object generated by serosv::mixture_model() -#' @param s - smoothing basis used to fit antibody level -#' @param sp - smoothing parameter -#' @param threshold_status - sero status using threshold approach in line listing (optional, for visualization and comparison only) -#' @param monotonize - whether to monotonize seroprevalence (default to TRUE) +#' @param age vector of age +#' @param antibody_level vector of the corresponding raw antibody level +#' @param mixture_model mixture_model object generated by serosv::mixture_model() +#' @param s smoothing basis used to fit antibody level +#' @param sp smoothing parameter +#' @param threshold_status sero status using threshold approach in line listing (optional, for visualization and comparison only) +#' @param monotonize whether to monotonize seroprevalence (default to TRUE) #' #' @importFrom mgcv gam #' @importFrom stats approx gaussian @@ -124,6 +195,7 @@ estimate_from_mixture <- function(age, antibody_level, threshold_status = NULL, # save fitted df model$df <- data.frame(age = age, antibody_level = antibody_level) + model$monotonize <- monotonize if (!is.null(threshold_status)){ model$df[["threshold_status"]] <- threshold_status diff --git a/R/mseir_model.R b/R/mseir_model.R deleted file mode 100644 index 0fcfa6c..0000000 --- a/R/mseir_model.R +++ /dev/null @@ -1,58 +0,0 @@ -#' MSEIR model -#' -#' Refers to section 3.4. -#' -#' @param a age sequence -#' @param gamma time in maternal class. -#' @param lambda time in susceptible class. -#' @param sigma time in latent class. -#' @param nu time in infected class. -#' -#' @examples -#' model <- mseir_model( -#' a=seq(from=1,to=20,length=500), # age range from 0 -> 20 yo -#' gamma=1/0.5, # 6 months in the maternal antibodies -#' lambda=0.2, # 5 years in the susceptible class -#' sigma=26.07, # 14 days in the latent class -#' nu=36.5 # 10 days in the infected class -#' ) -#' model -#' -#' @return list of class mseir_model with the following parameters -#' \item{parameters}{list of parameters used for fitting the model} -#' \item{output}{matrix of proportion for each compartment over time} -#' -#' @export -mseir_model <- function(a, gamma, lambda, sigma, nu) -{ - - - ma <- exp(-gamma*a) - sa <- (gamma/(gamma-lambda))*(exp(-lambda*a)-exp(-gamma*a)) - ea <- ((lambda*gamma)/(gamma-lambda))* - ( - ((exp(-sigma*a)-exp(-lambda*a))/(lambda-sigma)) - -((exp(-sigma*a)-exp(-gamma*a))/(gamma-sigma)) - ) - ia <- (sigma*lambda*gamma)* - ( - ((exp(-nu*a)-exp(-sigma*a))/((lambda-sigma)*(gamma-sigma)*(sigma-nu))) - +((exp(-nu*a)-exp(-lambda*a))/((lambda-gamma)*(lambda-sigma)*(lambda-nu))) - +((exp(-nu*a)-exp(-gamma*a))/((gamma-lambda)*(gamma-sigma)*(gamma-nu))) - ) - - - model <- list() - model$parameters <- list(gamma = gamma, lambda = lambda, sigma = sigma, nu = nu) - model$output <- data.frame( - a = c(0, a), - m = c(1, ma), - s = c(0, sa), - e = c(0, ea), - i = c(0, ia), - r = c(0, 1 - ma - sa - ea - ia) - ) - - class(model) <- "mseir_model" - model -} diff --git a/R/nonparametric.R b/R/nonparametric.R index 6f1ccf1..24fe1d6 100644 --- a/R/nonparametric.R +++ b/R/nonparametric.R @@ -1,16 +1,63 @@ #' A local polynomial model. #' -#' Refers to section 7.1. and 7.2. +#' @description Fit the age-specific seroprevalence to a local polynomial model, +#' where the linear predictor is approximated locally at one particular age. #' -#' @param data the input data frame, must either have `age`, `pos`, `tot` columns (for aggregated data) OR `age`, `status` for (linelisting data) +#' @details +#' Consider a linear predictor \eqn{\eta(a)} approximated locally at one particular value \eqn{a_0}. +#' +#' For a general degree \eqn{p}, the linear predictor for a neighbor of \eqn{a_0}, labeled \eqn{a_i} is equivalent to the Taylor approximation +#' \deqn{ +#' \eta(a_i) = \eta(a_0) + \eta^{(1)}(a_0)(a_i - a_0) + +#' \frac{\eta^{(2)}(a_0)}{2}(a_i - a_0)^2 + ... + \frac{\eta^{(p)}(a_0)}{p!}(a_i - a_0)^p +#' } +#' +#' \eqn{\eta(a_i)} can be estimated by maximizing +#' \deqn{ +#' \Sigma_{i=1}^{N} \ell_i \{Y_i, g^{-1} (\beta_0 + \beta_1(a_i-a_0)+ \beta_2(a_i-a_0)^2 ... + +#' \beta_p(a_i-a_0)^p) \} K_h(a_i - a_0) +#' } +#' +#' The estimator for the \eqn{k}-th derivative of \eqn{\eta(a_0)}, for \eqn{k = 0,1,…,p} +#' (degree of local polynomial) is thus: +#' \deqn{ +#' \hat{\eta}^{(k)}(a_0) = k!\hat{\beta}_k(a_0) +#' } +#' +#' The estimator for the prevalence at age \eqn{a_0} is then given by +#' \deqn{ +#' \hat{\pi}(a_0) = g^{-1}\{ \hat{\beta}_0(a_0) \} +#' } +#' Where \eqn{g} is the link function +#' +#' The estimator for the force of infection at age \eqn{a_0} by assuming \eqn{p \ge 1} is as followed +#' \deqn{ +#' \hat{\lambda}(a_0) = \hat{\beta}_1(a_0) \delta \{ \hat{\beta}_0 (a_0) \} +#' } +#' Where \eqn{\delta \{ \hat{\beta}_0(a_0) \} = \frac{dg^{-1} \{ \hat{\beta}_0(a_0) \} } {d\hat{\beta}_0(a_0)}} +#' +#' Refer to section 7.1 and 7.2. of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. +#' +#' @param data the input data frame, must either have columns for `age`, `pos`, `tot` (for aggregated data) OR `age`, `status` (for linelisting data) #' @param kern Weight function, default = "tcub". #' Other choices are "rect", "trwt", "tria", "epan", "bisq" and "gauss". #' Choices may be restricted when derivatives are required; #' e.g. for confidence bands and some bandwidth selectors. #' @param nn Nearest neighbor component of the smoothing parameter. #' Default value is 0.7, unless either h is provided, in which case the default is 0. -#' @param h The constant component of the smoothing parameter. Default: 0. +#' @param h The constant bandwidth of the smoothing parameter. Default: 0. #' @param deg Degree of polynomial to use. Default: 2. +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). #' #' @examples #' df <- mumps_uk_1986_1987 @@ -34,28 +81,106 @@ #' @seealso [locfit::locfit()] for more information on the fitted locfit object #' #' @export -lp_model <- function(data, kern="tcub", nn=0, h=0, deg=2) { - if (missing(nn) & missing(h)) { - nn <- 0.7 # default nn from lp() +lp_model <- function(data, kern="tcub", nn=0, h=0, deg=2, + age_col="age",pos_col="pos", tot_col="tot", status_col="status") { + if (all(nn==0) & all(h==0)) { + # default nn from lp() + nn <- 0.7 + h <- 0 } - model <- list() # check input whether it is line-listing or aggregated data - data <- check_input(data) + data <- check_input(data, stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) age <- data$age pos <- data$pos tot <- data$tot model$datatype <- data$type y <- pos/tot - estimator <- lp(age, deg=deg, nn=nn, h=h) - model$pi <- locfit(y~estimator, family="binomial", kern=kern) - model$eta <- locfit(y~estimator, family="binomial", kern=kern, deriv=1) - model$sp <- fitted(model$pi) - model$foi <- fitted(model$eta)*fitted(model$pi) # λ(a)=η′(a)π(a) + + if(length(nn) > 1 & length(h) > 1) stop("Tuning both `h` and `nn` may lead to parameter idenfifiability issues, please fix one of the parameters instead") + + if(length(nn) > 1 || length(h)>1){ + best_param <- best_lp_params( + data=data, + nn = nn, + h = h, + family="binomial", + kern=kern + ) + + nn <- best_param$nn + h <- best_param$h + } + + # print(paste0("nn: ", nn)) + # print(paste0("h: ", h)) + + model$info <- locfit(y~lp(age, deg=deg, nn=nn, h=h), family="binomial", kern=kern) + model$nn <- nn + model$h <- h + model$deg <- deg + model$kern <- kern + model$eta <- locfit(y~lp(age, deg=deg, nn=nn, h=h), family="binomial", kern=kern, deriv=1) + model$sp <- fitted(model$info) + model$foi <- fitted(model$eta)*fitted(model$info) # λ(a)=η′(a)π(a) model$df <- list(age=age, pos=pos, tot=tot) class(model) <- "lp_model" model } + +# function to return the best parameter of local polynomial model using GCV +# nn - range of values for nearest neighbor +# h - range of values for constant bandwidth +# if both nn and h are given, select either best nn or h, whichever gives the lowest GCV +#' @import tidyr dplyr locfit +best_lp_params <- function(data, nn=0, h=0, kern="tcub",deg=2, family="binomial"){ + # helper function to get df and GCV + summary.gcvplot <- function(object, ...){ + z <- cbind(object$df, object$values) + dimnames(z) <- list(NULL, c("df", object$cri)) + z + } + + # helper function which return parameter value which gives the lowest gcv + get_best_gcv <- function(nn_vals=0, h_vals=0){ + is_nn <- length(nn_vals) > 1 # check if we are tuning for nn or h + par_vals <- if (is_nn) nn_vals else h_vals + # generate parameters matrix + alpha <- + if (is_nn) + cbind(nn_vals, rep(h_vals, length(nn_vals))) + else + cbind(rep(nn_vals, length(h_vals)), h_vals) + + # compute gcv + gcv_out <- gcvplot( + pos / tot ~ age, + deg = deg, + kern = kern, + family = family, + alpha = alpha, + data=data + ) + + gcv_out <- cbind(par_vals, gcv_out$values) + best_idx <- which.min(gcv_out[,2]) + + c( + gcv_out[best_idx, 1], + gcv_out[best_idx, 2] + ) + } + + res <- list() + + if(length(nn) > 1 & length(h) > 1) stop("Tuning both `h` and `nn` may lead to parameter idenfifiability issues, please fix one of the parameters instead") + if(length(nn) > 1) res[c("h", "nn", "nn_gcv")] <- c(h, get_best_gcv(nn_vals=nn, h_vals=h)) + if(length(h) > 1) res[c("nn", "h", "h_gcv")] <- c(nn, get_best_gcv(nn_vals=nn, h_vals=h)) + + res +} + + diff --git a/R/plots.R b/R/plots.R index 78b4f89..b287811 100644 --- a/R/plots.R +++ b/R/plots.R @@ -28,8 +28,8 @@ set_plot_style <- function(sero = "blueviolet", ci = "royalblue1", foi = "#fc032 plot_data <- function(x){ if(x$datatype == "linelisting"){ # transform data before plotting - df_ <- transform_data(x$df$age, x$df$pos) - age <- df_$t + df_ <- transform_data(x$df, stratum_col="age", status_col="pos") + age <- df_$age pos <- df_$pos tot <- df_$tot # use pre-aggregated age for FOI @@ -45,7 +45,7 @@ plot_data <- function(x){ } #=== Helper function for plotting ===== -plot_util <- function(age, pos, tot, sero, foi, cex = 20){ +plot_util <- function(age, pos, tot, sero, foi, scale_foi=1, cex = 20){ # resolve no visible binding x <- y <- ymin <- ymax <- NULL @@ -55,7 +55,7 @@ plot_util <- function(age, pos, tot, sero, foi, cex = 20){ coord_cartesian(xlim=c(0,max(age)), ylim=c(0, 1)) + scale_y_continuous( name = "Seroprevalence", - sec.axis = sec_axis(~.*1, name = " Force of infection") + sec.axis = sec_axis(~.*scale_foi, name = " Force of infection") ) + set_plot_style() # === Add seroprevalence layer @@ -81,22 +81,51 @@ plot_util <- function(age, pos, tot, sero, foi, cex = 20){ # --- Handle cases where FOI is a data.frame (with CI) # plot <- plot + geom_smooth(aes_auto(foi, col = "foi", linetype="foi", fill="ci"), data=foi, # stat="identity",lwd=0.5) - plot <- plot + geom_smooth(aes(x = x, y=y, ymin =ymin, ymax = ymax, col = "foi", linetype="foi", fill="ci"), data=foi, - stat="identity",lwd=0.5) + plot <- plot + geom_smooth( + aes( + x = x, + # scale FOI if specified + y = y/scale_foi, + ymin = ymin/scale_foi, + ymax = ymax/scale_foi, + col = "foi", + linetype = "foi", + fill = "ci" + ), + data = foi, + stat = "identity", + lwd = 0.5 + ) }else{ # --- Handle cases where CI for FOI is not computable & length of age for foi differs from provided age vector - plot <- plot + geom_line(aes(x = x, y=y, col = "foi", linetype="foi"), data=foi, - stat="identity",lwd=0.5) + plot <- plot + geom_line( + aes( + x = x, + # scale FOI if specified + y = y/scale_foi, + col = "foi", + linetype = "foi" + ), + data = foi, + stat = "identity", + lwd = 0.5 + ) } }else if (length(age) != length(foi)){ # --- handle some cases when length of age differs from length of foi age <- age[c(-1,-length(age))] foi <- data.frame(x = age, y = foi) - plot <- plot + geom_line(aes(x = x, y=y, col = "foi", linetype="foi"), data = foi, - lwd = 0.5) + plot <- plot + geom_line(aes( + x = x, + y = y/scale_foi, + col = "foi", + linetype = "foi" + ), + data = foi, + lwd = 0.5) }else{ # --- Simply plot foi - plot <- plot + geom_line(aes(x = age, y = foi, col = "foi", linetype="foi"), + plot <- plot + geom_line(aes(x = age, y = foi/scale_foi, col = "foi", linetype="foi"), lwd = 0.5) } plot @@ -107,129 +136,6 @@ plot_util <- function(age, pos, tot, sero, foi, cex = 20){ # UseMethod("plot") # } -#### SIR model #### - -#' plot() overloading for SIR model -#' -#' @param x the sir_basic_model object. -#' @param ... arbitrary params. -#' @import ggplot2 -#' @importFrom methods is -#' @importFrom graphics plot -#' -#' @return ggplot object -#' @export -plot.sir_basic_model <- function(x, ...){ - comp_lvl <- c("S", "I", "R") - time <- S <- I <- R <- NULL - - ggplot(x$output) + - geom_line(aes(x = time, y = S, color = factor("S", levels = comp_lvl))) + - geom_line(aes(x = time, y = I, color = factor("I", levels = comp_lvl))) + - geom_line(aes(x = time, y = R, color = factor("R", levels = comp_lvl))) + - list( - scale_colour_manual( - values = c("S" = "blueviolet", "I" = "#fc0328", "R" = "royalblue1"), - labels = c("S"="susceptible", "I"="infected", "R"="recovered") - ), - labs( x = "Time", - y = "Count", - colour = "Compartment") - ) -} - -#' plot() overloading for SIR static model -#' -#' @param x the sir_static_model object. -#' @param ... arbitrary params. -#' @import ggplot2 -#' @importFrom methods is -#' @importFrom graphics plot -#' -#' @return ggplot object -#' @export -plot.sir_static_model <- function(x, ...){ - comp_lvl <- c("s", "i", "r") - time <- s <- i <- r <- NULL - - ggplot(x$output) + - geom_line(aes(x = time, y = s, color = factor("s", levels = comp_lvl))) + - geom_line(aes(x = time, y = i, color = factor("i", levels = comp_lvl))) + - geom_line(aes(x = time, y = r, color = factor("r", levels = comp_lvl))) + - list( - scale_colour_manual( - values = c("s" = "blueviolet", "i" = "#fc0328", "r" = "royalblue1"), - labels = c("s"="susceptible", "i"="infected", "r"="recovered") - ), - labs( - colour = "Compartment", - x = "Age", - y = "Fraction") - ) -} - - -#' plot() overloading for SIR sub populations model -#' -#' @param x the sir_subpops_models object. -#' @param ... arbitrary params. -#' @import ggplot2 -#' @importFrom methods is -#' @importFrom graphics plot -#' -#' @return list of ggplot objects, each object is the plot for the corresponding subpopulation -#' @export -plot.sir_subpops_model <- function(x, ...){ - time <- s <- i <- r <- NULL - comp_lvl <- c("s", "i", "r") - - # using for loop here would not work due to ggplot lazy eval - subpop_plots <- lapply(1:x$parameters$k, function(subpop) { - ggplot(x$output) + - geom_line(aes(x = time, y = get(paste0("s", subpop)), color = factor("s", levels = comp_lvl))) + - geom_line(aes(x = time, y = get(paste0("i", subpop)), color = factor("i", levels = comp_lvl))) + - geom_line(aes(x = time, y = get(paste0("r", subpop)), color = factor("r", levels = comp_lvl))) + - scale_colour_manual( - values = c("s" = "blueviolet", "i" = "#fc0328", "r" = "royalblue1"), - labels = c("s"="susceptible", "i"="infected", "r"="recovered") - ) + - labs(title= paste0("Plot for subpopulation ", subpop), - x = "Time", - y = "Fraction", - colour = "Compartment") - }) - - names(subpop_plots) <- paste0("subpop_", 1:x$parameters$k) - subpop_plots -} - -#' plot() overloading for MSEIR model -#' -#' @param x the mseir_model object. -#' @param ... arbitrary params. -#' @import ggplot2 -#' @importFrom methods is -#' @importFrom graphics plot -#' -#' @return ggplot object -#' @export -plot.mseir_model <- function(x, ...){ - a <- m <- s <- e <- i <- r <- NULL - # make leveled factor to force legend show color in order - comp_lvl <- c("m", "s", "e", "i", "r") - - ggplot(x$output) + - geom_line(aes(x = a, y = m, color = factor("m", levels = comp_lvl))) + - geom_line(aes(x = a, y = s, color = factor("s", levels = comp_lvl))) + - geom_line(aes(x = a, y = e, color = factor("e", levels = comp_lvl))) + - geom_line(aes(x = a, y = i, color = factor("i", levels = comp_lvl))) + - geom_line(aes(x = a, y = r, color = factor("r", levels = comp_lvl))) + - scale_colour_manual( - values = c("m"="#3ea379","s" = "blueviolet", "e"="#3e45a3", "i" = "#fc0328", "r" = "royalblue1"), - labels = c("m"="maternal immunity", "s"="susceptible", "e"="exposed", "i"="infected", "r"="recovered") - ) + - labs(x = "Age", y = "Fraction", color = "Compartment") -} #### Polynomial model #### #' plot() overloading for polynomial model @@ -511,12 +417,17 @@ plot.estimate_from_mixture <- function(x, ... ){ returned_plot <- ggplot() if(!is.null(x$df$threshold_status)){ - aggregated <- transform_data(round(x$df$age), x$df$threshold_status) + aggregated <- transform_data( + data.frame( + age = round(x$df$age), + status = x$df$threshold_status + ) + ) # resolve no visible binding note - t <- pos <- tot <- NULL + age <- pos <- tot <- NULL returned_plot <- returned_plot + - geom_point(aes( x = t, y = pos/tot, size = cex*(pos)/max(tot) ), data = aggregated, + geom_point(aes( x = age, y = pos/tot, size = cex*(pos)/max(tot) ), data = aggregated, shape = 1, show.legend = FALSE) } @@ -529,6 +440,7 @@ plot.estimate_from_mixture <- function(x, ... ){ returned_plot + set_plot_style() + labs(x = "Age", y="Seroprevalence") } +# ------- Plot age time varying seroprevalence ---------- #' Plot output for age_time_model #' #' @param x - a `age_time_model` object @@ -665,12 +577,13 @@ plot_gcv <- function(age, pos, tot, nn_seq, h_seq, kern="tcub", deg=2) { nn_plot + h_plot + plot_layout(ncol=2) } +# ----- Plot corrected prevalence ------- #' Plot output for corrected_prevalence #' #' @param x - the output of `correct_prevalence()` function #' @param y - another output of `correct_prevalence()` function (optional, for comparison only) #' @param facet - whether to plot as facets or on the same plot (only when y is provided) -#' @import ggplot2 tidyr patchwork +#' @import ggplot2 tidyr patchwork #' @importFrom magrittr %>% #' @importFrom assertthat assert_that #' @@ -712,7 +625,7 @@ plot_corrected_prev <- function(x, y=NULL, facet=FALSE){ ymax = sero_upr, color = label ), - alpha = 0.7, data = data + alpha = 0.5, data = data ) } # them add label layer @@ -739,9 +652,15 @@ plot_corrected_prev <- function(x, y=NULL, facet=FALSE){ plot <- ggplot() + geom_point(aes( - x = age, y = pos / tot, + x = age, y = sero, color = "apparent prevalence" - ), data = dat) + ), alpha = 0.8, data = dat)+ + geom_errorbar( + aes( + x = age, y = sero, ymin = sero_lwr, ymax = sero_upr, + color = "apparent prevalence" + ), alpha = 0.5,data = dat + ) if(!facet){ plot <- plot + generate_layers(corrected_dat, title = if(is.null(y)) paste0("Plot for ", x$method, " approach") else NULL) diff --git a/R/polynomial_models.R b/R/polynomial_models.R index c6113ef..6095f07 100644 --- a/R/polynomial_models.R +++ b/R/polynomial_models.R @@ -10,21 +10,67 @@ X <- function(t, degree) { #' Polynomial models #' -#' Refers to section 6.1.1 -#' @param data the input data frame, must either have `age`, `pos`, `tot` columns (for aggregated data) OR `age`, `status` for (linelisting data) -#' @param k degree of the model. -#' @param type name of method (Muench, Giffith, Grenfell). -#' @param link link function. +#' @description Fit age-stratified seroprevalence data to serocatalytic models formulated as polynomials. +#' +#' @details +#' The seroprevalence is assumed to follow the general format +#' \deqn{ +#' \pi(a) = 1 - e^{-\Sigma_{i=1}^k \beta_i a^i} +#' } +#' Which implies the force of infection to be \eqn{\lambda(a) = \Sigma_{i=1}^k \beta_i i a^{i-1}} +#' +#' Where: +#' +#' - \eqn{\pi} is the seroprevalence at age \eqn{a} +#' +#' - \eqn{a} is the variable age +#' +#' - \eqn{k} is the degree of the polynomial +#' +#' The seroprevalence \eqn{\pi(a)} is fitted using a GLM with log link with +#' the linear predictor \eqn{\eta(a) = \Sigma_{i=1}^k \beta_i a^{i}} +#' +#' Muench (1934) model is equivalent to a degree 1 (\eqn{k=1}) linear predictor +#' +#' Griffith model is equivalent to a degree 2 (\eqn{k=2}) linear predictor +#' +#' Grenfell & Anderson (1985) suggested a higher order polynomials (\eqn{k \geq 3}) +#' +#' Refer to section 6.1.1. of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. +#' +#' Grenfell, B. T., and R. M. Anderson. 1985. “The Estimation of +#' Age-Related Rates of Infection from Case Notifications and Serological Data.” +#' The Journal of Hygiene 95 (2): 419–36. \doi{https://doi.org/10.1017/s0022172400062859}. +#' +#' Muench, Hugo. 1934. “Derivation of Rates from Summation Data by the Catalytic Curve.” +#' Journal of the American Statistical Association 29 (185): +#' 25–38. \doi{https://doi.org/10.1080/01621459.1934.10502684}. +#' +#' @param data the input data frame, must either have columns for `age`, `pos`, `tot` (for aggregated data) OR `age`, `status` (for linelisting data) +#' @param k degree of the polynomial. (k=1 for Muench model, k=2 for Griffith model, k=3 for Grenfell model). +#' @param link link function (default link="log"). +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). #' #' @examples #' data <- parvob19_fi_1997_1998[order(parvob19_fi_1997_1998$age), ] -#' data$status <- data$seropositive -#' aggregated <- transform_data(data$age, data$seropositive, stratum_col = "age") +#' aggregated <- transform_data(data, stratum_col = "age", status_col="seropositive") #' #' # fit with aggregated data -#' model <- polynomial_model(aggregated, type = "Muench") +#' model <- polynomial_model(aggregated, k = 1) #' # fit with linelisting data -#' model <- polynomial_model(data, type = "Muench") +#' model <- polynomial_model(data, +#' status_col = "seropositive", +#' k = 1) #' plot(model) #' #' @return a list of class polynomial_model with 5 items @@ -35,97 +81,102 @@ X <- function(t, degree) { #' \item{foi}{force of infection} #' #' @export -polynomial_model <- function(data, k,type, link = "log"){ +polynomial_model <- function(data, k, link = "log", + age_col="age",pos_col="pos", tot_col="tot", status_col="status"){ model <- list() - data <- check_input(data) + data <- check_input(data, stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) model$datatype <- data$type - Age <- data$age - Pos <- data$pos - Neg <- data$tot - Pos + age <- data$age + pos <- data$pos + neg <- data$tot - pos + + df <- data.frame(cbind(age, pos,neg)) - df <- data.frame(cbind(Age, Pos,Neg)) - if(missing(k)){ - k <- switch(type, - "Muench" = 1 , - "Griffith" = 2, - "Grenfell" = 3)} - age <- function(k){ - if(k>1){ - formula<- paste0("I","(",paste("Age", 2:k,sep = "^"),")",collapse = "+") - paste0("cbind(Neg,Pos)"," ~","-1+Age+",formula) - } else { - paste0("cbind(Neg,Pos)"," ~","-1+Age") + # helper function to generate the polynomial given a k value + # to be used for parameter selection if multiple values for k are given + generate_polynomial <- function(k, df, link="log"){ + Age <- function(k){ + if(k>1){ + formula<- paste0("I","(",paste("age", 2:k,sep = "^"),")",collapse = "+") + paste0("cbind(neg,pos)"," ~","-1+age+",formula) + } else { + paste0("cbind(neg,pos)"," ~","-1+age") + } } + + glm(Age(k), family=binomial(link=link),df) + } + + # If a vector of values for k is provided -> select best value + if(length(k) > 1){ + out <- nested_mod_selection( + list("k" = k), + model_fn = \(k, df){ + generate_polynomial(k, df, link=link) + }, + dat = df + ) + + k <- out$best_par$k + model$info <- out$mod + }else{ + model$info <- generate_polynomial(k, df, link=link) } - model$info <- glm(age(k), family=binomial(link=link),df) - X <- X(Age, k) + + X <- X(age, k) model$sp <- 1 - model$info$fitted.values model$foi <- X%*%model$info$coefficients - model$df <- list(age=Age, pos=Pos, tot= Pos + Neg) + model$df <- list(age=age, pos=pos, tot= pos + neg) + model$k <- k class(model) <- "polynomial_model" model } -#' The Farrington (1990) model. -#' -#' Refers to section 6.1.2. -#' -#' @param data the input data frame, must either have `age`, `pos`, `tot` columns (for aggregated data) OR `age`, `status` for (linelisting data) -#' @param start Named list of vectors or single vector. -#' Initial values for optimizer. -#' @param fixed Named list of vectors or single vector. -#' Parameter values to keep fixed during optimization. -#' -#' @return a list of class farrington_model with 5 items -#' \item{datatype}{type of datatype used for model fitting (aggregated or linelisting)} -#' \item{df}{the dataframe used for fitting the model} -#' \item{info}{fitted "mle" object} -#' \item{sp}{seroprevalence} -#' \item{foi}{force of infection} -#' @seealso [stats4::mle()] for more information on the fitted mle object -#' -#' @examples -#' df <- rubella_uk_1986_1987 -#' model <- farrington_model( -#' df, -#' start=list(alpha=0.07,beta=0.1,gamma=0.03) -#' ) -#' plot(model) -#' -#' @importFrom stats4 mle -#' -#' @export -farrington_model <- function(data, start, fixed=list()) -{ - model <- list() +# TODO: check if this can be generalized to other functions as well (e.g. fractional polynomial) +# function to return the best parameter of nested glm models using LRT +# par_range - list of parameters and its possible values +# model_fn - function to fit and return a model, must takes 2 arguments: par, df +#' @import tidyr +#' @importFrom purrr pmap +nested_mod_selection <- function(par_range, model_fn, dat, method="LRT"){ + # generate all combinations of parameters values + par_combs <- tidyr::crossing(!!!par_range) - # check input whether it is line-listing or aggregated data - data <- check_input(data) - age <- data$age - pos <- data$pos - tot <- data$tot - model$datatype <- data$type + # fit the model using specified parameter values + mods_out <- par_combs %>% + purrr::pmap(\(...){model_fn(..., df=dat)}) + # perform LRT + lrt_out <- do.call( + anova, + c(mods_out, list(test="LRT")) + ) - farrington <- function(alpha,beta,gamma) { - p=1-exp((alpha/beta)*age*exp(-beta*age) - +(1/beta)*((alpha/beta)-gamma)*(exp(-beta*age)-1)-gamma*age) - ll=pos*log(p)+(tot-pos)*log(1-p) - return(-sum(ll)) - } + # get the best model + best_idx <- lrt_out %>% + as.data.frame() %>% + mutate( + idx = 1:n() + ) %>% + filter( + `Pr(>Chi)` < 0.05 + ) %>% + arrange( + Deviance + ) %>% + pull(idx) - model$info <- mle(farrington, fixed=fixed, start=start) - alpha <- model$info@coef[1] - beta <- model$info@coef[2] - gamma <- model$info@coef[3] - model$sp <- 1-exp( - (alpha/beta)*age*exp(-beta*age) - +(1/beta)*((alpha/beta)-gamma)*(exp(-beta*age)-1) - -gamma*age) - model$foi <- (alpha*age-gamma)*exp(-beta*age)+gamma - model$df <- list(age=age, pos=pos, tot=tot) + # handle scenario when the reference model (i.e., the first model) is in fact the best option + # i.e., when the other parameter combinations do not result in statistically significant improvement + if(length(best_idx)>1){ + best_idx <- best_idx[1] + }else{ + best_idx <- 1 + } - class(model) <- "farrington_model" - model + list( + best_par = par_combs[best_idx, ] %>% as.list(), + mod = mods_out[best_idx][[1]] + ) } diff --git a/R/predict.R b/R/predict.R new file mode 100644 index 0000000..5c65370 --- /dev/null +++ b/R/predict.R @@ -0,0 +1,253 @@ +# ====== Predict function for serosv models ====== +#' Prediction for serosv polynomial model +#' +#' A wrapper of predict.glm for direct prediction from polynomial_model object +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param ... arbitrary argument +#' +#' @importFrom stats predict.glm +#' @import dplyr +#' +#' @return prediction output +#' @seealso +#' [stats::predict.glm()] for more information on the predict function +#' @export +predict.polynomial_model <- function(object, newdata=NULL, ...){ + predict.glm(object$info, newdata, ...) +} + +#' Prediction for serosv fractional polynomial model +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param ... arbitrary argument +#' +#' @importFrom stats predict.glm +#' @return prediction output +#' @seealso +#' [stats::predict.glm()] for more information on the predict function +#' @export +predict.fp_model <- function(object, newdata=NULL, ...){ + + predict.glm(object$info,newdata=newdata, ...) +} + +#' Prediction for serosv Weibull model +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param ... arbitrary argument +#' +#' @importFrom stats predict.glm +#' @return prediction output +#' @seealso +#' [stats::predict.glm()] for more information on the predict function +#' @export +predict.weibull_model <- function(object, newdata=NULL, ...){ + predict.glm(object$info,data.frame("log(t)" = newdata$`log(t)`), ...) +} + + +#' Prediction for serosv local polynomial model +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param ... arbitrary argument +#' +#' @return prediction output +#' @export +predict.lp_model <- function(object, newdata=NULL,...){ + predict(object$info, data.frame(age = newdata[[1]]), ...) +} + + +#' Prediction for serosv penalized spline model +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param ... arbitrary argument +#' +#' @importFrom mgcv predict.gam +#' @return prediction output +#' @seealso +#' [mgcv::predict.gam()] for more information on the predict function +#' @export +predict.penalized_spline_model <- function(object, newdata=NULL,...){ + + # handle different output for different frameworks + if(object$framework == "pl"){ + gam_obj <- object$info + }else{ + gam_obj <- object$info$gam + } + + predict.gam(gam_obj, newdata, ...) +} + +#' Prediction for serosv Farrington model +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param ... arbitrary argument +#' +#' @return prediction output +#' @export +predict.farrington_model <- function(object, newdata=NULL,...){ + alpha <- object$info@coef[1] + beta <- object$info@coef[2] + gamma <- object$info@coef[3] + + 1-exp( + (alpha/beta)*newdata[[1]]*exp(-beta*newdata[[1]]) + +(1/beta)*((alpha/beta)-gamma)*(exp(-beta*newdata[[1]])-1) + -gamma*newdata[[1]]) +} + +#' Predict from an hierarchical bayesian model +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param ... arbitrary arguments +#' @import dplyr +#' +#' @return list of confidence interval for seroprevalence and foi. Each confidence interval dataframe with 4 variables, x and y for the fitted values and ymin and ymax for the confidence interval +#' @export +predict.hierarchical_bayesian_model <- function(object, newdata=NULL, ...){ + out_x <- object$df$age + out.DF <- NULL + + if (object$type == "far3"){ + alpha1 <- object$info["alpha1", "50%"] + alpha2 <- object$info["alpha2", "50%"] + alpha3 <- object$info["alpha3", "50%"] + + out.DF <- data.frame( + x = out_x, + y = object$sp_func(out_x, alpha1, alpha2, alpha3), + ) + }else if(object$type == "far2"){ + alpha1 <- object$info["alpha1", "50%"] + alpha2 <- object$info["alpha2", "50%"] + + out.DF <- data.frame( + x = out_x, + y = object$sp_func(out_x, alpha1, alpha2), + ) + + }else if(object$type == "log_logistic"){ + alpha1 <- object$info["alpha1", "50%"] + alpha2 <- object$info["alpha2", "50%"] + + out.DF <- data.frame( + x = out_x, + y = object$sp_func(out_x, alpha1, alpha2), + ) + }else{ + warning('Expect model type to be one of the following: "far3", "far2", "log_logistic"') + } + + out.DF +} + +#' Predict from the age_time_mdoel +#' +#' @param object serosv models +#' @param newdata data.frame with age column to generate prediction +#' @param modtype either "monotonized" (to predict using monotonized model) or "non-monotonized" +#' @param ... arbitrary argument +#' +#' @importFrom mgcv predict.gam +#' @import dplyr +#' +#' @return confidence interval dataframe with n_group x 3 cols, the columns are `group`, `sp_df`, `foi_df` +#' @export +predict.age_time_model <- function(object, newdata, modtype="monotonized", ...){ + # resolve no visible binding note + df <- monotonized_info <- monotonized_ci_mod <- age <- info <- fit <- se.fit <- sp_df <- foi_df <- NULL + + # check which type of model user wants to predict + modtype <- if (is.null(list(...)[["modtype"]])) "monotonized" else list(...)$modtype + assert_that( + modtype == "monotonized" | modtype == "non-monotonized", + msg = "modtype argument must be eithers 'monotonized' or 'non-monotonized'" + ) + + p <- (1 - object$ci) / 2 + + # check whether newdata match the requirement + if(!all(c(object$grouping_col, "age") %in% colnames(newdata)) ){ + stop(paste0( + "Data must have `", + object$grouping_col, + "`, `age` columns" + )) + } + + + # generate the newdata by survey time for prediction + out <- newdata %>% + group_by(.data[[object$grouping_col]]) %>% + nest() %>% + rename(age_df = data) %>% + left_join( + # can only predict for the survey time fitted to the model + object$out, + join_by(!!sym(object$grouping_col)) + ) + + # --- use the monotonized model for prediction and ci----- + if(modtype == "monotonized"){ + out <- out %>% + mutate( + sp_df = pmap(list(monotonized_info, monotonized_ci_mod, age_df), \(mod, ci_mod, grid){ + data.frame( + x = grid$age, + y = predict(mod, grid, type = "response"), + ymin = predict(ci_mod$ymin, grid, type = "response"), + ymax = predict(ci_mod$ymax, grid, type = "response") + ) + }) + ) + }else{ + # --- if user specify non-monotonized then simply compute CI from gam model----- + out <- out %>% + mutate( + sp_df = map2(info, age_df, \(mod, grid){ + link_inv <- mod$family$linkinv + dataset <- mod$model[,1:2] + n <- nrow(dataset) - length(mod$coefficients) + + predict(mod, grid, se.fit = TRUE) %>% + as_tibble() %>% + select(fit, se.fit) %>% + mutate( + x = grid, + ymin = link_inv(fit + qt( p, n) * se.fit), + ymax = link_inv(fit + qt(1 - p, n) * se.fit), + y = link_inv(fit) + ) %>% + select(- se.fit) + }) + ) + } + + # --- finally, compute FOI ----- + out <- out %>% + mutate( + foi_df = map2(age_df, sp_df, \(grid, sp){ + foi_x <- sort(unique(grid$age)) + foi_x <- foi_x[c(-1, -length(foi_x) )] + + data.frame( + x = foi_x, + y = est_foi(grid$age, sp$y) + ) + }) + ) %>% + select(!!sym(object$grouping_col), sp_df, foi_df) + + out +} + diff --git a/R/print.R b/R/print.R new file mode 100644 index 0000000..ca26328 --- /dev/null +++ b/R/print.R @@ -0,0 +1,143 @@ +#' @export +print.polynomial_model <- function(x, ...){ + cat("Polynomial model\n\n") + cat("Input type: ", x$datatype, "\n") + cat("Degree (k): ",x$k, "\n") + print(x$info) +} + +#' @export +print.fp_model <- function(x, ...){ + cat("Fractional polynomial model \n\n") + cat("Input type: ", x$datatype, "\n") + cat("Powers: ", paste(x$p, collapse=", "), "\n") + print(x$info) +} + +#' @export +print.age_time_model <- function(x, ...){ + cat("Age-time varying seroprevalence model \n\n") + cat("Input type: ", x$datatype, "\n") + cat("Grouping variable: ", x$grouping_col,"\n") + cat("Monotonization method: ", x$monotonize_method, "\n") + cat("Monotonize across: ", if(x$age_correct) "birth cohort" else "age group", "\n") + print(x$out) +} + +#' @export +print.farrington_model <- function(x, ...){ + cat("Farrington model \n\n") + cat("Input type: ", x$datatype, "\n") + # cat("Fitted parameters: ", + # paste(c("alpha", "beta", "gamma"), sprintf("%.4g", x$info@coef[1:3]), sep="=", collapse=", "), + # "\n\n") + print(x$info) +} + +#' @importFrom purrr compact +#' @export +print.hierarchical_bayesian_model <- function(x, ...){ + pars <- intersect(c("alpha1", "alpha2", "alpha3"), rownames(x$info)) + fitted_pars <- x$info[pars, "mean"] + # get the CrI + lower_pars <- x$info[pars, "2.5%"] + upper_pars <- x$info[pars, "97.5%"] + # get the sd + sd_pars <- purrr::compact(x$info[pars, "sd"]) + + cat("Hierarchical Bayesian model \n\n") + cat("Input type: ", x$datatype, "\n") + cat("Model: ", + switch(x$type, + far2 = "Farrington model with 2 parameters", + far3 = "Farrington model with 3 parameters", + log_logistic = "Log-logistic model"), + "\n\n") + cat("Fitted parameters:\n", + paste(names(fitted_pars), + # print fitted parameter in the format mean (95% CrI=lower-upper, sd=) + paste( + sprintf("%.4g", fitted_pars), + " (95% CrI [", + sprintf("%.4g", lower_pars), + ", ", + sprintf("%.4g", upper_pars), + "], sd = ", + sprintf("%.4g", sd_pars), + ")", + sep = "" + ), + sep=" = ", collapse="\n "), + "\n") +} + +#' @export +print.mixture_model <- function(x, ...){ + cat("Mixture model \n\n") + cat("Estimated proportion:\n", + paste( + c("Susceptible", "Infected"), + sprintf("%.4g", x$info$parameters$pi), + sep="=", collapse = ", " + ), + "\n\n" + ) + cat("Estimated mean Log(Antibody):\n", + paste( + c("Susceptible", "Infected"), + sprintf("%.4g", x$info$parameters$mu), + sep="=", collapse = ", " + ), + "\n") +} + +#' @export +print.estimate_from_mixture <- function(x, ...){ + cat("Age-varying seroprevalence estimated from mixture model \n\n") + cat("Monotonized seroprevalence: ", x$monotonize) + print(x$info) +} + +#' @export +print.lp_model <- function(x, ...){ + cat("Local polynomial model \n\n") + cat("Input type: ", x$datatype, "\n") + cat( + "Configs: ", + paste( + c("nn", "bandwidth(h)", "degree", "kernel"), + c(as.character(round( + as.numeric(x[c("nn", "h", "deg")]), digits = 4 + )), x$kern), + sep = "=", + collapse = ", " + ), + "\n\n" + ) + print(x$info) +} + +#' @export +print.penalized_spline_model <- function(x, ...){ + cat("Penalized spline model \n\n") + cat("Input type: ", x$datatype, "\n") + cat("Framework: ", if(x$framework=="pl") "Penalized likelihood" else "Mixed model", "\n") + print(x$info) +} + +#' @export +print.weibull_model <- function(x, ...){ + cat("Weibull model \n\n") + cat("Input type: ", x$datatype, "\n") + cat( + paste( + c("b0", "b1"), + sprintf("%.4g", coef(x$info)[1:2]), + sep="=", collapse=", " + ), + "\n" + ) + + print(x$info) +} + diff --git a/R/semiparametric_models.R b/R/semiparametric_models.R index 4fcd143..d867b8b 100644 --- a/R/semiparametric_models.R +++ b/R/semiparametric_models.R @@ -1,10 +1,99 @@ #' Penalized Spline model #' -#' @param data the input data frame, must either have `age`, `pos`, `tot` column for aggregated data OR `age`, `status` for linelisting data +#' @description Fit age-specific seroprevalence to a semi-parametric model where +#' predictor is modeled with penalized splines. The penalized splines can be estimated by either +#' (1) penalized likelihood framework or (2) mixed model framework +#' +#' @details +#' In the semi-parametric model, the predictor is formulated as a penalized spline +#' with truncated power basis functions of degree \eqn{p} +#' and fixed knots \eqn{\kappa_1,..., \kappa_k} as followed +#' +#' \deqn{ +#' \eta(a_i) = \beta_0 + \beta_1a_i + ... + \beta_p a_i^p + \Sigma_{k=1}^ku_k(a_i - \kappa_k)^p_+ +#' } +#' +#' - Where: +#' \deqn{ +#' (a_i - \kappa_k)^p_+ = \begin{cases} +#' 0, & a_i \le \kappa_k \\ +#' (a_i - \kappa_k)^p, & a_i > \kappa_k +#' \end{cases} +#' } +#' +#' FOI can then be derived by +#' +#' \deqn{\hat{\lambda}(a_i) = [\hat{\beta_1} , 2\hat{\beta_2}a_i, ..., +#' p \hat{\beta} a_i ^{p-1} + \Sigma^k_{k=1} p \hat{u}_k(a_i - \kappa_k)^{p-1}_+] \delta(\hat{\eta}(a_i)) +#' } +#' +#' Where \eqn{\delta(.)} is determined by the link function used in the model +#' +#' In matrix annotation, the mean structure model for \eqn{\eta(a_i)} becomes +#' \deqn{\eta = \textbf{X}\beta + \textbf{Zu}} +#' +#' Where \eqn{\eta = [\eta(a_i) ... \eta(a_N) ]^T}, \eqn{\beta = [\beta_0 \beta_1 .... \beta_p]^T}, +#' and \eqn{\textbf{u} = [u_1 u_2 ... u_k]^T} are the regression with corresponding design matrices +#' \deqn{ +#' +#' \textbf{X} = \begin{bmatrix} +#' 1 & a_1 & a_1^2 & ... & a_1^p \\ +#' 1 & a_2 & a_2^2 & ... & a_2^p \\ +#' \vdots & \vdots & \vdots & \dots & \vdots \\ +#' 1 & a_N & a_N^2 & ... & a_N^p +#' \end{bmatrix}, \textbf{Z} = \begin{bmatrix} +#' (a_1 - \kappa_1 )_+^p & (a_1 - \kappa_2 )_+^p & \dots & (a_1 - \kappa_k)_+^p \\ +#' (a_2 - \kappa_1 )_+^p & (a_2 - \kappa_2 )_+^p & \dots & (a_2 - \kappa_k)_+^p \\ +#' \vdots & \vdots & \dots & \vdots \\ +#' (a_N - \kappa_1 )_+^p & (a_N - \kappa_2 )_+^p & \dots & (a_N - \kappa_k)_+^p +#' \end{bmatrix} +#' } +#' +#' Under \bold{penalized likelihood framework}, the model is fitted by maximizing +#' the following likelihood +#' +#' \deqn{ +#' \phi^{-1}[y^T(\textbf{X}\beta + \textbf{Zu} ) - \textbf{1}^Tc(\textbf{X}\beta + \textbf{Zu} )] - \frac{1}{2}\lambda^2 +#' \begin{bmatrix} \beta \\ \textbf{u} \end{bmatrix}^T D\begin{bmatrix} \beta \\ \textbf{u} \end{bmatrix} +#' } +#' +#' Where: +#' +#' - \eqn{X\beta + Zu} is the predictor +#' +#' - \eqn{D} is a known semi-definite penalty matrix [@Wahba1978], [@Green1993] +#' +#' - \eqn{y} is the response vector +#' +#' - \eqn{\textbf{1}} the unit vector, \eqn{c(.)} is determined by the link function used +#' +#' - \eqn{\lambda} is the smoothing parameter (larger values –> smoother curves) +#' +#' - \eqn{\phi} is the overdispersion parameter and equals 1 if there is no overdispersion +#' +#' Under the \bold{mixed model} framework, +#' the model instead treats the coefficients \eqn{\textbf{u}} in the likelihood formulation +#' as random effects with \eqn{\textbf{u} \sim N(\textbf{0}, \mathbf{\sigma}^2_u \textbf{I})} +#' +#' Refer to section 8.1 and 8.2 of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. +#' +#' @param data the input data frame, must either have columns for `age`, `pos`, `tot` (for aggregated data) OR +#' columns for `age`, `status` (for linelisting data) #' @param s smoothing basis to use #' @param sp smoothing parameter #' @param link link function to use #' @param framework which approach to fit the model ("pl" for penalized likelihood framework, "glmm" for generalized linear mixed model framework) +#' @param age_col name of the `age` column (default age_col="age"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). #' #' @importFrom mgcv gam gamm #' @importFrom stats binomial @@ -27,10 +116,12 @@ #' model <- penalized_spline_model(data, framework="glmm") #' model$info$gam #' plot(model) -penalized_spline_model <- function(data, s = "bs", link = "logit", framework = "pl", sp = NULL){ +penalized_spline_model <- function(data, + age_col="age",pos_col="pos", tot_col="tot", status_col="status", + s = "bs", link = "logit", framework = "pl", sp = NULL){ model <- list() - data <- check_input(data) + data <- check_input(data, stratum_col=age_col,pos_col=pos_col, tot_col=tot_col, status_col=status_col) age <- data$age pos <- data$pos tot <- data$tot diff --git a/R/sir_basic_model.R b/R/sir_basic_model.R deleted file mode 100644 index 50cf76d..0000000 --- a/R/sir_basic_model.R +++ /dev/null @@ -1,72 +0,0 @@ -sir_basic <- function(t, state, parameters) -{ - with(as.list(c(state, parameters)), { - N <- S + I + R # total population - lambda <- beta*I # force of infection - dS <- N*mu*(1-p) - (lambda + mu)*S - dI <- lambda*S - (nu+alpha+mu)*I - dR <- N*mu*p + nu*I - mu*R - list(c(dS, dI, dR)) - }) -} - -#' Basic SIR model -#' -#' Refers to section 3.1.3. -#' -#' @details -#' In \code{state}: -#' -#' - \code{S}: number of susceptible -#' -#' - \code{I}: number of infected -#' -#' - \code{R}: number of recovered -#' -#' In \code{parameters}: -#' -#' - \code{alpha}: disease-related death rate -#' -#' - \code{mu}: natural death rate (= 1/life expectancy) -#' -#' - \code{beta}: transmission rate -#' -#' - \code{nu}: recovery rate -#' -#' - \code{p}: percent of population vaccinated at birth -#' -#' @param times time sequence. -#' -#' @param state the initial state of the model. -#' -#' @param parameters the parameters of the model. -#' -#' @examples -#' state <- c(S=4999, I=1, R=0) -#' parameters <- c( -#' mu=1/75, # 1 divided by life expectancy (75 years old) -#' alpha=0, # no disease-related death -#' beta=0.0005, # transmission rate -#' nu=1, # 1 year for infected to recover -#' p=0 # no vaccination at birth -#' ) -#' times <- seq(0, 250, by=0.1) -#' model <- sir_basic_model(times, state, parameters) -#' model -#' -#' @importFrom deSolve ode -#' @return list of class sir_basic_model with the following items -#' \item{parameters}{list of parameters used for fitting the model} -#' \item{output}{matrix of population for each compartment over time} -#' -#' -#' @export -sir_basic_model <- function(times, state, parameters) -{ - model <- list() - model$parameters <- parameters - model$output <- as.data.frame(ode(y=state,times=times,func=sir_basic,parms=parameters)) - - class(model) <- "sir_basic_model" - model -} diff --git a/R/sir_static_model.R b/R/sir_static_model.R deleted file mode 100644 index 2708abe..0000000 --- a/R/sir_static_model.R +++ /dev/null @@ -1,61 +0,0 @@ -sir_static <- function(t, state, parameters) -{ - with(as.list(c(state, parameters)), { - ds <- -lambda*s - di <- lambda*s - nu*i - dr <- nu*i - list(c(ds, di, dr)) - }) -} - -#' SIR static model (age-heterogeneous, endemic equilibrium) -#' -#' Refers to section 3.2.2. -#' -#' @details -#' In \code{state}: -#' -#' - \code{s}: proportion susceptible -#' -#' - \code{i}: proportion infected -#' -#' - \code{r}: proportion recovered -#' -#' In \code{parameters}: -#' -#' - \code{lambda}: natural death rate -#' -#' - \code{nu}: recovery rate -#' -#' @param a age sequence. -#' -#' @param state the initial state of the system. -#' -#' @param parameters the model's parameter. -#' -#' @examples -#' state <- c(s=0.99,i=0.01,r=0) -#' parameters <- c( -#' lambda = 0.05, -#' nu=1/(14/365) # 2 weeks to recover -#' ) -#' ages<-seq(0, 90, by=0.01) -#' model = sir_static_model(ages, state, parameters) -#' model -#' -#' @return list of class sir_static_model with the following items -#' \item{parameters}{list of parameters used for fitting the model} -#' \item{output}{matrix of proportion for each compartment over time} -#' -#' -#' @export -sir_static_model <- function(a, state, parameters) -{ - - model <- list() - model$parameters <- parameters - model$output <- as.data.frame(ode(y=state,times=a,func=sir_static,parms=parameters)) - - class(model) <- "sir_static_model" - model -} diff --git a/R/sir_subpops_model.R b/R/sir_subpops_model.R deleted file mode 100644 index 35a2dc8..0000000 --- a/R/sir_subpops_model.R +++ /dev/null @@ -1,114 +0,0 @@ -ds <- function(state, parameters, i) -{ - with(as.list(c(state, parameters)), { - sum_beta_i <- 0 - for (j in 1:k) { - sum_beta_i <- sum_beta_i + beta[i,j]*get(paste0("i", j)) - } - -sum_beta_i*get(paste0("s", i)) + mu - mu*get(paste0("s", i)) - }) -} - -di <- function(state, parameters, i) -{ - with(as.list(c(state, parameters)), { - sum_beta_i <- 0 - for (j in 1:k) { - sum_beta_i <- sum_beta_i + beta[i,j]*get(paste0("i", j)) - } - sum_beta_i*get(paste0("s", i)) - nu[i]*get(paste0("i", i)) - mu*get(paste0("i", i)) - }) -} - -dr <- function(state, parameters, i) -{ - with(as.list(c(state, parameters)), { - nu[i]*get(paste0("i", i)) - mu*get(paste0("r", i)) - }) -} - -sir_subpop <- function(t, state, parameters) { - with(as.list(c(state, parameters)), { - s_states <- c() - i_states <- c() - r_states <- c() - - for (i in 1:k) { - s_states <- c(s_states, ds(state, parameters, i)) - i_states <- c(i_states, di(state, parameters, i)) - r_states <- c(r_states, dr(state, parameters, i)) - } - - list(c(s_states, i_states, r_states)) - }) -} - -#' SIR Model with Interacting Subpopulations -#' -#' Refers to section 3.5.1. -#' -#' @details -#' In \code{state}: -#' -#' - \code{s}: Percent susceptible -#' -#' - \code{i}: Percent infected -#' -#' - \code{r}: Percent recovered -#' -#' In \code{parameters}: -#' -#' - \code{mu}: natural death rate (1/L). -#' -#' - \code{beta}: transmission rate w.r.t population (beta tilde) -#' -#' - \code{nu}: recovery rate -#' -#' - \code{k}: number of subpopulations -#' -#' @param times time sequence. -#' -#' @param state the initial state of the model. -#' -#' @param parameters the parameters of the model. -#' -#' @examples -#' \donttest{ -#' k <- 2 -#' state <- c( -#' s = c(0.8, 0.8), -#' i = c(0.2, 0.2), -#' r = c( 0, 0) -#' ) -#' beta_matrix <- c( -#' c(0.05, 0.00), -#' c(0.00, 0.05) -#' ) -#' parameters <- list( -#' beta = matrix(beta_matrix, nrow=k, ncol=k, byrow=TRUE), -#' nu = c(1/30, 1/30), -#' mu = 0.001, -#' k = k -#' ) -#' times<-seq(0,10000,by=0.5) -#' model <- sir_subpops_model(times, state, parameters) -#' model -#' } -#' -#' @return list of class sir_subpops_model with the following items -#' -#' \item{parameters}{list of parameters used for fitting the model} -#' \item{output}{matrix of proportion for each compartment over time} -#' -#' -#' @export -sir_subpops_model <- function(times, state, parameters) { - model <- list() - model$parameters <- parameters - model$output <- as.data.frame( - ode(y=state,times=times,func=sir_subpop,parms=parameters) - ) - - class(model) <- "sir_subpops_model" - model -} diff --git a/R/to_titer.R b/R/to_titer.R index 2735f89..5e98220 100644 --- a/R/to_titer.R +++ b/R/to_titer.R @@ -1,18 +1,28 @@ -#' Process assay test result to titer +#' Convert assay readings to titers +#' +#' to_titer() converts raw assay readings (e.g., OD, fluorescence intensity) to titer by +#' fitting a calibrating model +#' +#' @param df a standardized data.frame returned by`standardize_data()` +#' @param model either: +#' \enumerate{ +#' \item {A string naming a built-in model (currently supported: \code{"4PL"}), or} +#' \item {A named list with two functions: \code{$mod} for curve fitting and +#' \code{$quantify_ci} for titer estimation with confidence intervals.} +#' } +#' @param positive_threshold if not NULL, processed_data will have the serostatus labeled +#' @param ci confidence interval for the titer estimates (default is .95 i.e., 95\% CI) +#' @param negative_control if TRUE, output tibble will include the result for negative controls #' -#' @param df - a standardized data.frame generated by`standardize_data()` -#' @param model - either the name of a built-in model to fit standard curve or a named list of 2 functions for "mod" and "quantify_ci" -#' @param positive_threshold - if not NULL, processed_data will have the serostatus labeled -#' @param ci - confidence interval for the titer estimates -#' @param negative_control - if TRUE, output tibble will include the result for negative controls #' @importFrom magrittr %>% -#' @import dplyr purrr +#' @import dplyr +#' @importFrom purrr map_dfr #' #' @return a data.frame with 8 columns #' \item{plate_id}{id of the plate} -#' \item{data}{list of `data.frame`s containing the results from each plate} -#' \item{antitoxin_df}{list of `data.frame`s containing the results for antitoxins from each plate} -#' \item{standard_curve_func}{list of functions mapping from OD to titer for each plate} +#' \item{data}{list of `data.frame`s containing the raw sample results from each plate} +#' \item{antitoxin_df}{list of `data.frame`s containing the raw results for antitoxins from each plate} +#' \item{standard_curve_func}{list of functions mapping from assay reading to titer for each plate} #' \item{std_crv_midpoint}{midpoint of the standard curve, for qualitative analysis} #' \item{processed_data}{list of `tibble`s containing samples with titer estimates (lower, median, upper)} #' \item{negative_control}{list of `tibble`s containing negative control check results (if `negative_control=TRUE`)} @@ -73,13 +83,13 @@ to_titer <- function(df, model="4PL", positive_threshold=NULL, ci = .95, #' Preprocess data #' -#' @param df - data.frame with columns for plate id, sample id, result, dilution factor, and (optionally) negative controls -#' @param plate_id_col - name of the column storing plates id -#' @param id_col - name of the column storing sample id -#' @param result_col - name of the column storing result -#' @param dilution_fct_col - name of the column storing dilution factors -#' @param antitoxin_label - how antitoxin is label in the sample id column -#' @param negative_col - regex for columns for negative controls, assumed to be a label followed by the dilution factor (e.g. NEGATIVE_50, NEGATIVE_100) +#' @param df data.frame with columns for plate id, sample id, result, dilution factor, and (optionally) negative controls +#' @param plate_id_col name of the column storing plates id +#' @param id_col name of the column storing sample id +#' @param result_col name of the column storing result +#' @param dilution_fct_col name of the column storing dilution factors +#' @param antitoxin_label how antitoxin is label in the sample id column +#' @param negative_col regex for columns for negative controls, assumed to be a label followed by the dilution factor (e.g. NEGATIVE_50, NEGATIVE_100) #' #' @importFrom magrittr %>% #' @importFrom janitor clean_names @@ -157,9 +167,10 @@ data2function <- function(df) { } # Function to convert samples' OD to LC +#' @importFrom purrr map_dfr process_samples <- function(plate, std_crv, midpoint=2, positive_threshold = NULL) { out <- plate %>% - bind_cols(map_dfr(.$result, std_crv)) + bind_cols(purrr::map_dfr(.$result, std_crv)) if(is.numeric(positive_threshold)){ out <- out %>% mutate( @@ -205,6 +216,7 @@ nls4PL <- function(df) { # - sample different values for parameters (number of samples = nb) # - use model to compute OD for the new set of parameter values #' @importFrom mvtnorm rmvnorm +#' @importFrom purrr map_dfc simulate_nls_ci <- function(object, newdata, nb = 9999, alpha = .025) { rowsplit <- function(df) split(df, 1:nrow(df)) @@ -214,7 +226,7 @@ simulate_nls_ci <- function(object, newdata, nb = 9999, alpha = .025) { rowsplit() |> map(as.list) |> map(~ c(.x, newdata)) |> - map_dfc(eval, expr = parse(text = as.character(formula(object))[3])) %>% + purrr::map_dfc(eval, expr = parse(text = as.character(formula(object))[3])) %>% apply(1, quantile, c(alpha, .5, 1 - alpha)) %>% t() %>% as.data.frame() %>% setNames(c("lower", "median", "upper")) @@ -233,11 +245,11 @@ get_negative_controls <- function(plate, std_crv){ # ======== Plot functions ========== #' Visualize standard curves for each plate #' -#' @param x - output of `to_titer()` -#' @param facet - whether to faceted by plates or plot all standard curves on a single plot -#' @param xlab - label of the x axis -#' @param ylab - label of the y axis -#' @param datapoint_size - size of the data point (only applicable when `facet=TRUE`) +#' @param x output of `to_titer()` +#' @param facet whether to faceted by plates or plot all standard curves on a single plot +#' @param xlab label of the x axis +#' @param ylab label of the y axis +#' @param datapoint_size size of the data point (only applicable when `facet=TRUE`) #' #' @importFrom magrittr %>% #' @import ggplot2 dplyr @@ -293,7 +305,7 @@ plot_standard_curve <- function(x, facet=TRUE, #' @param shift_text adjust how much the text is shifted along the x-axis (relative to the threshold line) #' #' @importFrom magrittr %>% -#' @import ggplot2 dplyr purrr +#' @import ggplot2 dplyr #' #' @export add_thresholds <- function(dilution_factors, positive_threshold = 0.1, @@ -333,7 +345,8 @@ add_thresholds <- function(dilution_factors, positive_threshold = 0.1, #' @param n_dilutions - number of dilutions used for testing #' #' @importFrom magrittr %>% -#' @import ggplot2 dplyr purrr +#' @importFrom purrr walk +#' @import ggplot2 dplyr #' #' @export plot_titer_qc <- function(x, n_plates=18, n_samples=22, n_dilutions = 3){ @@ -371,7 +384,7 @@ plot_titer_qc <- function(x, n_plates=18, n_samples=22, n_dilutions = 3){ } opar <- par(mfrow = c(1, n_plates), mar = c(0,0.25,0,0.25)) - walk(head(x$processed_data, n=n_plates), plot_heatmap) + purrr::walk(head(x$processed_data, n=n_plates), plot_heatmap) par(opar) } diff --git a/R/utils.R b/R/utils.R index e288f77..af63e99 100644 --- a/R/utils.R +++ b/R/utils.R @@ -55,18 +55,19 @@ pava<- function(pos=pos,tot=rep(1,length(pos))) return(list(pai1=pai1,pai2=pai2)) } +# TODO: update aggregate func here to be "pipe"-able #' Aggregate data #' #' Generate a dataframe with `t`, `pos` and `tot` columns from #' `t` and `seropositive` vectors. #' -#' @param t the time vector (for stratification). -#' @param spos the seropositive vector. -#' @param stratum_col new name for the time vector (default to "t") +#' @param data a data frame with columns for age and serostatus +#' @param status_col name of the column for serostatus +#' @param stratum_col name of the column to stratify by (default to "age") #' #' @examples #' df <- hcv_be_2006 -#' hcv_df <- transform_data(df$dur, df$seropositive) +#' hcv_df <- transform_data(df, stratum_col="dur", status_col="seropositive") #' hcv_df #' #' @importFrom dplyr group_by @@ -76,15 +77,29 @@ pava<- function(pos=pos,tot=rep(1,length(pos))) #' #' @return dataframe in aggregated format #' @export -transform_data <- function(t, spos, stratum_col = "t") { - df <- data.frame(t, spos) +transform_data <- function(data, stratum_col="age", status_col="status") { + df <- NULL + + if( all(c(stratum_col, status_col) %in% names(data)) ) { + df <- data.frame( + age = data[[stratum_col]], + status = data[[status_col]] + ) + }else{ + stop(paste0( + "Data must have `", + stratum_col, + "`, `",status_col ,"` columns" + )) + } + + df_agg <- df %>% - group_by(t) %>% + group_by(age) %>% summarize( - pos = sum(spos), + pos = sum(status), tot = n() ) - colnames(df_agg) <- c(stratum_col, "pos", "tot") df_agg } @@ -94,7 +109,8 @@ transform_data <- function(t, spos, stratum_col = "t") { # - type of data (either linelisting or aggregated) # - preprocessed pos and tot columns #' @importFrom assertthat assert_that -check_input <- function(data, stratum_col = "age"){ +check_input <- function(data, pos_col="pos",tot_col="tot",status_col="status", + stratum_col = "age"){ assert_that( is.data.frame(data), msg = "Input must be a data.frame or tibble" @@ -106,23 +122,23 @@ check_input <- function(data, stratum_col = "age"){ type <- NULL - if( all(c(stratum_col, "pos", "tot") %in% colnames(data)) ){ + if( all(c(stratum_col, pos_col, tot_col) %in% colnames(data)) ){ age <- as.numeric(data[[stratum_col]]) - pos <- as.numeric(data$pos) - tot <- as.numeric(data$tot) + pos <- as.numeric(data[[pos_col]]) + tot <- as.numeric(data[[tot_col]]) type <- "aggregated" - }else if( all(c(stratum_col, "status") %in% colnames(data)) ){ + }else if( all(c(stratum_col, status_col) %in% colnames(data)) ){ age <- as.numeric(data[[stratum_col]]) - pos <- as.numeric(data$status) - tot <- rep(1, length(data$status)) + pos <- as.numeric(data[[status_col]]) + tot <- rep(1, length(data[[status_col]])) type <- "linelisting" }else{ stop(paste0( "Data must have `", stratum_col, - "`, `pos`, `tot` columns for aggregated data OR `", + "`, `", pos_col, "`, `", tot_col,"` columns for aggregated data OR `", stratum_col, - "`, `status` columns for linelisting data" + "`, `",status_col ,"` columns for linelisting data" )) } diff --git a/R/weibull_model.R b/R/weibull_model.R index 05631c6..9aac7d0 100644 --- a/R/weibull_model.R +++ b/R/weibull_model.R @@ -1,8 +1,35 @@ #' The Weibull model. #' -#' Refers to section 6.1.2. +#' @description Model seroprevalence as a function of duration since vaccination using the Weibull +#' model, where the force of infection is assumed to vary monotonically with duration. #' -#' @param data the input data frame, must either have `t`, `pos`, `tot` column for aggregated data OR `t`, `status` for linelisting data +#' @details +#' For a Weibull model, the prevalence is given by +#' \deqn{ +#' \pi (d) = 1 - e^{ - \beta_0 d ^ {\beta_1}} +#' } +#' Where \eqn{d} is exposure time (difference between age of vaccination and age at test) +#' +#' Which implies the force of infection to be the monotonic function +#' \deqn{ +#' \lambda(d) = \beta_0 \beta_1 d^{\beta_1 - 1} +#' } +#' +#' Refer to section 6.1.2. of the the book by Hens et al. (2012) for further details. +#' +#' @references +#' Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, +#' and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on +#' Serological and Social Contact Data: A Modern Statistical Perspective. +#' tatistics for Biology and Health. Springer New York. +#' \doi{https://doi.org/10.1007/978-1-4614-4072-7}. +#' +#' @param data the input data frame, must either have columns for `t`, `pos`, `tot` (for aggregated data) OR +#' `t`, `status` (for linelisting data) +#' @param t_lab name of the `t` column (default t_lab="t"). +#' @param pos_col name of the `pos` column (default pos_col="pos"). +#' @param tot_col name of the `tot` column (default tot_col="tot"). +#' @param status_col name of the `status` column (default status_col="status"). #' #' @importFrom stats coef #' @@ -10,7 +37,7 @@ #' df <- hcv_be_2006[order(hcv_be_2006$dur), ] #' df$t <- df$dur #' df$status <- df$seropositive -#' model <- weibull_model(df) +#' model <- weibull_model(df, t_lab="dur", status_col="seropositive") #' plot(model) #' #' @return list of class weibull_model with the following items @@ -23,12 +50,13 @@ #' @seealso [stats::glm()] for more information on the fitted "glm" object #' #' @export -weibull_model <- function(data) +weibull_model <- function(data, + t_lab="t",pos_col="pos", tot_col="tot", status_col="status") { model <- list() # check input whether it is line-listing or aggregated data - data <- check_input(data, stratum_col = "t") + data <- check_input(data, stratum_col = t_lab, pos_col=pos_col, tot_col=tot_col, status_col=status_col) t <- data$age pos <- data$pos tot <- data$tot diff --git a/README.Rmd b/README.Rmd index dadb13d..4095c80 100644 --- a/README.Rmd +++ b/README.Rmd @@ -49,14 +49,26 @@ The following methods are available to estimate seroprevalence and force of infe Parametric approaches: -- Polynomial models: - - Muench's model - - Griffiths' model - - Grenfell and Anderson's model -- Nonlinear models: - - Farrington's model - - Weibull model -- Fractional polynomial models +- Frequentist methods: + + - Polynomial models: + + - Muench's model + - Griffiths' model + - Grenfell and Anderson's model + + - Nonlinear models: + + - Farrington's model + - Weibull model + + - Fractional polynomial models + +- Bayesian methods: + + - Hierarchical Farrington model + + - Hierarchical log-logistic model Nonparametric approaches: @@ -70,12 +82,6 @@ Semiparametric approaches: - Generalized Linear Mixed Model framework -Hierarchical Bayesian approaches: - -- Hierarchical Farrington model - -- Hierarchical log-logistic model - ## Demo ### Fitting rubella data from the UK @@ -84,34 +90,33 @@ Load the rubella in UK dataset. ```{r} library(serosv) +rubella <- rubella_uk_1986_1987 ``` -Find the power for the best second degree fractional polynomial with monotonicity constraint and a logit link function. The power appears to be (-0.9,-0.9). +Fit the data using a fractional polynomial model via `fp_model()`. In this example, the model searches for the best combination of powers within a specified range. ```{r} -rubella <- rubella_uk_1986_1987 - -best_2d_mn <- find_best_fp_powers( +rubella_mod <- fp_model( rubella, - p=seq(-2,3,0.1), mc = T, degree=2, link="logit" + p=list( + p_range=seq(-2,3,0.1), # range of powers to search over + degree=2 # maximum degree for the search + ), + monotonic = T, # enforce model to be monotonic + link="logit" ) - -best_2d_mn +rubella_mod ``` -Finally, fit the second degree fractional polynomial. +Visualize the model ```{r} -fpmd <- fp_model( - rubella, - p=c(-0.9, -0.9), link="logit") - -plot(fpmd) +plot(rubella_mod) ``` ### Fitting Parvo B19 data from Finland -```{r} +```{r warning=FALSE, message=FALSE} library(dplyr) parvob19 <- parvob19_fi_1997_1998 @@ -120,12 +125,12 @@ transform_data( parvob19$age, parvob19$seropositive, stratum_col = "age") |> - polynomial_model(type = "Muench") |> + polynomial_model(k = 1) |> plot() # or fit data as is parvob19 |> rename(status = seropositive) |> - polynomial_model(type = "Muench") |> + polynomial_model(k = 1) |> plot() ``` diff --git a/README.md b/README.md index eaf53b6..5bb731c 100644 --- a/README.md +++ b/README.md @@ -47,14 +47,26 @@ of infection. Parametric approaches: -- Polynomial models: - - Muench’s model - - Griffiths’ model - - Grenfell and Anderson’s model -- Nonlinear models: - - Farrington’s model - - Weibull model -- Fractional polynomial models +- Frequentist methods: + + - Polynomial models: + + - Muench’s model + - Griffiths’ model + - Grenfell and Anderson’s model + + - Nonlinear models: + + - Farrington’s model + - Weibull model + + - Fractional polynomial models + +- Bayesian methods: + + - Hierarchical Farrington model + + - Hierarchical log-logistic model Nonparametric approaches: @@ -68,12 +80,6 @@ Semiparametric approaches: - Generalized Linear Mixed Model framework -Hierarchical Bayesian approaches: - -- Hierarchical Farrington model - -- Hierarchical log-logistic model - ## Demo ### Fitting rubella data from the UK @@ -82,30 +88,30 @@ Load the rubella in UK dataset. ``` r library(serosv) +rubella <- rubella_uk_1986_1987 ``` -Find the power for the best second degree fractional polynomial with -monotonicity constraint and a logit link function. The power appears to -be (-0.9,-0.9). +Fit the data using a fractional polynomial model via `fp_model()`. In +this example, the model searches for the best combination of powers +within a specified range. ``` r -rubella <- rubella_uk_1986_1987 - -best_2d_mn <- find_best_fp_powers( +rubella_mod <- fp_model( rubella, - p=seq(-2,3,0.1), mc = T, degree=2, link="logit" + p=list( + p_range=seq(-2,3,0.1), # range of powers to search over + degree=2 # maximum degree for the search + ), + monotonic = T, # enforce model to be monotonic + link="logit" ) - -best_2d_mn -#> $p -#> [1] -0.9 -0.9 -#> -#> $deviance -#> [1] 37.57966 +rubella_mod +#> Fractional polynomial model #> -#> $model +#> Input type: aggregated +#> Powers: -0.9, -0.9 #> -#> Call: glm(formula = as.formula(formulate(p_cur)), family = binomial(link = link)) +#> Call: glm(formula = as.formula(formulate(curr_p)), family = binomial(link = link)) #> #> Coefficients: #> (Intercept) I(age^-0.9) @@ -118,14 +124,10 @@ best_2d_mn #> Residual Deviance: 37.58 AIC: 210.1 ``` -Finally, fit the second degree fractional polynomial. +Visualize the model ``` r -fpmd <- fp_model( - rubella, - p=c(-0.9, -0.9), link="logit") - -plot(fpmd) +plot(rubella_mod) ``` @@ -134,15 +136,6 @@ plot(fpmd) ``` r library(dplyr) -#> Warning: package 'dplyr' was built under R version 4.3.1 -#> -#> Attaching package: 'dplyr' -#> The following objects are masked from 'package:stats': -#> -#> filter, lag -#> The following objects are masked from 'package:base': -#> -#> intersect, setdiff, setequal, union parvob19 <- parvob19_fi_1997_1998 # for linelisting data, either transform it to aggregated @@ -150,7 +143,7 @@ transform_data( parvob19$age, parvob19$seropositive, stratum_col = "age") |> - polynomial_model(type = "Muench") |> + polynomial_model(k = 1) |> plot() ``` @@ -161,7 +154,7 @@ transform_data( # or fit data as is parvob19 |> rename(status = seropositive) |> - polynomial_model(type = "Muench") |> + polynomial_model(k = 1) |> plot() ``` diff --git a/_pkgdown.yml b/_pkgdown.yml index 3e3609b..d634c66 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,6 +1,7 @@ url: https://oucru-modelling.github.io/serosv/ template: bootstrap: 5 + math-rendering: mathjax reference: - title: Dataset @@ -42,10 +43,6 @@ reference: - title: Models contents: - - sir_basic_model - - sir_static_model - - sir_subpops_model - - mseir_model - polynomial_model - farrington_model - weibull_model @@ -68,10 +65,6 @@ reference: - plot.weibull_model - plot.fp_model - plot.lp_model - - plot.mseir_model - - plot.sir_basic_model - - plot.sir_static_model - - plot.sir_subpops_model - plot.hierarchical_bayesian_model - plot.penalized_spline_model - plot.mixture_model @@ -90,7 +83,7 @@ reference: - transform_data - standardize_data - compare_models - - compute_ci + - compute_ci.default - compute_ci.fp_model - compute_ci.lp_model - compute_ci.weibull_model @@ -98,6 +91,14 @@ reference: - compute_ci.mixture_model - compute_ci.age_time_model - compute_ci.hierarchical_bayesian_model + - predict.age_time_model + - predict.farrington_model + - predict.fp_model + - predict.hierarchical_bayesian_model + - predict.lp_model + - predict.penalized_spline_model + - predict.weibull_model + - predict.polynomial_model - find_best_fp_powers articles: @@ -105,16 +106,15 @@ articles: navbar: ~ contents: - linelisting_vs_aggregated - - sir_model - parametric_model - nonparametric_model - semiparametric_model - - hierarchical_model - model_quantitative_data - repeated_cross_sectional - title: Utilities navbar: Utilities contents: + - reading_to_titer - data_transformation - imperfect_test - visualizing_model diff --git a/docs/404.html b/docs/404.html index 39e5eba..4d380e5 100644 --- a/docs/404.html +++ b/docs/404.html @@ -11,7 +11,13 @@ - + @@ -26,7 +32,7 @@ serosv - 1.2.0 + 1.2.0.9000 +

Penalized likelihood framework

-

Refer to Chapter 8.2.1

Proposed approach

-

A first approach to fit the model is by maximizing the following penalized likelihood

-\[\begin{equation} -\phi^{-1}[y^T(X\beta + Zu ) - 1^Tc(X\beta + Zu )] - \frac{1}{2}\lambda^2 -\begin{bmatrix} \beta \\u \end{bmatrix}^T D\begin{bmatrix} \beta \\u \end{bmatrix} +

The first approach to fit the model is by maximizing the following penalized likelihood

+

\[\begin{equation} +\phi^{-1}[y^T(\textbf{X}\beta + \textbf{Zu} ) - \textbf{1}^Tc(\textbf{X}\beta + \textbf{Zu} )] - \frac{1}{2}\lambda^2 +\begin{bmatrix} \beta \\ \textbf{u} \end{bmatrix}^T D\begin{bmatrix} \beta \\ \textbf{u} \end{bmatrix} \tag{1} -\end{equation}\] +\end{equation}\]

Where:

+

Refer to Chapter 8.2.1 of the book by Hens et al. (2012) for a more detailed explanation of the method.

Fitting data

To fit the data using the penalized likelihood framework, specify framework = "pl"

Basis function can be defined via the s parameter, some values for s includes:

-

For more options, refer to the mgcv documentation (Wood 2017)

+

For more options, refer to the mgcv documentation (Wood 2017)

 pl <- parvob19_be_2001_2003 %>% 
-  rename(status = seropositive) %>% 
-  penalized_spline_model(s = "tp", framework = "pl") 
-pl$info
+  penalized_spline_model(status_col = "seropositive", s = "tp", framework = "pl") 
+pl
+#> Penalized spline model 
+#> 
+#> Input type:  linelisting 
+#> Framework:  Penalized likelihood 
 #> 
 #> Family: binomial 
 #> Link function: logit 
@@ -216,11 +198,10 @@ 

Penalized likelihood framework

Generalized Linear Mixed Model framework

-

Refer to Chapter 8.2.2

Proposed approach

-

Looking back at (1), a constraint for \(u\) would be \(\Sigma_ku_k^2 < C\) for some positive value \(C\)

-

This is equivalent to choosing \((\beta, u)\) to maximise (1) with \(D = diag(0, 1)\) where \(0\) denotes zero vector length \(p+1\) and 1 denotes the unit vector of length \(K\)

-

For a fixed value for \(\lambda\) this is equivalent to fitting the following generalized linear mixed model Ngo and Wand (2004)

+

Looking back at (1), a constraint for \(\textbf{u}\) would be \(\Sigma_ku_k^2 < C\) for some positive value \(C\)

+

This is equivalent to choosing \((\beta, \textbf{u})\) to maximise (1) with \(D = diag(\textbf{0}, \textbf{1})\) where \(\textbf{0}\) denotes zero vector length \(p+1\) and \(\textbf{1}\) denotes the unit vector of length \(K\)

+

For a fixed value for \(\lambda\) this is equivalent to fitting the following generalized linear mixed model Ngo and Wand (2004)

\[ f(y|u) = exp\{ \phi^{-1} [y^T(X\beta + Zu) - c(X\beta + Zu)] + 1^Tc(y)\},\\ u \sim N(0, G) @@ -229,20 +210,47 @@

Generalized Linear Mixed Model
  • With similar notations as before and \(G = \sigma^2_uI_{K \times K}\)
  • -

    Thus \(Z\) is penalized by assuming the corresponding coefficients \(u\) are random effect with \(u \sim N(0, \sigma^2_uI)\).

    +

    Thus \(Z\) is penalized by assuming the corresponding coefficients \(\textbf{u}\) are random effect with \(\textbf{u} \sim N(\textbf{0}, \mathbf{\sigma}^2_u \textbf{I})\).

    +

    Refer to Chapter 8.2.2 of the book by Hens et al. (2012) for a more detailed explanation of the method.

    Fitting data

    To fit the data using the penalized likelihood framework, specify framework = "glmm"

     glmm <- parvob19_be_2001_2003 %>% 
    -  rename(status = seropositive) %>% 
    -  penalized_spline_model(s = "tp", framework = "glmm") 
    +  penalized_spline_model(status_col = "seropositive", s = "tp", framework = "glmm") 
     #> 
     #>  Maximum number of PQL iterations:  20
     #> iteration 1
     #> iteration 2
     #> iteration 3
     #> iteration 4
    -glmm$info$gam
    +glmm
    +#> Penalized spline model 
    +#> 
    +#> Input type:  linelisting 
    +#> Framework:  Mixed model 
    +#> $lme
    +#> Linear mixed-effects model fit by maximum likelihood
    +#>   Data: data 
    +#>   Log-likelihood: -6977.429
    +#>   Fixed: fixed 
    +#> X(Intercept)   Xs(age)Fx1 
    +#>    0.7122306    3.6123783 
    +#> 
    +#> Random effects:
    +#>  Formula: ~Xr - 1 | g
    +#>  Structure: pdIdnot
    +#>              Xr1      Xr2      Xr3      Xr4      Xr5      Xr6      Xr7      Xr8
    +#> StdDev: 6.020273 6.020273 6.020273 6.020273 6.020273 6.020273 6.020273 6.020273
    +#>         Residual
    +#> StdDev:        1
    +#> 
    +#> Variance function:
    +#>  Structure: fixed weights
    +#>  Formula: ~invwt 
    +#> Number of Observations: 3080
    +#> Number of Groups: 1 
    +#> 
    +#> $gam
     #> 
     #> Family: binomial 
     #> Link function: logit 
    @@ -251,17 +259,24 @@ 

    Generalized Linear Mixed Model #> spos ~ s(age, bs = s, sp = sp) #> #> Estimated degrees of freedom: -#> 6.45 total = 7.45

    +#> 6.45 total = 7.45 +#> +#> +#> attr(,"class") +#> [1] "gamm" "list"

     plot(glmm)

    -
    +
    Eilers, Paul H. C., and Brian D. Marx. 1996. “Flexible Smoothing with b-Splines and Penalties.” Statistical Science 11 (2). https://doi.org/10.1214/ss/1038425655.
    Green, P. J., and Bernard. W. Silverman. 1993. Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach. Chapman; Hall/CRC. https://doi.org/10.1201/b15710.
    +
    +Hens, Niel, Ziv Shkedy, Marc Aerts, Christel Faes, Pierre Van Damme, and Philippe Beutels. 2012. Modeling Infectious Disease Parameters Based on Serological and Social Contact Data: A Modern Statistical Perspective. Statistics for Biology and Health. Springer New York. https://doi.org/10.1007/978-1-4614-4072-7. +
    Ngo, Long, and Matthew P. Wand. 2004. “Smoothing with Mixed Model Software.” Journal of Statistical Software 9 (1). https://doi.org/10.18637/jss.v009.i01.
    @@ -280,28 +295,26 @@

    Generalized Linear Mixed Model

    -