Added ne wvignettes and COP modelling

Author: dchodge

R/SeroCOP.R

@@ -243,15 +243,54 @@ SeroCOP <- R6::R6Class(
         predictions <- matrix(NA, nrow = n_iter, ncol = n_new)
 
         for (i in 1:n_iter) {
-          predictions[i, ] <- params$floor[i] + 
-            (params$ceiling[i] - params$floor[i]) / 
-            (1 + exp(params$slope[i] * (newdata - params$ec50[i])))
+          # Using the new formula: ceiling * (floor + (1-floor) * inv_logit(-slope * (titre - ec50)))
+          logit_part <- 1 / (1 + exp(params$slope[i] * (newdata - params$ec50[i])))
+          predictions[i, ] <- params$ceiling[i] * (params$floor[i] + (1 - params$floor[i]) * logit_part)
         }
 
         return(predictions)
       }
     },
 
+    #' @description
+    #' Extract probability of protection from the fitted model
+    #' @param newdata Optional vector of new titre values for prediction
+    #' @return Matrix of protection probabilities (rows = MCMC samples, cols = observations)
+    #' @examples
+    #' sero <- SeroCOP$new()
+    #' sero$fit_model()
+    #' protection <- sero$predict_protection()
+    predict_protection = function(newdata = NULL) {
+      if (is.null(self$fit)) {
+        stop("Model has not been fitted yet. Run fit_model() first.")
+      }
+      
+      if (is.null(newdata)) {
+        # Extract fitted protection probabilities
+        prob_protection <- rstan::extract(self$fit, pars = "prob_protection")[[1]]
+        return(prob_protection)
+      } else {
+        # Predict protection for new data
+        # Get infection probabilities first
+        prob_infection <- self$predict(newdata = newdata)
+        
+        # Extract ceiling samples
+        params <- rstan::extract(self$fit)
+        ceiling_samples <- params$ceiling
+        
+        # Calculate protection: 1 - (prob_infection / ceiling)
+        n_iter <- nrow(prob_infection)
+        n_new <- ncol(prob_infection)
+        prob_protection <- matrix(NA, nrow = n_iter, ncol = n_new)
+        
+        for (i in 1:n_iter) {
+          prob_protection[i, ] <- 1 - (prob_infection[i, ] / ceiling_samples[i])
+        }
+        
+        return(prob_protection)
+      }
+    },
+    
     #' @description
     #' Get summary statistics for model parameters
     #' @return Data frame with parameter summaries
@@ -482,6 +521,29 @@ SeroCOP <- R6::R6Class(
         )
 
       return(p)
+    },
+    
+    #' @description
+    #' Extract the correlate of protection conditional on exposure.
+    #' @param correlate_of_risk Numeric vector of correlates of risk.
+    #' @param upper_bound Numeric value for the upper bound (default: 0.7).
+    #' @return Numeric vector of correlates of protection.
+    #' @examples
+    #' sero <- SeroCOP$new(titre = titre, infected = infected)
+    #' sero$fit_model()
+    #' cor <- sero$extract_cop(correlate_of_risk = c(0.1, 0.2), upper_bound = 0.7)
+    extract_cop = function(correlate_of_risk, upper_bound = 0.7) {
+      if (missing(correlate_of_risk)) {
+        stop("correlate_of_risk must be provided")
+      }
+      if (!is.numeric(correlate_of_risk) || any(correlate_of_risk < 0 | correlate_of_risk > 1)) {
+        stop("correlate_of_risk must be a numeric vector with values between 0 and 1")
+      }
+      if (!is.numeric(upper_bound) || upper_bound <= 0) {
+        stop("upper_bound must be a positive numeric value")
+      }
+      cop <- (1 - correlate_of_risk) / upper_bound
+      return(cop)
     }
   ),
 

R/SeroCOPMulti.R

@@ -300,6 +300,33 @@ SeroCOPMulti <- R6::R6Class(
         )
 
       return(p)
+    },
+    
+    #' @description
+    #' Extract the correlate of protection conditional on exposure for all biomarkers.
+    #' @param correlate_of_risk_list List of numeric vectors of correlates of risk for each biomarker.
+    #' @param upper_bound Numeric value for the upper bound (default: 0.7).
+    #' @return List of numeric vectors of correlates of protection for each biomarker.
+    #' @examples
+    #' multi_model <- SeroCOPMulti$new(titre = titre_matrix, infected = infected)
+    #' multi_model$fit_all()
+    #' cor_list <- list(IgG = c(0.1, 0.2), IgA = c(0.3, 0.4))
+    #' cop_list <- multi_model$extract_cop_multi(correlate_of_risk_list = cor_list, upper_bound = 0.7)
+    extract_cop_multi = function(correlate_of_risk_list, upper_bound = 0.7) {
+      if (!is.list(correlate_of_risk_list)) {
+        stop("correlate_of_risk_list must be a list of numeric vectors")
+      }
+      if (!is.numeric(upper_bound) || upper_bound <= 0) {
+        stop("upper_bound must be a positive numeric value")
+      }
+      cop_list <- lapply(correlate_of_risk_list, function(cor) {
+        if (!is.numeric(cor) || any(cor < 0 | cor > 1)) {
+          stop("Each correlate_of_risk must be a numeric vector with values between 0 and 1")
+        }
+        1 - (cor / upper_bound)
+      })
+      names(cop_list) <- names(correlate_of_risk_list)
+      return(cop_list)
     }
   )
 )

inst/stan/logistic_model.stan

@@ -1,10 +1,10 @@
 // Four-parameter logistic model for correlates of protection
 // Models the probability of infection as a function of antibody titre
 // Parameters:
-//   - floor: lower asymptote (minimum infection probability)
-//   - ceiling: upper asymptote (maximum infection probability)
-//   - ec50: titre at 50% between floor and ceiling (inflection point)
-//   - slope: steepness of the curve
+//   - floor: proportion of maximum risk remaining at high titre (relative protection at high titre)
+//   - ceiling: maximum infection probability (at low antibody titre)
+//   - ec50: titre at inflection point (50% reduction from ceiling to ceiling*floor)
+//   - slope: steepness of the protective curve (higher = steeper decline in risk with titre)
 
 data {
   int<lower=0> N;                    // number of observations
@@ -29,18 +29,21 @@ data {
 }
 
 parameters {
-  real<lower=0,upper=1> floor;       // lower asymptote
-  real<lower=0,upper=1> ceiling;     // upper asymptote
-  real ec50;                         // titre at midpoint
-  real<lower=0> slope;               // slope parameter (steepness)
+  real<lower=0,upper=1> floor;       // proportion of maximum risk at high titre
+  real<lower=0,upper=1> ceiling;     // maximum infection probability at low titre
+  real ec50;                         // titre at inflection point
+  real<lower=0> slope;               // steepness of protective curve
 }
 
 transformed parameters {
   vector[N] prob_infection;
-  
+  vector[N] prob_protection;
   // Four-parameter logistic function
+  // At low titre: prob → ceiling (maximum risk)
+  // At high titre: prob → ceiling * floor (minimum risk, as proportion of ceiling)
   for (n in 1:N) {
-    prob_infection[n] = floor + (ceiling - floor) / (1 + exp(slope * (titre[n] - ec50)));
+    prob_infection[n] = ceiling * (inv_logit(-slope * (titre[n] - ec50)) * (1 - floor) + floor);
+    prob_protection[n] = 1 - (prob_infection[n] / ceiling);
   }
 }