dsrTest - examples

Michael Nelson

2022-06-19

To demonstrate the functionality of dsrTest the various methods to construct confidence interval are demonstrated on the example data drawn from Table 14.4 in Fleiss (1981) (also used in Fay and Feuer (1997)).

The example data come from a study of Down Syndrome and maternal age (Stark and Mantell, 1966).

library(dsrTest)
# the data are in `downs.mi`
data("downs.mi", package = "dsrTest")
# Birth order 5 +
b5 <- downs.mi[downs.mi$BirthOrder == 5, ]
# Gamma Method
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "gamma"))

## 
##  Directly standardized rate: Gamma Method for Weighted Sum of Poissons
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##   67.70212 188.30017
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# Gamma Mid-p
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "gamma",
                 control = list(midp = TRUE)))

## 
##  Directly standardized rate: Confidence Distribution method for
##  Weighted Sum of Poissons
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##   59.71632 173.08170
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# Dobson (exact)
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "dobson"))

## 
##  Directly standardized rate: Dobson's method for Weighted Sum of
##  Poissons with Exact two-sided Poisson test (central method)
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  67.63284 83.86703
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# Dobson (Mid-p)
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "dobson",
                 control = list(midp = TRUE)))

## 
##  Directly standardized rate: Dobson's method for Weighted Sum of
##  Poissons with Exact two-sided Poisson test (central method), mid-p
##  version
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  67.70418 83.79030
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# Asymptotic (no transformation)
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "asymptotic"))

## 
##  Directly standardized rate: Asymptotic method for Weighted Sum of
##  Poissons normal approximation of MLE
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  67.48907 83.56895
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# Asymptotic (log transformation)
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "asymptotic",
                 control = list(trans = "log")))

## 
##  Directly standardized rate: Asymptotic method for Weighted Sum of
##  Poissons normal approximation of the log-transformed MLE
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  67.90220 84.01246
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# Beta Method
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "beta"))

## 
##  Directly standardized rate: Beta method for Weighted Sum of Poissons
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  67.70196 83.77768
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# Approximate Bootstrap Method
with(b5, dsrTest(Cases, Births, Standard, mult = 1e5, method = "bootstrap"))

## 
##  Directly standardized rate: Approximate bootstrap method for Weighted
##  Sum of Poissons
## 
## data:  x = Cases, time bases: n = Births, weights: w = Standard
## number of summands = 6, Multiplier = 1e+05, p-value = NA
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  68.35707 84.61566
## sample estimates:
## Directly standardized rate per 1e+05 
##                             75.52901

# A list of methods to implement
methods_list <- list(
  gamma = list(
    list(wmtype = "max"),
    list(midp = TRUE),
    list(wmtype = "tcz"),
    list(wmtype = "mean"),
    list(wmtype = "minmaxavg")),
  asymptotic = list(
    list(trans = "none"),
    list(trans = "log"),
    list(trans = "loglog"),
    list(trans = "logit")),
  dobson = list(
    list(midp = FALSE),
    list(midp = TRUE)),
  beta = list(
    list(wmtype = "none"),
    list(wmtype = "tcz"),
    list(wmtype = "mean"),
    list(wmtype = "minmaxavg"),
    list(wmtype = "max")),
  bootstrap = list(list())
)
# split out to allow call to mapply
methods <-rep(names(methods_list), times = lengths(methods_list))
controls <- do.call(c, unname(methods_list))
all_methods <- mapply(dsrTest,
  method = methods, control = controls,
  MoreArgs = list(mult = 1e5, x = b5$Cases, n = b5$Births, w = b5$Standard),
  SIMPLIFY = FALSE)
# create some "short" names
control_types <- unlist(controls)
control_names <- c(gsub("midp=FALSE", "Exact CI",
  gsub("=TRUE", "",
       sprintf("[%s=%s]", names(control_types), control_types))), "")
names(all_methods) <- paste(methods, control_names)
# combine CI into one data.frame
results <- do.call(rbind,lapply(all_methods,
  function(data) data.frame(
    estimate = data$estimate,
    lci = data$conf.int[1], 
    uci = data$conf.int[2])))
# and display
knitr::kable(results, digits = 3)