diff --git a/DESCRIPTION b/DESCRIPTION
index ac7845b..d7d01b0 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -33,10 +33,11 @@ Suggests:
     knitr,
     rmarkdown,
     spelling,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    survival
 VignetteBuilder:
     knitr
-Config/roxygen2/version: 8.1.0
+Config/roxygen2/version: 8.1.1
 Config/testthat/edition: 3
 Encoding: UTF-8
 Language: en-US
diff --git a/NAMESPACE b/NAMESPACE
index 7146dcf..a748a7e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,8 @@
 # Generated by roxygen2: do not edit by hand
-importFrom("stats", "complete.cases")
+
 export(acro_add_comments)
 export(acro_add_exception)
+export(acro_coxph)
 export(acro_crosstab)
 export(acro_custom_output)
 export(acro_disable_rounding)
@@ -25,3 +26,7 @@ export(create_factors)
 export(is_excluded)
 export(is_invalid)
 export(to_pandas_categorical)
+importFrom(stats,
+  complete.cases,
+  na.omit
+)
diff --git a/R/acro_regression.R b/R/acro_regression.R
index d0315da..efd6728 100644
--- a/R/acro_regression.R
+++ b/R/acro_regression.R
@@ -36,3 +36,126 @@ acro_glm <- function(formula, data, family) {
   }
   model$summary()
 }
