diff --git a/DESCRIPTION b/DESCRIPTION
index 7c9a3c2..afcbebb 100644
--- a/DESCRIPTION
+++ b/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'
diff --git a/R/integrate_estimates.R b/R/integrate_estimates.R
new file mode 100644
index 0000000..a6094df
--- /dev/null
+++ b/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))
+  )
+}
diff --git a/R/perform_MR.R b/R/perform_MR.R
index 0a26bc8..d7b81ed 100644
--- a/R/perform_MR.R
+++ b/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)
-  )
-
 }
diff --git a/man/calculate_mr_variance_trapz.Rd b/man/calculate_mr_variance_trapz.Rd
index 0d98f93..16f4570 100644
--- a/man/calculate_mr_variance_trapz.Rd
+++ b/man/calculate_mr_variance_trapz.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/perform_MR.R
+% Please edit documentation in R/integrate_estimates.R
 \name{calculate_mr_variance_trapz}
 \alias{calculate_mr_variance_trapz}
 \title{Calculate time-dependent effect variance of MR effects estimated by
diff --git a/man/integrate_estimates.Rd b/man/integrate_estimates.Rd
deleted file mode 100644
index 5e1e866..0000000
--- a/man/integrate_estimates.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/perform_MR.R
-\name{integrate_estimates}
-\alias{integrate_estimates}
-\title{Numerically integrate effect estimates by age}
-\usage{
-integrate_estimates(age, beta, method = c("trapezoidal", "midpoint"))
-}
-\arguments{
-\item{age}{Numeric vector containing ages}
-
-\item{beta}{Numeric vector containing effect estimates}
-
-\item{method}{The integration method to use. Either \code{trapezoidal} (default)
-or \code{midpoint}.}
-}
-\value{
-A numeric matrix with columns age and cumulative integral
-}
-\description{
-Numerically integrate effect estimates by age
-}
-\keyword{internal}
diff --git a/man/integrate_midpoint.Rd b/man/integrate_midpoint.Rd
new file mode 100644
index 0000000..30774d2
--- /dev/null
+++ b/man/integrate_midpoint.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/integrate_estimates.R
+\name{integrate_midpoint}
+\alias{integrate_midpoint}
+\title{Integrate time-dependent estimates based on midpoint rule and compute 95\% CIs}
+\usage{
+integrate_midpoint(age, wald_ratio, exposure_effect, outcome_variance)
+}
+\arguments{
+\item{age}{The ages at which to integrate}
+
+\item{wald_ratio}{The calculated Wald ratios}
+
+\item{exposure_effect}{The genetic effect on the exposure at every age}
+
+\item{outcome_variance}{The variance of the genetic effect on the outcome
+at every age}
+}
+\value{
+A data.frame with age, effect (Gamma), variance and 95\% CIs
+}
+\description{
+Integrate time-dependent estimates based on midpoint rule and compute 95\% CIs
+}
+\keyword{internal}
diff --git a/man/integrate_trapezoid.Rd b/man/integrate_trapezoid.Rd
new file mode 100644
index 0000000..813c5cb
--- /dev/null
+++ b/man/integrate_trapezoid.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/integrate_estimates.R
+\name{integrate_trapezoid}
+\alias{integrate_trapezoid}
+\title{Integrate time-dependent effects using the trapezoidal rule and compute 95\%
+CIs}
+\usage{
+integrate_trapezoid(age, wald_ratio, exposure_effect, outcome_variance)
+}
+\arguments{
+\item{age}{The ages at which to integrate}
+
+\item{wald_ratio}{The calculated Wald ratios}
+
+\item{exposure_effect}{The genetic effect on the exposure at every age}
+
+\item{outcome_variance}{The variance of the genetic effect on the outcome
+at every age}
+}
+\value{
+A data.frame with age, effect (Gamma), variance and 95\% CIs
+}
+\description{
+Integrate time-dependent effects using the trapezoidal rule and compute 95\%
+CIs
+}
+\keyword{internal}