Separated integration functions, fixed off-by-one error

DESCRIPTION

@@ -33,6 +33,7 @@ Collate:
     'LoessTimeDependentModel.R'
     'TimeResolvedMR-package.R'
     'calculate_genetic_effects.R'
+    'integrate_estimates.R'
     'perform_MR.R'
     'time_dependent_aalen.R'
     'time_dependent_glm.R'

R/integrate_estimates.R

@@ -0,0 +1,84 @@
+#' Integrate time-dependent estimates based on midpoint rule and compute 95% CIs
+#'
+#' @param age The ages at which to integrate
+#' @param wald_ratio The calculated Wald ratios
+#' @param exposure_effect The genetic effect on the exposure at every age
+#' @param outcome_variance The variance of the genetic effect on the outcome
+#'   at every age
+#'
+#' @return A data.frame with age, effect (Gamma), variance and 95% CIs
+#' @keywords internal
+integrate_midpoint <- function(age, wald_ratio, exposure_effect, outcome_variance) {
+  Gamma <- c(0, cumsum(diff(age) * wald_ratio))
+  Gamma_variance <- outcome_variance / exposure_effect ^ 2
+  data.frame(
+    age = age, Gamma = Gamma, Variance = Gamma_variance,
+    L95 = Gamma - 1.96 * sqrt(Gamma_variance),
+    H95 = Gamma + 1.96 * sqrt(Gamma_variance)
+  )
+}
+
+#' Integrate time-dependent effects using the trapezoidal rule and compute 95%
+#' CIs
+#'
+#' @param age The ages at which to integrate
+#' @param wald_ratio The calculated Wald ratios
+#' @param exposure_effect The genetic effect on the exposure at every age
+#' @param outcome_variance The variance of the genetic effect on the outcome
+#'   at every age
+#'
+#' @return A data.frame with age, effect (Gamma), variance and 95% CIs
+#' @keywords internal
+integrate_trapezoid <- function(age, wald_ratio, exposure_effect, outcome_variance) {
+  # Add (0,0) if it isn't there yet
+  if (age[1] != 0) {
+    age <- c(0, age)
+    wald_ratio <- c(0, wald_ratio)
+  }
+  Gamma <- pracma::cumtrapz(age, wald_ratio)
+  Gamma_variance <- calculate_mr_variance_trapz(
+    effect = exposure_effect,
+    variance = outcome_variance
+  )
+  data.frame(
+    age = age, Gamma = Gamma, Variance = Gamma_variance,
+    L95 = Gamma - 1.96 * sqrt(Gamma_variance),
+    H95 = Gamma + 1.96 * sqrt(Gamma_variance)
+  )
+}
+#' Calculate time-dependent effect variance of MR effects estimated by
+#' integration using the trapezoidal rules
+#'
+#' @param effect Numeric vector of total genetic effects on exposure for each
+#'   time point
+#' @param variance Numeric vector of variance of genetic effect on outcome for
+#'   each time point.
+#'
+#' @return A numeric vector of variance at each time point
+#' @keywords internal
+calculate_mr_variance_trapz <- function(effect, variance) {
+  stopifnot(
+    is.numeric(effect),
+    is.numeric(variance),
+    length(effect) == length(variance),
+    length(effect) > 0
+  )
+  n <- length(effect)
+  var1 <- variance[1] / (16 * effect[1] ^ 2)
+  if (n < 2) return(var1)
+
+  var2 <- var1 + variance[2] / (4 * effect[2] ^ 2)
+  if (n < 3) return (c(var1,var2))
+
+  # At this point, we know that n > 2, so the following is safe:
+  point_variances <- c(
+    var1,
+    vapply(2:n, \(k){ variance[k] / (4 * effect[k] ^2) }, numeric(1))
+  )
+
+  c(
+    var1,
+    var2,
+    vapply(3:n, \(k) { var1 + point_variances[k-1] + point_variances[k] }, numeric(1))
+  )
+}

R/perform_MR.R

