Corrected transposed matrix in reverse GP function

R/create_posterior_func.R

@@ -126,3 +126,35 @@ create_post_func_seeirr_combined <- function(parTab, data, ts, INCIDENCE_FUNC=de
   }
   f
 }
+
+
+#' @export
+create_posterior_func_detectable <- function(parTab,
+                                  data,
+                                  PRIOR_FUNC=NULL,
+                                  INCIDENCE_FUNC=NULL,
+                                  solve_likelihood=TRUE,
+                                  ...) {
+  par_names <- parTab$names
+  pars <- parTab$values
+  names(pars) <- par_names
+  times <- 0:max(data$t)
+  ages <- 1:max(data$t)
+  obs_times <- unique(data$t)
+
+  f <- function(pars){
+    names(pars) <- par_names
+    prob_infection_tmp <- INCIDENCE_FUNC(pars, times)
+    lik <- 0
+    if(solve_likelihood){
+      lik <- sum(likelihood_detectable(data, ages, pars, prob_infection_tmp))
+    }
+
+    if(!is.null(PRIOR_FUNC)){
+      prior <- PRIOR_FUNC(pars, ...)
+      lik <- lik+prior
+    }
+    return(lik)
+  }
+  f
+}

R/incidence_functions.R

@@ -101,7 +101,8 @@ reverse_gp_model <- function(desired_probs, pars, times){
   rho <- pars["rho"]
   K <- nu^2 * exp(-rho^2 * t_dist^2)
   diag(K) <- diag(K) + 0.01
-  L_K <- chol(K)
+  L_K <- t(chol(K))
+
 
   ps <- sum(desired_probs)*desired_probs + 0.000000001
   k1 <- -log((1/ps) -1)

R/likelihoods.R

@@ -48,6 +48,44 @@ likelihood_cpp_wrapper <- function(obs_dat, ages, times,
   liks_tj
 }
 
+#' @export
+prob_detectable_curve <- function(pars, ages){
+  mean <- pars["mean"]
+  var <- pars["var"]
+
+  scale <- var/mean
+  shape <- mean/scale
+
+  probs <- dgamma(ages,shape=shape,scale=scale)*pars["c"]
+  probs[probs > 1] <- 1
+
+  return(probs)
+
+}
+
+#' Likelihood for using proportion detectable only
+#'
+#' @export
+likelihood_detectable <- function(obs_dat, ages, pars, prob_infection){
+  liks <- 0
+  ## For each sampling times, calculate likelihood for Ct values measured at that time
+  for(i in 1:nrow(obs_dat)){
+    ## Only solve back until the earliest possible infection time
+    obs_time <- obs_dat$t[i]
+    ages1 <- ages[(obs_time - ages) > 0]
+    n_pos <- obs_dat$pos[i]
+    n_neg <- obs_dat$neg[i]
+    sum_detect_inf <- sum(vapply(ages1,
+                                 FUN=function(a) prob_infection[obs_time-a]*prob_detectable_curve(pars, a),
+                                 FUN.VALUE=numeric(1)))
+    liks_detect <- n_pos*log(sum_detect_inf)
+    liks_undetect <- n_neg*log(1-sum_detect_inf)
+    liks <- liks +  liks_detect+liks_undetect
+
+  }
+  liks
+}
+
 
 #' Function to give probability of observing x given age a and the viral kinetics curve
 #' @export

R/viral_load_functions.R

@@ -18,7 +18,6 @@ viral_load_func <- function(pars, obs_t, convert_vl=FALSE, infection_time=0){
   t_period2 <- obs_t > tshift & obs_t <= (desired_mode + tshift)
   t_period3 <- obs_t > (desired_mode + tshift) & obs_t <= (desired_mode + tshift + t_switch)
   t_period4 <- obs_t > (desired_mode + tshift + t_switch)
-
   y <- numeric(length(obs_t))
   y[t_period1] <- true_0
   y[t_period2] <- growth_rate * (obs_t[t_period2] - tshift) + true_0