diff --git a/DESCRIPTION b/DESCRIPTION
index d4637df..6138540 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,13 +1,13 @@
 Package: stabiot
 Type: Package
-Title: Common Statistical Analysis for Clinical Trials in Biotech
+Title: Common Statistical Analysis for Clinical Trials
 Version: 0.0.0.9000
 Authors@R: c(
     person("Kai", "Gu", , "gukai1212@163.com", role = c("aut", "cre", "cph"))
   )
 Maintainer: Kai Gu <gukai1212@163.com>
-Description: Providing common statistical analysis for oversighting process in 
-    clinical trials of biotech, making the outputs can be compared with SAS results.
+Description: Provides common statistical analysis for oversighting process in 
+    clinical trials of biotech, makes the outputs can be compared with SAS results.
 License: GPL (>= 3)
 Encoding: UTF-8
 LazyData: true
@@ -15,12 +15,17 @@ Depends:
     R (>= 3.6)
 Imports:
     checkmate,
+    DescTools,
     dplyr,
     emmeans,
     lifecycle,
     magrittr,
+    purrr,
+    rlang,
     stats,
-    tibble
+    survival,
+    tibble,
+    tidyr
 Suggests: 
     knitr,
     mmrm,
diff --git a/NAMESPACE b/NAMESPACE
index f8342eb..1e7e69b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,12 +1,29 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(print,or_ci)
+S3method(print,prop_ci)
 S3method(print,s_lsmeans)
 export("%>%")
+export(h_prep_prop)
 export(rrPostProb)
 export(s_get_lsmeans)
+export(s_odds_ratio)
+export(s_propci)
 import(checkmate)
+importFrom(dplyr,add_count)
+importFrom(dplyr,count)
+importFrom(dplyr,group_by)
 importFrom(lifecycle,deprecated)
 importFrom(magrittr,"%>%")
+importFrom(rlang,":=")
+importFrom(rlang,.data)
+importFrom(rlang,sym)
+importFrom(stats,as.formula)
+importFrom(stats,coef)
 importFrom(stats,confint)
 importFrom(stats,pbeta)
 importFrom(stats,rbinom)
+importFrom(stats,setNames)
+importFrom(survival,Surv)
+importFrom(survival,coxph)
+importFrom(survival,strata)
diff --git a/R/package.R b/R/package.R
index 8c96bd6..ffb77f9 100644
--- a/R/package.R
+++ b/R/package.R
@@ -6,7 +6,10 @@
 
 #' @import checkmate
 #' @importFrom lifecycle deprecated
