In mnel/dsrci: Confidence Intervals on Directly Standardized Rates for Several Parameterizations
Notation
Let $x_1, \ldots, x_k$ be the cell counts with associated random variables $X_1, \ldots, X_n$, $n_1, \ldots, n_k$ be the number of person-years for each cell and assume $X_i~\mathcal{P}\left(\theta_i\right)$.
Let $c_i, \ldots, c_k$ be the number of person-years in a standard population.
Let the weight for the $i$th cell be
$$
w_i = \frac{c_i}{n_i\left(\sum_{j=1}^k{c_j}\right)}
$$
The directly standardized rate (DSR) can be calculated as
$$
y = \sum_{i=1}^k w_ix_i
$$
An estimate of the variance of the DSR is
$$
\upsilon = \sum_{i=1}^k w_i^2 x_i
$$
Let $Y = \sum_{i=1}^k w_i X_i$, and $E\left[Y\right] = \mu = \sum_{i=1}^k w_i \theta_i$.
Constructing confidence intervals
Fay & Feuer (1997), Tiwari et al (2006) and Ng et al (2008) describe and compare a number of confidence interval methods.
Throughout, let $L\left(y\right), U\left(y\right)$ be the lower and upper bounds of a $100(1-\alpha)\%$ confidence interval.
The example data from example(wspoissonTest, package = 'asht')
library(dsrci)
## from example(wspoissonTest, package = 'asht')
## birth data on Down's syndrome from Michigan, 1950-1964
## see Table II of Fay and Feuer (1997)
## xfive= counts for mothers who have had 5 or more children
## nfive and ntotal are number of live births
xfive <- c(0, 8, 63, 112, 262, 295)
nfive <- c(327, 30666, 123419, 149919, 104088, 34392)
ntotal <- c(319933, 931318, 786511, 488235, 237863, 61313)
Asymptotic method ci.method = "asymptotic"
Normal approxmation (trans = "none")
# ci.method = "asymptotic", trans = "none"
# are the default values
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95)
# you can also call ci.asymptotic directly,
# only when weights = weights not standard population counts
# this returns the results without applying a multiplicative
# factor
ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95)
Log transformation (trans = "log")
# ci.method = "asymptotic", trans = "log"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, trans = "log")
# using ci.asymptotic
ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, trans = "log")
Cube root (trans = "cube.root")
Using the delta method,
# ci.method = "asymptotic", trans = "cube.root"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, trans = "cube.root")
# using ci.asymptotic
ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, trans = "cube.root")
Edgeworth correction for skew (trans = "skew")
# ci.method = "asymptotic", trans = "skew"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, trans = "skew")
# using ci.asymptotic
ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, trans = "skew")
Moment matching (ci.method = "moments")
Following Dobson et al (1991), an approximate confidence interval can be obtained as a linear function of the confidence interval for a single Poisson paramenter ($X = \sum_{i=1}^k X_i$) where the confidence interval for this unweighted sum of poisson parameters $\sum_{i=1}^k \theta_i$ is $\left(X_L,X_U \right)$.
An approximate confidence interval for the weighted sum of $\theta$ is
$$
T_L = Y + \sqrt{\frac{\upsilon}{y}}\left(X_L - X\right)
$$
$$
T_U = Y + \sqrt{\frac{\upsilon}{y}}\left(X_U - X\right)
$$
A number of methods for estimating $\left(X_L,X_L\right)$ are implemented in dsrci
Dobson (exact) (type = "dobson")
# ci.method = "moments", type = "dobson"
# note type = "dobson" is the default
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "dobson")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95)
Boise-Monson (type = "boise.monson")
# ci.method = "moments", type = "boise.monson"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "boise.monson")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "boise.monson")
Normal approximation (type = "normal")
# ci.method = "moments", type = "normal"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "normal")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "normal")
Wilson-Hilferty (type = "wilson.hilferty")
# ci.method = "moments", type = "wilson.hilferty"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "wilson.hilferty")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "wilson.hilferty")
Byar's method (type = "byar")
# ci.method = "moments", type = "byar"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "byar")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "byar")
Exact mid-p confidence interval
# ci.method = "moments", type = "midp"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "midp")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "midp")
Approximation to mid-p confidence interval
# ci.method = "moments", type = "approx.midp"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "approx.midp")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "approx.midp")
Simple Approximation to mid-p confidence interval
# ci.method = "moments", type = "approx.midp"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "simple.midp")
# using ci.moments
ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "simple.midp")
Gamma distribution
All Gamma-distribution methods call asht::wspoissonTest
Fay and Feuer (`type = "max")
# ci.method = "gamma", type = "max" (this is the default type)
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "max")
# using ci.gamma
ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "max")
Mid-p (type = "midp")
# ci.method = "gamma", type = "midp"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "midp")
# using ci.gamma
ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "midp")
Tiwari et al 2006 mean (type = "mean")
# ci.method = "gamma", type = "mean"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "mean")
# using ci.gamma
ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "mean")
Tiwari et al 2006 minmaxavg (type = "minmaxavg")
# ci.method = "gamma", type = "minmaxavg"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "minmaxavg")
# using ci.gamma
ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "minmaxavg")
Tiwari et al 2006 (type = "tcz")
# ci.method = "gamma", type = "tcz"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "tcz")
# using ci.gamma
ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "tcz")
Beta Distribution
Tiwari et al 2006 (type = "tcz")
# ci.method = "beta", type = "tcz"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "beta", type = "tcz")
# using ci.beta
ci.beta(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "tcz")
Tiwari et al 2006 Continuity Correction (type = "cc")
# ci.method = "beta", type = "cc"
dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "beta", type = "cc")
# using ci.beta
ci.beta(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "cc")
Comparison of Methods
# a function to loop through methods and return a data.frame
loop_methods <- function(list, level = 0.95, x = xfive, n = nfive, w= ntotal, mult = 1e5, ci.method){
# loop through "types"
results <- do.call(rbind,lapply(list, dsr, level = level, x = x, w = w, n = n, mult = mult, ci.method = ci.method))
# combine into a single data.frame
}
# transformations implemented for "asymptotic"
trans <- c("none", "log", "cube.root", "skew")
asymp.results <- loop_methods(trans, ci.method = "asymptotic")
# types implemented for "moments"
m.types <- c("dobson", "boise.monson", "normal", "wilson.hilferty", "byar",
"midp", "approx.midp", "simple.midp")
moments.results <- loop_methods(m.types, ci.method = "moments")
# types implemented for "gamma"
g.types <- c("max", "midp", "mean", "minmaxavg", "tcz")
gamma.results <- loop_methods(g.types, ci.method = "gamma")
# types implemented for "beta"
b.types <- c("tcz", "cc")
beta.results <- loop_methods(b.types, ci.method = "beta")
# there are no "types" for the approximate bootstrap
# method - no need to loop
boot.results <- dsr(xfive, nfive, ntotal, mult = 1e5, ci.method = 'bootstrap')
# combine into a single data.frame
all.results <- rbind(asymp.results, moments.results,
gamma.results, beta.results,
boot.results)
# show results
knitr::kable(all.results[,c("ci.method", "method.arg", "estimate", "lower", "upper")])
References
Dobson, AJ, Kuulasmaa, K, Eberle, E and Scherer, J (1991)
'Confidence intervals for weighted sums of Poisson parameters',Statistics in Medicine, 10: 457—462.
Swift, MB (1995). 'Simple confidence intervals for standardized rates based on the approximate bootstrap method', Statistics in Medicine, 14, 1875—1888.
Fay & Feuer (1997). 'Confidence intervals for directly standardized rates: a method based on the gamma distribution.' Statistics in Medicine. 16: 791-801.
Tiwari, Clegg, & Zou (2006). 'Efficient interval estimation for age-adjusted cancer rates.' Statistical Methods in Medical Research 15: 547-569.
Ng, Filardo, & Zheng (2008). 'Confidence interval estimating procedures for standardized incidence rates.' Computational Statistics and Data Analysis 52: 3501-3516.
Fay MP (2010). 'Two-sided Exact Tests and Matching Confidence Intervals for Discrete Data'. R Journal 2(1):53-58.
Fay & Kim (2017). 'Confidence intervals for directly standardized rates using mid-p gamma intervals.' Biometrical Journal. doi:10.1002/bimj.201600111.
Fay MP (2017). 'asht: Applied Statistical Hypothesis Tests'. R package version 0.9. https://CRAN.R-project.org/package=asht
mnel/dsrci documentation built on May 23, 2019, 5:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Notation
Let $x_1, \ldots, x_k$ be the cell counts with associated random variables $X_1, \ldots, X_n$, $n_1, \ldots, n_k$ be the number of person-years for each cell and assume $X_i~\mathcal{P}\left(\theta_i\right)$.
Let $c_i, \ldots, c_k$ be the number of person-years in a standard population.
Let the weight for the $i$th cell be
$$ w_i = \frac{c_i}{n_i\left(\sum_{j=1}^k{c_j}\right)} $$ The directly standardized rate (DSR) can be calculated as
$$ y = \sum_{i=1}^k w_ix_i $$ An estimate of the variance of the DSR is
$$ \upsilon = \sum_{i=1}^k w_i^2 x_i $$ Let $Y = \sum_{i=1}^k w_i X_i$, and $E\left[Y\right] = \mu = \sum_{i=1}^k w_i \theta_i$.
Constructing confidence intervals
Fay & Feuer (1997), Tiwari et al (2006) and Ng et al (2008) describe and compare a number of confidence interval methods.
Throughout, let $L\left(y\right), U\left(y\right)$ be the lower and upper bounds of a $100(1-\alpha)\%$ confidence interval.
The example data from example(wspoissonTest, package = 'asht')
library(dsrci) ## from example(wspoissonTest, package = 'asht') ## birth data on Down's syndrome from Michigan, 1950-1964 ## see Table II of Fay and Feuer (1997) ## xfive= counts for mothers who have had 5 or more children ## nfive and ntotal are number of live births xfive <- c(0, 8, 63, 112, 262, 295) nfive <- c(327, 30666, 123419, 149919, 104088, 34392) ntotal <- c(319933, 931318, 786511, 488235, 237863, 61313)
Asymptotic method ci.method = "asymptotic"
Normal approxmation (trans = "none")
# ci.method = "asymptotic", trans = "none" # are the default values dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95) # you can also call ci.asymptotic directly, # only when weights = weights not standard population counts # this returns the results without applying a multiplicative # factor ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95)
Log transformation (trans = "log")
# ci.method = "asymptotic", trans = "log" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, trans = "log") # using ci.asymptotic ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, trans = "log")
Cube root (trans = "cube.root")
Using the delta method,
# ci.method = "asymptotic", trans = "cube.root" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, trans = "cube.root") # using ci.asymptotic ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, trans = "cube.root")
Edgeworth correction for skew (trans = "skew")
# ci.method = "asymptotic", trans = "skew" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, trans = "skew") # using ci.asymptotic ci.asymptotic(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, trans = "skew")
Moment matching (ci.method = "moments")
Following Dobson et al (1991), an approximate confidence interval can be obtained as a linear function of the confidence interval for a single Poisson paramenter ($X = \sum_{i=1}^k X_i$) where the confidence interval for this unweighted sum of poisson parameters $\sum_{i=1}^k \theta_i$ is $\left(X_L,X_U \right)$.
An approximate confidence interval for the weighted sum of $\theta$ is
$$ T_L = Y + \sqrt{\frac{\upsilon}{y}}\left(X_L - X\right) $$ $$ T_U = Y + \sqrt{\frac{\upsilon}{y}}\left(X_U - X\right) $$
A number of methods for estimating $\left(X_L,X_L\right)$ are implemented in dsrci
Dobson (exact) (type = "dobson")
# ci.method = "moments", type = "dobson" # note type = "dobson" is the default dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "dobson") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95)
Boise-Monson (type = "boise.monson")
# ci.method = "moments", type = "boise.monson" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "boise.monson") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "boise.monson")
Normal approximation (type = "normal")
# ci.method = "moments", type = "normal" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "normal") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "normal")
Wilson-Hilferty (type = "wilson.hilferty")
# ci.method = "moments", type = "wilson.hilferty" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "wilson.hilferty") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "wilson.hilferty")
Byar's method (type = "byar")
# ci.method = "moments", type = "byar" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "byar") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "byar")
Exact mid-p confidence interval
# ci.method = "moments", type = "midp" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "midp") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "midp")
Approximation to mid-p confidence interval
# ci.method = "moments", type = "approx.midp" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "approx.midp") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "approx.midp")
Simple Approximation to mid-p confidence interval
# ci.method = "moments", type = "approx.midp" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "moments", type = "simple.midp") # using ci.moments ci.moments(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "simple.midp")
Gamma distribution
All Gamma-distribution methods call asht::wspoissonTest
Fay and Feuer (`type = "max")
# ci.method = "gamma", type = "max" (this is the default type) dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "max") # using ci.gamma ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "max")
Mid-p (type = "midp")
# ci.method = "gamma", type = "midp" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "midp") # using ci.gamma ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "midp")
Tiwari et al 2006 mean (type = "mean")
# ci.method = "gamma", type = "mean" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "mean") # using ci.gamma ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "mean")
Tiwari et al 2006 minmaxavg (type = "minmaxavg")
# ci.method = "gamma", type = "minmaxavg" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "minmaxavg") # using ci.gamma ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "minmaxavg")
Tiwari et al 2006 (type = "tcz")
# ci.method = "gamma", type = "tcz" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "gamma", type = "tcz") # using ci.gamma ci.gamma(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "tcz")
Beta Distribution
Tiwari et al 2006 (type = "tcz")
# ci.method = "beta", type = "tcz" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "beta", type = "tcz") # using ci.beta ci.beta(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "tcz")
Tiwari et al 2006 Continuity Correction (type = "cc")
# ci.method = "beta", type = "cc" dsr(x = xfive, n = nfive, w = ntotal, mult = 1e5, level = 0.95, ci.method = "beta", type = "cc") # using ci.beta ci.beta(xfive, ntotal/(nfive*sum(ntotal)), level = 0.95, type = "cc")
Comparison of Methods
# a function to loop through methods and return a data.frame loop_methods <- function(list, level = 0.95, x = xfive, n = nfive, w= ntotal, mult = 1e5, ci.method){ # loop through "types" results <- do.call(rbind,lapply(list, dsr, level = level, x = x, w = w, n = n, mult = mult, ci.method = ci.method)) # combine into a single data.frame } # transformations implemented for "asymptotic" trans <- c("none", "log", "cube.root", "skew") asymp.results <- loop_methods(trans, ci.method = "asymptotic") # types implemented for "moments" m.types <- c("dobson", "boise.monson", "normal", "wilson.hilferty", "byar", "midp", "approx.midp", "simple.midp") moments.results <- loop_methods(m.types, ci.method = "moments") # types implemented for "gamma" g.types <- c("max", "midp", "mean", "minmaxavg", "tcz") gamma.results <- loop_methods(g.types, ci.method = "gamma") # types implemented for "beta" b.types <- c("tcz", "cc") beta.results <- loop_methods(b.types, ci.method = "beta") # there are no "types" for the approximate bootstrap # method - no need to loop boot.results <- dsr(xfive, nfive, ntotal, mult = 1e5, ci.method = 'bootstrap') # combine into a single data.frame all.results <- rbind(asymp.results, moments.results, gamma.results, beta.results, boot.results)
# show results knitr::kable(all.results[,c("ci.method", "method.arg", "estimate", "lower", "upper")])
References
Dobson, AJ, Kuulasmaa, K, Eberle, E and Scherer, J (1991) 'Confidence intervals for weighted sums of Poisson parameters',Statistics in Medicine, 10: 457—462.
Swift, MB (1995). 'Simple confidence intervals for standardized rates based on the approximate bootstrap method', Statistics in Medicine, 14, 1875—1888.
Fay & Feuer (1997). 'Confidence intervals for directly standardized rates: a method based on the gamma distribution.' Statistics in Medicine. 16: 791-801.
Tiwari, Clegg, & Zou (2006). 'Efficient interval estimation for age-adjusted cancer rates.' Statistical Methods in Medical Research 15: 547-569.
Ng, Filardo, & Zheng (2008). 'Confidence interval estimating procedures for standardized incidence rates.' Computational Statistics and Data Analysis 52: 3501-3516.
Fay MP (2010). 'Two-sided Exact Tests and Matching Confidence Intervals for Discrete Data'. R Journal 2(1):53-58.
Fay & Kim (2017). 'Confidence intervals for directly standardized rates using mid-p gamma intervals.' Biometrical Journal. doi:10.1002/bimj.201600111.
Fay MP (2017). 'asht: Applied Statistical Hypothesis Tests'. R package version 0.9. https://CRAN.R-project.org/package=asht
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.