Miettinen-Nurminen (MN) method has increasingly been requested by regulatory agencies, rather than the traditional Wald method that is based on the asymptotic normal distribution. This is particularly relevant for non-inferiority trials where it's appropriate for the variance to be estimated under the null hypothesis. Calculate MN Test Statistics within PROC FREQ.
The formula for compute the MN test statistic can be found in many articles, such as a paper in lexjansen Constructing Confidence Intervals for the Differences of Binomial Proportions in SAS.
As can be seen from the above description, if we would like to get the CI of M&N method, we should compute the root of the maximum likelihood estimations with the closed-form solution provided by M&N. Referring to the following method from metalite.ae rate-compare vignettes, we can create an equation p2 - p1 = delta and the chi-square statistic to find out the root.
As can be seen from the above explanation, if we would like to get the confidence interval (CI) of M&N method, we should compute the root of the maximum likelihood estimations with the closed-form solution provided by M&N. Referring to the following method, we can construct an equation to find out the root from the chi-square distribution with 1 degree of freedom. So I roughly understand that the delta can be regarded as the CI when the equation of p1 - p2 - delta divided by variance is equivalent to the chi-square statistic. Therefore, I create a function so that I can obtain the root using a solving process.
mn_root <- function(delta, alpha = 0.05) {
diff <- p1 - p2 - delta
theta <- n2 / n1
a <- 1 + theta
b <- -(1 + theta + p1 + theta * p2 + delta * (theta + 2))
c <- delta^2 + delta * (2 * p1 + theta + 1) + p1 + theta * p2
d <- -p1 * delta * (1 + delta)
v <- b^3 / (27 * a^3) - b * c / (6 * a^2) + d / (a * 2)
# here should be '-' instead of '+' inside the square root.
u <- ifelse(v > 0, 1, -1) * sqrt(b^2 / (9 * a^2) - c / (a * 3))
w <- (pi + acos(v / u^3)) / 3
p1d <- 2 * u * cos(w) - b / (a * 3)
p2d <- p1d - delta
# n <- n1 + n2
var <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2) * (n1 + n2) / (n1 + n2 - 1)
chisq <- diff^2 / var
return(chisq - qchisq(1 - alpha, 1))
}Now I also create an example data as shown below to compute the risk difference (RD) of response rate between two treatments.
library(tidyverse)
ana <- data.frame(
treatment = factor(c(rep(2, 100), rep(1, 100)), level = c(2, 1)),
response = c(rep(0, 80), rep(1, 20), rep(0, 40), rep(1, 60)),
stratum = c(rep(1:4, 12), 1, 3, 3, 1, rep(1:4, 12), rep(1:4, 25))
)By this example data, I can know the sample size and response rate for each treatment group. And then in order to obtain the root, someone will use the bisection method but I feel that we can just use the uniroot() to solve the single root, or rootSolve::uniroot.all() to get multiple roots.
p1 <- 0.6
p2 <- 0.2
n1 <- 100
n2 <- 100
rootSolve::uniroot.all(mn_root, interval = c(-0.999, 0.999), n = 100, tol = 1e-6)
## [1] 0.2696618 0.5165744 Then I use the metalite.ae::rate_compare() to check the result from the above self-defined function. Both of the upper and lower values of CI are consistent.
metalite.ae::rate_compare(response ~ treatment, data = ana)
## est z_score p lower upper
## 1 0.4 5.759051 4.229411e-09 0.269662 0.5165743According to the stratified M&N test, I will just talk about it briefly here. If you would like to use the Miettinen-Nurminen stratified method as described in the 1985 paper by SAS, unfortunately the current SAS version cannot support it yet. Although I can use COMMONRISKDIFF(CL=score) to compute the CI of risk difference, the results are not the MN 1985 method that is referred to as Agresti Score at the summary level. More discussions can be found in Confidence limits using the Miettinen-Nurminen stratified method as described in the 1985 paper, and Calculate MN Test Statistics within PROC FREQ.
The basic SAS and R codes for example data used to compute the unstratified and stratified Miettinen and Nurminen test are provided below for reference. It should be emphasized that the results of stratified M&N results are not consistent.
R code:
# Unstratified M&N test
metalite.ae::rate_compare(
formula = response ~ treatment, data = ana
)
# Stratified M&N test
metalite.ae::rate_compare(
formula = response ~ treatment, data = ana,
strata = stratum, weight = "ss"
)
## est z_score p lower upper
## 1 0.3998397 5.712797 5.556727e-09 0.2684383 0.5172779SAS code:
data ana;
input id trtpn stratum response count;
datalines;
1 0 1 0 21
2 0 1 1 5
3 0 2 0 19
4 0 2 1 5
5 0 3 0 21
6 0 3 1 5
7 0 4 0 19
8 0 4 1 5
9 1 1 0 10
10 1 1 1 15
11 1 2 0 10
12 1 2 1 15
13 1 3 0 10
14 1 3 1 15
15 1 4 0 10
16 1 4 1 15
;
run;
ods listing close;
proc freq data=ana;
weight count/zeros;
tables trtpn*response/ riskdiff(CL=mn) alpha=0.05;
ods output PdiffCLs=respdiff;
run;
ods listing;
ods listing close;
proc freq data=ana;
weight count/zeros;
tables stratum*trtpn*response/ riskdiff(CL=mn)
COMMONRISKDIFF(CL=(score) test=(score)) alpha=0.05;
ods output PdiffCLs=respdiff CommonPdiff=compdiff;
run;
ods listing;