-#' @importFrom stats pbeta rbinom confint
+#' @importFrom stats pbeta rbinom confint as.formula setNames coef
+#' @importFrom dplyr count add_count group_by
+#' @importFrom rlang sym := .data
+#' @importFrom survival coxph Surv strata
 NULL
 
 .onLoad <- function(libname, pkgname) {
diff --git a/R/pkg-methods.R b/R/pkg-methods.R
index 9e085b2..ee4e7ee 100644
--- a/R/pkg-methods.R
+++ b/R/pkg-methods.R
@@ -33,3 +33,44 @@ print.s_lsmeans <- function(x, ...) {
 
   invisible(x)
 }
+
+#' @describeIn prop_odds_ratio prints proportion and confidence interval.
+#' @exportS3Method
+#' @keywords internal
+print.prop_ci <- function(x, ...) {
+  cat(sprintf("Proportion and %s confidence interval:", x$params$method))
+  cat("\n")
+  print(x$prop_est)
+
+  if (!is.null(x$params$by.level)) {
+    cat(sprintf("\nProportion Difference and %s confidence interval:", x$params$diff.method))
+    cat("\n")
+    print(x$prop_diff)
+  }
+
+  invisible(x)
+}
+
+#' @describeIn prop_odds_ratio prints odds ratio and confidence interval.
+#' @exportS3Method
+#' @keywords internal
+print.or_ci <- function(x, ...) {
+  comp <- paste0(rev(x$params$by.level), collapse = "/")
+  cat(sprintf("Common Odds Ratio (%s) and %s confidence interval:", comp, x$params$or.method))
+  cat("\n")
+  print(x$or)
+
+  if (!is.null(x$params$strata)) {
+    cat(sprintf(
+      "\nStratified Odds Ratio (%s) using %s", comp,
+      ifelse(x$params$strata.method == "CMH",
+        "Cochran-Mantel-Haenszel Chi-Squared Test:",
+        "Conditional logistic regression:"
+      )
+    ))
+    cat("\n")
+    print(x$strata_or)
+  }
+
+  invisible(x)
+}
diff --git a/R/proportion.R b/R/proportion.R
new file mode 100644
index 0000000..84af805
--- /dev/null
+++ b/R/proportion.R
@@ -0,0 +1,388 @@
+# s_propci ----
+
+#' Computing Proportion and Odds Ratio
+#'
+#' @description `r lifecycle::badge("experimental")`
+#'
+#' Compute confidence interval for proportion and difference of binomial response,
+#' and odds ratio using the `DescTools` package. As an alternative, use `stats::glm`
+#' with `logit` link to estimate the odds ratio. Regarding to stratified odds ratio,
+#' use Cochran-Mantel-Haenszel test (`stats::mantelhaen.test`) or conditional logistic
+#' regression (`survival::clogit`) to handle with. However this CMH test function is
+#' limited to the 2 by 2 by K scenario, whereas conditional logistic regression
+#' can handle 2 by N by K case.
+#'
+#' @note
+#' * The proportion difference can only support the two-levels for `by` variable.
+#' * Compute odds ratio when `or.glm` is FALSE, 0.5 will be added to 2 x 2 table
+#' in case of zero entries. More details can be found in [DescTools::OddsRatio].
+#' * The arguments of this function can only match few parameters of SAS `proc freq`
+#' to calculate the odds ratio. For example, `mantelhaen.test` corresponds the
+#' Mantel-Haenszel Estimator with the same reference article, and using `correct = FALSE`
+#' can produce the Wald confidence interval same as the default parameters of SAS.
+#'
+#' @name prop_odds_ratio
+#' @order 1
+NULL
+
+
+#' @describeIn prop_odds_ratio Computes the confidence interval of proportion/response
+#' rate and difference of two proportions/rates.
+#'
+#' @param data (`data.frame`)\cr a data frame as input.
+#' @param var (`string`)\cr target variable name for estimation.
+#' @param by (`string`)\cr a optional variable to group by. If null, use the whole data.
+#' @param by.level (`vector`)\cr an optional vector for encoding `var` as a factor
+#'  and the second level will be as the reference group. If null, use the default
+#'  order to encode.
+#' @param event (`numeric` or `character`)\cr an option to define which one as the
+#'  event in the elements of `var`. By default, the positive and maximal one if
+#'  the `var` variable is numeric, or the first one of if the `var` variable is
+#'  character/factor.
+#' @param conf.level (`numeric`)\cr significance level for the returned confidence
+#'  interval.
+#' @param method (`string`)\cr string specifying which method to calculate the
+#'  confidence interval for binomial proportions, default is Clopper-Pearson. More
+#'  details see [DescTools::BinomCI].
+#' @param diff.method (`string`)\cr string specifying which method to calculate the
+#'  confidence interval for difference between binomial proportions, default is
+#'  Wald. More details see [DescTools::BinomDiffCI].
+#' @param alternative (`string`)\cr string specifying the alternative hypothesis,
+#'    must be one of "two.sided" (default), "greater" or "less".
+#' @param ... other arguments to be passed to [DescTools::BinomCI].
+#'
+#'
+#' @return
+#' * `s_propci` returns an object of class `prop_ci` that is a list contains
+#' proportion rate, proportion difference and input arguments.
+#'
+#' @export
+#'
+#' @examples
+#' set.seed(12)
+#' dta <- data.frame(
+#'   orr = sample(c(1, 0), 100, TRUE),
+#'   trtp = factor(rep(c("TRT", "PBO"), each = 50)),
+#'   strata1 = factor(sample(c("A", "B"), 100, TRUE)),
+#'   strata2 = factor(sample(c("C", "D"), 100, TRUE)),
+#'   strata3 = factor(sample(c("E", "F"), 100, TRUE))
+#' )
+#'
+#' # not having by variable:
+#' s_propci(dta, var = "orr")
+#'
+#' # two levels of by variable:
+#' s_propci(dta, var = "orr", by = "trtp")
+#' # opposite order of trtp with above:
+#' s_propci(dta, var = "orr", by = "trtp", by.level = c("TRT", "PBO"), event = 1)
+s_propci <- function(data,
+                     var,
+                     by = NULL,
+                     by.level = NULL,
+                     event = NULL,
+                     conf.level = 0.95,
+                     method = c(
+                       "clopper-pearson", "wald", "waldcc", "wilson",
+                       "wilsoncc", "modified wilson", "jeffreys",
+                       "modified jeffreys", "agresti-coull"
+                     ),
+                     diff.method = c(
+                       "wald", "waldcc", "score", "scorecc", "mn", "mee",
+                       "blj", "ha", "hal", "jp"
+                     ),
+                     alternative = c("two.sided", "less", "greater"),
+                     ...) {
+  assert_class(data, "data.frame")
+  assert_subset(var, names(data), empty.ok = FALSE)
+  assert_subset(by, names(data))
+  assert_subset(event, data[[var]])
+  assert_number(conf.level, lower = 0, upper = 1)
+  method <- match.arg(method, c(
+    "clopper-pearson", "wald", "waldcc", "wilson", "wilsoncc",
+    "modified wilson", "jeffreys", "modified jeffreys", "agresti-coull"
+  ), several.ok = FALSE)
+  diff.method <- match.arg(diff.method, c(
+    "wald", "waldcc", "score", "scorecc", "mn", "mee",
+    "blj", "ha", "hal", "jp"
+  ), several.ok = FALSE)
+  alternative <- match.arg(alternative, c(
+    "two.sided", "less", "greater"
+  ), several.ok = FALSE)
+
+  object <- h_prep_prop(data, var = var, by = by, by.level = by.level, event = event)
+  by <- object$by
+
+  est_res <- object$data %>%
+    count(!!sym(by), !!sym(var)) %>%
+    add_count(!!sym(by), wt = .data$n, name = "tot") %>%
+    dplyr::filter(!!sym(var) == object$event) %>%
+    split(as.formula(paste("~", by))) %>%
+    purrr::map(function(x) {
+      tibble::tibble(
+        group = x[[by]],
+        tibble::as_tibble(
+          DescTools::BinomCI(
+            x = x$n, n = x$tot, method = method, conf.level = conf.level,
+            sides = alternative, ...
+          )
+        )
+      )
+    }) %>%
+    purrr::list_rbind()
+
+  diff_res <- if (!is.null(object$by.level)) {
+    object$data %>%
+      count(!!sym(by), !!sym(var)) %>%
+      add_count(!!sym(by), wt = .data$n, name = "tot") %>%
+      dplyr::filter(!!sym(var) == object$event) %>%
+      tidyr::pivot_wider(names_from = by, values_from = c("n", "tot")) %>%
+      split(as.formula(paste("~", var))) %>%
+      purrr::map(function(x) {
+        tibble::tibble(
+          group = paste0(object$by.level, collapse = " - "),
+          tibble::as_tibble(
+            DescTools::BinomDiffCI(
+              x1 = x[[2]], n1 = x[[4]], x2 = x[[3]], n2 = x[[5]],
+              method = diff.method, sides = alternative, conf.level = conf.level
+            )
+          )
+        )
+      }) %>%
+      purrr::list_rbind()
+  }
+
+  structure(
+    list(
+      prop_est = est_res,
+      prop_diff = diff_res,
+      params = list(
+        var = var,
+        by = by,
+        by.level = object$by.level,
+        event = object$event,
+        conf.level = conf.level,
+        method = method,
+        diff.method = diff.method,
+        alternative = alternative
+      )
+    ),
+    class = "prop_ci"
+  )
+}
+
+
+# s_odds_ratio ----
+
+#' @describeIn prop_odds_ratio Computes the confidence interval of common odds ratio.
+#' And also provides stratified odds ratio with Cochran-Mantel-Haenszel chi-squared
+#' test and conditional logistic regression.
+#'
+#' @param strata (`vector`)\cr a character vector is used for stratification.
+#' @param or.glm (`logical`)\cr a logical indicating whether to use logit method to
+#'  calculate the odds ratio. If TRUE, the `glm` with logit is used, otherwise the
+#'  common method from `DescTools::OddsRatio`.
+#' @param or.method (`string`)\cr string indicating the method from `DescTools::OddsRatio`
+#'  used to calculate odds ratio when the `or.glm` is FALSE.
+#' @param strata.method (`string`)\cr string indicating the method used to calculate
+#'  stratified odds ratio, Cochran-Mantel-Haenszel Chi-Squared test (CMH) or conditional
+#'  logistic regression (clogit).
+#' @param correct (`logical`)\cr a logical indicating whether to apply continuity
+#'  correction when computing the test statistic.
+#' @param exact (`logical`)\cr a logical indicating whether the Mantel-Haenszel
+#'  test or the exact conditional test (given the strata margins) should be computed.
+#'
+#' @return
+#' * `s_odds_ratio` returns an object of class `or_ci` that is a list contains
+#' odds ratio with or without stratification and input arguments.
+#'
+#' @export
+#'
+#' @examples
+#'
+#' # without stratification:
+#' s_odds_ratio(dta, var = "orr", by = "trtp", or.method = "wald")
+#'
+#' # with three stratifications, strata1 - strata3:
+#' s_odds_ratio(
+#'   dta,
+#'   var = "orr", by = "trtp",
+#'   strata = c("strata1", "strata2", "strata3"),
+#'   strata.method = "CMH",
+#'   correct = FALSE
+#' )
+#'
+#' # using condition logistic regression:
+#' s_odds_ratio(
+#'   dta,
+#'   var = "orr", by = "trtp",
+#'   strata = c("strata1", "strata2", "strata3"),
+#'   strata.method = "clogit"
+#' )
+s_odds_ratio <- function(data,
+                         var,
+                         by,
+                         by.level = NULL,
+                         event = NULL,
+                         strata = NULL,
+                         conf.level = 0.95,
+                         or.glm = FALSE,
+                         or.method = c("wald", "mle", "midp"),
+                         strata.method = c("CMH", "clogit"),
+                         correct = FALSE,
+                         exact = FALSE) {
+  assert_class(data, "data.frame")
+  assert_subset(var, names(data), empty.ok = FALSE)
+  assert_subset(by, names(data), empty.ok = FALSE)
+  assert_subset(event, data[[var]])
+  assert_subset(strata, names(data))
+  assert_number(conf.level, lower = 0, upper = 1)
+  assert_logical(or.glm)
+  assert_logical(correct)
+  assert_logical(exact)
+  or.method <- match.arg(or.method, c("wald", "mle", "midp"), several.ok = FALSE)
+  strata.method <- match.arg(strata.method, c("CMH", "clogit"), several.ok = FALSE)
+
+  object <- h_prep_prop(data, var = var, by = by, by.level = by.level, event = event)
+  if (length(object$by.level) != 2) {
+    stop("The by.level should have two levels as input.")
+  }
+
+  or_res <- if (or.glm) {
+    mod <- stats::glm(
+      as.formula(paste(var, "~", by)),
+      data = object$data,
+      family = stats::binomial(link = "logit")
+    )
+    exp(c(coef(mod)[-1], confint(mod, level = conf.level)[-1, ]))
+  } else {
+    object$data %>%
+      count(!!sym(by), !!sym(var)) %>%
+      dplyr::arrange(!!sym(var) == object$event, !!sym(by)) %>%
+      tidyr::pivot_wider(names_from = var, values_from = "n") %>%
+      tibble::column_to_rownames(var = by) %>%
+      as.matrix() %>%
+      # follow the preferred 2x2 table structure
+      DescTools::OddsRatio(method = or.method, conf.level = conf.level)
+  }
+  or_res <- tibble::tibble(
+    !!!setNames(or_res, c("or.est", "lwr.ci", "upr.ci"))
+  )
+
+  stra_or_res <- if (!is.null(strata)) {
+    if (strata.method == "CMH") {
+      assert_set_equal(length(unique(data[[var]])), 2)
+      tab <- stats::xtabs(
+        as.formula(paste("~", paste(c(by, var, strata), collapse = "+"))),
+        data = data
+      )
+      # get the number of factors for each stratification
+      grpn <- strata %>%
+        purrr::map(~ count(data, data[[.x]])) %>%
+        purrr::map_int(nrow)
+      tb <- as.table(array(c(tab), dim = c(2, 2, prod(grpn))))
+      mod <- stats::mantelhaen.test(
+        tb,
+        conf.level = conf.level,
+        correct = correct, exact = exact
+      )
+      tibble::tibble(
+        !!!setNames(
+          c(mod$estimate, mod$conf.int, mod$p.value),
+          c("or.est", "lwr.ci", "upr.ci", "pval")
+        )
+      )
+    } else {
+      mod <- survival::clogit(
+        formula = as.formula(paste0(
+          var, " ~ ", by, " + strata(",
+          paste(strata, collapse = ", "), ")"
+        )),
+        data = data
+      )
+
+      # defaultly the first level is regarded as the reference
+      names(coef(mod)) %>%
+        purrr::map(
+          ~ data.frame(
+            or.est = exp(coef(mod)[.x]),
+            lwr.ci = exp(confint(mod, level = conf.level)[.x, 1]),
+            upr.ci = exp(confint(mod, level = conf.level)[.x, 2]),
+            row.names = gsub(pattern = paste0("^", by), x = .x, "")
+          )
+        ) %>%
+        purrr::list_rbind() %>%
+        tibble::tibble()
+    }
+  } else {
+    NULL
+  }
+
+  structure(
+    list(
+      or = or_res,
+      strata_or = stra_or_res,
+      params = list(
+        var = var,
+        by = by,
+        by.level = object$by.level,
+        event = object$event,
+        strata = strata,
+        conf.level = conf.level,
+        or.method = ifelse(or.glm, "logit", or.method),
+        strata.method = ifelse(
+          !is.null(strata),
+          strata.method,
+          NA
+        )
+      )
+    ),
+    class = "or_ci"
+  )
+}
+
+
+# h_prep_prop ----
+
+#' @describeIn prop_odds_ratio Helper Function for Pre-processing Proportion Data.
+#'
+#' @export
+#'
+h_prep_prop <- function(data,
+                        var,
+                        by,
+                        by.level,
+                        event) {
+  if (is.null(by)) {
+    by <- "Total"
+    data[[by]] <- "Total"
+  } else {
+    if (!is.null(by.level)) {
+      assert_set_equal(by.level, as.character(unique(data[[by]])))
+      data[[by]] <- factor(data[[by]], levels = by.level)
+    } else {
+      if (is.factor(data[[by]])) {
+        by.level <- levels(data[[by]])
+      } else {
+        data[[by]] <- factor(data[[by]], levels = unique(data[[by]]))
+        by.level <- levels(data[[by]])
+      }
+    }
+  }
+
+  if (is.null(event)) {
+    event <- if (is.numeric(data[[var]])) {
+      max(unique(data[[var]])[unique(data[[var]]) > 0])
+    } else if (is.character(event)) {
+      unique(data[[var]])[1]
+    } else if (is.factor(data[[var]])) {
+      levels(data[[var]])[1]
+    }
+  }
+
+  list(
+    data = data,
+    by = by,
+    by.level = by.level,
+    event = event
+  )
+}
diff --git a/R/summarize-lsmeans.R b/R/summarize-lsmeans.R
index 0cf3cbd..6381ba9 100644
--- a/R/summarize-lsmeans.R
+++ b/R/summarize-lsmeans.R
@@ -55,6 +55,20 @@
 #' when the null hypothesis is θ not superiority to `2`, the `side` should be
 #' `greater`, and `null` is defined as `2`.
 #'
+#' @references
+#' SAS code for your reference with consistent results.
+#' ```
+#' proc mixed data=fev_data;
+#'   class ARMCD(ref='PBO') AVISIT RACE SEX USUBJID;
+#'   model FEV1 = RACE SEX ARMCD ARMCD*AVISIT / ddfm=KR;
+#'   repeated AVISIT / subject=USUBJID type=UN r rcorr;
+#'   lsmeans ARMCD*AVISIT / cl alpha=0.05 diff slice=AVISIT;
+#'   lsmeans ARMCD / cl alpha=0.05 diff;
+#'   lsmestimate ARMCD*AVISIT [1,1 4] [-1,2 4] / cl upper alpha=0.025 testvalue=2;
+#'   ods output lsmeans=lsm diffs=diff LSMEstimates=est;
+#' run;
+#' ```
+#'
 #' @export
 #'
 #' @examples
diff --git a/README.Rmd b/README.Rmd
index 856a5d2..d85990f 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -19,10 +19,19 @@ knitr::opts_chunk$set(
 [![R-CMD-check](https://github.com/kaigu1990/stabiot/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/kaigu1990/stabiot/actions/workflows/R-CMD-check.yaml)
 <!-- badges: end -->
 
-The goal of `stabiot` is to help statistician and programmer using R functions 
-and method to oversight the outputs usually produced by SAS from vendors in clinical 
-trials. The data sets should be ADaM format perfectly, but not must follow the CDISC 
-standard.
+The goal of `stabiot` is to assist statisticians and programmers in using R functions 
+and methods to oversee the outputs often produced by SAS from vendors in clinical 
+trials. The data sets would be ADaM format preferably, but they do not have to 
+follow the CDISC standard.
+
+In order to ensure the quality and accuracy of the results, I prefer to wrap mature 
+R package rather than rebuild the statistical methods. For now the completed sections 
+are listed as shown below.
+
+- Simulation of sample size determination by Bayesian.
+- Summarize Least-squares Means from models, such as ANCOVA and MMRM.
+- Computing response rate, odds ratio with or without stratification, and corresponding 
+confidence interval.
 
 ## Installation
 
diff --git a/README.md b/README.md
index d6b9472..5931c9a 100644
--- a/README.md
+++ b/README.md
@@ -8,10 +8,19 @@
 [![R-CMD-check](https://github.com/kaigu1990/stabiot/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/kaigu1990/stabiot/actions/workflows/R-CMD-check.yaml)
 <!-- badges: end -->
 
-The goal of `stabiot` is to help statistician and programmer using R
-functions and method to oversight the outputs usually produced by SAS
-from vendors in clinical trials. The data sets should be ADaM format
-perfectly, but not must follow the CDISC standard.
+The goal of `stabiot` is to assist statisticians and programmers in
+using R functions and methods to oversee the outputs often produced by
+SAS from vendors in clinical trials. The data sets would be ADaM format
+preferably, but they do not have to follow the CDISC standard.
+
+In order to ensure the quality and accuracy of the results, I prefer to
+wrap mature R package rather than rebuild the statistical methods. For
+now the completed sections are listed as shown below.
+
+- Simulation of sample size determination by Bayesian.
+- Summarize Least-squares Means from models, such as ANCOVA and MMRM.
+- Computing response rate, odds ratio with or without stratification,
+  and corresponding confidence interval.
 
 ## Installation
 
diff --git a/inst/WORDLIST b/inst/WORDLIST
index f2971c8..85127db 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -1,10 +1,19 @@
 ADaM
 ANCOVA
+BinomCI
+BinomDiffCI
 CDISC
 CMD
+CMH
+Clopper
+DescTools
+Haenszel
 MMRM
+OddsRatio
 PN
+Pre
 TRT
+clogit
 magrittr
 oversighting
 θ
diff --git a/man/prop_odds_ratio.Rd b/man/prop_odds_ratio.Rd
new file mode 100644
index 0000000..a423f44
--- /dev/null
+++ b/man/prop_odds_ratio.Rd
@@ -0,0 +1,186 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/proportion.R, R/pkg-methods.R
+\name{prop_odds_ratio}
+\alias{prop_odds_ratio}
+\alias{print.prop_ci}
+\alias{print.or_ci}
+\alias{s_propci}
+\alias{s_odds_ratio}
+\alias{h_prep_prop}
+\title{Computing Proportion and Odds Ratio}
+\usage{
+\method{print}{prop_ci}(x, ...)
+
+\method{print}{or_ci}(x, ...)
+
+s_propci(
+  data,
+  var,
+  by = NULL,
+  by.level = NULL,
+  event = NULL,
+  conf.level = 0.95,
+  method = c("clopper-pearson", "wald", "waldcc", "wilson", "wilsoncc",
+    "modified wilson", "jeffreys", "modified jeffreys", "agresti-coull"),
+  diff.method = c("wald", "waldcc", "score", "scorecc", "mn", "mee", "blj", "ha", "hal",
+    "jp"),
+  alternative = c("two.sided", "less", "greater"),
+  ...
+)
+
+s_odds_ratio(
+  data,
+  var,
+  by,
+  by.level = NULL,
+  event = NULL,
+  strata = NULL,
+  conf.level = 0.95,
+  or.glm = FALSE,
+  or.method = c("wald", "mle", "midp"),
+  strata.method = c("CMH", "clogit"),
+  correct = FALSE,
+  exact = FALSE
+)
+
+h_prep_prop(data, var, by, by.level, event)
+}
+\arguments{
+\item{...}{other arguments to be passed to \link[DescTools:BinomCI]{DescTools::BinomCI}.}
+
+\item{data}{(\code{data.frame})\cr a data frame as input.}
+
+\item{var}{(\code{string})\cr target variable name for estimation.}
+
+\item{by}{(\code{string})\cr a optional variable to group by. If null, use the whole data.}
+
+\item{by.level}{(\code{vector})\cr an optional vector for encoding \code{var} as a factor
+and the second level will be as the reference group. If null, use the default
+order to encode.}
+
+\item{event}{(\code{numeric} or \code{character})\cr an option to define which one as the
+event in the elements of \code{var}. By default, the positive and maximal one if
+the \code{var} variable is numeric, or the first one of if the \code{var} variable is
+character/factor.}
+
+\item{conf.level}{(\code{numeric})\cr significance level for the returned confidence
+interval.}
+
+\item{method}{(\code{string})\cr string specifying which method to calculate the
+confidence interval for binomial proportions, default is Clopper-Pearson. More
+details see \link[DescTools:BinomCI]{DescTools::BinomCI}.}
+
+\item{diff.method}{(\code{string})\cr string specifying which method to calculate the
+confidence interval for difference between binomial proportions, default is
+Wald. More details see \link[DescTools:BinomDiffCI]{DescTools::BinomDiffCI}.}
+
+\item{alternative}{(\code{string})\cr string specifying the alternative hypothesis,
+must be one of "two.sided" (default), "greater" or "less".}
+
+\item{strata}{(\code{vector})\cr a character vector is used for stratification.}
+
+\item{or.glm}{(\code{logical})\cr a logical indicating whether to use logit method to
+calculate the odds ratio. If TRUE, the \code{glm} with logit is used, otherwise the
+common method from \code{DescTools::OddsRatio}.}
+
+\item{or.method}{(\code{string})\cr string indicating the method from \code{DescTools::OddsRatio}
+used to calculate odds ratio when the \code{or.glm} is FALSE.}
+
+\item{strata.method}{(\code{string})\cr string indicating the method used to calculate
+stratified odds ratio, Cochran-Mantel-Haenszel Chi-Squared test (CMH) or conditional
+logistic regression (clogit).}
+
+\item{correct}{(\code{logical})\cr a logical indicating whether to apply continuity
+correction when computing the test statistic.}
+
+\item{exact}{(\code{logical})\cr a logical indicating whether the Mantel-Haenszel
+test or the exact conditional test (given the strata margins) should be computed.}
+}
+\value{
+\itemize{
+\item \code{s_propci} returns an object of class \code{prop_ci} that is a list contains
+proportion rate, proportion difference and input arguments.
+}
+
+\itemize{
+\item \code{s_odds_ratio} returns an object of class \code{or_ci} that is a list contains
+odds ratio with or without stratification and input arguments.
+}
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+Compute confidence interval for proportion and difference of binomial response,
+and odds ratio using the \code{DescTools} package. As an alternative, use \code{stats::glm}
+with \code{logit} link to estimate the odds ratio. Regarding to stratified odds ratio,
+use Cochran-Mantel-Haenszel test (\code{stats::mantelhaen.test}) or conditional logistic
+regression (\code{survival::clogit}) to handle with. However this CMH test function is
+limited to the 2 by 2 by K scenario, whereas conditional logistic regression
+can handle 2 by N by K case.
+}
+\section{Functions}{
+\itemize{
+\item \code{print(prop_ci)}: prints proportion and confidence interval.
+
+\item \code{print(or_ci)}: prints odds ratio and confidence interval.
+
+\item \code{s_propci()}: Computes the confidence interval of proportion/response
+rate and difference of two proportions/rates.
+
+\item \code{s_odds_ratio()}: Computes the confidence interval of common odds ratio.
+And also provides stratified odds ratio with Cochran-Mantel-Haenszel chi-squared
+test and conditional logistic regression.
+
+\item \code{h_prep_prop()}: Helper Function for Pre-processing Proportion Data.
+
+}}
+\note{
+\itemize{
+\item The proportion difference can only support the two-levels for \code{by} variable.
+\item Compute odds ratio when \code{or.glm} is FALSE, 0.5 will be added to 2 x 2 table
+in case of zero entries. More details can be found in \link[DescTools:OddsRatio]{DescTools::OddsRatio}.
+\item The arguments of this function can only match few parameters of SAS \verb{proc freq}
+to calculate the odds ratio. For example, \code{mantelhaen.test} corresponds the
+Mantel-Haenszel Estimator with the same reference article, and using \code{correct = FALSE}
+can produce the Wald confidence interval same as the default parameters of SAS.
+}
+}
+\examples{
+set.seed(12)
+dta <- data.frame(
+  orr = sample(c(1, 0), 100, TRUE),
+  trtp = factor(rep(c("TRT", "PBO"), each = 50)),
+  strata1 = factor(sample(c("A", "B"), 100, TRUE)),
+  strata2 = factor(sample(c("C", "D"), 100, TRUE)),
+  strata3 = factor(sample(c("E", "F"), 100, TRUE))
+)
+
+# not having by variable:
+s_propci(dta, var = "orr")
+
+# two levels of by variable:
+s_propci(dta, var = "orr", by = "trtp")
+# opposite order of trtp with above:
+s_propci(dta, var = "orr", by = "trtp", by.level = c("TRT", "PBO"), event = 1)
+
+# without stratification:
+s_odds_ratio(dta, var = "orr", by = "trtp", or.method = "wald")
+
+# with three stratifications, strata1 - strata3:
+s_odds_ratio(
+  dta,
+  var = "orr", by = "trtp",
+  strata = c("strata1", "strata2", "strata3"),
+  strata.method = "CMH",
+  correct = FALSE
+)
+
+# using condition logistic regression:
+s_odds_ratio(
+  dta,
+  var = "orr", by = "trtp",
+  strata = c("strata1", "strata2", "strata3"),
+  strata.method = "clogit"
+)
+}
+\keyword{internal}
diff --git a/man/s_get_lsmeans.Rd b/man/s_get_lsmeans.Rd
index 971b095..9b41ce4 100644
--- a/man/s_get_lsmeans.Rd
+++ b/man/s_get_lsmeans.Rd
@@ -109,4 +109,18 @@ fit2 <- fev_data \%>\%
 
 s_get_lsmeans(fit2, "ARMCD")
 }
+\references{
+SAS code for your reference with consistent results.
+
+\if{html}{\out{<div class="sourceCode">}}\preformatted{proc mixed data=fev_data;
+  class ARMCD(ref='PBO') AVISIT RACE SEX USUBJID;
+  model FEV1 = RACE SEX ARMCD ARMCD*AVISIT / ddfm=KR;
+  repeated AVISIT / subject=USUBJID type=UN r rcorr;
+  lsmeans ARMCD*AVISIT / cl alpha=0.05 diff slice=AVISIT;
+  lsmeans ARMCD / cl alpha=0.05 diff;
+  lsmestimate ARMCD*AVISIT [1,1 4] [-1,2 4] / cl upper alpha=0.025 testvalue=2;
+  ods output lsmeans=lsm diffs=diff LSMEstimates=est;
+run;
+}\if{html}{\out{</div>}}
+}
 \keyword{internal}
diff --git a/tests/testthat/test-proportion.R b/tests/testthat/test-proportion.R
new file mode 100644
index 0000000..e14145b
--- /dev/null
+++ b/tests/testthat/test-proportion.R
@@ -0,0 +1,217 @@
+# s_propci ----
+
+test_that("s_propci works as expected with default arguments", {
+  set.seed(12)
+  dta <- data.frame(
+    orr = sample(c(1, 0), 100, TRUE),
+    trtp = factor(rep(c("TRT", "PBO"), each = 50))
+  )
+
+  res <- s_propci(dta, var = "orr", by = "trtp")
+  expect_class(res, "prop_ci")
+  expect_equal(
+    res$prop_est,
+    tibble::tibble(
+      group = factor(c("PBO", "TRT"), levels = c("PBO", "TRT")),
+      est = c(0.44, 0.40),
+      lwr.ci = c(0.2999072, 0.2640784),
+      upr.ci = c(0.5874559, 0.5482060)
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res$prop_diff,
+    tibble::tibble(
+      group = "PBO - TRT",
+      est = 0.04,
+      lwr.ci = -0.1533125,
+      upr.ci = 0.2333125
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res$params,
+    list(
+      var = "orr",
+      by = "trtp",
+      by.level = c("PBO", "TRT"),
+      event = 1,
+      conf.level = 0.95,
+      method = "clopper-pearson",
+      diff.method = "wald",
+      alternative = "two.sided"
+    )
+  )
+})
+
+test_that("s_propci works as expected with specific group level and event", {
+  set.seed(12)
+  dta <- data.frame(
+    orr = sample(c(1, 0), 100, TRUE),
+    trtp = factor(rep(c("TRT", "PBO"), each = 50))
+  )
+
+  res <- s_propci(dta, var = "orr", by = "trtp", by.level = c("TRT", "PBO"), event = 0)
+  expect_equal(
+    res$prop_est,
+    tibble::tibble(
+      group = factor(c("TRT", "PBO"), levels = c("TRT", "PBO")),
+      est = c(0.6, 0.56),
+      lwr.ci = c(0.4517940, 0.4125441),
+      upr.ci = c(0.7359216, 0.7000928)
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res$prop_diff,
+    tibble::tibble(
+      group = "TRT - PBO",
+      est = 0.04,
+      lwr.ci = -0.1533125,
+      upr.ci = 0.2333125
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res$params[c("by.level", "event")],
+    list(
+      by.level = c("TRT", "PBO"),
+      event = 0
+    )
+  )
+})
+
+test_that("s_propci works as expected with different method for binomial CI", {
+  set.seed(12)
+  dta <- data.frame(
+    orr = sample(c(1, 0), 100, TRUE),
+    trtp = factor(rep(c("TRT", "PBO"), each = 50))
+  )
+
+  res <- s_propci(dta, var = "orr", by = "trtp", method = "wald")
+  expect_equal(
+    res$prop_est,
+    tibble::tibble(
+      group = factor(c("PBO", "TRT"), levels = c("PBO", "TRT")),
+      est = c(0.44, 0.40),
+      lwr.ci = c(0.3024111, 0.2642097),
+      upr.ci = c(0.5775889, 0.5357903)
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(res$params$method, "wald")
+})
+
+# s_odds_ratio ----
+
+test_that("s_odds_ratio works as expected with default arguments, no stratification", {
+  set.seed(12)
+  dta <- data.frame(
+    orr = sample(c(1, 0), 100, TRUE),
+    trtp = factor(rep(c("TRT", "PBO"), each = 50)),
+    strata1 = factor(sample(c("A", "B"), 100, TRUE)),
+    strata2 = factor(sample(c("C", "D"), 100, TRUE)),
+    strata3 = factor(sample(c("E", "F"), 100, TRUE))
+  )
+
+  res <- s_odds_ratio(dta, var = "orr", by = "trtp")
+  expect_class(res, "or_ci")
+  expect_null(res$strata_or)
+  expect_equal(
+    res$or,
+    tibble::tibble(
+      or.est = 0.8484848,
+      lwr.ci = 0.3831831,
+      upr.ci = 1.8788054
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res$params,
+    list(
+      var = "orr",
+      by = "trtp",
+      by.level = c("PBO", "TRT"),
+      event = 1,
+      strata = NULL,
+      conf.level = 0.95,
+      or.method = "wald",
+      strata.method = NA
+    )
+  )
+})
+
+test_that("s_odds_ratio works as expected with stratification", {
+  set.seed(12)
+  dta <- data.frame(
+    orr = sample(c(1, 0), 100, TRUE),
+    trtp = factor(rep(c("TRT", "PBO"), each = 50)),
+    strata1 = factor(sample(c("A", "B"), 100, TRUE)),
+    strata2 = factor(sample(c("C", "D"), 100, TRUE)),
+    strata3 = factor(sample(c("E", "F"), 100, TRUE))
+  )
+
+  res <- s_odds_ratio(
+    dta,
+    var = "orr", by = "trtp",
+    strata = c("strata1", "strata2", "strata3")
+  )
+  expect_equal(
+    res$strata_or,
+    tibble::tibble(
+      or.est = 0.7499121,
+      lwr.ci = 0.3248143,
+      upr.ci = 1.7313529,
+      pval = 0.5089342
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(res$params$strata.method, "CMH")
+})
+
+test_that("s_odds_ratio works as expected using or.glm and clogit", {
+  set.seed(12)
+  dta <- data.frame(
+    orr = sample(c(1, 0), 100, TRUE),
+    trtp = factor(rep(c("TRT", "PBO"), each = 50)),
+    strata1 = factor(sample(c("A", "B"), 100, TRUE)),
+    strata2 = factor(sample(c("C", "D"), 100, TRUE)),
+    strata3 = factor(sample(c("E", "F"), 100, TRUE))
+  )
+
+  res <- s_odds_ratio(
+    dta,
+    var = "orr", by = "trtp",
+    or.glm = TRUE,
+    strata = c("strata1", "strata2", "strata3"),
+    strata.method = "clogit"
+  )
+  expect_equal(
+    res$or,
+    tibble::tibble(
+      or.est = 0.8484848,
+      lwr.ci = 0.3811997,
+      upr.ci = 1.8797355
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res$strata_or,
+    tibble::tibble(
+      data.frame(
+        or.est = 0.7592608,
+        lwr.ci = 0.335024,
+        upr.ci = 1.720704,
+        row.names = "TRT"
+      )
+    ),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res$params[c("or.method", "strata.method")],
+    list(
+      or.method = "logit",
+      strata.method = "clogit"
+    )
+  )
+})