From 14ccc4b41c3edd64e2c0f988f7e0c905eab41264 Mon Sep 17 00:00:00 2001 From: Morgan Kain Date: Thu, 5 Sep 2024 11:47:52 -0400 Subject: [PATCH] Updating some code for manuscript figures and supplemental code walkthrough --- analysis_walkthrough.Rmd | 165 +++++++++++++++++-------------- scr_R/comparae_cluster_methods.R | 3 +- scr_R/main_approach_coverage.R | 44 ++++++++- scr_R/main_approach_prop_sero.R | 20 +++- scr_R/scr_manuscript_figures.R | 139 +++++++++++++++++++------- 5 files changed, 258 insertions(+), 113 deletions(-) diff --git a/analysis_walkthrough.Rmd b/analysis_walkthrough.Rmd index 5ae1a79..57e5c4b 100644 --- a/analysis_walkthrough.Rmd +++ b/analysis_walkthrough.Rmd @@ -10,8 +10,10 @@ library(tidyverse) library(ggplot2) library(magrittr) library(cmdstanr) +library(formatR) register_knitr_engine(override = FALSE) knitr::opts_chunk$set(echo = TRUE) +knitr::opts_chunk$set(tidy = FALSE) build_stan_priors <- function(simulated_data, skew_fit, fit_attempt) { @@ -175,22 +177,20 @@ dataset3 <- sim.data %>% filter(param_set == 150) ### Overview -Here we present the analysis of three serological datasets. We analyze each with a single method well-suited for that dataset (the consequences of using a poorly-suited analysis are described in the main text). +Here we present the analysis of three serological datasets (3 of the 500 analyzed in the main text). These three datasets represent three alternatives in regards to percent population seroposivity, differences in mean MFI between seropositive and seronegative groups, and distributional skew. As described in the main text, seropositivity was simulated to be a function of 2 categorical covariates with two levels each (that could capture, for example, sex and species) and 1 continuous covariate (e.g., age in years). -The three datasets we focus on here are three of the 500 analyzed in the main text. These datasets were specifically chosen because they provide three different examples in terms of percent population seroposivity, the difference in mean MFI between seropositive and seronegative groups, and distributional skew. As described in the main text, seropositivity was simulated to be a function of 2 categorical covariates with two levels each (that could capture, for example, sex and species) and 1 continuous covariate (e.g., age in years). +We analyze each dataset here with a single method well-suited for that dataset (the consequences of using a poorly-suited analysis are summarized in the main text). Specifically, we consider both a two step extreme value + logistic regression approach as well as a single step Bayesian latent class regression (LCR) approach. Though we generally advocate against a two step approach when seeking inference about risk factors of seropositivity (what we call here an "epidemiological analysis"), an extreme value approach (e.g., 3sd method) is well suited for serostatus assignment and population-level seropositivity estimates when seropositivity is low (e.g., under 3%) and the positive MFI values have high variance--which together results in a unimodal right-skewed MFI distribution. Though using dichotomous seropositive assignments can result in over-confidence in logistic regression-estimated parameter values (too-narrow confidence intervals), attempting to estimate risk factors with very low seroposivity may nevertheless result in a lack of clear inference on risk factors (i.e., CI overlapping zero) because of low power. We show a two step extreme value approach with one of the three datasets here. -We conduct serostatus assignment and population-level seropositivity as well as derive inference about risk factors determining seropositivity using both a two step extreme value + logistic regression approach as well as a single step Bayesian latent class regression (LCR) approach. Though we generally advocate against a two step approach for epidemiological analysis (here an analysis of risk factors determining seropositivity), the 3sd method (i.e., step one of the two step approach) can be better suited for serostatus assignment and population-level seropositivity estimates when seropositivity is low (e.g., under 3%) and the positive MFI values have high variance, which together results in a unimodal right-skewed MFI distribution. Though using dichotomous seropositive assignments can result in over-confidence in logistic regression-estimated parameter values (too-narrow confidence intervals), attempting to estimate risk factors with very low seroposivity may nevertheless result in a lack of clear inference on risk factors (i.e., CI overlapping zero) because of low power. We illustrate this with one of the three datasets here. +We work through the analysis from raw data visualization through diagnostics prior to comparing estimated values and simulated "truth", keeping the true simulated coefficients hidden until this point. -We work through the analysis process from visualization through diagnostics and finally a comparison between estimated values and truth, keeping truth hidden from the reader until this point. +### Visualization -#### Visualization +For the analysis of any new serological (MFI: Median/Mean Fluorescence Intensity) dataset, the first step is a visual inspection of the MFI distribution to check for: a) bimodality, and b) if bimodality is present, evidence of right-skew in the right mode. -For the analysis of any new serological (MFI: Median/Mean Fluorescence Intensity) dataset, the first step should be a visual inspection of the MFI distribution to check for: a) bimodality, and b) if bimodality is present, evidence of right-skew in the right mode. +Given our findings (see main text) that false-negative and false-positive serostatus assignment error rates--and thus bias in population-level seropositivity estimates--are lower for the 3sd extreme value approach when MFI is analyzed on a linear scale but for clustering methods when MFI is analyzed on a log scale, we advise visualization of data on both scales. -Given our findings (see main text) that false-negative and false-positive serostatus assignment error rates, and thus bias in population-level seropositivity estimates, are lower for the 3sd extreme value when MFI is analyzed on a linear scale but for clustering methods when MFI is analyzed on a log scale, we advise visualization of of data on both scales. - -The three datasets we will analyze here are visualized in Figure 1. +Our three datasets chosen for analysis here are visualized in Figure 1. ```{r mfi_hist, echo=FALSE, warning=FALSE, message=FALSE, fig.cap = "Three example MFI datasets (columns), visualized on the linear scale (top row) and log scale (bottom row)."} gg.1 <- sim.data %>% filter(log_mfi == "mfi") %>% { @@ -224,23 +224,23 @@ gg.2 <- sim.data %>% filter(log_mfi == "mfi") %>% { gridExtra::grid.arrange(gg.1, gg.2, ncol = 1) ``` -A few things stand out from these visualizations. First, the dataset pictured in the left column has a long right tail on the linear scale but no strong bimodality when visualized on either scale. This observation hints that seropositivity is likely quite low in this dataset (or that seronegative and seropositive individuals have very similar MFI values; however, this is unlikely in any empirical dataset). Second, both of the datasets in the right two columns have clear bimodality, which is especially prominent on the log scale. Third, the right mode of the right-most dataset has some clear left skew due to the MFI values getting "stacked up" at the machine maximum. +A few things stand out from these visualizations. First, the dataset pictured in the left column has a long right tail on the linear scale and no strong bimodality when visualized on either scale. This observation hints that seropositivity is likely quite low in this dataset (or that seronegative and seropositive individuals have very similar MFI values; however, this is both unlikely in an empirical dataset and would be very difficult to analyze regardless: see Figure S19 in manuscript's supplemental material). Second, both of the datasets in the right two columns have clear bimodality, which is especially prominent on the log scale. Third, the right mode of the right-most dataset shows clear left skew due to the MFI values getting "stacked up" at the machine maximum. -#### An analysis plan +### An analysis plan -This visual inspection leads us to the following analysis plan for the three datasets for assigning serostatus. Dataset 1 (Figure 1, left column): 3sd with a cutoff derived from a robust mean and sd calculation using the R functions dplR::tbrm and jointseg::estimateSd, respectively; MFI analyzed on a linear scale. Estimation of population-level percent seropositivity using a binomial confidence interval. Analysis of risk factors of seropositivity using logistic regression with dichotomous serostatus assignments. +This visual inspection leads us to the following analysis plan for the three datasets for assigning serostatus (see the main text for further exposition). Dataset 1 (Figure 1, left column): 3sd with a cutoff derived from a robust mean and sd calculation using the R functions dplR::tbrm and jointseg::estimateSd, respectively; MFI analyzed on a linear scale. Estimation of population-level percent seropositivity using a binomial confidence interval. Analysis of risk factors of seropositivity using logistic regression with dichotomous serostatus assignments. Dataset 2 (Figure 1, center column): Bayesian LCR fitting normal distributions to log MFI. Estimation of population seroposivity and risk factors within the Bayesian model. Dataset 3 (Figure 1, right column): Bayesian LCR fitting a normal and a skewnormal distribution to log MFI. Estimation of population seroposivity and risk factors within the Bayesian model. -#### Analysis +### Analysis -##### Dataset 1 +### Dataset 1 -Assign serostatus based on a 3sd cutoff +A) Assign serostatus based on a 3sd cutoff ```{r set1_a, echo=TRUE, warning=FALSE, message=FALSE} dataset1.linear_mfi <- dataset1 %>% @@ -257,12 +257,14 @@ thresh <- mean1 + 3 * sd1 dataset1.linear_mfi %<>% mutate(seropositive = ifelse(mfi > thresh, 1, 0)) ``` -Estimate population level percent seropositivity +B) Estimate population level percent seropositivity ```{r set1_b, echo=TRUE, warning=FALSE, message=FALSE} -dataset1.pop_mod <- glm(seropositive ~ 1, family = "binomial", data = dataset1.linear_mfi) +dataset1.pop_mod <- glm(seropositive ~ 1, family = "binomial" + , data = dataset1.linear_mfi) -(dataset1.pop_pos <- c(confint(dataset1.pop_mod) %>% plogis(), coef(dataset1.pop_mod) %>% plogis()) %>% +(dataset1.pop_pos <- c(confint(dataset1.pop_mod) %>% plogis() + , coef(dataset1.pop_mod) %>% plogis()) %>% t() %>% as.data.frame() %>% mutate( @@ -276,10 +278,11 @@ dataset1.pop_mod <- glm(seropositive ~ 1, family = "binomial", data = dataset1. ``` -Predict impact of risk factors on seropositivity +C) Estimate the impacts of risk factors on seropositivity ```{r set1_c, echo=TRUE, warning=FALSE, message=FALSE} -dataset1.risk_mod <- glm(seropositive ~ cat1f + cat2f + con1f, family = "binomial", data = dataset1.linear_mfi) +dataset1.risk_mod <- glm(seropositive ~ cat1f + cat2f + con1f + , family = "binomial", data = dataset1.linear_mfi) (dataset1.risk_coef <- confint(dataset1.risk_mod) %>% as.data.frame() %>% @@ -301,13 +304,9 @@ dataset1.risk_mod <- glm(seropositive ~ cat1f + cat2f + con1f, family = "binom ) ``` -##### Dataset 2 - -Assign serostatus, predict population-level seropositivity and risk factors +### Dataset 2 -```{cmdstan, output.var="stan_fit2"} - -``` +A) Assign serostatus, predict population-level seropositivity and risk factors ```{r set2_a, echo=TRUE, warning=FALSE, message=FALSE} @@ -323,7 +322,9 @@ Assign serostatus, predict population-level seropositivity and risk factors ) ## Compile model - stan_model <- cmdstanr::cmdstan_model("stan_models/publication_model_normal_2.stan", pedantic = FALSE) + stan_model <- cmdstanr::cmdstan_model( + "stan_models/publication_model_normal_2.stan" + , pedantic = FALSE) ## Fit model. For model definition see ".stan" in the online supplemental material stan_fit2 <- try(R.utils::withTimeout(stan_model$sample( @@ -333,16 +334,26 @@ Assign serostatus, predict population-level seropositivity and risk factors , cat1f = dataset2.log$cat1f , cat2f = dataset2.log$cat2f , con1f = dataset2.log$con1f - , beta_base_prior_m = stan_priors$priors %>% filter(param == "beta_base_prior_m") %>% pull(prior) - , beta_base_prior_v = stan_priors$priors %>% filter(param == "beta_base_prior_v") %>% pull(prior) - , mu_base_prior_m = stan_priors$priors %>% filter(param == "mu_base_prior_m") %>% pull(prior) - , mu_diff_prior_m = stan_priors$priors %>% filter(param == "mu_diff_prior_m") %>% pull(prior) - , sigma_base_prior_m = stan_priors$priors %>% filter(param == "sigma_base_prior_m") %>% pull(prior) - , sigma_diff_prior_m = stan_priors$priors %>% filter(param == "sigma_diff_prior_m") %>% pull(prior) - , mu_base_prior_v = stan_priors$priors %>% filter(param == "mu_base_prior_v") %>% pull(prior) - , mu_pos_prior_v = stan_priors$priors %>% filter(param == "mu_pos_prior_v") %>% pull(prior) - , sigma_base_prior_v = stan_priors$priors %>% filter(param == "sigma_base_prior_v") %>% pull(prior) - , sigma_diff_prior_v = stan_priors$priors %>% filter(param == "sigma_diff_prior_v") %>% pull(prior) + , beta_base_prior_m = stan_priors$priors %>% + filter(param == "beta_base_prior_m") %>% pull(prior) + , beta_base_prior_v = stan_priors$priors %>% + filter(param == "beta_base_prior_v") %>% pull(prior) + , mu_base_prior_m = stan_priors$priors %>% + filter(param == "mu_base_prior_m") %>% pull(prior) + , mu_diff_prior_m = stan_priors$priors %>% + filter(param == "mu_diff_prior_m") %>% pull(prior) + , sigma_base_prior_m = stan_priors$priors %>% + filter(param == "sigma_base_prior_m") %>% pull(prior) + , sigma_diff_prior_m = stan_priors$priors %>% + filter(param == "sigma_diff_prior_m") %>% pull(prior) + , mu_base_prior_v = stan_priors$priors %>% + filter(param == "mu_base_prior_v") %>% pull(prior) + , mu_pos_prior_v = stan_priors$priors %>% + filter(param == "mu_pos_prior_v") %>% pull(prior) + , sigma_base_prior_v = stan_priors$priors %>% + filter(param == "sigma_base_prior_v") %>% pull(prior) + , sigma_diff_prior_v = stan_priors$priors %>% + filter(param == "sigma_diff_prior_v") %>% pull(prior) ) , init = stan_priors$starting_conditions , chains = 4 @@ -360,17 +371,14 @@ Assign serostatus, predict population-level seropositivity and risk factors stanfit <- rstan::read_stan_csv(stan_fit2$output_files()) samps.mat2 <- rstan::extract(stanfit, permuted = FALSE) samps.simp2 <- rstan::extract(stanfit, permuted = TRUE) -#stansummary2 <- summary(stanfit) ``` -Some quick diagnostics +B) Run some diagnostics (for more extensive diagnostics and a broader discussion of model convergence see: Bayesian Workflow. Gelman et al. 2020, arXiv and the Stan manual) ```{r set2_b, echo=TRUE, warning=FALSE, message=FALSE} stan_fit2$diagnostic_summary() -#paste("Max Rhat = ", stansummary2$summary %>% as.data.frame() %>% filter(!is.na(Rhat)) %>% pull(Rhat) %>% max() %>% round(4), sep = "") - ``` ```{r set2_c, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "Traceplot for baseline seropositivity (`Intercept`) parameter"} @@ -397,6 +405,8 @@ samps.mat2[ , , "beta_base"] %>% ``` +C) Summarize model output for downstream model comparisons + ```{r set2_e, echo=TRUE, warning=FALSE, message=FALSE} quants <- c(0.025, 0.975, 0.500) ## get CI for predictions of interest @@ -427,15 +437,11 @@ dataset2.pop_pos <- (samps.simp2$pop_sero/nrow(dataset2.log)) %>% ) ``` -##### Dataset 3 +### Dataset 3 -Assign serostatus, predict population-level seropositivity and risk factors +A) Assign serostatus, predict population-level seropositivity and risk factors -```{cmdstan, output.var="stan_fit3"} - -``` - -```{r set3_a, echo=TRUE, warning=FALSE, message=FALSE} +```{r set3_a, echo=TRUE, message=FALSE, warning=FALSE} ## just log mfi dataset3.log <- dataset3 %>% filter(log_mfi == "log_mfi") @@ -449,7 +455,9 @@ Assign serostatus, predict population-level seropositivity and risk factors ) ## Compile model - stan_model <- cmdstanr::cmdstan_model("stan_models/publication_model_skew_normal_wf_2.stan", pedantic = FALSE) + stan_model <- cmdstanr::cmdstan_model( + "stan_models/publication_model_skew_normal_wf_2.stan" + , pedantic = FALSE) ## Fit model. For model definition see ".stan" in the online supplemental material stan_fit3 <- try(R.utils::withTimeout(stan_model$sample( @@ -459,17 +467,29 @@ Assign serostatus, predict population-level seropositivity and risk factors , cat1f = dataset3.log$cat1f , cat2f = dataset3.log$cat2f , con1f = dataset3.log$con1f - , beta_base_prior_m = stan_priors$priors %>% filter(param == "beta_base_prior_m") %>% pull(prior) - , beta_base_prior_v = stan_priors$priors %>% filter(param == "beta_base_prior_v") %>% pull(prior) - , mu_diff_prior_m = stan_priors$priors %>% filter(param == "mu_diff_prior_m") %>% pull(prior) - , mu_diff_prior_v = stan_priors$priors %>% filter(param == "mu_diff_prior_v") %>% pull(prior) - , mu_pos_prior_m = stan_priors$priors %>% filter(param == "mu_pos_prior_m") %>% pull(prior) - , mu_pos_prior_v = stan_priors$priors %>% filter(param == "mu_pos_prior_v") %>% pull(prior) - , sigma_base_prior_m = stan_priors$priors %>% filter(param == "sigma_base_prior_m") %>% pull(prior) - , sigma_diff_prior_m = stan_priors$priors %>% filter(param == "sigma_diff_prior_m") %>% pull(prior) - , sigma_base_prior_v = stan_priors$priors %>% filter(param == "sigma_base_prior_v") %>% pull(prior) - , sigma_diff_prior_v = stan_priors$priors %>% filter(param == "sigma_diff_prior_v") %>% pull(prior) - , skew_pos_prior_m = min(dataset3.log %>% filter(group == 2) %>% pull(mfi) %>% moments::skewness(), -1) + , beta_base_prior_m = stan_priors$priors %>% + filter(param == "beta_base_prior_m") %>% pull(prior) + , beta_base_prior_v = stan_priors$priors %>% + filter(param == "beta_base_prior_v") %>% pull(prior) + , mu_diff_prior_m = stan_priors$priors %>% + filter(param == "mu_diff_prior_m") %>% pull(prior) + , mu_diff_prior_v = stan_priors$priors %>% + filter(param == "mu_diff_prior_v") %>% pull(prior) + , mu_pos_prior_m = stan_priors$priors %>% + filter(param == "mu_pos_prior_m") %>% pull(prior) + , mu_pos_prior_v = stan_priors$priors %>% + filter(param == "mu_pos_prior_v") %>% pull(prior) + , sigma_base_prior_m = stan_priors$priors %>% + filter(param == "sigma_base_prior_m") %>% pull(prior) + , sigma_diff_prior_m = stan_priors$priors %>% + filter(param == "sigma_diff_prior_m") %>% pull(prior) + , sigma_base_prior_v = stan_priors$priors %>% + filter(param == "sigma_base_prior_v") %>% pull(prior) + , sigma_diff_prior_v = stan_priors$priors %>% + filter(param == "sigma_diff_prior_v") %>% pull(prior) + , skew_pos_prior_m = min(dataset3.log %>% + filter(group == 2) %>% pull(mfi) %>% + moments::skewness(), -1) , skew_pos_prior_v = 3 , skew_neg_prior_m = 0 , skew_neg_prior_v = 0.5 @@ -490,18 +510,15 @@ Assign serostatus, predict population-level seropositivity and risk factors stanfit <- rstan::read_stan_csv(stan_fit3$output_files()) samps.mat3 <- rstan::extract(stanfit, permuted = FALSE) samps.simp3 <- rstan::extract(stanfit, permuted = TRUE) -#stansummary3 <- summary(stanfit) ``` -Some quick diagnostics +B) Some diagnostics for this second fit ```{r set3_b, echo=TRUE, warning=FALSE, message=FALSE} stan_fit3$diagnostic_summary() -#paste("Max Rhat = ", stansummary3$summary %>% as.data.frame() %>% filter(!is.na(Rhat)) %>% pull(Rhat) %>% max() %>% round(4), sep = "") - ``` ```{r set3_c, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "Traceplot for baseline seropositivity (`Intercept`) parameter"} @@ -528,6 +545,8 @@ samps.mat3[ , , "beta_base"] %>% ``` +C) Summarize model output for downstream model comparisons + ```{r set3_e, echo=TRUE, warning=FALSE, message=FALSE} ## get CI for predictions of interest quants <- c(0.025, 0.975, 0.500) @@ -560,14 +579,14 @@ dataset3.pop_pos <- (samps.simp3$pop_sero/nrow(dataset3.log)) %>% ``` -#### Comparison to truth +### Comparison to truth ```{r comb_a, echo=TRUE, warning=FALSE, message=FALSE} pop_pos_pred <- rbind(dataset1.pop_pos, dataset2.pop_pos, dataset3.pop_pos) risk_coef_pred <- rbind(dataset1.risk_coef, dataset2.risk_coef, dataset3.risk_coef) ``` -Individual seropositivity, dataset 1: 3sd. Five false-positives (true positive assigned negative; positive-negative), 13 false-negatives (negative-positive), and many correct negative-negative and positive-positive values. +Individual seropositivity, dataset 1, analyzed with a 3sd approach: visualized in Figure 6. This analysis resulted in five false-positives (true positive assigned negative; positive-negative), 13 false-negatives (negative-positive), and many correct negative-negative and positive-positive values. ```{r comb_b, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "Serostatus assignments. Dashed vertical line gives 3sd cutoff; values to the left are assigned a negative status, values to the right a positive status."} @@ -583,7 +602,7 @@ dataset1.linear_mfi %>% mutate(group = as.numeric(group) - 1) %>% ``` -Individual seropositivity, dataset 2: Normal-Normal LCR. Raw probability estimates (point values) vs truth shown in Figure X. With these clearly separated distributions the mixture model is highly accurate. +Individual seropositivity, dataset 2, analyzed with a Normal-Normal LCR approach: visualized in Figure 7. With these clearly separated distributions the mixture model is highly accurate; that is, few true seropositive individuals are given a high probability of being seronegative and vice versa. ```{r comb_c, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "Seropositive probability estimate (median) shown in blue; true serostatus shown as jittered black points (0 = seronegative, 1 = seropositive)."} @@ -600,9 +619,10 @@ dataset2.log %>% mutate( ``` -Individual seropositivity, dataset 3: Normal-Skew Normal LCR. Raw probability estimates (point values) vs truth shown in Figure X. The sizeable overlap in distributions and large left-skew in the seropositive MFI distribution leads to a number of false-negatives (see Figure X). +Individual seropositivity, dataset 3, analyzed with a Normal-Skew Normal LCR approach: visualized in Figure 8. Because of the sizeable overlap in distributions and large left-skew in the seropositive MFI distribution, a large number of false-negatives arise (i.e., small positive probabilities assigned to true seropositivies). We note however that analyzing such a dataset with a 3sd approach using a robust mean and sd calculation would lead to 100% of seropositives being assigned a seronegative status. + -```{r comb_d, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "Serostatus assignments. Dashed vertical line gives 3sd cutoff; values to the left are assigned a negative status, values to the right a positive status."} +```{r comb_d, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "Seropositive probability estimate (median) shown in blue; true serostatus shown as jittered black points (0 = seronegative, 1 = seropositive)."} dataset3.log %>% mutate( pred_prob = samps.simp3$membership_p[,1,] %>% colMeans() @@ -617,7 +637,7 @@ dataset3.log %>% mutate( ``` -Population seropositivity, all datasets. These two examples for population seropositivity for the two LCR models show high confidence (narrow CI) because of high confidence in the shapes of the two distributions. This leads to over-confidence (and thus CI that do not cover the true value) for the Normal-Skew Normal model. +Population seropositivity, all datasets. Estimates for the two LCR models show high confidence (narrow CI) because of high confidence in the shapes of the two distributions. This leads to over-confidence (and thus CI that do not cover the true value) for the Normal-Skew Normal model. ```{r comb_e, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "True vs estimated population-level seropositivity. Intervals are 95% CI."} @@ -633,7 +653,8 @@ data.frame( ) %>% mutate( dataset = paste("Dataset", dataset) , app_set = interaction(dataset, approach, sep = ": ") - , app_set = factor(app_set, levels = rev(unique(app_set))) + , app_set = factor(app_set + , levels = rev(unique(app_set))) ) %>% { ggplot(., aes(true_pos, app_set)) + geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0.2) + @@ -643,7 +664,7 @@ data.frame( ``` -Risk factors (regression coefficients). As expected from the low seropositivity, CI on most of the coefficients for the 3sd + logistic regression are very wide. As was found for a number of LCR models (see main text for more details), for the Normal-Skew Normal example fit here, one level of the categorical covariates (the "Intercept") was estimated to be too large and another too small. +Risk factors (regression coefficients). Here, because of the the low seropositivity in dataset 1, power is very low to estimate the effects of risk factors on seropositivity and thus CI on most of the coefficients for the 3sd + logistic regression are very wide. This is in contrast to using the two step 3sd + logistic regression approach in general. CI on these regression coefficients tend to be much too narrow (undercover) when seropositivity is higher (see main text). As was found for a number of other LCR model fits (see main text for more details), for the Normal-Skew Normal example fit here, one level of the categorical covariates (the "Intercept") was estimated to be too large and another too small. ```{r comb_f, echo=TRUE, warning=FALSE, message=FALSE, fig.cap = "True vs estimated risk factors (regression coefficients). Intervals are 95% CI."} risk_coef_pred %>% left_join(., sim.params.t, by = c("dataset", "coef")) %>% mutate( diff --git a/scr_R/comparae_cluster_methods.R b/scr_R/comparae_cluster_methods.R index 5e4e9f6..fa7174b 100644 --- a/scr_R/comparae_cluster_methods.R +++ b/scr_R/comparae_cluster_methods.R @@ -189,7 +189,8 @@ titer.out.just_groups %<>% rbind(., pred_pos) pop_seropos <- (samps$pop_sero / nrow(titer.out)) %>% quantile(c(0.025, 0.200, 0.500, 0.800, 0.975)) %>% - t() %>% as.data.frame() + t() %>% + as.data.frame() names(pop_seropos) <- c("lwr", "lwr_n", "mid", "upr_n", "upr") diff --git a/scr_R/main_approach_coverage.R b/scr_R/main_approach_coverage.R index d0e9ddf..3a60794 100644 --- a/scr_R/main_approach_coverage.R +++ b/scr_R/main_approach_coverage.R @@ -109,7 +109,13 @@ coef.cover <- best_fits %>% main_approach = factor(main_approach, levels = unique(main_approach)) ) -gg.1 <- coef.cover %>% mutate(name = plyr::mapvalues( +gg.1 <- coef.cover %>% + mutate(main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM")) %>% + mutate(main_approach = factor( + main_approach + , levels = c("3sd", "Bayesian LCR", "GMM")) + ) %>% + mutate(name = plyr::mapvalues( name , from = c("beta_base", "beta_cat1f_delta", "beta_cat2f_delta", "beta_con1f_delta") , to = c("Baseline (Intercept)", "Categorical Covariate 1" @@ -137,6 +143,7 @@ Parameter", values = c(16, 6, 2, 12)) + } gg.2 <- best_fits %>% + mutate(main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM")) %>% group_by(main_approach, param_set, name) %>% slice(1) %>% ungroup() %>% @@ -151,14 +158,14 @@ gg.2 <- best_fits %>% , "Categorical Covariate 2", "Continuous Covariate 1") )) %>% mutate(main_approach = factor( main_approach - , levels = c("3sd", "Mclust", "Bayesian LCR")) + , levels = c("3sd", "Bayesian LCR", "GMM")) ) %>% filter( m_diff > -20 ) %>% { ggplot(., aes(main_approach, m_diff)) + #geom_violin(aes(colour = main_approach, fill = main_approach)) + - geom_boxplot(aes(fill = main_approach), alpha = 0.4) + + geom_boxplot(aes(fill = main_approach), alpha = 0.4, linewidth = 0.3) + scale_shape_manual(name = "Regression Parameter", values = c(16, 6, 2, 12)) + scale_colour_manual(values = c( @@ -248,11 +255,38 @@ coef.cover %>% , from = c("beta_base", "beta_cat1f_delta", "beta_cat2f_delta", "beta_con1f_delta") , to = c("Baseline (Intercept)", "Categorical Covariate 1" , "Categorical Covariate 2", "Continuous Covariate 1") - )) %>% { + )) %>% + mutate(main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM")) %>% + mutate(main_approach = factor( + main_approach + , levels = c("3sd", "Bayesian LCR", "GMM")) + ) %>% mutate( + + m.s = plyr::mapvalues( + m.s + , from = c( + "Mclust -- (Unconstrained; sum of 2-n) + unweighted regression" + , "Mclust -- (Unconstrained; sum of 2-n) + probability weighted regression" + , "Mclust -- (2 group constrained) + unweighted regression" + , "Mclust -- (2 group constrained) + probability weighted regression" + , "Mclust -- (Unconstrained; BIC collapsed) + unweighted regression" + , "Mclust -- (Unconstrained; BIC collapsed) + probability weighted regression" + ) + , to = c( + "GMM -- (Unconstrained; sum of 2-n) + unweighted regression" + , "GMM -- (Unconstrained; sum of 2-n) + probability weighted regression" + , "GMM -- (2 group constrained) + unweighted regression" + , "GMM -- (2 group constrained) + probability weighted regression" + , "GMM -- (Unconstrained; BIC collapsed) + unweighted regression" + , "GMM -- (Unconstrained; BIC collapsed) + probability weighted regression" + ) + ) + + ) %>% { ggplot(., aes(perc_cov, m.s)) + geom_point(aes(colour = main_approach, shape = name), size = 3) + scale_colour_manual(values = c( - "#1b9e77" + "#1b9e77" , "#7570b3" , "#e6ab02" ), name = "Approach") + diff --git a/scr_R/main_approach_prop_sero.R b/scr_R/main_approach_prop_sero.R index 4e98ad0..314be87 100644 --- a/scr_R/main_approach_prop_sero.R +++ b/scr_R/main_approach_prop_sero.R @@ -176,6 +176,21 @@ all_coefs_for_gg.2 %<>% ungroup() %>% mutate( ) gg.1 <- param_coverage_by_model.all.adj %>% + mutate( + main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM") + , m.s = plyr::mapvalues( + m.s + , from = c( + "Mclust (Unconstrained; sum of 2-n)" + , "Mclust (Unconstrained; BIC collapsed)" + , "Mclust (2 group constrained)" + ) + , to = c( + "GMM (Unconstrained; sum of 2-n)" + , "GMM (Unconstrained; BIC collapsed)" + , "GMM (2 group constrained)" + ) + )) %>% rename(Approach = m.s) %>% { ggplot(., aes(cover, Approach)) + geom_point(aes(colour = main_approach), size = 3) + @@ -194,14 +209,15 @@ gg.1 <- param_coverage_by_model.all.adj %>% , legend.position = "none" ) + scale_x_continuous( - breaks = c(0.1, 0.35, 0.60, 0.85, 0.95) - , labels = c("10%", "35%", "60%", "85%", "95%")) + + breaks = c(0.1, 0.35, 0.60, 0.80, 0.95) + , labels = c("10%", "35%", "60%", "80%", "95%")) + xlab("Coverage") + ylab("Model") + geom_vline(xintercept = 0.95, linetype = "dashed") } gg.2 <- all_coefs_for_gg.2 %>% + mutate(main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM")) %>% filter(quantile == "mid") %>% mutate(m.s = factor(m.s, levels = levvs)) %>% rename(Approach = m.s) %>% { diff --git a/scr_R/scr_manuscript_figures.R b/scr_R/scr_manuscript_figures.R index 04cc651..b68fb0f 100644 --- a/scr_R/scr_manuscript_figures.R +++ b/scr_R/scr_manuscript_figures.R @@ -30,7 +30,6 @@ all.out$coefficient_ests <- all.out$coefficient_ests %>% all.out$coverage <- all.out$coverage %>% left_join(., sim.data.summaries.j) - ##### ## Figure 1: 3sd approaches, log and linear ##### @@ -187,7 +186,7 @@ all.out$pop_seropositivity %>% , "publication_model_skew_normal_wf_2.stan.Bayesian LCR" ) , to = c( - "Mclust -- (Unconstrained -- BIC collapsed)" + "GMM -- (Unconstrained -- BIC collapsed)" , "Bayesian LCR -- normal-normal" , "Bayesian LCR -- normal-skewnormal" )) @@ -465,10 +464,28 @@ gg.2 <- all_compare1 %>% { gridExtra::grid.arrange(gg.1, gg.2, ncol = 1) +## !!NOTE: Need to run main_approach_prop_sero.R to build "param_coverage_by_model.all.adj" for + ## use in this figure + gg.1 <- all_compare1 %>% mutate(Approach = factor(Approach, levels = c( unique(param_coverage_by_model.all.adj$m.s) - ))) %>% { + ))) %>% + mutate( + main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM") + , Approach = plyr::mapvalues( + Approach + , from = c( + "Mclust (Unconstrained; sum of 2-n)" + , "Mclust (Unconstrained; BIC collapsed)" + , "Mclust (2 group constrained)" + ) + , to = c( + "GMM (Unconstrained; sum of 2-n)" + , "GMM (Unconstrained; BIC collapsed)" + , "GMM (2 group constrained)" + ) + )) %>% { ggplot(., aes(tot_cov, Approach)) + geom_point(aes(colour = main_approach, shape = log_mfi), size = 3) + ylab("Model") + @@ -486,8 +503,8 @@ gg.1 <- all_compare1 %>% , legend.position = "none" ) + scale_x_continuous( - breaks = c(0.1, 0.35, 0.60, 0.85, 0.95) - , labels = c("10%", "35%", "60%", "85%", "95%") + breaks = c(0.1, 0.35, 0.60, 0.80, 0.95) + , labels = c("10%", "35%", "60%", "80%", "95%") , lim = c(0, 1)) + scale_shape_discrete(name = "MFI Scale") + geom_vline(xintercept = 0.95, linetype = "dashed") @@ -521,6 +538,21 @@ gg.2 <- all_compare1 %>% gg.2 <- all_coefs_for_gg.2 %>% mutate(m.s = factor(m.s, levels = levvs)) %>% rename(Approach = m.s) %>% + mutate( + main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM") + , Approach = plyr::mapvalues( + Approach + , from = c( + "Mclust (Unconstrained; sum of 2-n)" + , "Mclust (Unconstrained; BIC collapsed)" + , "Mclust (2 group constrained)" + ) + , to = c( + "GMM (Unconstrained; sum of 2-n)" + , "GMM (Unconstrained; BIC collapsed)" + , "GMM (2 group constrained)" + ) + )) %>% filter(quantile == "mid") %>% filter(true < 0.05, prop_overlap_quant < 0.2) %>% { ggplot(., aes(prop_pos_diff, Approach)) + @@ -551,6 +583,8 @@ gg.2 <- all_coefs_for_gg.2 %>% gridExtra::grid.arrange(gg.1, gg.2, ncol = 2) +## !!NOTE: Need to run main_approach_coverage.R first + gg.1 <- all_compare2 %>% { ggplot(., aes(Approach, tot_cov)) + geom_point(aes(colour = Approach, shape = log_mfi), size = 3) + @@ -641,6 +675,20 @@ parm_sets <- all_coefs_for_gg.2 %>% pull(param_set) %>% unique() +to_vals <- c( + "3sd -- (Approximated control)" + , "3sd -- (Robust Mean and SD)" + , "Bayesian LCR -- lognormal-lognormal" + , "Bayesian LCR -- normal-normal" + , "Bayesian LCR -- normal-skewnormal" + , "GMM -- (2 group constrained) + unweighted regression" + , "GMM -- (2 group constrained) + probability weighted regression" + , "GMM -- (Unconstrained; sum of 2-n) + unweighted regression" + , "GMM -- (Unconstrained; sum of 2-n) + probability weighted regression" + , "GMM -- (Unconstrained; BIC collapsed) + unweighted regression" + , "GMM -- (Unconstrained; BIC collapsed) + probability weighted regression" +) + these_to_vals <- to_vals[c( 4, 5, 3, 6, 7 , 8, 9, 10, 11 @@ -651,11 +699,11 @@ gg.2 <- all_coef_ests %>% filter(param_set %in% parm_sets) %>% mutate( Approach = plyr::mapvalues(mod.meth, from = unique(mod.meth), to = these_to_vals) - ) %>% mutate( - Approach = factor(Approach, levels = c( - coef.cover$m.s %>% levels() - )) - ) %>% filter( + ) %>% + mutate( + Approach = factor(Approach, levels = c( + these_to_vals[c(10, 11, 6, 7, 4, 5, 8, 9, 2, 3, 1)] + ))) %>% filter( m_diff > -5, m_diff < 5 ) %>% { ggplot(., aes(Approach, m_diff)) + @@ -673,8 +721,8 @@ gg.2 <- all_coef_ests %>% theme( strip.text.x = element_blank() , plot.margin = unit(c(0.0,0.2,0.2,0.15), "cm") - # , axis.text.x = element_text(size = 8, angle = 330, hjust = 0, vjust = 0.5) - , axis.text.x = element_blank() + , axis.text.x = element_text(size = 8, angle = 330, hjust = 0, vjust = 0.5) + # , axis.text.x = element_blank() , axis.text.y = element_text(size = 10) , legend.position = "none" ) + @@ -714,28 +762,35 @@ all_coefs_for_gg.2 %>% ## See main_approach_coverage for required code ci_wid_by_model %>% + left_join(., sim.params %>% dplyr::select(param_set, n_samps)) %>% mutate(m.s = interaction(main_approach, sub_approach, sep = " -- ")) %>% mutate(m.s = plyr::mapvalues(m.s, from = c( - "Mclust -- mclust (2 group constrained) + unweighted regression" - , "Mclust -- mclust (2 group constrained) + probability weighted regression" - , "Mclust -- mclust (Unconstrained) + unweighted regression" - , "Mclust -- mclust (Unconstrained) + probability weighted regression" - , "Mclust -- mclust (Unconstrained reduced) + unweighted regression" - , "Mclust -- mclust (Unconstrained reduced) + probability weighted regression" + "Mclust -- (2 group constrained) + unweighted regression" + , "Mclust -- (2 group constrained) + probability weighted regression" + , "Mclust -- (Unconstrained; BIC collapsed) + unweighted regression" + , "Mclust -- (Unconstrained; BIC collapsed) + probability weighted regression" + , "Mclust -- (Unconstrained; sum of 2-n) + unweighted regression" + , "Mclust -- (Unconstrained; sum of 2-n) + probability weighted regression" + , "3sd -- (Approximated control)" + , "3sd -- (Robust Mean and SD)" ), to = c( - "Mclust -- mclust (2 group constrained) + + "GMM -- 2 group constrained + unweighted regression" - , "Mclust -- mclust (2 group constrained) + + , "GMM -- 2 group constrained + probability weighted regression" - , "Mclust -- mclust (Unconstrained) + + , "GMM -- Unconstrained; BIC collapsed + unweighted regression" - , "Mclust -- mclust (Unconstrained) + + , "GMM -- Unconstrained; BIC collapsed + probability weighted regression" - , "Mclust -- mclust (Unconstrained reduced) + + , "GMM -- Unconstrained; sum of 2-n + unweighted regression" - , "Mclust -- mclust (Unconstrained reduced) + + , "GMM -- Unconstrained; sum of 2-n + probability weighted regression" + , "3sd -- Approximated control" + , "3sd -- Robust Mean and SD" ))) %>% + mutate(main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM")) %>% + mutate(main_approach = factor(main_approach, levels = c("3sd", "Bayesian LCR", "GMM"))) %>% mutate(log_mfi = plyr::mapvalues(log_mfi, from = c( "log_mfi", "mfi" ), to = c("Log MFI", "Linear MFI"))) %>% @@ -747,8 +802,14 @@ probability weighted regression" filter(ci_wid_mean < 200) %>% { ggplot(., aes(n_samps, ci_wid_mean)) + geom_point(aes(colour = `Main Approach`, shape = `MFI Scale`), alpha = 0.5) + + scale_colour_manual(values = c( + "#1b9e77" + , "#7570b3" + , "#e6ab02" + ), name = "Approach") + theme( strip.text.x = element_text(size = 9) + , axis.text.y = element_text(size = 11) , plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm") ) + scale_y_log10( @@ -776,26 +837,32 @@ param_coverage_by_model %>% ungroup() %>% mutate(m.s = interaction(main_approach, sub_approach, sep = " -- ")) %>% mutate(m.s = plyr::mapvalues(m.s, from = c( - "Mclust -- mclust (2 group constrained) + unweighted regression" + "Mclust -- mclust (2 group constrained) + unweighted regression" , "Mclust -- mclust (2 group constrained) + probability weighted regression" - , "Mclust -- mclust (Unconstrained) + unweighted regression" - , "Mclust -- mclust (Unconstrained) + probability weighted regression" , "Mclust -- mclust (Unconstrained reduced) + unweighted regression" , "Mclust -- mclust (Unconstrained reduced) + probability weighted regression" + , "Mclust -- mclust (Unconstrained) + unweighted regression" + , "Mclust -- mclust (Unconstrained) + probability weighted regression" + , "3sd -- 3sd (Approximated control)" + , "3sd -- 3sd (Robust Mean and SD)" ), to = c( - "Mclust -- mclust (2 group constrained) + + "GMM -- 2 group constrained + unweighted regression" - , "Mclust -- mclust (2 group constrained) + + , "GMM -- 2 group constrained + probability weighted regression" - , "Mclust -- mclust (Unconstrained) + + , "GMM -- Unconstrained; BIC collapsed + unweighted regression" - , "Mclust -- mclust (Unconstrained) + + , "GMM -- Unconstrained; BIC collapsed + probability weighted regression" - , "Mclust -- mclust (Unconstrained reduced) + + , "GMM -- Unconstrained; sum of 2-n + unweighted regression" - , "Mclust -- mclust (Unconstrained reduced) + + , "GMM -- Unconstrained; sum of 2-n + probability weighted regression" + , "3sd -- Approximated control" + , "3sd -- Robust Mean and SD" ))) %>% + mutate(main_approach = plyr::mapvalues(main_approach, from = "Mclust", to = "GMM")) %>% + mutate(main_approach = factor(main_approach, levels = c("3sd", "Bayesian LCR", "GMM"))) %>% mutate(log_mfi = plyr::mapvalues(log_mfi, from = c( "log_mfi", "mfi" ), to = c("Log MFI", "Linear MFI"))) %>% @@ -807,8 +874,14 @@ probability weighted regression" mutate(m_cov = m_cov / 4) %>% { ggplot(., aes(n_samps, m_cov)) + geom_point(aes(colour = `Main Approach`, shape = `MFI Scale`)) + + scale_colour_manual(values = c( + "#1b9e77" + , "#7570b3" + , "#e6ab02" + ), name = "Approach") + theme( strip.text.x = element_text(size = 9) + , axis.text.y = element_text(size = 11) , plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm") ) + xlab("Sample Size") +