+
+#' Fits a Cox proportional hazards regression model
+#'
+#' @param formula a formula object, with the response on the left of a ~ operator, and the terms on the right. The response must be a survival object as returned by the Surv function.
+#' @param data a data.frame in which to interpret the variables named in the formula, or in the subset and the weights argument.
+#' @param weights vector of case weights
+#' @param subset expression indicating which subset of the rows of data should be used in the fit. All observations are included by default.
+#' @param na.action a missing-data filter function.
+#' @param robust should a robust variance be computed.
+#' @param init vector of initial values of the iteration.
+#' @param cluster optional variable which clusters the observations, for the purposes of a robust variance. If present, it implies robust. This variable will normally be found in data.
+#' @param istate optional variable giving the current state at the start each interval.
+#' @param statedata optional data set used to describe multistate models.
+#' @param ... Any other parameters.
+#'
+#' @returns Fitted CoxPHFitter results object
+#' @export
+
+#' @importFrom stats na.omit
+
+acro_coxph <- function(formula, data, weights = NULL, subset = NULL, na.action = NULL, robust = TRUE, init = NULL, cluster = NULL, istate = NULL, statedata = NULL, ...) {
+  if (is.null(acroEnv$ac)) {
+    stop("ACRO has not been initialised. Please first call acro_init()")
+  }
+
+  # Check for multi-state model arguments (not supported by lifelines)
+  if (!is.null(istate) || !is.null(statedata)) {
+    stop("Multi-state models utilizing 'istate' or 'statedata' are not supported at the moment.")
+  }
+
+  # Handle subsetting if provided
+  if (!missing(subset)) {
+    subset_expr <- substitute(subset)
+    subsetted_columns <- eval(subset_expr, data)
+    data <- data[subsetted_columns & !is.na(subsetted_columns), ]
+  }
+
+  # Handle na.action
+  if (is.null(na.action)) {
+    na.action <- getOption("na.action")
+  }
+
+  na_function <- match.fun(na.action)
+  data <- na_function(data)
+
+  # Parse the formula to extract duration and event columns from the left-hand side
+  surv_func <- formula[[2]]
+  if (!inherits(surv_func, "call") || as.character(surv_func[[1]]) != "Surv") {
+    stop("The left-hand side of the formula must use Surv(), e.g., Surv(time, event) ~ x1 + x2")
+  }
+  duration_col <- as.character(surv_func[[2]])
+  event_col <- as.character(surv_func[[3]])
+
+  # Validate and transform event column for python lifelines
+  event_values <- unique(na.omit(data[[event_col]]))
+
+  if (all(event_values %in% c(1, 2))) {
+    data[[event_col]][data[[event_col]] == 1] <- 0
+    data[[event_col]][data[[event_col]] == 2] <- 1
+  } else if (!all(event_values %in% c(0, 1))) {
+    # If it's neither {1, 2} nor {0, 1}, throw an error
+    stop(paste0("In Surv(", duration_col, ", ", event_col, ") : Invalid status values"))
+  }
+
+  # Extract the right-hand side (independent variables) as a formula string
+  independent_vars <- deparse(formula[[3]])
+
+  # Handle weights if provided as a column or a vector
+  weights_col <- NULL
+  if (!missing(weights)) {
+    weights_expr <- substitute(weights)
+    if (!is.null(weights_expr)) {
+      # If weights passed as a quoted character string matching a column name
+      if (is.character(weights_expr) && weights_expr %in% names(data)) {
+        weights_col <- weights_expr
+      }
+      # If weights passed as an unquoted column name
+      else if (deparse(weights_expr) %in% names(data)) {
+        weights_col <- deparse(weights_expr)
+      }
+      # If weights passed as an external vector, expression, or function
+      else {
+        data$acro_weight_col__ <- eval(weights_expr, data, parent.frame())
+        weights_col <- "acro_weight_col__"
+      }
+    }
+  }
+
+  # Handle cluster
+  cluster_col <- NULL
+  if (!missing(cluster)) {
+    cluster_expr <- substitute(cluster)
+    if (!is.null(cluster_expr)) {
+      if (is.character(cluster_expr) && cluster_expr %in% names(data)) {
+        cluster_col <- cluster_expr
+      } else if (deparse(cluster_expr) %in% names(data)) {
+        cluster_col <- deparse(cluster_expr)
+      }
+      # If weights passed as an external vector, expression, or function
+      else {
+        data$acro_cluster_col__ <- eval(cluster_expr, data, parent.frame())
+        cluster_col <- "acro_cluster_col__"
+      }
+      robust <- TRUE # Automatically trigger robust variance when clustering is provided
+    }
+  }
+
+  # Check for any unused arguments
+  if (length(list(...)) > 0) {
+    warning("Unused arguments were provided: ", paste0(names(list(...)), collapse = ", "), "\n", "Please use the help command to learn more about the function.")
+  }
+
+  model <- acroEnv$ac$coxph(
+    df = data,
+    duration_col = duration_col,
+    event_col = event_col,
+    weights_col = weights_col,
+    robust = robust,
+    formula = independent_vars,
+    initial_point = init,
+    cluster_col = cluster_col
+  )
+}
diff --git a/R/acro_tables.R b/R/acro_tables.R
index 57b58cb..bf0d5e8 100644
--- a/R/acro_tables.R
+++ b/R/acro_tables.R
@@ -177,6 +177,8 @@ acro_table <- function(index, columns, dnn = NULL, deparse.level = 0, useNA = "n
 #' @returns Summary of the data
 #' @export
 
+#' @importFrom stats complete.cases
+
 acro_summarise <- function(.data, ..., .groups = NULL, .by = NULL) {
   if (is.null(acroEnv$ac)) {
     stop("ACRO has not been initialised. Please first call acro_init()")
diff --git a/man/acro_coxph.Rd b/man/acro_coxph.Rd
new file mode 100644
index 0000000..87e5a4c
--- /dev/null
+++ b/man/acro_coxph.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/acro_regression.R
+\name{acro_coxph}
+\alias{acro_coxph}
+\title{Fits a Cox proportional hazards regression model}
+\usage{
+acro_coxph(
+  formula,
+  data,
+  weights = NULL,
+  subset = NULL,
+  na.action = NULL,
+  robust = TRUE,
+  init = NULL,
+  cluster = NULL,
+  istate = NULL,
+  statedata = NULL,
+  ...
+)
+}
+\arguments{
+\item{formula}{a formula object, with the response on the left of a ~ operator, and the terms on the right. The response must be a survival object as returned by the Surv function.}
+
+\item{data}{a data.frame in which to interpret the variables named in the formula, or in the subset and the weights argument.}
+
+\item{weights}{vector of case weights}
+
+\item{subset}{expression indicating which subset of the rows of data should be used in the fit. All observations are included by default.}
+
+\item{na.action}{a missing-data filter function.}
+
+\item{robust}{should a robust variance be computed.}
+
+\item{init}{vector of initial values of the iteration.}
+
+\item{cluster}{optional variable which clusters the observations, for the purposes of a robust variance. If present, it implies robust. This variable will normally be found in data.}
+
+\item{istate}{optional variable giving the current state at the start each interval.}
+
+\item{statedata}{optional data set used to describe multistate models.}
+
+\item{...}{Any other parameters.}
+}
+\value{
+Fitted CoxPHFitter results object
+}
+\description{
+Fits a Cox proportional hazards regression model
+}
diff --git a/tests/testthat/test-acro_coxph.R b/tests/testthat/test-acro_coxph.R
new file mode 100644
index 0000000..5be1006
--- /dev/null
+++ b/tests/testthat/test-acro_coxph.R
@@ -0,0 +1,258 @@
+test_data <-  lung
+
+test_that("acro_coxph without initialising ACRO object first", {
+  acroEnv$ac <- NULL
+  expect_error(acro_coxph(Surv(time, status) ~ age + sex, data = test_data), "ACRO has not been initialised. Please first call acro_init()")
+})
+
+test_that("acro_coxph throws an error with multi-state arguments", {
+  acro_init()
+  expect_error(
+    acro_coxph(Surv(time, status) ~ age + sex, data = test_data, istate = 1),
+    "Multi-state models utilizing 'istate' or 'statedata' are not supported at the moment."
+  )
+  expect_error(
+    acro_coxph(Surv(time, status) ~ age + sex, data = test_data, statedata = test_data),
+    "Multi-state models utilizing 'istate' or 'statedata' are not supported at the moment."
+  )
+})
+
+test_that("acro_coxph validates the left side of the formula uses Surv() function", {
+  acro_init()
+  expect_error(
+    acro_coxph(time ~ age + sex, data = test_data),
+    "The left-hand side of the formula must use Surv\\(\\)"
+  )
+})
+
+test_that("acro_coxph works", {
+  acro_init()
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data)
+  r_coefs <- stats::coef(r_model)
+
+  # ACRO model
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data)
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with subsetting", {
+  acro_init()
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, subset = age > 50)
+  r_coefs <- stats::coef(r_model)
+
+  # ACRO model
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, subset = age > 50)
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with na.action", {
+  acro_init()
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, na.action = na.pass)
+  r_coefs <- stats::coef(r_model)
+
+  # ACRO model
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, na.action = na.pass)
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with unvalid status options(eg. having a status as 4)", {
+  acro_init()
+  test_data$status[test_data$status == 1] <- 4
+
+  expect_error(
+    acro_coxph(Surv(time, status) ~ age + sex, data = test_data),
+    "In Surv\\(time, status\\) : Invalid status values"
+  )
+})
+
+test_that("acro_coxph with weights as a column from the dataset", {
+  acro_init()
+
+  #add a weight column in the dataset
+  test_data$my_weights <- runif(nrow(test_data), 0.8, 1.2)
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, weights = my_weights)
+  r_coefs <- stats::coef(r_model)
+
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, weights = my_weights)
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with weights as a string column name from the dataset", {
+  acro_init()
+
+  #add a weight column in the dataset
+  test_data$my_weights <- runif(nrow(test_data), 0.8, 1.2)
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, weights = my_weights)
+  r_coefs <- stats::coef(r_model)
+
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, weights = "my_weights")
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with weights provided by the user", {
+  acro_init()
+
+  #define a vector weight
+  my_weights <-  runif(nrow(test_data), 0.8, 1.2)
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, weights = my_weights)
+  r_coefs <- stats::coef(r_model)
+
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, weights = my_weights)
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+   tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with cluster as a column in the dataset", {
+  acro_init()
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, cluster = inst)
+  r_coefs <- stats::coef(r_model)
+
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, cluster = inst)
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with cluster as a quoted column name in the dataset", {
+  acro_init()
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, cluster = inst)
+  r_coefs <- stats::coef(r_model)
+
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, cluster = "inst")
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph with cluster provided by the user", {
+  acro_init()
+
+  #define a cluster vector
+  my_cluster <- sample(1:10, size = nrow(test_data), replace = TRUE)
+
+  # R survival model
+  r_model <- survival::coxph(Surv(time, status) ~ age + sex, data = test_data, cluster = my_cluster)
+  r_coefs <- stats::coef(r_model)
+
+  acro_model <- acro_coxph(Surv(time, status) ~ age + sex, data = test_data, cluster = my_cluster)
+
+  # Compare the coefficients
+  expect_equal(
+    py_to_r(unname(acro_model$params_["age"])),
+    unname(r_coefs["age"]),
+    tolerance = 0.05
+  )
+  expect_equal(
+    py_to_r(unname(acro_model$params_["sex"])),
+    unname(r_coefs["sex"]),
+    tolerance = 0.05
+  )
+})
+
+test_that("acro_coxph warns on unused arguments in ...", {
+  acro_init()
+
+  expect_warning(
+    acro_coxph(Surv(time, status) ~ age, data = test_data, dump_arg = 999),
+    "Unused arguments were provided: dump_arg"
+  )
+})