diff --git a/DESCRIPTION b/DESCRIPTION
index d6693f21..b0bbbc79 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Type: Package
 Package: RobinCar2
 Title: ROBust INference for Covariate Adjustment in Randomized Clinical Trials
-Version: 0.2.3
-Date: 2026-07-02
+Version: 0.2.3.9000
+Date: 2026-07-04
 Authors@R:
     c(
       person("Liming", "Li", , "liming.li1@astrazeneca.com", role = c("aut", "cre"), comment = c(ORCID = "0009-0008-6870-0878")),
@@ -44,6 +44,7 @@ Imports:
 Suggests:
     knitr,
     rmarkdown,
+    speff2trial,
     testthat (>= 3.0)
 VignetteBuilder:
     knitr
diff --git a/NAMESPACE b/NAMESPACE
index 4fedc5dd..ad8ed7fa 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -35,6 +35,7 @@ export(predict_counterfactual)
 export(robin_glm)
 export(robin_lm)
 export(robin_surv)
+export(surv_control)
 export(table)
 export(treatment_effect)
 export(vcovG)
diff --git a/NEWS.md b/NEWS.md
index 541d519b..667ff61e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,13 @@
+# RobinCar2 0.2.3.9000
+
+### New Features
+
+* The new `surv_control` argument in `robin_surv()` allows to fine-control the root finding algorithm used for the hazard ratio estimation.
+
+### Bug Fixes
+
+* Previously `robin_surv()` gave small numerical differences to `survival::coxph()` for the hazard ratio estimate. This is now fixed.
+
 # RobinCar2 0.2.3
 
 ### Bug Fixes
diff --git a/R/RobinCar2-package.R b/R/RobinCar2-package.R
index 95c18c98..1ed32fbd 100644
--- a/R/RobinCar2-package.R
+++ b/R/RobinCar2-package.R
@@ -6,7 +6,7 @@
 #' @import checkmate
 #' @importFrom numDeriv grad
 #' @importFrom stats predict residuals fitted model.response model.matrix coefficients family
-#' gaussian terms glm var family pnorm var as.formula qnorm lm confint
+#' @importFrom stats gaussian terms glm var family pnorm var as.formula qnorm lm confint
 #' @importFrom sandwich vcovHC
 #' @importFrom MASS negative.binomial
 #' @importFrom utils combn tail
diff --git a/R/survival.R b/R/survival.R
index a20bd1ea..50bac240 100644
--- a/R/survival.R
+++ b/R/survival.R
@@ -10,6 +10,7 @@
 #'
 #' @param score_fun (`function`) The log-rank score function to be used for estimation.
 #' @param interval (`numeric`) A numeric vector of length 2 specifying the interval in which to search for the root.
+#' @param control (`list`) Control parameters from [surv_control()] for the root-finding algorithm.
 #' @param ... Additional arguments passed to `score_fun`.
 #' @return A list containing:
 #' - `theta`: The estimated log hazard ratio.
@@ -18,10 +19,11 @@
 #' - `n`: The number of observations used in the calculation.
 #'
 #' @keywords internal
-h_log_hr_est_via_score <- function(score_fun, interval = c(-5, 5), ...) {
+h_log_hr_est_via_score <- function(score_fun, interval = c(-5, 5), control = surv_control(), ...) {
   assert_function(score_fun, args = c("theta", "use_ties_factor", "calculate_variance"))
   assert_numeric(interval, len = 2L, finite = TRUE)
   assert_true(interval[1] < interval[2])
+  h_assert_surv_control(control)
 
   score_solution <- stats::uniroot(
     score_fun,
@@ -30,7 +32,9 @@ h_log_hr_est_via_score <- function(score_fun, interval = c(-5, 5), ...) {
     check.conv = TRUE, # If the root cannot be found, an error is thrown.
     use_ties_factor = FALSE,
     calculate_variance = FALSE, # We will only do this at the solution.
-    tol = .Machine$double.eps^0.1, # Use a small tolerance for convergence.
+    tol = control$tol,
+    maxiter = control$maxiter,
+    trace = control$trace,
     ...
   )
   score_root <- score_solution$root
@@ -125,6 +129,7 @@ robin_surv_comparison <- function(
   exp_level,
   control_level,
   contrast,
+  control = surv_control(),
   unadj_score_fun = NULL,
   ...
 ) {
@@ -138,6 +143,7 @@ robin_surv_comparison <- function(
   assert_count(control_level)
   assert_true(exp_level != control_level)
   assert_string(contrast)
+  h_assert_surv_control(control)
 
   # Subset data to the two treatment arms of interest.
   trt_levels <- vars$levels[c(control_level, exp_level)]
@@ -147,34 +153,35 @@ robin_surv_comparison <- function(
   data[[vars$treatment]] <- stats::relevel(data[[vars$treatment]], ref = trt_levels[1L])
 
   # Prepare arguments for the test and estimation calls below.
-  args <- list(
+  score_args <- list(
     score_fun = score_fun,
     df = data,
     ...
   )
 
   # Perform the log-rank test via the score function.
-  test_result <- do.call(h_lr_test_via_score, args)
+  test_result <- do.call(h_lr_test_via_score, score_args)
 
   # Estimate the log hazard ratio via the score function, if requested.
   hr_result <- if (contrast == "hazardratio") {
+    hr_args <- c(score_args, list(control = control))
     # If an unadjusted score function is provided, use it to estimate the log hazard ratio first.
     if (!is.null(unadj_score_fun)) {
       assert_function(unadj_score_fun)
       assert_true(length(vars$covariates) > 0)
       args_to_drop <- c("model", "hr_se_plugin_adjusted", "check_rand_strat_warning")
-      unadj_args <- args[!(names(args) %in% args_to_drop)]
+      unadj_args <- hr_args[!(names(hr_args) %in% args_to_drop)]
       unadj_args$score_fun <- unadj_score_fun
       # Get theta_hat from the unadjusted score function.
       unadj_hr_result <- do.call(h_log_hr_est_via_score, unadj_args)
       # Add this to the arguments for the adjusted score function call below.
-      args$theta_hat <- unadj_hr_result$theta
+      hr_args$theta_hat <- unadj_hr_result$theta
     } else {
       # We enforce to have no covariates in this case.
       assert_true(length(vars$covariates) == 0L)
     }
     # Estimate the log hazard ratio via the score function.
-    do.call(h_log_hr_est_via_score, args)
+    do.call(h_log_hr_est_via_score, hr_args)
   } else {
     list(
       theta = NA_real_,
@@ -218,6 +225,7 @@ robin_surv_no_strata_no_cov <- function(
   exp_level,
   control_level,
   contrast,
+  control = surv_control(),
   check_rand_strat_warning = FALSE
 ) {
   robin_surv_comparison(
@@ -227,6 +235,7 @@ robin_surv_no_strata_no_cov <- function(
     exp_level = exp_level,
     control_level = control_level,
     contrast = contrast,
+    control = control,
     treatment = vars$treatment,
     time = vars$time,
     status = vars$status,
@@ -243,6 +252,7 @@ robin_surv_strata <- function(
   exp_level,
   control_level,
   contrast,
+  control = surv_control(),
   check_rand_strat_warning = FALSE
 ) {
   robin_surv_comparison(
@@ -252,6 +262,7 @@ robin_surv_strata <- function(
     exp_level = exp_level,
     control_level = control_level,
     contrast = contrast,
+    control = control,
     treatment = vars$treatment,
     time = vars$time,
     status = vars$status,
@@ -264,7 +275,7 @@ robin_surv_strata <- function(
 #' @describeIn survival_comparison_functions without strata and without covariates, based on
 #'   [h_lr_score_cov()] and [h_lr_score_no_strata_no_cov()] (which is used to find the unadjusted
 #'   log hazard ratio estimate).
-robin_surv_cov <- function(vars, data, exp_level, control_level, contrast, ...) {
+robin_surv_cov <- function(vars, data, exp_level, control_level, contrast, control = surv_control(), ...) {
   robin_surv_comparison(
     score_fun = h_lr_score_cov,
     unadj_score_fun = h_lr_score_no_strata_no_cov,
@@ -273,6 +284,7 @@ robin_surv_cov <- function(vars, data, exp_level, control_level, contrast, ...)
     exp_level = exp_level,
     control_level = control_level,
     contrast = contrast,
+    control = control,
     treatment = vars$treatment,
     time = vars$time,
     status = vars$status,
@@ -285,7 +297,7 @@ robin_surv_cov <- function(vars, data, exp_level, control_level, contrast, ...)
 #' @describeIn survival_comparison_functions with strata and covariates, based on
 #'   [h_lr_score_strat_cov()] and [h_lr_score_strat()] (which is used to find the unadjusted
 #'   log hazard ratio estimate).
-robin_surv_strata_cov <- function(vars, data, exp_level, control_level, contrast, ...) {
+robin_surv_strata_cov <- function(vars, data, exp_level, control_level, contrast, control = surv_control(), ...) {
   robin_surv_comparison(
     score_fun = h_lr_score_strat_cov,
     unadj_score_fun = h_lr_score_strat,
@@ -294,6 +306,7 @@ robin_surv_strata_cov <- function(vars, data, exp_level, control_level, contrast
     exp_level = exp_level,
     control_level = control_level,
     contrast = contrast,
+    control = control,
     treatment = vars$treatment,
     time = vars$time,
     status = vars$status,
@@ -397,6 +410,43 @@ h_events_table <- function(data, vars) {
   )
 }
 
+#' Control Survival Analysis Root Finding
+#'
+#' Control parameters for the root-finding algorithm used by [robin_surv()] to estimate
+#' the log hazard ratio.
+#'
+#' @param tol (`number`) The desired accuracy, passed to [stats::uniroot()].
+#' @param maxiter (`count`) The maximum number of iterations, passed to [stats::uniroot()].
+#' @param trace (`count`) Tracing level, passed to [stats::uniroot()].
+#'
+#' @return A named list with elements `tol`, `maxiter`, and `trace`.
+#'
+#' @export
+surv_control <- function(tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0) {
+  control <- list(
+    tol = tol,
+    maxiter = maxiter,
+    trace = trace
+  )
+  h_assert_surv_control(control)
+
+  list(
+    tol = tol,
+    maxiter = as.integer(maxiter),
+    trace = as.integer(trace)
+  )
+}
+
+h_assert_surv_control <- function(control) {
+  assert_list(control, len = 3L, names = "unique")
+  assert_names(names(control), must.include = c("tol", "maxiter", "trace"))
+  assert_number(control$tol, lower = 0, finite = TRUE)
+  assert_integerish(control$maxiter, len = 1L, lower = 1L, any.missing = FALSE)
+  assert_integerish(control$trace, len = 1L, lower = 0L, any.missing = FALSE)
+
+  invisible(control)
+}
+
 #' Covariate Adjusted and Stratified Survival Analysis
 #'
 #' Calculate log-rank test as well as hazard ratio estimates for survival data, optionally adjusted
@@ -416,6 +466,9 @@ h_events_table <- function(data, vars) {
 #' @param contrast (`character(1)`) The contrast statistic to be used, currently only `"hazardratio"`
 #'   is supported. Can be disabled by specifying `"none"`, in which case only the log-rank test is performed.
 #' @param test (`character(1)`) The test to be used, currently only `"logrank"` is supported.
+#' @param control (`list`) Control parameters for the root-finding algorithm used to estimate the
+#'   log hazard ratio, usually created with [surv_control()]. The arguments correspond to
+#'   `tol`, `maxiter`, and `trace` in [stats::uniroot()].
 #' @param ... Additional arguments passed to the survival analysis functions, in particular `hr_se_plugin_adjusted`
 #'   (please see [here][survival_score_functions] for details).
 #' @return A `surv_effect` object containing the results of the survival analysis.
@@ -462,6 +515,7 @@ robin_surv <- function(
   comparisons,
   contrast = c("hazardratio", "none"),
   test = "logrank",
+  control = surv_control(),
   ...
 ) {
   attr(formula, ".Environment") <- environment()
@@ -472,6 +526,7 @@ robin_surv <- function(
   assert_subset(all.vars(treatment), names(data))
   contrast <- match.arg(contrast)
   test <- match.arg(test)
+  h_assert_surv_control(control)
 
   input <- h_prep_survival_input(formula, data, treatment)
 
@@ -529,6 +584,7 @@ robin_surv <- function(
         exp_level = exp_level,
         control_level = control_level,
         contrast = contrast,
+        control = control,
         check_rand_strat_warning = !give_rand_strat_warning,
         ...
       )
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 9e9eee24..53880bf0 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -67,3 +67,4 @@ reference:
     contents:
     - find_data
     - confint
+    - surv_control
diff --git a/man/h_log_hr_est_via_score.Rd b/man/h_log_hr_est_via_score.Rd
index 287feb2e..a23ebe07 100644
--- a/man/h_log_hr_est_via_score.Rd
+++ b/man/h_log_hr_est_via_score.Rd
@@ -4,13 +4,20 @@
 \alias{h_log_hr_est_via_score}
 \title{Estimate Log Hazard Ratio via Score Function}
 \usage{
-h_log_hr_est_via_score(score_fun, interval = c(-5, 5), ...)
+h_log_hr_est_via_score(
+  score_fun,
+  interval = c(-5, 5),
+  control = surv_control(),
+  ...
+)
 }
 \arguments{
 \item{score_fun}{(\code{function}) The log-rank score function to be used for estimation.}
 
 \item{interval}{(\code{numeric}) A numeric vector of length 2 specifying the interval in which to search for the root.}
 
+\item{control}{(\code{list}) Control parameters from \code{\link[=surv_control]{surv_control()}} for the root-finding algorithm.}
+
 \item{...}{Additional arguments passed to \code{score_fun}.}
 }
 \value{
diff --git a/man/robin_surv.Rd b/man/robin_surv.Rd
index 7f30b8fc..2caebc36 100644
--- a/man/robin_surv.Rd
+++ b/man/robin_surv.Rd
@@ -11,6 +11,7 @@ robin_surv(
   comparisons,
   contrast = c("hazardratio", "none"),
   test = "logrank",
+  control = surv_control(),
   ...
 )
 }
@@ -35,6 +36,10 @@ is supported. Can be disabled by specifying \code{"none"}, in which case only th
 
 \item{test}{(\code{character(1)}) The test to be used, currently only \code{"logrank"} is supported.}
 
+\item{control}{(\code{list}) Control parameters for the root-finding algorithm used to estimate the
+log hazard ratio, usually created with \code{\link[=surv_control]{surv_control()}}. The arguments correspond to
+\code{tol}, \code{maxiter}, and \code{trace} in \code{\link[stats:uniroot]{stats::uniroot()}}.}
+
 \item{...}{Additional arguments passed to the survival analysis functions, in particular \code{hr_se_plugin_adjusted}
 (please see \link[=survival_score_functions]{here} for details).}
 }
diff --git a/man/robin_surv_comparison.Rd b/man/robin_surv_comparison.Rd
index 1be19df7..923d801d 100644
--- a/man/robin_surv_comparison.Rd
+++ b/man/robin_surv_comparison.Rd
@@ -11,6 +11,7 @@ robin_surv_comparison(
   exp_level,
   control_level,
   contrast,
+  control = surv_control(),
   unadj_score_fun = NULL,
   ...
 )
diff --git a/man/surv_control.Rd b/man/surv_control.Rd
new file mode 100644
index 00000000..4b474547
--- /dev/null
+++ b/man/surv_control.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/survival.R
+\name{surv_control}
+\alias{surv_control}
+\title{Control Survival Analysis Root Finding}
+\usage{
+surv_control(tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)
+}
+\arguments{
+\item{tol}{(\code{number}) The desired accuracy, passed to \code{\link[stats:uniroot]{stats::uniroot()}}.}
+
+\item{maxiter}{(\code{count}) The maximum number of iterations, passed to \code{\link[stats:uniroot]{stats::uniroot()}}.}
+
+\item{trace}{(\code{count}) Tracing level, passed to \code{\link[stats:uniroot]{stats::uniroot()}}.}
+}
+\value{
+A named list with elements \code{tol}, \code{maxiter}, and \code{trace}.
+}
+\description{
+Control parameters for the root-finding algorithm used by \code{\link[=robin_surv]{robin_surv()}} to estimate
+the log hazard ratio.
+}
diff --git a/man/survival_comparison_functions.Rd b/man/survival_comparison_functions.Rd
index 02ac97be..490f030f 100644
--- a/man/survival_comparison_functions.Rd
+++ b/man/survival_comparison_functions.Rd
@@ -14,6 +14,7 @@ robin_surv_no_strata_no_cov(
   exp_level,
   control_level,
   contrast,
+  control = surv_control(),
   check_rand_strat_warning = FALSE
 )
 
@@ -23,12 +24,29 @@ robin_surv_strata(
   exp_level,
   control_level,
   contrast,
+  control = surv_control(),
   check_rand_strat_warning = FALSE
 )
 
-robin_surv_cov(vars, data, exp_level, control_level, contrast, ...)
+robin_surv_cov(
+  vars,
+  data,
+  exp_level,
+  control_level,
+  contrast,
+  control = surv_control(),
+  ...
+)
 
-robin_surv_strata_cov(vars, data, exp_level, control_level, contrast, ...)
+robin_surv_strata_cov(
+  vars,
+  data,
+  exp_level,
+  control_level,
+  contrast,
+  control = surv_control(),
+  ...
+)
 }
 \arguments{
 \item{vars}{(\code{list}) A list containing \code{levels}, \code{treatment}, and \code{covariates}.}
diff --git a/tests/testthat/_snaps/surv_effect.md b/tests/testthat/_snaps/surv_effect.md
index f3e37644..cbe399f8 100644
--- a/tests/testthat/_snaps/surv_effect.md
+++ b/tests/testthat/_snaps/surv_effect.md
@@ -9,7 +9,7 @@
       Contrast     : Log Hazard Ratio
       
                        Estimate Std.Err Z Value Pr(>|z|)   
-      Male v.s. Female  0.53343 0.16727   3.189 0.001428 **
+      Male v.s. Female  0.53037 0.16718  3.1724 0.001512 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
@@ -32,7 +32,7 @@
       Contrast     : Stratified Log Hazard Ratio
       
                        Estimate Std.Err Z Value Pr(>|z|)   
-      Male v.s. Female  0.55482 0.17063  3.2516 0.001147 **
+      Male v.s. Female   0.5536  0.1706  3.2451 0.001174 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
@@ -54,8 +54,8 @@
       
       Contrast     : Covariate-adjusted Log Hazard Ratio
       
-                       Estimate Std.Err Z Value Pr(>|z|)  
-      Male v.s. Female  0.47686 0.18608  2.5626  0.01039 *
+                       Estimate Std.Err Z Value Pr(>|z|)   
+      Male v.s. Female  0.48335 0.18631  2.5944 0.009477 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
@@ -79,7 +79,7 @@
       Contrast     : Covariate-adjusted Stratified Log Hazard Ratio
       
                        Estimate Std.Err Z Value Pr(>|z|)   
-      Male v.s. Female  0.55219 0.19133  2.8861   0.0039 **
+      Male v.s. Female  0.54791 0.19118   2.866 0.004157 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
diff --git a/tests/testthat/_snaps/survival.md b/tests/testthat/_snaps/survival.md
index f8c0188d..ca3c6218 100644
--- a/tests/testthat/_snaps/survival.md
+++ b/tests/testthat/_snaps/survival.md
@@ -1,6 +1,6 @@
 # h_log_hr_est_via_score works as expected
 
-    list(theta = 0.53342923961288, se = 0.167270504465807, sigma_l2 = 0.156756810948865, 
+    list(theta = 0.530368384142102, se = 0.167180005660338, sigma_l2 = 0.156926569811932, 
         n = 228L)
 
 # h_log_hr_est_via_score does not give spurious warning
@@ -10,7 +10,7 @@
 
 # h_log_hr_est_via_score extends the search interval as needed
 
-    list(theta = 0.529999478789383, se = 0.167169124542081, sigma_l2 = 0.156946999327738, 
+    list(theta = 0.530397994806822, se = 0.167180879290955, sigma_l2 = 0.156924929726076, 
         n = 228L)
 
 # h_lr_test_via_score works as expected
@@ -20,8 +20,8 @@
 
 # robin_surv_comparison works as expected without covariate adjustment
 
-    list(estimate = 0.53342923961288, se = 0.167270504465807, hr_n = 228L, 
-        hr_sigma_l2 = 0.156756810948865, test_stat = 3.21352484896035, 
+    list(estimate = 0.530368384142102, se = 0.167180005660338, hr_n = 228L, 
+        hr_sigma_l2 = 0.156926569811932, test_stat = 3.21352484896035, 
         p_value = 0.0013111645203555, test_score = 0.0895537761860842, 
         test_n = 228L, test_sigma_l2 = 0.17706769289317, give_rand_strat_warning = FALSE)
 
@@ -34,22 +34,22 @@
 
 # robin_surv_no_strata_no_cov works as expected
 
-    list(estimate = -0.53342923961288, se = 0.167270504465807, hr_n = 228L, 
-        hr_sigma_l2 = 0.156756810948865, test_stat = -3.21352484896035, 
+    list(estimate = -0.530368384142102, se = 0.167180005660338, hr_n = 228L, 
+        hr_sigma_l2 = 0.156926569811932, test_stat = -3.21352484896035, 
         p_value = 0.0013111645203555, test_score = -0.0895537761860842, 
         test_n = 228L, test_sigma_l2 = 0.17706769289317, give_rand_strat_warning = FALSE)
 
 # robin_surv_strata works as expected
 
-    list(estimate = -0.554820697243091, se = 0.170628735365298, hr_n = 227L, 
-        hr_sigma_l2 = 0.151310755205329, test_stat = -3.2855836062253, 
+    list(estimate = -0.553595870294275, se = 0.17059512525311, hr_n = 227L, 
+        hr_sigma_l2 = 0.151370382611702, test_stat = -3.2855836062253, 
         p_value = 0.00101771334472415, test_score = -0.0896871248297144, 
         test_n = 227L, test_sigma_l2 = 0.169145720705825, give_rand_strat_warning = FALSE)
 
 # robin_surv_cov works as expected
 
-    list(estimate = -0.503821677843015, se = 0.1656799127954, hr_n = 228L, 
-        hr_sigma_l2 = 0.156990532256188, test_stat = -3.06610317830763, 
+    list(estimate = -0.5016113861505, se = 0.16561813992359, hr_n = 228L, 
+        hr_sigma_l2 = 0.157110113177543, test_stat = -3.06610317830763, 
         p_value = 0.0021686846372504, test_score = -0.0851266162073163, 
         test_n = 228L, test_sigma_l2 = 0.175748674116587, give_rand_strat_warning = FALSE)
 
@@ -104,8 +104,8 @@
 
 # robin_surv works as expected with strata
 
-    structure(c(0.554820697243091, 0.170628735365298, 3.25162520870402, 
-    0.00114747222174709), dim = c(1L, 4L), dimnames = list("Male v.s. Female", 
+    structure(c(0.553595870294275, 0.17059512525311, 3.24508610356194, 
+    0.00117415114967595), dim = c(1L, 4L), dimnames = list("Male v.s. Female", 
         c("Estimate", "Std.Err", "Z Value", "Pr(>|z|)")))
 
 ---
@@ -130,8 +130,8 @@
 
 # robin_surv works as expected with covariates
 
-    structure(c(-0.00438916915405738, 0.155227542623415, -0.0282757111262503, 
-    0.977442252566086), dim = c(1L, 4L), dimnames = list("1 v.s. 0", 
+    structure(c(-0.00445020431030736, 0.155227462755351, -0.0286689238573794, 
+    0.977128641344646), dim = c(1L, 4L), dimnames = list("1 v.s. 0", 
         c("Estimate", "Std.Err", "Z Value", "Pr(>|z|)")))
 
 ---
@@ -152,8 +152,8 @@
 
 # robin_surv works as expected with strata and covariates
 
-    structure(c(-0.0122310580846328, 0.159010010937409, -0.0769200505837796, 
-    0.938687146743022), dim = c(1L, 4L), dimnames = list("1 v.s. 0", 
+    structure(c(-0.0122920932408829, 0.159009954919017, -0.0773039225571954, 
+    0.938381770498044), dim = c(1L, 4L), dimnames = list("1 v.s. 0", 
         c("Estimate", "Std.Err", "Z Value", "Pr(>|z|)")))
 
 ---
@@ -177,8 +177,8 @@
 
 # robin_surv works also with character variable in the correlation case
 
-    structure(c(0.553358244734497, 0.236938504567525, 2.33545090421045, 
-    0.0195198778286913), dim = c(1L, 4L), dimnames = list("Male v.s. Female", 
+    structure(c(0.546676517528469, 0.237207661991659, 2.30463262838319, 
+    0.0211871569907112), dim = c(1L, 4L), dimnames = list("Male v.s. Female", 
         c("Estimate", "Std.Err", "Z Value", "Pr(>|z|)")))
 
 # robin_surv gives a warning if stratified randomization was specified but simple log rank test used
@@ -200,7 +200,7 @@
       Contrast     : Log Hazard Ratio
       
                        Estimate Std.Err Z Value Pr(>|z|)   
-      Male v.s. Female   0.5258  0.1675  3.1391 0.001695 **
+      Male v.s. Female  0.52309 0.16742  3.1244 0.001782 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
@@ -232,7 +232,7 @@
       Contrast     : Stratified Log Hazard Ratio
       
                        Estimate Std.Err Z Value Pr(>|z|)   
-      Male v.s. Female  0.55178 0.16997  3.2463 0.001169 **
+      Male v.s. Female  0.54989 0.16992  3.2362 0.001211 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
@@ -264,7 +264,7 @@
       Contrast     : Covariate-adjusted Log Hazard Ratio
       
                        Estimate Std.Err Z Value Pr(>|z|)   
-      Male v.s. Female  0.51573 0.16373  3.1499 0.001633 **
+      Male v.s. Female  0.51335 0.16366  3.1366 0.001709 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
@@ -297,7 +297,7 @@
       Contrast     : Covariate-adjusted Stratified Log Hazard Ratio
       
                        Estimate Std.Err Z Value Pr(>|z|)   
-      Male v.s. Female  0.54011 0.16667  3.2406 0.001193 **
+      Male v.s. Female  0.53849 0.16663  3.2317 0.001231 **
       ---
       Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
       
diff --git a/tests/testthat/test-survival.R b/tests/testthat/test-survival.R
index c3232525..d9e5eeb2 100644
--- a/tests/testthat/test-survival.R
+++ b/tests/testthat/test-survival.R
@@ -37,6 +37,33 @@ test_that("h_log_hr_est_via_score extends the search interval as needed", {
   expect_snapshot_value(result, tolerance = 1e-4, style = "deparse")
 })
 
+test_that("surv_control works as expected", {
+  expect_identical(
+    surv_control(tol = 1e-6, maxiter = 10, trace = 1),
+    list(tol = 1e-6, maxiter = 10L, trace = 1L)
+  )
+  expect_error(surv_control(tol = -1), "Assertion")
+  expect_error(surv_control(maxiter = 0), "Assertion")
+  expect_error(surv_control(trace = -1), "Assertion")
+})
+
+test_that("h_log_hr_est_via_score passes control to uniroot", {
+  expect_error(
+    h_log_hr_est_via_score(
+      h_lr_score_no_strata_no_cov,
+      interval = c(-0.2, 0.2),
+      control = surv_control(maxiter = 1),
+      df = surv_data,
+      treatment = "sex",
+      time = "time",
+      status = "status",
+      randomization_strata = character()
+    ),
+    "no sign change found in 1 iteration",
+    fixed = TRUE
+  )
+})
+
 test_that("h_lr_test_via_score works as expected", {
   result <- h_lr_test_via_score(
     h_lr_score_no_strata_no_cov,
@@ -418,7 +445,8 @@ test_that("robin_surv works as expected without strata or covariates", {
   result <- robin_surv(
     Surv(time, status) ~ 1,
     data = surv_data,
-    treatment = ecog ~ sr(1)
+    treatment = ecog ~ sr(1),
+    control = surv_control()
   )
   expect_s3_class(result, "surv_effect")
   expect_snapshot_value(result$log_hr_coef_mat, tolerance = 1e-4, style = "deparse")
@@ -447,6 +475,44 @@ test_that("robin_surv gives the same results as RobinCar functions without strat
   expect_equal(result$log_hr_coef_mat[, "Std.Err"], robincar_result$se, tolerance = 1e-4)
 })
 
+test_that("robin_surv gives the same estimate as coxph", {
+  skip_if_not_installed("speff2trial")
+
+  data("ACTG175", package = "speff2trial", envir = environment())
+  trial_dat <- ACTG175
+  trial_dat <- trial_dat[trial_dat$arms %in% c(0, 3), ]
+  trial_dat <- data.frame(
+    id = trial_dat$pidnum,
+    drugs = trial_dat$drugs,
+    cd40 = trial_dat$cd40,
+    days = trial_dat$days,
+    obs = trial_dat$cens,
+    didanosine = factor(trial_dat$arms == 3),
+    strat = factor(trial_dat$strat)
+  )
+
+  coxres <- survival::coxph(
+    survival::Surv(days, obs) ~ didanosine,
+    data = trial_dat,
+    ties = "breslow"
+  )
+  expect_warning(
+    robinres <- robin_surv(
+      survival::Surv(days, obs) ~ 1,
+      treatment = didanosine ~ pb(strat),
+      data = trial_dat
+    ),
+    "It looks like you have not included all of the variables that were used during randomization",
+    fixed = TRUE
+  )
+  expect_equal(
+    robinres$estimate,
+    coxres$coefficients,
+    ignore_attr = TRUE,
+    tolerance = 1e-6
+  )
+})
+
 test_that("robin_surv works as expected with strata", {
   result <- robin_surv(
     Surv(time, status) ~ 1 + strata(strata),