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$.
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)
ci.method = "asymptotic"
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)
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")
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")
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")
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
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)
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")
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")
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")
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")
# 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")
# 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")
# 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")
All Gamma-distribution methods call asht::wspoissonTest
# 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")
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")
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")
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")
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")
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")
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")
# 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")])
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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.