PF Vignette
David Siev, Christopher H. Tong
January 4, 2012
PFvignette.rmd1 Introduction
The PF package is a collection of functions related to estimating prevented fraction, \(PF=1-RR\) ,
where \(RR={{{\pi}_{2}}}/{{{\pi}_{1}}}\;\) .
1.1 Technical notes
Optimization. Unless otherwise stated, optimization is by the DUD algorithm (Ralston and Jennrich 1978).
Level tested. The help files indicate the level of testing undergone by each function. In some cases that is a subjective judgement, since most of these functions were originally tested in SPlus and have been ported to R more recently.
2 Score based methods
2.1 The score statistic
Confidence intervals for the risk ratio may be based on the score statistic Miettinen and Nurminen (1985), \[\frac{{{{\hat{\pi }}}_{2}}-{{\rho }_{0}}{{{\hat{\pi }}}_{1}}}{\sqrt{\left( {{{{\tilde{\pi }}}_{2}}(1-{{{\tilde{\pi }}}_{2}})}/{{{n}_{2}}}\; \right)+\rho _{0}^{2}\left( {{{{\tilde{\pi }}}_{1}}(1-{{{\tilde{\pi }}}_{1}})}/{{{n}_{1}}}\; \right)}}\] where hat indicates the MLE and tilde indicates the MLE under the restriction that \(\rho ={{\rho }_{0}}\).
2.2 Asymptotic intervals
RRsc() estimates asymptotic confidence intervals for the risk ratio or prevented fraction based on the score statistic.
Interval estimates are returned for three estimators. The score method was originally introduced by (Koopman 1984). Gart and Nam’s modification includes a skewness correction (Gart and Nam 1988). The method of (Miettinen and Nurminen 1985) is a version made slightly more conservative than Koopman’s by including a factor of \({(N-1)}/{N}\;\).
## Loading required package: PF##
## PF
## 95% interval estimates
##
## PF LL UL
## MN method 0.611 0.0251 0.857
## score method 0.611 0.0328 0.855
## skew corr 0.611 0.0380 0.876Starting estimates for the algorithm are obtained by the modified Katz method (log method with 0.5 added to each cell).
Both forms of the Katz estimate may be retrieved from the returned object using RRsc()\$estimate.
2.3 Exact intervals
These methods give intervals that are `exact’ in the sense that they are based on the actual sampling distribution rather than an approximation to it. The score statistic is used to select the \(2 \times 2\) tables that would fall in the tail area, and the binomial probability is estimated over the tail area by taking the maximum over the nuisance parameter. The search over the tail area is made more efficient by the Berger-Boos correction (Berger and Boos 1994).
RRtosst() gives intervals obtained by inverting two one-sided score tests; RRotsst() gives intervals obtained by inverting one two-sided score test. RRtosst() is thus more conservative, preserving at least \({\alpha }/{2}\;\) in each tail.
(Agresti and Min 2001) discuss the properties and relative benefits of the two approaches. The price of exactnesss is conservatism, due to the discreteness of the binomial distribution (Agresti 2001). This means that the actual coverage of the confidence interval does not exactly conform to the nominal coverage, but it will not be less than it. (See also (Agresti 2003).) Both functions use a simple step search algorithm.
##
## PF
## 95% interval estimates
##
## PF LL UL
## 0.6111 0.0148 0.8519##
## PF
## 95% interval estimates
##
## PF LL UL
## 0.611 0.012 0.9023 Stratified designs
Methods for estimating a common RR from stratified or clustered designs depend on homogeneity with respect to the common parameter.
3.1 Gart-Nam score method
(Gart 1985) and (Gart and Nam 1988) derived a score statistic for a common estimate of RR from designs with multiple independent strata, and they showed that it is identical to one proposed by (Radhakrishnan 1965) from a different perspective.
RRstr() provides confidence intervals and a homogeneity test based on Gart’s statistic.
Data may be input two ways, either using a formula and data frame, or as a matrix.
##
## Test of homogeneity across clusters
##
## stat 0.954
## df 3
## p 0.812
##
## RR
## 95% interval estimates
##
## RR LL UL
## starting 2.66 1.37 5.18
## mle 2.65 1.39 5.03
## skew corr 2.65 1.31 5.08
# Data matrix input: RRstr(Y = table6, pf = F)3.2 Mantel-Haenszel estimator
A widely-used heuristic method for sparse frequency tables is the weighted average approach of (Mantel and Haenszel 1959)1. MH interval estimates are based on the asymptotic normality of the log of the risk ratio.
RRmh() utilizes the variance estimator given by (Greenland and Robins 1985) for sparse strata.
The resulting asymptotic estimator is consistent for both the case of sparse strata where the number of strata is assumed increasing,
and the case of limited number of strata where the stratum size is assumed increasing. In the latter case, however, it is less
efficient than maximum likelihood Greenland and Robins (1985).
Additional discussion may be found in Section 4.3.1 (Lachin 2000), (Landis et al. 2005), and (Somes and O’Brien 2006)
2.
##
## RR
## 95% interval estimates
##
## RR LL UL
## 2.67 1.37 5.23
# Data matrix input: RRmh(Y = table6, pf = F)4 Model based intervals
4.1 Logistic regression estimates
Intervals may be estimated from logistic regression models with RRor(). It takes the fit of a glm() object and estimates the intervals by the delta method.
##
## 95% t intervals on 4 df
##
## PF
## PF LL UL
## 0.5000 -0.0406 0.7598
##
## mu.hat LL UL
## txcon 0.733 0.896 0.466
## txvac 0.367 0.624 0.1684.2 Estimating the dispersion parameter
The binomial GLM weights are \[\frac{\hat{\pi }(1-\hat{\pi })}{a(\hat{\varphi })/n\;}\] where \(a(\hat{\varphi })\) is a function of the dispersion parameter.
4.2.1 Dispersion parameter \(\varphi\)
A simple estimator of the dispersion parameter, \(\varphi\), may be estimated by the method of moments [Wed74]. It is given by phiWt().
This form of the dispersion parameter has \(a(\varphi)=\varphi\), and \(\varphi\) is estimated by \({X^2}/{df}\;\), the Pearson statistic divided by the degrees of freedom.
Note that \(\varphi\) is the same estimator as may be obtained by the quasibinomial family in glm() which is, in fact, what is used by phiWt() to reweight the original fit:
phiWt(bird.fit, fit.only = FALSE)$phi## all
## 2.471592## [1] 2.471592phiWt() makes it easy to estimate PF intervals with a single command.
##
## 95% t intervals on 4 df
##
## PF
## PF LL UL
## 0.500 -0.583 0.842
##
## mu.hat LL UL
## txcon 0.733 0.943 0.3121
## txvac 0.367 0.752 0.0997It also allows different estimates of \({\hat{\varphi}}\) for specified subsets of the data.
# model with separate phi for vaccinates and controls
RRor(phiWt(bird.fit, subset.factor = bird$tx))##
## 95% t intervals on 4 df
##
## PF
## PF LL UL
## 0.500 -0.645 0.848
##
## mu.hat LL UL
## txcon 0.733 0.938 0.3330
## txvac 0.367 0.767 0.0925If you want to subtract a degree of freedom for each additional parameter, you can do that by entering the degrees of freedom as an argument to RRor().
##
## 95% t intervals on 2 df
##
## PF
## PF LL UL
## 0.500 -2.164 0.921
##
## mu.hat LL UL
## txcon 0.733 0.975 0.1635
## txvac 0.367 0.895 0.03774.2.2 Dispersion parameter \(\tau\)
When overdispersion is due to intra-cluster correlation, it may make sense to estimate the dispersion as a function of the intra-cluster correlation parameter \(\tau\). In other words, \({a({\varphi }_{ij})}=1+{{\tau }_{j}}({{n}_{ij}}-1)\). tauWt() does this using the Williams procedure (Williams 1982).
##
## 95% t intervals on 4 df
##
## PF
## PF LL UL
## 0.500 -0.645 0.848
##
## mu.hat LL UL
## txcon 0.733 0.938 0.3330
## txvac 0.367 0.767 0.0925In this example the tauWt() estimates are the same as the phiWt() estimates. That is because the cluster sizes are all the same. Let’s see what happens if we modify the bird data set. The birdm data set has the same cluster fractions but differing cluster sizes.
# different cluster sizes, same cluster fractions
birdm.fit <- glm(cbind(y, n - y) ~ tx - 1, binomial, birdm)
RRor(tauWt(birdm.fit, subset.factor = birdm$tx))##
## 95% t intervals on 4 df
##
## PF
## PF LL UL
## 0.490 -0.605 0.838
##
## mu.hat LL UL
## txcon 0.737 0.942 0.328
## txvac 0.376 0.764 0.101Note that increasing cluster size can make things worse when there is intra-cluster correlation.
Now let’s compare the weights from phiWt() and tauWt() with unequal cluster sizes.
In the output below, w represents \(1/a(\hat{\varphi })\;\) and nw is \(n/a(\hat{\varphi })\;\)
# Compare phi and tau weights
#
phi.wts <- phiWt(birdm.fit, fit.only = FALSE, subset.factor = birdm$tx)$weights
tau.wts <- tauWt(birdm.fit, fit.only = FALSE, subset.factor = birdm$tx)$weights
w <- cbind(w.phi = phi.wts,
w.tau = tau.wts,
nw.phi = phi.wts * birdm$n,
nw.tau = tau.wts * birdm$n)
print(cbind(birdm[, c(3, 1, 2)], round(w, 2)), row.names = FALSE)## tx y n w.phi w.tau nw.phi nw.tau
## vac 1 10 0.32 0.35 3.20 3.55
## vac 8 20 0.32 0.21 6.40 4.13
## vac 9 15 0.32 0.26 4.80 3.92
## con 8 16 0.21 0.27 3.39 4.33
## con 8 10 0.21 0.38 2.12 3.82
## con 27 30 0.21 0.16 6.36 4.84Look at the last two rows. Note that the nw.phi are directly proportional to n within treatment group, while the nw.tau are not. With intra-cluster correlation, increasing cluster size does not give a corresponding increase in information.
4.3 Rao-Scott weights
(Rao and Scott 1992) give a method of weighting clustered binomial observations based on the variance inflation due to clustering.
They relate their approach to the concepts of design effect and effective sample size familiar in survey sampling, and they illustrate its use in a variety of contexts.
rsbWt() implements it in the same manner as phiWt() and tauWt(). For more general use, the function rsb() just returns the design effect estimates and the weights.
##
## 95% t intervals on 4 df
##
## PF
## PF LL UL
## 0.479 -0.314 0.793
##
## mu.hat LL UL
## txcon 0.768 0.960 0.311
## txvac 0.400 0.717 0.149
# just the design effect estimates
rsb(birdm$y, birdm$n)$d## all
## 6.3345145 Likelihood based intervals
The RRlsi() function estimates likelihood support intervals for RR by the profile likelihood Section 7.6 (Royall 1997).
Likelihood support intervals are usually formed based on the desired likelihood ratio, \({1}/{k}\;\), often \({1}/{8}\;\) or \({1}/{32}\;\). Under some conditions the log likelihood ratio may follow the chi-square distribution. If so, then \(\alpha =1-{{F}_{{{\chi }^{2}}}}\left( 2\log (k),1 \right)\). RRsc() will make the conversion from \(\alpha\) to \(k\) with the argument use.alpha = T.
##
## 1/8 likelihood support interval for PF
##
## corresponds to 95.858% confidence
## (under certain assumptions)
##
## PF
## PF LL UL
## 0.6111 0.0168 0.8859##
## 1/6.826 likelihood support interval for PF
##
## corresponds to 95% confidence
## (under certain assumptions)
##
## PF
## PF LL UL
## 0.6111 0.0495 0.87926 Incidence ratio
The incidence is the number of cases per subject-time. Its distribution is assumed Poisson. Under certain designs, the incidence ratio (\(IR\)) is used as a measure of treatment effect. Correspondingly, \({{PF}_{IR}}=1-IR\) would be used as a measure of effect for an intervention that is preventive, such as vaccination. \(IR\) is also called incidence density ratio (\(IDR\)), and that is the term used in the following functions.
6.1 Score method
IDRsc() estimates a confidence interval for the incidence density ratio using Siev’s formula Appendix 1 (Siev 1994) based on the Poisson score statistic3.
\[IDR=\widehat{IDR}\left\{ 1+\left( \frac{1}{{{y}_{1}}}+\frac{1}{{{y}_{2}}} \right)\frac{z_{\alpha /2}^{2}}{2}\ \ \pm \ \ \frac{z_{\alpha /2}^{2}}{2{{y}_{1}}{{y}_{2}}}\sqrt{{{y}_{\bullet }}\left( {{y}_{\bullet }}z_{\alpha /2}^{2}+4{{y}_{1}}{{y}_{2}} \right)} \right\}\]
##
## IDR
## 95% interval estimates
##
## IDR LL UL
## 2.61 1.28 5.346.2 Likelihood method
A likelihood support interval for \(IDR\) may be estimated based on orthogonal factoring of the reparameterized likelihood. Section 7.2 (Royall 1997) IDRlsi() implements this method.
Likelihood support intervals are usually formed based on the desired likelihood ratio, \({1}/{k}\;\), often \({1}/{8}\;\) or \({1}/{32}\;\). Under some conditions the log likelihood ratio may follow the chi square distribution. If so, then \(\alpha =1-{{F}_{{{\chi }^{2}}}}\left( 2\log (k),1 \right)\). IDRlsi() will make the conversion from \(\alpha\) tp \(k\) with the argument use.alpha = T.
##
## 1/8 likelihood support interval for IDR
##
## corresponds to 95.858% confidence
## (under certain assumptions)
##
## IDR
## IDR LL UL
## 2.61 1.26 5.88##
## 1/6.826 likelihood support interval for IDR
##
## corresponds to 95% confidence
## (under certain assumptions)
##
## IDR
## IDR LL UL
## 2.61 1.30 5.69