@@ -1,63 +1,3 @@
-#' Numerically integrate effect estimates by age
-#'
-#' @param age Numeric vector containing ages
-#' @param beta Numeric vector containing effect estimates
-#' @param method The integration method to use. Either `trapezoidal` (default)
-#'   or `midpoint`.
-#'
-#' @return A numeric matrix with columns age and cumulative integral
-#' @keywords internal
-integrate_estimates <- function(age, beta, method = c("trapezoidal", "midpoint")){
-  stopifnot(
-    is.numeric(age),
-    is.numeric(beta)
-  )
-  method <- match.arg(method)
-  effects <- switch(method,
-    trapezoidal = pracma::cumtrapz(age, beta),
-    midpoint =  c(0, cumsum(diff(age) * zoo::rollmean(beta, 2))),
-    "You should not be here"
-  )
-  effects
-}
-
-#' Calculate time-dependent effect variance of MR effects estimated by
-#' integration using the trapezoidal rules
-#'
-#' @param effect Numeric vector of total genetic effects on exposure for each
-#'   time point
-#' @param variance Numeric vector of variance of genetic effect on outcome for
-#'   each time point.
-#'
-#' @return A numeric vector of variance at each time point
-#' @keywords internal
-calculate_mr_variance_trapz <- function(effect, variance) {
-  stopifnot(
-    is.numeric(effect),
-    is.numeric(variance),
-    length(effect) == length(variance),
-    length(effect) > 0
-  )
-  n <- length(effect)
-  var1 <- variance[1] / (16 * effect[1] ^ 2)
-  if (n < 2) return(var1)
-
-  var2 <- var1 + variance[2] / (4 * effect[2] ^ 2)
-  if (n < 3) return (c(var1,var2))
-
-  # At this point, we know that n > 2, so the following is safe:
-  point_variances <- c(
-    var1,
-    vapply(2:n, \(k){ variance[k] / (4 * effect[k] ^2) }, numeric(1))
-  )
-
-  c(
-    var1,
-    var2,
-    vapply(3:n, \(k) { var1 + point_variances[k-1] + point_variances[k] }, numeric(1))
-  )
-}
-
 #' Perform time-dependent MR analysis
 #'
 #' @param age_seq Ages at which to estimate exposure-outcome effect
@@ -116,28 +56,28 @@ time_dependent_MR <- function(age_seq, exposure_model, outcome_model,
     rlang::abort(paste0("age_seq is an invalid range: ", age_seq))
   }
 
-  midpoints <-  zoo::rollmean(age_seq, 2)
+  midpoints <- zoo::rollmean(age_seq, 2)
+
   total_exposure_effects <- totalEffect(exposure_model, age_seq)
   midpoint_exposure_effects <- totalEffect(exposure_model, midpoints)
+
   outcome_variances <- variance(outcome_model, age_seq)
   outcome_effects <- totalEffect(outcome_model, age_seq)
   dBeta <- diff(outcome_effects) / diff(age_seq)
+
   wald_ratio <- dBeta / midpoint_exposure_effects
-  Gamma <- c(0, integrate_estimates(midpoints, wald_ratio))
-  Gamma_variance <- switch(method,
-    trapezoidal = calculate_mr_variance_trapz(
-      effect = total_exposure_effects,
-      variance = outcome_variances
+
+  switch(method,
+    trapezoidal = integrate_trapezoid(
+      age = midpoints, wald_ratio = wald_ratio,
+      outcome_variance = outcome_variances,
+      exposure_effect = total_exposure_effects
+    ),
+    midpoint = integrate_midpoint(
+      age = age_seq, wald_ratio = wald_ratio,
+      outcome_variance = outcome_variances,
+      exposure_effect = total_exposure_effects
     ),
-    midpoint =  outcome_variances / total_exposure_effects ^ 2,
     "Method is neither trapezoidal nor midpoint. There should've been an error earlier."
   )
-  data.frame(
-    age = c(0, midpoints),
-    Gamma = Gamma,
-    var = Gamma_variance,
-    L95 = Gamma - 1.96 * sqrt(Gamma_variance),
-    H95 = Gamma + 1.96 * sqrt(Gamma_variance)
-  )
-
 }