Package 'epiR'

Summary measures for count data presented in a 2 by 2 table

Description

Computes summary measures of risk and a chi-squared test for difference in the observed proportions from count data presented in a 2 by 2 table. With multiple strata the function returns crude and Mantel-Haenszel adjusted measures of association and chi-squared tests of homogeneity.

Usage

epi.2by2(dat, method = "cohort.count", digits = 2, conf.level = 0.95, 
   units = 100, interpret = FALSE, outcome = "as.columns")

## S3 method for class 'epi.2by2'
print(x, ...)

## S3 method for class 'epi.2by2'
summary(object, ...)

Arguments

dat	a vector of length four, an object of class table or an object of class grouped_df from package dplyr containing the individual cell frequencies (see below).
method	a character string indicating the study design on which the tabular data has been based. Options are cohort.count, cohort.time, case.control, or cross.sectional. Based on the study design specified by the user, appropriate measures of association, measures of effect in the exposed and measures of effect in the population are returned by the function.
digits	scalar, number of digits to be reported for print output. Must be an integer of either 2, 3 or 4.
conf.level	magnitude of the returned confidence intervals. Must be a single number between 0 and 1.
units	multiplier for prevalence and incidence (risk or rate) estimates.
interpret	logical. If TRUE interpretive statements are appended to the printepi.2by2 object.
outcome	a character string indicating how the outcome variable is represented in the contingency table. Options are as.columns (outcome as columns) or as.rows (outcome as rows).
x, object	an object of class epi.2by2.
...	Ignored.

Details

Where method is cohort.count, case.control, or cross.sectional and outcome = as.columns the required 2 by 2 table format is:

-----------	----------	----------	----------
	Disease +	Disease -	Total
-----------	----------	----------	----------
Expose +	a	b	a+b
Expose -	c	d	c+d
-----------	----------	----------	----------
Total	a+c	b+d	a+b+c+d
-----------	----------	----------	----------

Where method is cohort.time and outcome = as.columns the required 2 by 2 table format is:

-----------	----------	-------------
	Disease +	Time at risk
-----------	----------	-------------
Expose +	a	b
Expose -	c	d
-----------	----------	-------------
Total	a+c	b+d
-----------	----------	-------------

A summary of the methods used for each of the confidence interval calculations in this function is as follows:

Value

An object of class epi.2by2 comprised of:

method	character string returning the study design specified by the user.
n.strata	number of strata.
conf.level	magnitude of the returned confidence intervals.
interp	logical. Are interpretative statements included?
units	character string listing the outcome measure units.
tab	a data frame comprised of of the contingency table data.
massoc.summary	a data frame listing the computed measures of association, measures of effect in the exposed and measures of effect in the population and their confidence intervals.
massoc.interp	a data frame listing the interpretive statements for each computed measure of association.
massoc.detail	a list comprised of the computed measures of association, measures of effect in the exposed and measures of effect in the population. See below for details.

When method equals cohort.count the following measures of association, measures of effect in the exposed and measures of effect in the population are returned:

RR	Wald, Taylor and score confidence intervals for the incidence risk ratios for each strata. Wald, Taylort and score confidence intervals for the crude incidence risk ratio. Wald confidence interval for the Mantel-Haenszel adjusted incidence risk ratio.
OR	Wald, score, Cornfield and maximum likelihood confidence intervals for the odds ratios for each strata. Wald, score, Cornfield and maximum likelihood confidence intervals for the crude odds ratio. Wald confidence interval for the Mantel-Haenszel adjusted odds ratio.
ARisk	Wald and score confidence intervals for the attributable risk (risk difference) for each strata. Wald and score confidence intervals for the crude attributable risk. Wald, Sato and Greenland-Robins confidence intervals for the Mantel-Haenszel adjusted attributable risk.
NNT	Wald and score confidence intervals for the number needed to treat for benefit (NNTB) or number needed to treat for harm (NNTH).
PARisk	Wald and Pirikahu confidence intervals for the population attributable risk for each strata. Wald and Pirikahu confidence intervals for the crude population attributable risk. The Pirikahu confidence intervals are calculated using the delta method.
AFRisk	Wald confidence intervals for the attributable fraction for each strata. Wald confidence intervals for the crude attributable fraction.
PAFRisk	Wald confidence intervals for the population attributable fraction for each strata. Wald confidence intervals for the crude population attributable fraction.
chisq.strata	chi-squared test for difference in exposed and non-exposed proportions for each strata.
chisq.crude	chi-squared test for difference in exposed and non-exposed proportions across all strata.
chisq.mh	Mantel-Haenszel chi-squared test that the combined odds ratio estimate is equal to 1.
RR.homog	Mantel-Haenszel (Woolf) test of homogeneity of the individual strata incidence risk ratios.
OR.homog	Mantel-Haenszel (Woolf) test of homogeneity of the individual strata odds ratios.

When method equals cohort.time the following measures of association and effect are returned:

IRR	Wald confidence interval for the incidence rate ratios for each strata. Wald confidence interval for the crude incidence rate ratio. Wald confidence interval for the Mantel-Haenszel adjusted incidence rate ratio.
ARate	Wald confidence interval for the attributable rate for each strata. Wald confidence interval for the crude attributable rate. Wald confidence interval for the Mantel-Haenszel adjusted attributable rate.
PARate	Wald confidence interval for the population attributable rate for each strata. Wald confidence intervals for the crude population attributable rate.
AFRate	Wald confidence interval for the attributable fraction for each strata. Wald confidence interval for the crude attributable fraction.
PAFRate	Wald confidence interval for the population attributable fraction for each strata. Wald confidence interval for the crude poulation attributable fraction.
chisq.strata	chi-squared test for difference in exposed and non-exposed proportions for each strata.
chisq.crude	chi-squared test for difference in exposed and non-exposed proportions across all strata.
chisq.mh	Mantel-Haenszel chi-squared test that the combined odds ratio estimate is equal to 1.

When method equals case.control the following measures of association and effect are returned:

OR	Wald, score, Cornfield and maximum likelihood confidence intervals for the odds ratios for each strata. Wald, score, Cornfield and maximum likelihood confidence intervals for the crude odds ratio. Wald confidence interval for the Mantel-Haenszel adjusted odds ratio.
ARisk	Wald and score confidence intervals for the attributable risk for each strata. Wald and score confidence intervals for the crude attributable risk. Wald, Sato and Greenland-Robins confidence intervals for the Mantel-Haenszel adjusted attributable risk.
PARisk	Wald and Pirikahu confidence intervals for the population attributable risk for each strata. Wald and Pirikahu confidence intervals for the crude population attributable risk.
AFest	Wald confidence intervals for the estimated attributable fraction for each strata. Wald confidence intervals for the crude estimated attributable fraction.
PAFest	Wald confidence intervals for the population estimated attributable fraction for each strata. Wald confidence intervals for the crude population estimated attributable fraction.
chisq.strata	chi-squared test for difference in exposed and non-exposed proportions for each strata.
chisq.crude	chi-squared test for difference in exposed and non-exposed proportions across all strata.
chisq.mh	Mantel-Haenszel chi-squared test that the combined odds ratio estimate is equal to 1.
OR.homog	Mantel-Haenszel (Woolf) test of homogeneity of the individual strata odds ratios.

When method equals cross.sectional the following measures of association and effect are returned:

PR	Wald, Taylor and score confidence intervals for the prevalence ratios for each strata. Wald, Taylor and score confidence intervals for the crude prevalence ratio. Wald confidence interval for the Mantel-Haenszel adjusted prevalence ratio.
OR	Wald, score, Cornfield and maximum likelihood confidence intervals for the odds ratios for each strata. Wald, score, Cornfield and maximum likelihood confidence intervals for the crude odds ratio. Wald confidence interval for the Mantel-Haenszel adjusted odds ratio.
ARisk	Wald and score confidence intervals for the attributable risk for each strata. Wald and score confidence intervals for the crude attributable risk. Wald, Sato and Greenland-Robins confidence intervals for the Mantel-Haenszel adjusted attributable risk.
NNT	Wald and score confidence intervals for the number needed to treat for benefit (NNTB) or number needed to treat for harm (NNTH).
PARisk	Wald and Pirikahu confidence intervals for the population attributable risk for each strata. Wald and Pirikahu confidence intervals for the crude population attributable risk.
AFRisk	Wald confidence intervals for the attributable fraction for each strata. Wald confidence intervals for the crude attributable fraction.
PAFRisk	Wald confidence intervals for the population attributable fraction for each strata. Wald confidence intervals for the crude population attributable fraction.
chisq.strata	chi-squared test for difference in exposed and non-exposed proportions for each strata.
chisq.crude	chi-squared test for difference in exposed and non-exposed proportions across all strata.
chisq.mh	Mantel-Haenszel chi-squared test that the combined odds ratio estimate is equal to 1.
PR.homog	Mantel-Haenszel (Woolf) test of homogeneity of the individual strata prevalence ratios.
OR.homog	Mantel-Haenszel (Woolf) test of homogeneity of the individual strata odds ratios.

The point estimates of the wald, score and cfield odds ratios are calculated using the cross product method. Method mle computes the conditional maximum likelihood estimate of the odds ratio.

Confidence intervals for the Cornfield (cfield) odds ratios are computed using the hypergeometric distribution and computation times are slow when the cell frequencies are large. For this reason, Cornfield confidence intervals are only calculated if the total number of event frequencies is less than 500. Maximum likelihood estimates of the odds ratio and Fisher's exact test are only calculated when the total number of observations is less than 2E09.

If the Haldane-Anscombe (Haldane 1940, Anscombe 1956) correction is applied (i.e., addition of 0.5 to each cell of the 2 by 2 table when at least one of the cell frequencies is zero) Cornfield (cfield) odds ratios are not computed.

Variable phi.coef equals the phi coefficient (Fleiss et al. 2003, Equation 6.2, p. 98) and is included with the output for each of the uncorrected chi-squared tests. This value can be used for argument rho.cc in epi.sscc. Refer to the documentation for epi.sscc for details.

The Mantel-Haenszel chi-squared test that the combined odds ratio estimate is equal to 1 uses a two-sided test without continuity correction.

Interpretive statements for NNTB and NNTH follow the approach described by Altman (1998). See the examples for details. Note that number needed to treat to benefit (NNTB) and number needed to treat to harm (NNTH) estimates are not computed when method = "cohort.time" or method = "case.control" because attributable risk can't be calculated using these study designs.

Note

For cohort studies where the outcome of interest is expressed as an incidence risk, measures of association include the incidence risk ratio and the incidence odds ratio. When the incidence of the outcome in the study population is low (less than 5%) the incidence odds ratio will provide a reliable estimate of the incidence risk ratio. The more frequent the outcome becomes, the more the incidence odds ratio will overestimate the incidence risk ratio when it is greater than than 1 or understimate the incidence risk ratio when it is less than 1.

For cohort studies where outcome of interest is expressed as an incidence rate, the incidence rate ratio is used as a measure of association. For case-control studies the exposure odds ratio is used as a measure of association.

For cross-sectional studies where the outcome of interest is expressed as a prevalence, measures of association include the prevalence risk ratio and the prevalence odds ratio. When the prevalence of the outcome in the study population is low (less than 5%) the prevalence odds ratio will provide a reliable estimate of the incidence risk ratio. The more frequent the outcome becomes, the more the prevalence odds ratio will overestimate the prevalence risk ratio when it is greater than than 1 or understimate the prevalence risk ratio when it is less than 1.

Measures of effect in the exposed include the attributable risk (for cohort studies where the outcome of interest is expressed as an incidence risk), attributable rate (for cohort studies where the outcome of interest is expressed as an incidence rate) and prevalence risk (for cross-sectional studies where the outcome of interest is expressed as a prevalence risk). The attributable risk is the incidence risk of the outcome in the exposed minus the incidence risk of the outcome in the unexposed. The attributable risk provides a measure of the absolute increase or decrease in outcome risk that is associated with exposure. The attributable fraction is the proportion of study outcomes in the exposed group that is attributable to exposure.

The number needed to treat (NNT) is the inverse of the attributable risk. Depending on the outcome of interest we use different labels for NNT. When dealing with an outcome that is desirable (e.g., treatment success) we call NNT the number needed to treat for benefit, NNTB. NNTB equals the number of subjects who would have to be exposed to result in a single (desirable) outcome. When dealing with an outcome that is undesirable (e.g., death) we call NNT the number needed to treat for harm, NNTH. NNTH equals the number of subjects who would have to be exposed to result in a single (undesirable) outcome.

Measures of effect in the population include the attributable risk in the population (for cohort studies where the outcome of interest is expressed as an incidence risk) and the attributable prevalence in the population (for cross-sectional studies where the outcome of interest is expressed as a prevalence risk). The population attributable risk provides a measure of the absolute increase or decrease in outcome risk in the population that is associated with exposure. The population attributable fraction is the proportion of study outcomes in the populaton that is attributable to exposure.

Point estimates and confidence intervals for the prevalence ratio and incidence risk ratio are calculated using the Wald (1943) and score methods (Miettinen and Nurminen 1985). Point estimates and confidence intervals for the incidence rate ratio are calculated using the exact method described by Kirkwood and Sterne (2003) and Juul (2004). Point estimates and confidence intervals for the odds ratio are calculated using the Wald (1943), score (Miettinen and Nurminen 1985) and maximum likelihood methods (Fleiss et al. 2003). Point estimates and confidence intervals for the population attributable risk are calculated using formulae provided by Lash et al (2021) and Pirikahu (2014). Point estimates and confidence intervals for the population attributable fraction are calculated using formulae provided by Jewell (2004, p 84 - 85). Point estimates and confidence intervals for the Mantel-Haenszel adjusted attributable risk are calculated using formulae provided by Klingenberg (2014).

Wald confidence intervals are provided in the summary table simply because they are widely used and would be familiar to most users.

The Mantel-Haenszel adjusted measures of association are valid when the measures of association across the different strata are similar (homogenous), that is when the test of homogeneity of the odds (risk) ratios is not statistically significant.

The Mantel-Haenszel (Woolf) test of homogeneity of the odds ratio are based on Jewell (2004, p 152 - 158). Thanks to Jim Robison-Cox for sharing his implementation of these functions.

Author(s)

Mark Stevenson (Melbourne Veterinary School, Faculty of Science, The University of Melbourne, Australia), Cord Heuer (EpiCentre, IVABS, Massey University, Palmerston North, New Zealand), Jim Robison-Cox (Department of Math Sciences, Montana State University, Montana, USA), Kazuki Yoshida (Brigham and Women's Hospital, Boston Massachusetts, USA) and Simon Firestone (Melbourne Veterinary School, Faculty of Science, The University of Melbourne, Australia). Thanks to Ian Dohoo for numerous helpful suggestions to improve the documentation for this function.

References

Altman D (1998). British Medical Journal 317, 1309 - 1312.

Altman D, Machin D, Bryant T, Gardner M (2000). Statistics with Confidence. British Medical Journal, London, pp. 69.

Anscombe F (1956). On estimating binomial response relations. Biometrika 43, 461 - 464.

Cornfield, J (1956). A statistical problem arising from retrospective studies. In: Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, University of California Press, Berkeley California 4: 135 - 148.

Elwood JM (2007). Critical Appraisal of Epidemiological Studies and Clinical Trials. Oxford University Press, London.

Feinstein AR (2002). Principles of Medical Statistics. Chapman Hall/CRC, London, pp. 332 - 336.

Fisher RA (1962). Confidence limits for a cross-product ratio. Australian Journal of Statistics 4: 41.

Feychting M, Osterlund B, Ahlbom A (1998). Reduced cancer incidence among the blind. Epidemiology 9: 490 - 494.

Fleiss JL, Levin B, Paik MC (2003). Statistical Methods for Rates and Proportions. John Wiley and Sons, New York.

Haldane J (1940). The mean and variance of the moments of chi square, when used as a test of homogeneity, when expectations are small. Biometrika 29, 133 - 143.

Hanley JA (2001). A heuristic approach to the formulas for population attributable fraction. Journal of Epidemiology and Community Health 55: 508 - 514.

Hightower AW, Orenstein WA, Martin SM (1988) Recommendations for the use of Taylor series confidence intervals for estimates of vaccine efficacy. Bulletin of the World Health Organization 66: 99 - 105.

Jewell NP (2004). Statistics for Epidemiology. Chapman & Hall/CRC, London, pp. 84 - 85.

Juul S (2004). Epidemiologi og evidens. Munksgaard, Copenhagen.

Kirkwood BR, Sterne JAC (2003). Essential Medical Statistics. Blackwell Science, Malden, MA, USA.

Klingenberg B (2014). A new and improved confidence interval for the Mantel-Haenszel risk difference. Statistics in Medicine 33: 2968 - 2983.

Lancaster H (1961) Significance tests in discrete distributions. Journal of the American Statistical Association 56: 223 - 234.

Lash TL, VanderWeele TJ, Haneuse S, Rothman KJ (2021). Modern Epidemiology. Lippincott - Raven Philadelphia, USA, pp. 79 - 103.

Lawson R (2004). Small sample confidence intervals for the odds ratio. Communications in Statistics Simulation and Computation 33: 1095 - 1113.

Martin SW, Meek AH, Willeberg P (1987). Veterinary Epidemiology Principles and Methods. Iowa State University Press, Ames, Iowa, pp. 130.

McNutt L, Wu C, Xue X, Hafner JP (2003). Estimating the relative risk in cohort studies and clinical trials of common outcomes. American Journal of Epidemiology 157: 940 - 943.

Miettinen OS, Nurminen M (1985). Comparative analysis of two rates. Statistics in Medicine 4: 213 - 226.

Pirikahu S (2014). Confidence Intervals for Population Attributable Risk. Unpublished MSc thesis. Massey University, Palmerston North, New Zealand.

Robbins AS, Chao SY, Fonesca VP (2002). What's the relative risk? A method to directly estimate risk ratios in cohort studies of common outcomes. Annals of Epidemiology 12: 452 - 454.

Sullivan KM, Dean A, Soe MM (2009). OpenEpi: A Web-based Epidemiologic and Statistical Calculator for Public Health. Public Health Reports 124: 471 - 474.

Wald A (1943). Tests of statistical hypotheses concerning several parameters when the number of observations is large. Transactions of the American Mathematical Society 54: 426 - 482.

Willeberg P (1977). Animal disease information processing: Epidemiologic analyses of the feline urologic syndrome. Acta Veterinaria Scandinavica. Suppl. 64: 1 - 48.

Woodward M (2014). Epidemiology Study Design and Data Analysis. Chapman & Hall/CRC, New York, pp. 89 - 124.

Zhang J, Yu KF (1998). What's the relative risk? A method for correcting the odds ratio in cohort studies of common outcomes. Journal of the American Medical Association 280: 1690 - 1691.

Examples

## EXAMPLE 1:
## A cross sectional study investigating the relationship between dry cat 
## food (DCF) and feline urologic syndrome (FUS) was conducted (Willeberg 
## 1977). Counts of individuals in each group were as follows:

## DCF-exposed cats (cases, non-cases) 13, 2163
## Non DCF-exposed cats (cases, non-cases) 5, 3349

## Outcome variable (FUS) as columns:
dat.v01 <- c(13,2163,5,3349); dat.v01

epi.2by2(dat = dat.v01, method = "cross.sectional", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")

## Outcome variable (FUS) as rows:
dat.v01 <- c(13,5,2163,3349); dat.v01

epi.2by2(dat =  dat.v01, method = "cross.sectional", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.rows")

## The prevalence of FUS in DCF exposed cats was 4.01 (95% CI 1.43 to 11.23) 
## times greater than the prevalence of FUS in non-DCF exposed cats.

## In DCF exposed cats, 75% (95% CI 30% to 91%) of the FUS cases were 
## attributable to DCF.

## Fifty-four percent of FUS cases in the population was attributable 
## to DCF (95% CI 4% to 78%).


## EXAMPLE 2:
## This example shows how the table function in base R can be used to pass 
## data to epi.2by2. Here we use the birthwt data set from the MASS package.

library(MASS)
dat.df02 <- birthwt; head(dat.df02)

## Generate a table of cell frequencies. First, set the outcome and exposure 
## as factors and set their levels appropriately so the frequencies in the 
## 2 by 2 table come out in the conventional format:
dat.df02$low <- factor(dat.df02$low, levels = c(1,0))
dat.df02$smoke <- factor(dat.df02$smoke, levels = c(1,0))
dat.df02$race <- factor(dat.df02$race, levels = c(1,2,3))
dat.tab02 <- table(dat.df02$smoke, dat.df02$low, dnn = c("Smoke", "Low BW"))
print(dat.tab02)

## Compute the odds ratio and other measures of association:
epi.2by2(dat = dat.tab02, method = "cohort.count", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")

## The odds of having a low birth weight child for smokers was 2.02 
## (95% CI 1.08 to 3.78) times greater than the odds of having a low birth 
## weight child for non-smokers.

## Stratify by race:
dat.tab02 <- table(dat.df02$smoke, dat.df02$low, dat.df02$race, 
   dnn = c("Smoke", "Low BW", "Race"))
print(dat.tab02)

## Compute the crude odds ratio, the Mantel-Haenszel adjusted odds ratio 
## and other measures of association:
dat.epi02 <- epi.2by2(dat = dat.tab02, method = "cohort.count", digits = 2,
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
print(dat.epi02)

## The Mantel-Haenszel test of homogeneity of the strata odds ratios is not 
## statistically significant (chi square test statistic 2.800; df 2; 
## p-value = 0.25). We accept the null hypothesis and conclude that the 
## odds ratios for each strata of race are the same. 

## After accounting for the confounding effect of race, the odds of 
## having a low birth weight child for smokers was 3.09 (95% CI 1.49 to 6.39)
## times that of non-smokers.

## Compare the Greenland-Robins confidence intervals for the Mantel-Haenszel
## adjusted attributable risk with the Wald confidence intervals for the 
## Mantel-Haenszel adjusted attributable risk:

dat.epi02$massoc.detail$ARisk.mh.green
dat.epi02$massoc.detail$ARisk.mh.wald

## How many mothers need to stop smoking to avoid one low birth weight baby? 
dat.epi02$massoc.interp$text[dat.epi02$massoc.interp$var == 
   "NNTB NNTH (crude)"]
   
## Exposure (smoking) increased the crude outcome incidence risk in the exposed. 
## The number needed to treat (i.e., expose) to increase the outcome 
## frequency by one was 7 (95% CI 3 to 62).

dat.epi02$massoc.interp$text[dat.epi02$massoc.interp$var == "NNTB NNTH (M-H)"]

## After adjusting for confounding exposure (smoking) increased the outcome 
## incidence risk in the exposed. The number needed to treat (i.e., expose) 
## to increase the outcome frequency by one was 5 (95% CI 2 to 71).

## Now turn dat.tab02 into a data frame where the frequencies of individuals in 
## each exposure-outcome category are provided. Often your data will be 
## presented in this summary format:
dat.df02 <- data.frame(dat.tab02); head(dat.df02)

## Re-format dat.df02 (a summary count data frame) into tabular format using 
## the xtabs function:
dat.tab02 <- xtabs(Freq ~ Smoke + Low.BW + Race, data = dat.df02)
print(dat.tab02)

# dat02.tab can now be passed to epi.2by2:
dat.epi02 <- epi.2by2(dat = dat.tab02, method = "cohort.count", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
print(dat.epi02)

## The Mantel-Haenszel adjusted incidence odds ratio is 3.09 (95% CI 1.49 to 
## 6.39). The ratio of the crude incidene odds ratio to the Mantel-Haensel 
## adjusted incidence odds ratio is 0.66.

## What are the Cornfield confidence limits, the maximum likelihood 
## confidence limits and the score confidence limits for the crude 
## incidence odds ratio?
dat.epi02$massoc.detail$OR.crude.cfield
dat.epi02$massoc.detail$OR.crude.mle
dat.epi02$massoc.detail$OR.crude.score

## Cornfield: 2.02 (95% CI 1.07 to 3.79)
## Maximum likelihood: 2.01 (1.03 to 3.96)
## Score: 2.02 (95% CI 1.08 to 3.77)

## Plot the individual strata-level incidence odds ratios and compare them 
## with the Mantel-Haenszel adjusted incidence odds ratio.

## Not run: 
library(ggplot2); library(scales)

nstrata <- 1:dim(dat.tab02)[3]
strata.lab <- paste("Strata ", nstrata, sep = "")
y.at <- c(nstrata, max(nstrata) + 1)
y.lab <- c("M-H", strata.lab)
x.at <- c(0.25,0.5,1,2,4,8,16,32)

or.p <- c(dat.epi02$massoc.detail$OR.mh$est, 
   dat.epi02$massoc.detail$OR.strata.cfield$est)
or.l <- c(dat.epi02$massoc.detail$OR.mh$lower, 
   dat.epi02$massoc.detail$OR.strata.cfield$lower)
or.u <- c(dat.epi02$massoc.detail$OR.mh$upper, 
   dat.epi02$massoc.detail$OR.strata.cfield$upper)
dat.df02 <- data.frame(y.at, y.lab, or.p, or.l, or.u)

ggplot(data = dat.df02, aes(x = or.p, y = y.at)) +
  theme_bw() +
  geom_point() + 
  geom_errorbarh(aes(xmax = or.l, xmin = or.u, height = 0.2)) + 
  labs(x = "Odds ratio", y = "Strata") + 
  scale_x_continuous(trans = log2_trans(), breaks = x.at, 
     limits = c(0.25,32)) + 
  scale_y_continuous(breaks = y.at, labels = y.lab) + 
  geom_vline(xintercept = 1, lwd = 1) + 
  coord_fixed(ratio = 0.75 / 1) + 
  theme(axis.title.y = element_text(vjust = 0))

## End(Not run)


## EXAMPLE 3:
## Same as Example 2 but showing how a 2 by 2 contingency table can be prepared 
## using functions from the tidyverse package:

## Not run: 
library(MASS); library(tidyverse)

dat.df03 <- birthwt; head(dat.df03)

dat.df03 <- dat.df03 %>%
  mutate(low = factor(low, levels = c(1,0), labels = c("yes","no"))) %>%
  mutate(smoke = factor(smoke, levels = c(1,0), labels = c("yes","no"))) %>%
  group_by(smoke, low) %>%
  summarise(n = n()) 
dat.df03

## View the data in conventional 2 by 2 table format:
pivot_wider(dat.df03, id_cols = c(smoke), names_from = low, values_from = n)

dat.epi03 <- epi.2by2(dat = dat.df03, method = "cohort.count", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
dat.epi03

## End(Not run)

## The incidence odds of having a low birth weight child for smokers was 2.02 
## (95% CI 1.08 to 3.78) times greater than the incidence odds of having a 
## low birth weight child for non-smokers.

## Stratify by race:
## Not run: 
library(MASS); library(tidyverse)

dat.df04 <- birthwt; head(dat.df04)

dat.df04 <- dat.df04 %>%
  mutate(low = factor(low, levels = c(1,0), labels = c("yes","no"))) %>%
  mutate(smoke = factor(smoke, levels = c(1,0), labels = c("yes","no"))) %>%
  mutate(race = factor(race)) %>%
  group_by(race, smoke, low) %>%
  summarise(n = n())
dat.df04

## View the data in conventional 2 by 2 table format:
pivot_wider(dat.df04, id_cols = c(race, smoke),
   names_from = low, values_from = n)

dat.epi04 <- epi.2by2(dat = dat.df04, method = "cohort.count", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
dat.epi04

## End(Not run)

## The Mantel-Haenszel test of homogeneity of the strata odds ratios is not 
## statistically significant (chi square test statistic 2.800; df 2; 
## p-value = 0.25). We accept the null hypothesis and conclude that the 
## incidence odds ratios for each strata of race are the same. 

## After accounting for the confounding effect of race, the incidence odds of 
## having a low birth weight child for smokers was 3.09 (95% CI 1.49 to 6.39)
## times that of non-smokers.


## EXAMPLE 4:
## Sometimes you'll have only event count data for a stratified analysis. This
## example shows how to coerce a three column matrix listing (in order) counts
## of outcome positive individuals, counts of outcome negative individuals (or
## total time at risk, as in the example below) and strata into a three 
## dimensional array.

## We assume that two rows are recorded for each strata. The first for those 
## exposed and the second for those unexposed. Note that the strata variable
## needs to be numeric (not a factor).

dat.m04 <- matrix(c(1308,884,200,190,4325264,13142619,1530342,5586741,1,1,2,2), 
   nrow = 4, ncol = 3, byrow = FALSE)
colnames(dat.m04) <- c("obs","tar","grp")
dat.df04 <- data.frame(dat.m04)

## Here we use the apply function to coerce the two rows for each strata into
## tabular format. An array is created of with the length of the third 
## dimension of the array equal to the number of strata: 
dat.tab04 <- sapply(1:length(unique(dat.df04$grp)), function(x) 
   as.matrix(dat.df04[dat.df04$grp == x,1:2], ncol = 2, byrow = TRUE), 
   simplify = "array")
dat.tab04

epi.2by2(dat = dat.tab04, method = "cohort.time", digits = 2, 
   conf.level = 0.95, units = 1000 * 365.25, interpret = FALSE, 
   outcome = "as.columns")

## The Mantel-Haenszel adjusted incidence rate ratio was 4.39 (95% CI 4.06 
## to 4.75).


## EXAMPLE 5:
## A study was conducted by Feychting et al (1998) comparing cancer occurrence
## among the blind with occurrence among those who were not blind but had 
## severe visual impairment. From these data we calculate a cancer rate of
## 136/22050 person-years among the blind compared with 1709/127650 person-
## years among those who were visually impaired but not blind.

dat.v05 <- c(136,22050,1709,127650)

dat.epi05 <- epi.2by2(dat = dat.v05, method = "cohort.time", digits = 2,
   conf.level = 0.95, units = 1000, interpret = FALSE, outcome = "as.columns")
summary(dat.epi05)$massoc.detail$ARate.strata.wald

## The incidence rate of cancer was 7.22 (95% CI 6.00 to 8.43) cases per 
## 1000 person-years less in the blind, compared with those who were not 
## blind but had severe visual impairment.

round(summary(dat.epi05)$massoc.detail$IRR.strata.wald, digits = 2)

## The incidence rate of cancer in the blind group was less than half that 
## of the comparison group (incidence rate ratio 0.46, 95% CI 0.38 to 0.55).


## EXAMPLE 6:
## A study has been conducted to assess the effect of a new treatment for 
## mastitis in dairy cows. Eight herds took part in the study. The following 
## data were obtained. The vectors ai, bi, ci and di list (for each herd) the 
## number of cows in the E+D+, E+D-, E-D+ and E-D- groups, respectively.

## Not run: 
hid <- 1:8
ai <- c(23,10,20,5,14,6,10,3)
bi <- c(10,2,1,2,2,2,3,0)
ci <- c(3,2,3,2,1,3,3,2)
di <- c(6,4,3,2,6,3,1,1)
dat.df06 <- data.frame(hid, ai, bi, ci, di)
head(dat.df06)

## Re-format data into a format suitable for epi.2by2:
hid <- rep(1:8, times = 4)
exp <- factor(rep(c(1,1,0,0), each = 8), levels = c(1,0))
out <- factor(rep(c(1,0,1,0), each = 8), levels = c(1,0))
dat.df06 <- data.frame(hid, exp, out, n = c(ai,bi,ci,di)) 
dat.tab06 <- xtabs(n ~ exp + out + hid, data = dat.df06)
print(dat.tab06)

epi.2by2(dat = dat.tab06, method = "cohort.count", digits = 2, 
   conf.level = 0.95, units = 1000, interpret = FALSE, outcome = "as.columns")

## The Mantel-Haenszel test of homogeneity of the strata incidence odds ratios 
## is not statistically significant (chi square test statistic 5.276; df 7; 
## p-value = 0.63). We accept the null hypothesis and conclude that the 
## incidence odds ratios for each strata of herd are the same. 

## After adjusting for the effect of herd, compared to untreated cows, treatment
## increased the incidence odds of recovery by a factor of 5.97 
## (95% CI 2.72 to 13.13).
 
## End(Not run)

---

The library epiR: summary information

Description

Tools for the analysis of epidemiological data.

Usage

epi.about()

Details

More information about the epiR package can be found at https://mvs.unimelb.edu.au/research/groups/veterinary-epidemiology-melbourne and https://www.ausvet.com.au/.

FUNCTIONS AND DATASETS

The following is a summary of the main functions and datasets in the epiR package. An alphabetical list of all functions and datasets is available by typing library(help = epiR).

For further information on any of these functions, type help(name) or ?name where name is the name of the function or dataset.

For details on how to use epiR for routine epidemiological work start R, type help.start() to open the help browser and navigate to Packages > epiR > Vignettes.

CONTENTS:

The functions in epiR can be categorised into two main groups: tools for epidemiological analysis and tools for the analysis of surveillance data. A summary of the package functions is as follows:

I. EPIDEMIOLOGY

1. Descriptive statistics

epi.conf	Confidence intervals.
epi.descriptives	Descriptive statistics.

2. Measures of health and measures of association

epi.directadj	Directly adjusted incidence rate estimates.
epi.edr	Compute estimated dissemination ratios from outbreak event data.
epi.empbayes	Empirical Bayes estimates of observed event counts.
epi.indirectadj	Indirectly adjusted incidence risk estimates.
epi.insthaz	Instantaneous hazard estimates based on Kaplan-Meier survival estimates.
epi.2by2	Measures of association from data presented in a 2 by 2 table.

3. Diagnostic tests

epi.betabuster	An R version of Wes Johnson and Chun-Lung Su's Betabuster.
epi.herdtest	Estimate the characteristics of diagnostic tests applied at the herd (group) level.
epi.nomogram	Compute the post-test probability of disease given characteristics of a diagnostic test.
epi.pooled	Estimate herd test characteristics when samples are pooled.
epi.prev	Compute the true prevalence of a disease in a population on the basis of an imperfect test.
epi.tests	Sensitivity, specificity and predictive value of a diagnostic test.

4. Meta-analysis

epi.dsl	Mixed-effects meta-analysis of binary outcome data using the DerSimonian and Laird method.
epi.iv	Fixed-effects meta-analysis of binary outcome data using the inverse variance method.
epi.mh	Fixed-effects meta-analysis of binary outcome data using the Mantel-Haenszel method.
epi.smd	Fixed-effects meta-analysis of continuous outcome data using the standardised mean difference method.

5. Regression analysis tools

epi.cp	Extract unique covariate patterns from a data set.
epi.cpresids	Compute covariate pattern residuals from a logistic regression model.
epi.interaction	Relative excess risk due to interaction in a case-control study.

6. Data manipulation tools

epi.asc	Write matrix to an ASCII raster file.
epi.convgrid	Convert British National Grid georeferences to easting and northing coordinates.
epi.dms	Convert decimal degrees to degrees, minutes and seconds and vice versa.
epi.ltd	Calculate lactation to date and standard lactation (that is, 305 or 270 day) milk yields.
epi.offset	Create an offset vector based on a list suitable for WinBUGS.
epi.RtoBUGS	Write data from an R list to a text file in WinBUGS-compatible format.

7. Sample size calculations

The naming convention for the sample size functions in epiR is: epi.ss (sample size) + an abbreviation to represent the sampling design (e.g., simple, strata, clus1, clus2) + an abbreviation of the objectives of the study (est when you want to estimate a population parameter or comp when you want to compare two groups) + a single letter defining the outcome variable type (b for binary, c for continuous and s for survival data).

epi.sssimpleestb	Sample size to estimate a binary outcome using simple random sampling.
epi.sssimpleestc	Sample size to estimate a continuous outcome using simple random sampling.
epi.ssstrataestb	Sample size to estimate a binary outcome using stratified random sampling.
epi.ssstrataestc	Sample size to estimate a continuous outcome using stratified random sampling.
epi.ssclus1estb	Sample size to estimate a binary outcome using one-stage cluster sampling.
epi.ssclus1estc	Sample size to estimate a continuous outcome using one-stage cluster sampling.
epi.ssclus2estb	Sample size to estimate a binary outcome using two-stage cluster sampling.
epi.ssclus2estc	Sample size to estimate a continuous outcome using two-stage cluster sampling.
epi.ssxsectn	Sample size, power or detectable prevalence ratio for a cross-sectional study.
epi.sscohortc	Sample size, power or detectable risk ratio for a cohort study using count data.
epi.sscohortt	Sample size, power or detectable risk ratio for a cohort study using time at risk data.
epi.sscc	Sample size, power or detectable odds ratio for case-control studies.
epi.sscompb	Sample size, power and detectable risk ratio when comparing binary outcomes.
epi.sscompc	Sample size, power and detectable risk ratio when comparing continuous outcomes.
epi.sscomps	Sample size, power and detectable hazard when comparing time to event.
epi.ssequb	Sample size for a parallel equivalence or equality trial, binary outcome.
epi.ssequc	Sample size for a parallel equivalence or equality trial, continuous outcome.
epi.sssupb	Sample size for a parallel superiority trial, binary outcome.
epi.sssupc	Sample size for a parallel superiority trial, continuous outcome.
epi.ssninfb	Sample size for a non-inferiority trial, binary outcome.
epi.ssninfc	Sample size for a non-inferiority trial, continuous outcome.
epi.ssdetect	Sample size to detect an event.
epi.ssdxsesp	Sample size to estimate the sensitivity or specificity of a diagnostic test.
epi.ssdxtest	Sample size to validate a diagnostic test in the absence of a gold standard.

8. Miscellaneous functions

epi.prcc	Compute partial rank correlation coefficients.
epi.psi	Compute proportional similarity indices.
epi.realrisk	Return absolute risks from odds, incidence risk and hazard ratios.

9. Data sets

epi.epidural	Rates of use of epidural anaesthesia in trials of caregiver support.
epi.incin	Laryngeal and lung cancer cases in Lancashire 1974 - 1983.
epi.SClip	Lip cancer in Scotland 1975 - 1980.

II. SURVEILLANCE

Below, SSe stands for surveillance system sensitivity. That is, the average probability that a surveillance system (as a whole) will return a positive surveillance outcome, given disease is present in the population at a level equal to or greater than a specified design prevalence.

1. Representative sampling — sample size

rsu.sspfree.rs	Defined probability of disease freedom.
rsu.sssep.rs	SSe, perfect test specificity.
rsu.sssep.rs2st	SSe, two stage sampling.
rsu.sssep.rsfreecalc	SSe, imperfect test specificity.
rsu.sssep.rspool	SSe, pooled sampling.

2. Representative sampling — surveillance system sensitivity and specificity

rsu.sep.rs	SSe, representative sampling.
rsu.sep.rs2st	SSe, representative two-stage sampling.
rsu.sep.rsmult	SSe, representative multiple surveillance components.
rsu.sep.rsfreecalc	SSe, imperfect test specificity.
rsu.sep.rspool	SSe, representative pooled sampling.
rsu.sep.rsvarse	SSe, varying surveillance unit sensitivity.
rsu.spp.rs	Surveillance system specificity.

3. Representative sampling — probability of disease freedom

rsu.pfree.rs	Probability of disease freedom for a single or multiple time periods.
rsu.pfree.equ	Equilibrium probability of disease freedom.

4. Risk-based sampling — sample size

rsu.sssep.rbsrg	SSe, single sensitivity for each risk group.
rsu.sssep.rbmrg	SSe, multiple sensitivities within risk groups.
rsu.sssep.rb2st1rf	SSe, 2 stage sampling, 1 risk factor.
rsu.sssep.rb2st2rf	SSe, 2 stage sampling, 2 risk factors.

5. Risk-based sampling — surveillance system sensitivity and specificity

rsu.sep.rb	SSe, risk-based sampling.
rsu.sep.rb1rf	SSe, risk-based sampling, 1 risk factor.
rsu.sep.rb2rf	SSe, risk-based sampling, 2 risk factors.
rsu.sep.rbvarse	SSe, risk-based sampling, varying unit sensitivity.
rsu.sep.rb2st	SSe, 2-stage risk-based sampling.

6. Risk-based sampling — probability of disease freedom

rsu.pfree.equ	Equilibrium probability of disease freedom.

7. Census sampling — surveillance system sensitivity

rsu.sep.cens	SSe, census sampling.

8. Passive surveillance — surveillance system sensitivity

rsu.sep.pass	SSe, passive surveillance.

9. Miscellaneous functions

rsu.adjrisk	Adjusted risk values.
rsu.dxtest	Series and parallel diagnostic test interpretation.
rsu.epinf	Effective probability of disease.
rsu.pstar	Design prevalence back calculation.
rsu.sep	Probability disease is less than specified design prevalence.

Author(s)

Mark Stevenson ([email protected]), Melbourne Veterinary School, Faculty of Science, The University of Melbourne, Parkville Victoria 3010, Australia.

Evan Sergeant ([email protected]), Ausvet Pty Ltd, Level 1 34 Thynne St, Bruce ACT 2617, Australia.

Simon Firestone, Melbourne Veterinary School, Faculty of Science, The University of Melbourne, Parkville Victoria 3010, Australia.

Telmo Nunes, UISEE/DETSA, Faculdade de Medicina Veterinaria — UTL, Rua Prof. Cid dos Santos, 1300 - 477 Lisboa Portugal.

Javier Sanchez, Atlantic Veterinary College, University of Prince Edward Island, Charlottetown Prince Edward Island, C1A 4P3, Canada.

Ron Thornton, Ministry for Primary Industries New Zealand, PO Box 2526 Wellington, New Zealand.

With contributions from: Cord Heuer, Jonathon Marshall, Jeno Reiczigel, Jim Robison-Cox, Paola Sebastiani, Peter Solymos, Yoshida Kazuki, Geoff Jones, Sarah Pirikahu, Ryan Kyle, Johann Popp, Methew Jay, Allison Cheung, Nagendra Singanallur, Aniko Szabo and Ahmad Rabiee.

---

Write matrix to an ASCII raster file

Description

Writes a data frame to an ASCII raster file, suitable for display in a Geographic Information System.

Usage

epi.asc(dat, file, xllcorner, yllcorner, cellsize, na = -9999)

Arguments

dat	a matrix with data suitable for plotting using the image function.
file	character string specifying the name and path of the ASCII raster output file.
xllcorner	the easting coordinate corresponding to the lower left hand corner of the matrix.
yllcorner	the northing coordinate corresponding to the lower left hand corner of the matrix.
cellsize	number, defining the size of each matrix cell.
na	scalar, defines null values in the matrix. NAs are converted to this value.

Value

Writes an ASCII raster file (typically with *.asc extension), suitable for display in a Geographic Information System.

Note

The image function in R rotates tabular data counter clockwise by 90 degrees for display. A matrix of the form:

1	3
2	4

is displayed (using image) as:

3	4
1	2

It is recommended that the source data for this function is a matrix. Replacement of NAs in a data frame extends processing time for this function.

---

An R version of Wes Johnson and Chun-Lung Su's Betabuster

Description

A function to return shape1 and shape2 parameters for a beta distribution, based on expert elicitation.

Usage

epi.betabuster(mode, conf, imsure, x, conf.level = 0.95, max.shape1 = 100, 
   step = 0.001)

Arguments

mode	scalar, the mode of the variable of interest. Must be a number between 0 and 1.
conf	level of confidence (expressed on a 0 to 1 scale) that the true value of the variable of interest is greater or less than argument x.
imsure	a character string, if "greater than" you are making the statement that you are conf confident that the true value of the variable of interest is greater than x. If "less than" you are making the statement that you are conf confident that the true value of the variable of interest is less than x.
x	scalar, value of the variable of interest (see above).
conf.level	magnitude of the returned confidence interval for the estimated beta distribution. Must be a single number between 0 and 1.
max.shape1	scalar, maximum value of the shape1 parameter for the beta distribution.
step	scalar, step value for the shape1 parameter. See details.

Details

The beta distribution has two parameters: shape1 and shape2, corresponding to a and b in the original version of BetaBuster. If r equals the number of times an event has occurred after n trials, shape1 = (r + 1) and shape2 = (n - r + 1).

Take care when you're parameterising probability estimates that are at the extremes of the 0 to 1 bounds. If the returned shape1 parameter is equal to the value of max.shape1 (which, by default is 100) consider increasing the value of the max.shape1 argument. The epi.betabuster functions issues a warning if these conditions are met.

Value

A list containing the following:

shape1	the shape1 parameter for the estimated beta distribution.
shape2	the shape2 parameter for the estimated beta distribution.
mode	the mode of the estimated beta distribution.
mean	the mean of the estimated beta distribution.
median	the median of the estimated beta distribution.
lower	the lower bound of the confidence interval of the estimated beta distribution.
upper	the upper bound of the confidence interval of the estimated beta distribution.
variance	the variance of the estimated beta distribution.
exp	a statement of the arguments used for this instance of the function.

Author(s)

Simon Firestone (Melbourne Veterinary School, Faculty of Science, The University of Melbourne, Parkville Victoria 3010, Australia) with acknowledgements to Wes Johnson and Chun-Lung Su for the original standalone software.

References

Christensen R, Johnson W, Branscum A, Hanson TE (2010). Bayesian Ideas and Data Analysis: An Introduction for Scientists and Statisticians. Chapman and Hall, Boca Raton.

Su C-L, Johnson W (2014) Beta Buster. Software for obtaining parameters for the Beta distribution based on expert elicitation. URL: https://cadms.vetmed.ucdavis.edu/diagnostic/software.

Examples

## EXAMPLE 1:
## If a scientist is asked for their best guess for the diagnostic sensitivity
## of a particular test and the answer is 0.90, and if they are also willing 
## to assert that they are 80% certain that the sensitivity is greater than 
## 0.75, what are the shape1 and shape2 parameters for a beta distribution
## satisfying these constraints? 

rval.beta01 <- epi.betabuster(mode = 0.90, conf = 0.80, imsure = "greater than", 
   x = 0.75, conf.level = 0.95, max.shape1 = 100, step = 0.001)
rval.beta01$shape1; rval.beta01$shape2

## The shape1 and shape2 parameters for the beta distribution that satisfy the
## constraints listed above are 9.875 and 1.986, respectively.

## This beta distribution reflects the probability distribution obtained if 
## there were 9 successes, r:
r <- rval.beta01$shape1 - 1; r

## from 10 trials, n:
n <- rval.beta01$shape2 + rval.beta01$shape1 - 2; n

dat.df01 <- data.frame(x = seq(from = 0, to = 1, by = 0.001), 
   y = dbeta(x = seq(from = 0, to = 1,by = 0.001), 
   shape1 = rval.beta01$shape1, shape2 = rval.beta01$shape2))

## Density plot of the estimated beta distribution:

## Not run: 
library(ggplot2)

ggplot(data = dat.df01, aes(x = x, y = y)) +
  theme_bw() +
  geom_line() + 
  scale_x_continuous(name = "Test sensitivity") +
  scale_y_continuous(name = "Density")

## End(Not run)


## EXAMPLE 2:
## The most likely value of the specificity of a PCR for coxiellosis in 
## small ruminants is 1.00 and we're 97.5% certain that this estimate is 
## greater than 0.99. What are the shape1 and shape2 parameters for a beta 
## distribution satisfying these constraints?

epi.betabuster(mode = 1.00, conf = 0.975, imsure = "greater than", x = 0.99, 
   conf.level = 0.95, max.shape1 = 100, step = 0.001)

## The shape1 and shape2 parameters for the beta distribution that satisfy the
## constraints listed above are 100 and 1, respectively. epi.betabuster 
## issues a warning that the value of shape1 equals max.shape1. Increase
## max.shape1 to 500:

epi.betabuster(mode = 1.00, conf = 0.975, imsure = "greater than", x = 0.99, 
   conf.level = 0.95, max.shape1 = 500, step = 0.001)

## The shape1 and shape2 parameters for the beta distribution that satisfy the
## constraints listed above are 367.04 and 1, respectively.

---

Number of parameters to be inferred and number of informative priors required for a Bayesian latent class model

Description

Returns the number of unknown parameters to be inferred and the number of informative priors likely to be needed for an identifiable Bayesian latent class model to estimate diagnostic sensitivity and specificity in the absence of a gold standard.

Usage

epi.blcm.paras(ntest.dep = 2, ntest.indep = 1, npop = 2)

Arguments

ntest.dep	scalar, the number of conditionally dependent tests evaluated.
ntest.indep	scalar, the number of conditionally independent tests evaluated.
npop	scalar, the number of populations with a distinct prevalence investigated.

Value

A list containing the following:

df	scalar, the degrees of freedom in the available data.
npars	scalar, the number of unknown parameters to be inferred.
ninf.priors	scalar, the number of informative priors likely to be needed for an identifiable model.

Note

A model may still be useful for inference if it has less informative priors, though cautious interpretation is warranted, typically with a sensitivity analysis of the influence of the priors on the findings.

Author(s)

Simon Firestone and Allison Cheung (Faculty of Veterinary and Agricultural Sciences, The University of Melbourne, Parkville Victoria, Australia), Nagendra Singanallur (Australian Centre for Disease Preparedness, Newcomb Victoria, Australia).

Examples

## EXAMPLE 1 --- Two conditionally dependent tests, 1 population:
epi.blcm.paras(ntest.dep = 2, ntest.indep = 0, npop = 1)

## This model has 3 degrees of freedom. The model has 7 unknown parameters to
## be inferred. At least 4 informative priors are required.


## EXAMPLE 2 --- Two conditionally dependent tests, 2 populations:
epi.blcm.paras(ntest.dep = 2, ntest.indep = 0, npop = 2)

## This model has 6 degrees of freedom. The model has 8 unknown parameters to
## be inferred. At least 2 informative priors are required.


## EXAMPLE 3 --- Two conditionally dependent tests, 3 populations:
epi.blcm.paras(ntest.dep = 2, ntest.indep = 0, npop = 3)

## This model has 9 degrees of freedom. The model has 9 unknown parameters to
## be inferred. This model may be able to proceed without informative priors.


## EXAMPLE 4 --- Two conditionally dependent tests, 1 independent test, 1 
## population:
epi.blcm.paras(ntest.dep = 2, ntest.indep = 1, npop = 1)

## This model has 7 degrees of freedom. The model has 9 unknown parameters to
## be inferred. At least 2 informative priors are required.


## EXAMPLE 5 --- Two conditionally dependent tests, 1 independent test, 2 
## populations:
epi.blcm.paras(ntest.dep = 2, ntest.indep = 1, npop = 2)

## This model has 14 degrees of freedom. The model has 10 unknown parameters to
## be inferred. This model may be able to proceed without informative priors.


## EXAMPLE 6 --- Three conditionally dependent tests, 1 population:
epi.blcm.paras(ntest.dep = 3, ntest.indep = 0, npop = 1)

## This model has 7 degrees of freedom. The model has 13 unknown parameters to
## be inferred. At least 6 informative priors are required.


## EXAMPLE 7 --- Three conditionally dependent tests, 2 populations:
epi.blcm.paras(ntest.dep = 3, ntest.indep = 0, npop = 2)

## This model has 14 degrees of freedom. The model has 14 unknown parameters to
## be inferred. This model may be able to proceed without informative priors.

---

Bohning's test for overdispersion of Poisson data

Description

A test for overdispersion of Poisson data.

Usage

epi.bohning(obs, exp, alpha = 0.05)

Arguments

obs	the observed number of cases in each area.
exp	the expected number of cases in each area.
alpha	alpha level to be used for the test of significance. Must be a single number between 0 and 1.

Value

A data frame with two elements: test.statistic, Bohning's test statistic and p.value the associated P-value.

References

Bohning D (2000). Computer-assisted Analysis of Mixtures and Applications. Chapman and Hall, Boca Raton.

Ugarte MD, Ibanez B, Militino AF (2006). Modelling risks in disease mapping. Statistical Methods in Medical Research 15: 21 - 35.

Examples

## EXAMPLE 1:
data(epi.SClip)
obs <- epi.SClip$cases
pop <- epi.SClip$population
exp <- (sum(obs) / sum(pop)) * pop

epi.bohning(obs, exp, alpha = 0.05)

## Bohning's test was used to determine if there was statistically significant
## overdispersion in lip cancer cases across 56 Scottish districts for the 
## period 1975 to 1980. 

## The test statistic was 53.33. The associated P value was <0.01. We reject 
## the null hypothesis of no over dispersion and accept the null hypothesis
## concluding that the lip cancer data are over dispersed.

---

Concordance correlation coefficient

Description

Calculates Lin's (1989, 2000) concordance correlation coefficient for agreement on a continuous measure.

Usage

epi.ccc(x, y, ci = "z-transform", conf.level = 0.95, rep.measure = FALSE, 
   subjectid)

Arguments

x	a vector, representing the first set of measurements.
y	a vector, representing the second set of measurements.
ci	a character string, indicating the method to be used. Options are z-transform or asymptotic.
conf.level	magnitude of the returned confidence interval. Must be a single number between 0 and 1.
rep.measure	logical. If TRUE there are repeated observations across subject.
subjectid	a factor providing details of the observer identifier if rep.measure == TRUE.

Details

Computes Lin's (1989, 2000) concordance correlation coefficient for agreement on a continuous measure obtained by two methods. The concordance correlation coefficient combines measures of both precision and accuracy to determine how far the observed data deviate from the line of perfect concordance (that is, the line at 45 degrees on a square scatter plot). Lin's coefficient increases in value as a function of the nearness of the data's reduced major axis to the line of perfect concordance (the accuracy of the data) and of the tightness of the data about its reduced major axis (the precision of the data).

Both x and y values need to be present for a measurement pair to be included in the analysis. If either or both values are missing (i.e., coded NA) then the measurement pair is deleted before analysis.

Value

A list containing the following:

rho.c	the concordance correlation coefficient.
s.shift	the scale shift.
l.shift	the location shift.
C.b	a bias correction factor that measures how far the best-fit line deviates from a line at 45 degrees. No deviation from the 45 degree line occurs when C.b = 1. See Lin (1989, page 258).
blalt	a data frame with two columns: mean the mean of each pair of measurements, delta vector y minus vector x.
sblalt	a data frame listing the average difference between the two sets of measurements, the standard deviation of the difference between the two sets of measurements and the lower and upper confidence limits of the difference between the two sets of measurements. If rep.measure == TRUE the confidence interval of the difference is adjusted to account for repeated observations across individual subjects.
nmissing	a count of the number of measurement pairs ignored due to missingness.

References

Bland J, Altman D (1986). Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet 327: 307 - 310.

Bland J, Altman D (1999). Measuring agreement in method comparison studies. Statistical Methods in Medical Research 8: 135 - 160.

Bland J, Altman D (2007). Agreement between methods of measurement with multiple observations per individual. Journal of Biopharmaceutical Statistics 17: 571 - 582. (Corrects the formula quoted in the 1999 paper).

Bradley E, Blackwood L (1989). Comparing paired data: a simultaneous test for means and variances. American Statistician 43: 234 - 235.

Burdick RK, Graybill FA (1992). Confidence Intervals on Variance Components. New York: Dekker.

Dunn G (2004). Statistical Evaluation of Measurement Errors: Design and Analysis of Reliability Studies. London: Arnold.

Euser AM, Dekker FW, le Cessie S (2008). A practical approach to Bland-Altman plots and variation coefficients for log transformed variables. Journal of Clinical Epidemiology 61: 978 - 982.

Hsu C (1940). On samples from a normal bivariate population. Annals of Mathematical Statistics 11: 410 - 426.

Krippendorff K (1970). Bivariate agreement coefficients for reliability of data. In: Borgatta E, Bohrnstedt G (eds) Sociological Methodology. San Francisco: Jossey-Bass, pp. 139 - 150.

Lin L (1989). A concordance correlation coefficient to evaluate reproducibility. Biometrics 45: 255 - 268.

Lin L (2000). A note on the concordance correlation coefficient. Biometrics 56: 324 - 325.

Pitman E (1939). A note on normal correlation. Biometrika 31: 9 - 12.

Rashid M, Stevenson M, Waenga S, Mirams G, Campbell A, Vaughan J, Jabbar A (2018). Comparison of McMaster and FECPAK methods for counting nematode eggs in the faeces of alpacas. Parasites & Vectors 11, 278. DOI: 10.1186/s13071-018-2861-1.

Reynolds M, Gregoire T (1991). Comment on Bradley and Blackwood. American Statistician 45: 163 - 164.

Snedecor G, Cochran W (1989). Statistical Methods. Ames: Iowa State University Press.

See Also

epi.occc

Examples

## EXAMPLE 1:
set.seed(seed = 1234)
method1 <- rnorm(n = 100, mean = 0, sd = 1)
method2 <- method1 + runif(n = 100, min = -0.25, max = 0.25)

## Add some missing values:
method1[50] <- NA
method2[75] <- NA

dat.df01 <- data.frame(method1, method2)
rval.ccc01 <- epi.ccc(method1, method2, ci = "z-transform", conf.level = 0.95, 
   rep.measure = FALSE)

rval.lab01 <- data.frame(lab = paste("CCC: ", 
   round(rval.ccc01$rho.c[,1], digits = 2), " (95% CI ", 
   round(rval.ccc01$rho.c[,2], digits = 2), " - ",
   round(rval.ccc01$rho.c[,3], digits = 2), ")", sep = ""))

z <- lm(method2 ~ method1)
alpha <- summary(z)$coefficients[1,1]
beta <-  summary(z)$coefficients[2,1]
rval.lm01 <- data.frame(alpha, beta)

## Concordance correlation plot:
## Not run: 
library(ggplot2)

ggplot(data = dat.df01, aes(x = method1, y = method2)) +
  theme_bw() +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  geom_abline(data = rval.lm01, aes(intercept = alpha, slope = beta), 
     linetype = "dashed") +
  scale_x_continuous(limits = c(0,3), name = "Method 1") +
  scale_y_continuous(limits = c(0,3), name = "Method 2") +
  geom_text(data = rval.lab01, x = 0.5, y = 2.95, label = rval.lab01$lab) + 
  coord_fixed(ratio = 1 / 1)

## In this plot the dashed line represents the line of perfect concordance. 
## The solid line represents the reduced major axis.  

## End(Not run)


## EXAMPLE 2:
## Bland and Altman plot (Figure 2 from Bland and Altman 1986):
x <- c(494,395,516,434,476,557,413,442,650,433,417,656,267,
   478,178,423,427)

y <- c(512,430,520,428,500,600,364,380,658,445,432,626,260,
   477,259,350,451)

rval.ccc02 <- epi.ccc(x, y, ci = "z-transform", conf.level = 0.95, 
   rep.measure = FALSE)
rval.df02 <- data.frame(mean = rval.ccc02$blalt[,1], 
   delta = rval.ccc02$blalt[,2])

## Not run: 
library(ggplot2)

ggplot(data = rval.ccc02$blalt, aes(x = mean, y = delta)) +
  theme_bw() + 
  geom_point() +
  geom_hline(data = rval.ccc02$sblalt, aes(yintercept = lower), linetype = 2) +
  geom_hline(data = rval.ccc02$sblalt, aes(yintercept = upper), linetype = 2) +
  geom_hline(data = rval.ccc02$sblalt, aes(yintercept = est), linetype = 1) +
  scale_x_continuous(limits = c(0,800), 
     name = "Average PEFR by two meters (L/min)") +
  scale_y_continuous(limits = c(-150,150), 
     name = "Difference in PEFR (L/min)")   

## End(Not run)


## EXAMPLE 3:
## Setting limits of agreement when your data are skewed. See Euser et al. 
## (2008) for details and Rashid et al. (2018) for an applied example.
x <- c(0,210,15,90,0,0,15,0,0,0,15,0,15,0,0,0,0,15,0,0,15,135,0,0,15,
   120,30,15,30,0,0,5235,780,1275,10515,1635,1905,1830,720,450,225,420,
   300,15,285,0,225,525,675,5280,465,270,0,1485,15,420,0,60,0,0,0,750,
   570,0)
y <- c(0,70,0,0,0,0,35,0,0,0,0,0,0,35,0,0,0,0,0,0,35,35,70,0,0,140,35,
   105,0,0,0,1190,385,1190,6930,560,1260,700,840,0,105,385,245,35,105,
   0,140,350,350,3640,385,350,0,1505,0,630,70,0,0,140,0,420,490,0)   
   
rval.ccc03 <- epi.ccc(x, y, ci = "z-transform", 
   conf.level = 0.95, rep.measure = FALSE)

## Not run: 
library(ggplot2)

ggplot(data = rval.ccc03$blalt, aes(x = mean, y = delta)) +
  theme_bw() + 
  geom_point() +
  geom_hline(data = rval.ccc03$sblalt, aes(yintercept = lower), linetype = 2) +
  geom_hline(data = rval.ccc03$sblalt, aes(yintercept = upper), linetype = 2) +
  geom_hline(data = rval.ccc03$sblalt, aes(yintercept = est), linetype = 1) +
  scale_x_continuous(limits = c(0,8000), 
     name = "Average of the two measurements") +
  scale_y_continuous(limits = c(-8000,8000), 
     name = "Difference in the two measurements")  

## End(Not run)

## In the above plot the spread of the differences increases with increasing 
## mean of the observations. The Bland Altman limits of agreement should be 
## calculated on a log scale.

logx <- log(x + 50, base = 10)
logy <- log(y + 50, base = 10)

log10.ccc03 <- epi.ccc(x = logx, y = logy, ci = "z-transform", 
   conf.level = 0.95, rep.measure = FALSE)

## Transform the limits of agreement back to the original scale by taking 
## anti-logs. If the limits of agreement for Z = log10(x) are between -a 
## and +a, with a = 1.96 * s, the ratio between two measures on the original 
## scale is between 10^-a and 10^a. See page 979 of Euser et al. (2008). 

a <- 1.96 * log10.ccc03$sblalt$delta.sd

## For a given value for the mean Xbar, it can be shown that x - y is between 
## -2Xbar(10^a - 1) / (10^a + 1) and +2Xbar(10^a - 1) / (10^a + 1):

Xbar = seq(from = 0, to = 8000, by = 100)
Xbar.low <- (-2 * Xbar * (10^a - 1)) / (10^a + 1)
Xbar.upp <- (+2 * Xbar * (10^a - 1)) / (10^a + 1)
limits <- data.frame(mean = Xbar, lower = Xbar.low, upper = Xbar.upp)

## Not run: 
library(ggplot2)ggplot(

data = rval.ccc03$blalt, aes(x = mean, y = delta)) +
  theme_bw() +
  geom_point() +
  geom_line(data = limits, aes(x = mean, y = lower), linetype = 2) +
  geom_line(data = limits, aes(x = mean, y = upper), linetype = 2) +
  geom_line(data = limits, aes(x = mean, y = 0), linetype = 1) +
  scale_x_continuous(limits = c(0,8000), 
     name = "Average of the two measurements") +
  scale_y_continuous(limits = c(-8000,8000), 
     name = "Difference in the two measurements") 

## End(Not run)

---

Confidence intervals for means, proportions, incidence, and standardised mortality ratios

Description

Computes confidence intervals for means, proportions, incidence, and standardised mortality ratios.

Usage

epi.conf(dat, ctype = "mean.single", method, N, design = 1,
   conf.level = 0.95)

Arguments

dat	the data, either a vector or a matrix depending on the method chosen.
ctype	a character string indicating the type of confidence interval to calculate. Options are mean.single, mean.unpaired, mean.paired, prop.single, prop.unpaired, prop.paired, prevalence, inc.risk, inc.rate, odds, ratio and smr.
method	a character string indicating the method to use. Where ctype = "inc.risk" or ctype = "prevalence" options are exact, wilson, fleiss, agresti, clopper-pearson and jeffreys. Where ctype = "inc.rate" options are exact and byar.
N	scalar, representing the population size.
design	scalar, representing the design effect.
conf.level	magnitude of the returned confidence interval. Must be a single number between 0 and 1.

Details

Method mean.single requires a vector as input. Method mean.unpaired requires a two-column data frame; the first column defining the groups must be a factor. Method mean.paired requires a two-column data frame; one column for each group. Method prop.single requires a two-column matrix; the first column specifies the number of positives, the second column specifies the number of negatives. Methods prop.unpaired and prop.paired require a four-column matrix; columns 1 and 2 specify the number of positives and negatives for the first group, columns 3 and 4 specify the number of positives and negatives for the second group. Method prevalence and inc.risk require a two-column matrix; the first column specifies the number of positives, the second column specifies the total number tested. Method inc.rate requires a two-column matrix; the first column specifies the number of positives, the second column specifies individual time at risk. Method odds requires a two-column matrix; the first column specifies the number of positives, the second column specifies the number of negatives. Method ratio requires a two-column matrix; the first column specifies the numerator, the second column specifies the denominator. Method smr requires a two-colum matrix; the first column specifies the total number of positives, the second column specifies the total number tested.

The methodology implemented here follows Altman, Machin, Bryant, and Gardner (2000). Where method is inc.risk or prevalence if the numerator equals zero the lower bound of the confidence interval estimate is set to zero. Where method is smr the method of Dobson et al. (1991) is used. A summary of the methods used for each of the confidence interval calculations in this function is as follows:

-----------------------------------	------------------------
ctype-method	Reference
-----------------------------------	------------------------
mean.single	Altman et al. (2000)
mean.unpaired	Altman et al. (2000)
mean.paired	Altman et al. (2000)
prop.single	Altman et al. (2000)
prop.unpaired	Altman et al. (2000)
prop.paired	Altman et al. (2000)
inc.risk, exact	Collett (1999)
inc.risk, wilson	Rothman (2012)
inc.risk, fleiss	Fleiss (1981)
inc.risk, agresti	Agresti and Coull (1998)
inc.risk, clopper-pearson	Clopper and Pearson (1934)
inc.risk, jeffreys	Brown et al. (2001)
prevalence, exact	Collett (1999)
prevalence, wilson	Wilson (1927)
prevalence, fleiss	Fleiss (1981)
prevalence, agresti	Agresti and Coull (1998)
prevalence, clopper-pearson	Clopper and Pearson (1934)
prevalence, jeffreys	Brown et al. (2001)
inc.rate, exact	Ulm (1990)
inc.rate, byar	Rothman (2012)
odds	Ederer and Mantel (1974)
ratio	Ederer and Mantel (1974)
smr	Dobson et al. (1991)
-----------------------------------	------------------------

The Wald interval often has inadequate coverage, particularly for small sample sizes and proportion estimates close to 0 or 1. Conversely, the Clopper-Pearson Exact method is conservative and tends to produce wider intervals than necessary. Brown et al. recommends the Wilson or Jeffreys methods when sample sizes are small and Agresti-Coull, Wilson, or Jeffreys methods for larger sample sizes.

The Clopper-Pearson interval is an early and very common method for calculating binomial confidence intervals. The Clopper-Pearson interval is sometimes called an 'exact' method because it is based on the cumulative probabilities of the binomial distribution (i.e., exactly the correct distribution rather than an approximation).

The design effect is used to adjust the confidence interval around a prevalence or incidence risk estimate in the presence of clustering. The design effect is a measure of the variability between clusters and is calculated as the ratio of the variance calculated assuming a complex sample design divided by the variance calculated assuming simple random sampling. Adjustment for the effect of clustering can only be made on those prevalence and incidence risk methods that return a standard error (i.e., method = "wilson" or method = "fleiss").

References

Agresti A, Coull B (1998). Approximate is better than 'exact' for interval estimation of binomial proportions. The American Statistician 52. DOI: 10.2307/2685469.

Altman DG, Machin D, Bryant TN, and Gardner MJ (2000). Statistics with Confidence, second edition. British Medical Journal, London, pp. 28 - 29 and pp. 45 - 56.

Brown L, Cai T, Dasgupta A (2001). Interval estimation for a binomial proportion. Statistical Science 16: 101 - 133.

Clopper C, Pearson E (1934) The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika 26: 404 - 413. DOI: 10.1093/biomet/26.4.404.

Collett D (1999). Modelling Binary Data. Chapman & Hall/CRC, Boca Raton Florida, pp. 24.

Dobson AJ, Kuulasmaa K, Eberle E, and Scherer J (1991). Confidence intervals for weighted sums of Poisson parameters. Statistics in Medicine 10: 457 - 462.

Ederer F, and Mantel N (1974). Confidence limits on the ratio of two Poisson variables. American Journal of Epidemiology 100: 165 - 167

Fleiss JL (1981). Statistical Methods for Rates and Proportions. 2nd edition. John Wiley & Sons, New York.

Killip S, Mahfoud Z, Pearce K (2004). What is an intracluster correlation coefficient? Crucial concepts for primary care researchers. Annals of Family Medicine 2: 204 - 208.

Otte J, Gumm I (1997). Intra-cluster correlation coefficients of 20 infections calculated from the results of cluster-sample surveys. Preventive Veterinary Medicine 31: 147 - 150.

Rothman KJ (2012). Epidemiology An Introduction. Oxford University Press, London, pp. 164 - 175.

Ulm K (1990). A simple method to calculate the confidence interval of a standardized mortality ratio. American Journal of Epidemiology 131: 373 - 375.

Wilson EB (1927) Probable inference, the law of succession, and statistical inference. Journal of the American Statistical Association 22: 209 - 212.

Examples

## EXAMPLE 1:
dat.v01 <- rnorm(n = 100, mean = 0, sd = 1)
epi.conf(dat = dat.v01, ctype = "mean.single")


## EXAMPLE 2:
group <- c(rep("A", times = 5), rep("B", times = 5))
val = round(c(rnorm(n = 5, mean = 10, sd = 5), 
   rnorm(n = 5, mean = 7, sd = 5)), digits = 0)
dat.df02 <- data.frame(group = group, val = val)
epi.conf(dat = dat.df02, ctype = "mean.unpaired")


## EXAMPLE 3:
## Two paired samples (Altman et al. 2000, page 31):
## Systolic blood pressure levels were measured in 16 middle-aged men
## before and after a standard exercise test. The mean rise in systolic 
## blood pressure was 6.6 mmHg. The standard deviation of the difference
## was 6.0 mm Hg. The standard error of the mean difference was 1.49 mm Hg.

before <- c(148,142,136,134,138,140,132,144,128,170,162,150,138,154,126,116)
after <- c(152,152,134,148,144,136,144,150,146,174,162,162,146,156,132,126)
dat.df03 <- data.frame(before, after)
epi.conf(dat = dat.df03, ctype = "mean.paired", conf.level = 0.95)

## The 95% confidence interval for the population value of the mean 
## systolic blood pressure increase after standard exercise was 3.4 to 9.8
## mm Hg.


## EXAMPLE 4:
## Single sample (Altman et al. 2000, page 47):
## Out of 263 giving their views on the use of personal computers in 
## general practice, 81 thought that the privacy of their medical file
## had been reduced.

pos <- 81
neg <- (263 - 81)
dat.m04 <- as.matrix(cbind(pos, neg))
round(epi.conf(dat = dat.m04, ctype = "prop.single"), digits = 3)

## The 95% confidence interval for the population value of the proportion
## of patients thinking their privacy was reduced was from 0.255 to 0.366. 


## EXAMPLE 5:
## Two samples, unpaired (Altman et al. 2000, page 49):
## Goodfield et al. report adverse effects in 85 patients receiving either
## terbinafine or placebo treatment for dermatophyte onchomychois.
## Out of 56 patients receiving terbinafine, 5 patients experienced
## adverse effects. Out of 29 patients receiving a placebo, none experienced
## adverse effects.

grp1 <- matrix(cbind(5, 51), ncol = 2)
grp2 <- matrix(cbind(0, 29), ncol = 2)
dat.m05 <- as.matrix(cbind(grp1, grp2))
round(epi.conf(dat = dat.m05, ctype = "prop.unpaired"), digits = 3)

## The 95% confidence interval for the difference between the two groups is
## from -0.038 to +0.193.


## EXAMPLE 6:
## Two samples, paired (Altman et al. 2000, page 53):
## In a reliability exercise, 41 patients were randomly selected from those
## who had undergone a thalium-201 stress test. The 41 sets of images were
## classified as normal or not by the core thalium laboratory and, 
## independently, by clinical investigators from different centres.
## Of the 19 samples identified as ischaemic by clinical investigators 
## 5 were identified as ischaemic by the laboratory. Of the 22 samples 
## identified as normal by clinical investigators 0 were identified as 
## ischaemic by the laboratory. 

## Clinic       | Laboratory    |             |   
##              | Ischaemic     | Normal      | Total
## ---------------------------------------------------------
##  Ischaemic   | 14            | 5           | 19
##  Normal      | 0             | 22          | 22
## ---------------------------------------------------------
##  Total       | 14            | 27          | 41
## ---------------------------------------------------------

dat.m06 <- as.matrix(cbind(14, 5, 0, 22))
round(epi.conf(dat = dat.m06, ctype = "prop.paired", conf.level = 0.95), 
   digits = 3)

## The 95% confidence interval for the population difference in 
## proportions is 0.011 to 0.226 or approximately +1% to +23%.


## EXAMPLE 7:
## A herd of 1000 cattle were tested for brucellosis. Four samples out of 200
## test returned a positive result. Assuming 100% test sensitivity and 
## specificity, what is the estimated prevalence of brucellosis in this 
## group of animals?

pos <- 4; pop <- 200
dat.m07 <- as.matrix(cbind(pos, pop))
epi.conf(dat = dat.m07, ctype = "prevalence", method = "exact", N = 1000, 
   design = 1, conf.level = 0.95) * 100

## The estimated prevalence of brucellosis in this herd is 2.0 (95% CI 0.54 to
## 5.0) cases per 100 cattle at risk.


## EXAMPLE 8:
## The observed disease counts and population size in four areas are provided
## below. What are the the standardised morbidity ratios of disease for each
## area and their 95% confidence intervals?

obs <- c(5, 10, 12, 18); pop <- c(234, 189, 432, 812)
dat.m08 <- as.matrix(cbind(obs, pop))
round(epi.conf(dat = dat.m08, ctype = "smr"), digits = 2)


## EXAMPLE 9:
## A survey has been conducted to determine the proportion of broilers
## protected from a given disease following vaccination. We assume that 
## the intra-cluster correlation coefficient for protection (also known as the 
## rate of homogeneity, rho) is 0.4 and the average number of birds per
## flock is 30. A total of 5898 birds from a total of 10363 were identified
## as protected. What proportion of birds are protected and what is the 95%
## confidence interval for this estimate?

## Calculate the design effect, given rho = (design - 1) / (nbar - 1), where
## nbar equals the average number of individuals sampled per cluster:

D <- 0.4 * (30 - 1) + 1; D

## The design effect is 12.6. Now calculate the proportion protected. We set
## N to large number.

dat.m09 <- as.matrix(cbind(5898, 10363))
epi.conf(dat = dat.m09, ctype = "prevalence", method = "fleiss", N = 1000000, 
   design = D, conf.level = 0.95)

## The estimated proportion of the population protected is 0.57 (95% CI
## 0.53 to 0.60). Recalculate this estimate assuming the data were from a 
## simple random sample (i.e., where the design effect is one):

epi.conf(dat = dat.m09, ctype = "prevalence", method = "fleiss", N = 1000000, 
   design = 1, conf.level = 0.95)

## If we had mistakenly assumed that data were a simple random sample the 
## confidence interval for the proportion of birds protect would have been 
## 0.56 -- 0.58.

---

Convert British National Grid georeferences to easting and northing coordinates

Description

Convert British National Grid georeferences to Brtish National Grid (EPSG 27700) easting and northing coordinates.

Usage

epi.convgrid(osref)

Arguments

osref

a vector of character strings listing the British National Grid georeferences to be converted.

Note

If an invalid georeference is encountered in the vector os.ref the method returns a NA.

Examples

## EXAMPLE 1:
os.ref <- c("SJ505585","SJ488573","SJ652636")
epi.convgrid(os.ref)

os.ref <- c("SJ505585","SJ488573","ZZ123456")
epi.convgrid(os.ref)

---

Extract unique covariate patterns from a data set

Description

Extract the set of unique patterns from a set of covariates (explanatory variables).

Usage

epi.cp(dat)

Arguments

dat	an i row by j column data frame where the i rows represent individual observations and the m columns represent a set of m covariates. The function allows for one or more covariates for each observation.

Details

This function extracts the k unique covariate patterns in a data set comprised of i observations, labelling them from 1 to k. The frequency of occurrence of each covariate pattern is listed. A vector of length i is also returned, listing the 1:k covariate pattern identifier for each observation.

Value

A list containing the following:

cov.pattern	a data frame with columns: id the unique covariate pattern identifier (labelled 1 to k), n the number of occasions each of the listed covariate pattern appears in the data, and the unique covariate combinations.
id	a vector of length i listing the 1:k covariate pattern identifier for each observation.

Author(s)

Thanks to Johann Popp and Mathew Jay for providing code and suggestions to enhance the utility of this function.

References

Dohoo I, Martin W, Stryhn H (2003). Veterinary Epidemiologic Research. AVC Inc, Charlottetown, Prince Edward Island, Canada.

Examples

## EXAMPLE 1:

## Generate a set of covariates:
set.seed(seed = 1234)
obs <- round(runif(n = 100, min = 0, max = 1), digits = 0)
v1 <- round(runif(n = 100, min = 0, max = 4), digits = 0)
v2 <- round(runif(n = 100, min = 0, max = 4), digits = 0)
dat.df01 <- data.frame(obs, v1, v2)

dat.glm01 <- glm(obs ~ v1 + v2, family = binomial, data = dat.df01)
dat.mf01 <- model.frame(dat.glm01)

## Covariate pattern. Drop the first column of dat.mf01 (since column 1 is the
## outcome variable:
epi.cp(dat.mf01[,2:3])

## There are 25 covariate patterns in this data set. Subject 100 has
## covariate pattern 21.

---

Covariate pattern residuals from a logistic regression model

Description

Returns covariate pattern residuals and delta betas from a logistic regression model.

Usage

epi.cpresids(obs, fit, covpattern)

Arguments

obs	a vector of observed values (i.e., counts of ‘successes’) for each covariate pattern).
fit	a vector defining the predicted (i.e., fitted) probability of success for each covariate pattern.
covpattern	a epi.cp object.

Value

A data frame with 13 elements: cpid the covariate pattern identifier, n the number of subjects in this covariate pattern, obs the observed number of successes, pred the predicted number of successes, raw the raw residuals, sraw the standardised raw residuals, pearson the Pearson residuals, spearson the standardised Pearson residuals, deviance the deviance residuals, leverage leverage, deltabeta the delta-betas, sdeltabeta the standardised delta-betas, and deltachi delta chi statistics.

References

Hosmer DW, Lemeshow S (1989). Applied Logistic Regression. John Wiley & Sons, New York, USA, pp. 137 - 138.

See Also

epi.cp

Examples

## EXAMPLE 1:
dat.glm01 <- glm(case ~ spontaneous + induced, family = binomial(), 
   data = infert)

## Covariate patterns:   
dat.mf01 <- model.frame(dat.glm01)
dat.cp01 <- epi.cp(dat.mf01[-1])

dat.obs01 <- as.vector(by(infert$case, as.factor(dat.cp01$id), 
   FUN = sum))
dat.fit01 <- as.vector(by(fitted(dat.glm01), as.factor(dat.cp01$id), 
   FUN = min))
dat.cpr01 <- epi.cpresids(obs = dat.obs01, fit = dat.fit01, 
   covpattern = dat.cp01)
head(dat.cpr01)

---

Descriptive statistics

Description

Computes descriptive statistics for a numeric vector or a table of frequencies for a factor.

Usage

epi.descriptives(dat, conf.level = 0.95)

Arguments

dat	either a numeric vector or a factor.
conf.level	magnitude of the returned confidence intervals. Must be a single number between 0 and 1.

Value

If dat is numeric a list containing the following:

arithmetic	n number of observations, mean arithmetic mean, sd arithmetic standard deviation, q25 25th quantile, q50 50th quantile, q75 75th quantile, lower lower bound of the confidence interval, upper upper bound of the confidence interval, min minimum value, max maximum value, and na number of missing values.
geometric	n number of observations, mean geometric mean, sd geometric standard deviation, q25 25th quantile, q50 50th quantile, q75 75th quantile, lower lower bound of the confidence interval, upper upper bound of the confidence interval, min minimum value, max maximum value, and na number of missing values.
symmetry	skewness and kurtosis.

If dat is a factor a data frame listing:

level	The levels of the factor
n	The frequency of the respective factor level, including the column totals.

Examples

## EXAMPLE 1:
## Generate some data:
id <- 1:100
n <- rnorm(100, mean = 0, sd = 1)
dat.df01 <- data.frame(id, n)

# Add missing values:
missing <- dat.df01$id %in% sample(dat.df01$id, size = 20)
dat.df01$n[missing] <- NA

epi.descriptives(dat.df01$n, conf.level = 0.95)


## EXAMPLE 2:
## Generate some data:
n <- 1000; p.exp <- 0.50; p.dis <- 0.75
strata <- c(rep("A", times = n / 2), rep("B", times = n / 2))
exp <- rbinom(n = n, size = 1, prob = p.exp)
dis <- rbinom(n = n, size = 1, prob = p.dis)
dat.df02 <- data.frame(strata, exp, dis)

dat.df02$strata <- factor(dat.df02$strata)
dat.df02$exp <- factor(dat.df02$exp, levels = c("1", "0"))
head(dat.df02)

epi.descriptives(dat.df02$exp, conf.level = 0.95)

---

Estimate the precision of a [structured] heterogeneity term

Description

Returns the precision of a [structured] heterogeneity term after one has specified the amount of variation a priori.

Usage

epi.dgamma(rr, quantiles = c(0.05, 0.95))

Arguments

rr	the lower and upper limits of relative risk, estimated a priori.
quantiles	a vector of length two defining the quantiles of the lower and upper relative risk estimates.

Value

Returns the precision (the inverse variance) of the heterogeneity term.

References

Best, NG. WinBUGS 1.3.1 Short Course, Brisbane Australia, November 2000.

Examples

## EXAMPLE 1:
## Suppose we are expecting the lower 5% and upper 95% confidence interval 
## of relative risk in a data set to be 0.5 and 3.0, respectively. 
## A prior estimate of the precision of the heterogeneity term would be:

tau <- epi.dgamma(rr = c(0.5, 3.0), quantiles = c(0.05, 0.95))
tau

## The estimate of the precision of the heterogeneity term (tau) is 3.37. 
## This can be re-expressed using the gamma distribution. We set the mean of the 
## distribution as tau and specify a large variance (that is, we are not 
## certain about tau).

mean <- tau; var <- 1000
shape <- mean^2 / var
inv.scale <- mean / var

## In WinBUGS the precision of the heterogeneity term is parameterised 
## as tau ~ dgamma(shape, inv.scale). Plot the probability density function
## of tau:

z <- seq(0.01, 10, by = 0.01)
fz <- dgamma(z, shape = shape, scale = 1 / inv.scale)
plot(x = z, y = fz, type = "l", ylab = "Probability density of tau")

---

Directly adjusted incidence rate estimates

Description

Compute directly adjusted incidence rate estimates.

Usage

epi.directadj(obs, tar, std, units = 1, conf.level = 0.95)

Arguments

obs	a matrix providing counts of the observed number of events. Rows represent strata (e.g., region); columns represent the different levels of the variable to be adjusted for (e.g., age class) — the covariate. The sum of each row will equal the total number of observed events for each stratum. The row names of matrix obs are the appropriate strata names. The column names of obs are the appropriate level identifiers for the adjustment covariate. See the examples, below.
tar	a matrix representing population time at risk. Rows represent strata (e.g., region); columns represent the different levels of the variable to be adjusted for (e.g., age class) — the covariate. The sum of each row will equal the total population time at risk for each stratum. The row names of matrix pop must be the same as the row names used for matrix obs. The column names of matrix pop must be the same as the column names used for matrix pop. See the examples, below.
std	a matrix representing the standard population size for the different levels of the covariate to be adjusted for. The column names of matrix std must be the same as the column names used for matrices obs and pop.
units	multiplier for the incidence rate estimates.
conf.level	magnitude of the returned confidence interval. Must be a single number between 0 and 1.

Details

This function returns unadjusted (crude) and directly adjusted incidence rate estimates for each of the specified population strata. The term ‘covariate’ is used here to refer to the factors we want to control (i.e., adjust) for when calculating the directly adjusted incidence rate estimates.

When the outcome of interest is rare, the confidence intervals for the adjusted incidence rates returned by this function (based on Fay and Feuer, 1997) will be appropriate for incidence risk data. In this situation the argument tar is assumed to represent the size of the population at risk (instead of population time at risk). Example 3 (below) provides an approach if you are working with incidence risk data and the outcome of interest is not rare.

Value

A list containing the following:

crude	the crude incidence rate estimates for each stratum-covariate combination.
crude.strata	the crude incidence rate estimates for each stratum.
adj.strata	the directly adjusted incidence rate estimates for each stratum.

Author(s)

Thanks to Karl Ove Hufthammer for helpful suggestions to improve the execution and documentation of this function.

References

Fay M, Feuer E (1997). Confidence intervals for directly standardized rates: A method based on the gamma distribution. Statistics in Medicine 16: 791 - 801.

Fleiss JL (1981). Statistical Methods for Rates and Proportions, Wiley, New York, USA, pp. 240.

Frome E, Checkoway H (1985). Use of Poisson regression models in estimating incidence rates and ratios. American Journal of Epidemiology 121: 309 - 323.

Haneuse S, Rothman KJ. Stratification and Standardization. In: Lash TL, VanderWeele TJ, Haneuse S, Rothman KJ (2021). Modern Epidemiology. Lippincott - Raven Philadelphia, USA, pp. 415 - 445.

Thrusfield M (2007). Veterinary Epidemiology, Blackwell Publishing, London, UK, pp. 63 - 64.

Wilcosky T, Chambless L (1985). A comparison of direct adjustment and regression adjustment of epidemiologic measures. Journal of Chronic Diseases 38: 849 - 956.

See Also

epi.indirectadj

Examples

## EXAMPLE 1 (from Thrusfield 2007 pp. 63 - 64):
## A study was conducted to estimate the seroprevalence of leptospirosis in 
## dogs in Glasgow and Edinburgh, Scotland. Data frame dat.df lists counts 
## of leptospirosis cases and the number of dog years at risk for male and 
## female dogs. In this example we stratify the data by city and the covariate 
## is gender (with levels male and female). The data are presented as a matrix
## with two rows and two columns with city (Edinburgh, Glasgow) as the row 
## names and levels of gender (male, female) as column names: 

dat.df01 <- data.frame(obs = c(15,46,53,16), tar = c(48,212,180,71), 
   sex = c("M","F","M","F"), city = c("ED","ED","GL","GL"))

obs01 <- matrix(dat.df01$obs, nrow = 2, byrow = TRUE, 
   dimnames = list(c("ED","GL"), c("M","F")))
tar01 <- matrix(dat.df01$tar, nrow = 2, byrow = TRUE, 
   dimnames = list(c("ED","GL"), c("M","F")))

## Create a standard population with equal numbers of male and female dogs:

std01 <- matrix(data = c(250,250), nrow = 1, byrow = TRUE, 
   dimnames = list("", c("M","F")))

## Directly adjusted incidence rates:
epi.directadj(obs01, tar01, std01, units = 1, conf.level = 0.95)

## $crude
##  strata cov obs tar       est     lower     upper
##      ED   M  15  48 0.3125000 0.1749039 0.5154212
##      GL   M  53 180 0.2944444 0.2205591 0.3851406
##      ED   F  46 212 0.2169811 0.1588575 0.2894224
##      GL   F  16  71 0.2253521 0.1288082 0.3659577

## $crude.strata
##  strata obs tar       est     lower     upper
##      ED  61 260 0.2346154 0.1794622 0.3013733
##      GL  69 251 0.2749004 0.2138889 0.3479040

## $adj.strata
##  strata obs tar       est     lower     upper
##      ED  61 260 0.2647406 0.1866047 0.3692766
##      GL  69 251 0.2598983 0.1964162 0.3406224

## The adjusted incidence rate of leptospirosis in Glasgow dogs is 26 (95%
## CI 20 to 34) cases per 100 dog-years at risk. The confounding effect of 
## gender has been removed by the adjusted incidence rate estimates.


## EXAMPLE 2:
## Here we provide a more flexible approach for calculating 
## adjusted incidence rate estimates using Poisson regression. See Frome and 
## Checkoway (1985) for details.

dat.glm02 <- glm(obs ~ city, offset = log(tar), family = poisson, 
   data = dat.df01)
summary(dat.glm02)

## To obtain adjusted incidence rate estimates, use the predict method on a 
## new data set with the time at risk (tar) variable set to 1 (which means 
## log(tar) = 0). This will return the predicted number of cases per one unit 
## of individual time, i.e., the incidence rate.

dat.pred02 <- predict(object = dat.glm02, newdata = 
   data.frame(city = c("ED","GL"), tar = c(1,1)), 
   type  = "link", se = TRUE)

conf.level <- 0.95
critval <- qnorm(p = 1 - ((1 - conf.level) / 2), mean = 0, sd = 1)
est <- dat.glm02$family$linkinv(dat.pred02$fit)
lower <- dat.glm02$family$linkinv(dat.pred02$fit - 
   (critval * dat.pred02$se.fit))
upper <- dat.glm02$family$linkinv(dat.pred02$fit + 
   (critval * dat.pred02$se.fit))
round(x = data.frame(est, lower, upper), digits = 3)

##   est lower upper
## 0.235 0.183 0.302
## 0.275 0.217 0.348

## Results identical to the crude incidence rate estimates from epi.directadj.


## EXAMPLE 3:
## Now adjust for the effect of gender and city and report the adjusted 
## incidence rate estimates for each city:

dat.glm03 <- glm(obs ~ city + sex, offset = log(tar), 
   family = poisson, data = dat.df01)
dat.pred03 <- predict(object = dat.glm03, newdata = 
   data.frame(sex = c("F","F"), city = c("ED","GL"), tar = c(1,1)), 
   type  = "link", se.fit = TRUE)

conf.level <- 0.95
critval <- qnorm(p = 1 - ((1 - conf.level) / 2), mean = 0, sd = 1)
est <- dat.glm03$family$linkinv(dat.pred03$fit)
lower <- dat.glm03$family$linkinv(dat.pred03$fit - 
   (critval * dat.pred03$se.fit))
upper <- dat.glm03$family$linkinv(dat.pred03$fit + 
   (critval * dat.pred03$se.fit))
round(x = data.frame(est, lower, upper), digits = 3)

##   est lower upper
## 0.220 0.168 0.287
## 0.217 0.146 0.323

## Using Poisson regression the gender adjusted incidence rate of leptospirosis 
## in Glasgow dogs was 22 (95% CI 15 to 32) cases per 100 dog-years at risk.
## These results won't be the same as those using direct adjustment because
## for direct adjustment we use a contrived standard population.


## EXAMPLE 4 --- Logistic regression to return adjusted incidence risk 
## estimates:

## Say, for argument's sake, that we are now working with incidence risk data. 
## Here we'll re-label the variable 'tar' (time at risk) as 'pop' 
## (population size). We adjust for the effect of gender and city and 
## report the adjusted incidence risk of canine leptospirosis estimates for 
## each city:

dat.df01$pop <- dat.df01$tar

dat.glm04 <- glm(cbind(obs, pop - obs) ~ city + sex, 
   family = "binomial", data = dat.df01)
dat.pred04 <- predict(object = dat.glm04, newdata = 
   data.frame(sex = c("F","F"), city = c("ED","GL")), 
   type  = "link", se.fit = TRUE)

conf.level <- 0.95
critval <- qnorm(p = 1 - ((1 - conf.level) / 2), mean = 0, sd = 1)
est <- dat.glm04$family$linkinv(dat.pred04$fit)
lower <- dat.glm04$family$linkinv(dat.pred04$fit - 
   (critval * dat.pred04$se.fit))
upper <- dat.glm04$family$linkinv(dat.pred04$fit + 
   (critval * dat.pred04$se.fit))
round(x = data.frame(est, lower, upper), digits = 3)

##   est lower upper
## 0.220 0.172 0.276
## 0.217 0.150 0.304

## The adjusted incidence risk of leptospirosis in Glasgow dogs is 22 (95%
## CI 15 to 30) cases per 100 dogs at risk.

---

Decimal degrees and degrees, minutes and seconds conversion

Description

Converts decimal degrees to degrees, minutes and seconds. Converts degrees, minutes and seconds to decimal degrees.

Usage

epi.dms(dat)

Arguments

dat	the data. A one-column matrix is assumed when converting decimal degrees to degrees, minutes, and seconds. A two-column matrix is assumed when converting degrees and decimal minutes to decimal degrees. A three-column matrix is assumed when converting degrees, minutes and seconds to decimal degrees.

Examples

## EXAMPLE 1:
## Degrees, minutes, seconds to decimal degrees:
dat.m01 <- matrix(c(41, 38, 7.836, -40, 40, 27.921), 
   byrow = TRUE, nrow = 2)
epi.dms(dat.m01)


## EXAMPLE 2: 
## Decimal degrees to degrees, minutes, seconds:
dat.m02 <- matrix(c(41.63551, -40.67442), nrow = 2)
epi.dms(dat.m02)

---

Mixed-effects meta-analysis of binary outcomes using the DerSimonian and Laird method

Description

Computes individual study odds or risk ratios for binary outcome data. Computes the summary odds or risk ratio using the DerSimonian and Laird method. Performs a test of heterogeneity among trials. Performs a test for the overall difference between groups (that is, after pooling the studies, do treated groups differ significantly from controls?).

Usage

epi.dsl(ev.trt, n.trt, ev.ctrl, n.ctrl, names, method = "odds.ratio", 
   alternative = c("two.sided", "less", "greater"), conf.level = 0.95)

Arguments

ev.trt	observed number of events in the treatment group.
n.trt	number in the treatment group.
ev.ctrl	observed number of events in the control group.
n.ctrl	number in the control group.
names	character string identifying each trial.
method	a character string indicating the method to be used. Options are odds.ratio or risk.ratio.
alternative	a character string specifying the alternative hypothesis, must be one of two.sided, greater or less.
conf.level	magnitude of the returned confidence interval. Must be a single number between 0 and 1.

Details

alternative = "greater" tests the hypothesis that the DerSimonian and Laird summary measure of association is greater than 1.

Value

A list containing the following:

OR	the odds ratio for each trial and the lower and upper bounds of the confidence interval of the odds ratio for each trial.
RR	the risk ratio for each trial and the lower and upper bounds of the confidence interval of the risk ratio for each trial.
OR.summary	the DerSimonian and Laird summary odds ratio and the lower and upper bounds of the confidence interval of the DerSimonian and Laird summary odds ratio.
RR.summary	the DerSimonian and Laird summary risk ratio and the lower and upper bounds of the confidence interval of the DerSimonian and Laird summary risk ratio.
weights	the inverse variance and DerSimonian and Laird weights for each trial.
heterogeneity	a vector containing Q the heterogeneity test statistic, df the degrees of freedom and its associated P-value.
Hsq	the relative excess of the heterogeneity test statistic Q over the degrees of freedom df.
Isq	the percentage of total variation in study estimates that is due to heterogeneity rather than chance.
tau.sq	the variance of the treatment effect among trials.
effect	a vector containing z the test statistic for overall treatment effect and its associated P-value.

Note

Under the random-effects model, the assumption of a common treatment effect is relaxed, and the effect sizes are assumed to have a normal distribution with variance tau.sq.

Using this method, the DerSimonian and Laird weights are used to compute the pooled odds ratio.

The function checks each strata for cells with zero frequencies. If a zero frequency is found in any cell, 0.5 is added to all cells within the strata.

References

Deeks JJ, Altman DG, Bradburn MJ (2001). Statistical methods for examining heterogeneity and combining results from several studies in meta-analysis. In: Egger M, Davey Smith G, Altman D (eds). Systematic Review in Health Care Meta-Analysis in Context. British Medical Journal, London, 2001, pp. 291 - 299.

DerSimonian R, Laird N (1986). Meta-analysis in clinical trials. Controlled Clinical Trials 7: 177 - 188.

Higgins J, Thompson S (2002). Quantifying heterogeneity in a meta-analysis. Statistics in Medicine 21: 1539 - 1558.

See Also

epi.iv, epi.mh, epi.smd

Examples

## EXAMPLE 1:
data(epi.epidural)
epi.dsl(ev.trt = epi.epidural$ev.trt, n.trt = epi.epidural$n.trt, 
   ev.ctrl = epi.epidural$ev.ctrl, n.ctrl = epi.epidural$n.ctrl, 
   names = as.character(epi.epidural$trial), method = "odds.ratio", 
   alternative = "two.sided", conf.level = 0.95)

---

Estimated dissemination ratio

Description

Computes estimated dissemination ratios on the basis of a vector of count data (usually incident disease cases identified on each day of an epidemic).

Usage

epi.edr(dat, n = 4, conf.level = 0.95, nsim = 99, na.zero = TRUE)

Arguments

dat	a numeric vector listing the number of incident cases for each day of an epidemic.
n	scalar, defining the number of days to be used when computing the estimated dissemination ratio.
conf.level	magnitude of the returned confidence interval. Must be a single number between 0 and 1.
nsim	scalar, defining the number of simulations to be used for the confidence interval calculations.
na.zero	logical, replace NaN or Inf values with zeros?

Details

In infectious disease outbreaks the n-day estimated dissemination ratio (EDR) at day i equals the total number of incident disease cases between day i and day [i - (n - 1)] (inclusive) divided by the total number of incident disease cases between day (i - n) and day (i - 2n) + 1 (inclusive). EDR values are often calculated for each day of an outbreak and presented as a time series analysis. If the EDR is consistently less than one, the epidemic is said to be ‘under control’.

A simulation approach is used to calculate confidence intervals for each daily EDR estimate. The numerator and denominator of the EDR estimate for each day is taken in turn and a random number drawn from a Poisson distribution, using the calculated numerator and denominator value as the mean. EDR is then calculated for these simulated values and the process repeated nsim times. Confidence intervals are then derived from the vector of simulated values for each day.

Value

Returns the point estimate of the EDR and the lower and upper bounds of the confidence interval of the EDR.

References

Miller W (1976). A state-transition model of epidemic foot-and-mouth disease. In: Proceedings of an International Symposium: New Techniques in Veterinary Epidemiology and Economics, University of Reading, Reading, pp. 56 - 72.

Morris R, Sanson R, Stern M, Stevenson M, Wilesmith J (2002). Decision-support tools for foot-and-mouth disease control. Revue Scientifique et Technique de l'Office International des Epizooties 21, 557 - 567.

Perez-Reche FJ, Taylor N, McGuigan C, Conaglen P, Forbes K, Strachan N, Honhold N (2021) Estimated Dissemination Ratio — A practical alternative to the reproduction number for infectious diseases. Frontiers in Public Health 9. DOI: 10.3389/fpubh.2021.675065.

Examples

## EXAMPLE 1 Foot and mouth disease, Cumbria 2001:

## Counts of incident FMD positive farms by date:
edate <- seq(from = as.Date("2001-02-28", format = "%Y-%m-%d"), 
   to = as.Date("2001-06-15", format = "%Y-%m-%d"), by = 1)
ncas <- c(1,2,0,0,1,1,0,0,4,2,3,3,5,2,8,2,5,0,5,7,15,13,6,7,11,8,7,11,
   6,5,10,7,8,8,7,5,6,3,3,3,3,4,1,4,6,2,1,4,3,3,1,1,1,2,2,2,2,0,4,1,1,
   0,0,2,1,2,1,0,1,2,2,4,0,1,0,1,0,0,0,1,0,3,1,1,3,0,0,1,1,2,0,0,1,1,1,
   3,3,1,1,1,0,0,0,1,1,1,1,1)
dat.df01 <- data.frame(edate = edate, ncas = ncas)


## Seven day EDR:
edr <- epi.edr(dat.df01$ncas, n = 7, conf.level = 0.95, nsim = 99, 
   na.zero = FALSE)

dat.df01$edr.500 <- edr$est
dat.df01$edr.025 <- edr$lower
dat.df01$edr.975 <- edr$upper

## log2 transformed EDRs:
dat.df01$ledr.500 <- log(edr$est, base = 2)
dat.df01$ledr.025 <- log(edr$lower, base = 2)
dat.df01$ledr.975 <- log(edr$upper, base = 2)

## You may want to smooth the EDRs:
dat.df01$sedr.500 <- lowess(dat.df01$edate, dat.df01$est, f = 0.15)$y
dat.df01$sedr.025 <- lowess(dat.df01$edate, dat.df01$lower, f = 0.15)$y
dat.df01$sedr.975 <- lowess(dat.df01$edate, dat.df01$upper, f = 0.15)$y

## Not run: 
library(ggplot2); library(scales)

## Frequency histogram of the number of incident FMD positive farms as a
## function of date:
ggplot() +
  theme_bw() +
  geom_histogram(data = dat.df01, aes(x = edate, weight = ncas), 
     binwidth = 1, fill = "#008080", alpha = 0.45) +
  scale_x_date(breaks = date_breaks("7 days"), name = "Date") +
  scale_y_continuous(limits = c(0,15), 
     name = "Number of events") + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1))


## EDR line plot (and its 95% confidence interval) as a function of 
## date:
ybreaks <- -5:5
ylabels <- 2^(-5:5)


## Drop EDR estimates outside of the 0.03125 and 32 range since distracting on
## plot:
dat.df01$ledr.025[dat.df01$ledr.025^2 < 0.03125 | 
   dat.df01$ledr.025^2 > 32] <- NA
dat.df01$ledr.975[dat.df01$ledr.975^2 < 0.03125 | 
   dat.df01$ledr.975^2 > 32] <- NA

ggplot() +
  theme_bw() +
  geom_line(data = dat.df01, aes(x = edate, y = ledr.500), 
     col = "black", linewidth = 1) +
  geom_line(data = dat.df01, aes(x = edate, y = ledr.025), 
     col = "black", linetype = "dashed", linewidth = 0.5) +
  geom_line(data = dat.df01, aes(x = edate, y = ledr.975), 
     col = "black", linetype = "dashed", linewidth = 0.5) +
  scale_x_date(breaks = date_breaks("7 days"), name = "Date") +
  scale_y_continuous(name = "Estimated dissemination ratio (EDR)",
     breaks = ybreaks, labels = ylabels, limits = range(ybreaks)) + 
  geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))  


## EDR line plot (and its 95% confidence interval) as a function of 
## date superimposed on the epidemic curve. Set the upper limit for the 
## vertical axis of the histogram as ymax. The vertical position of 
## the EDR curve will be shifted by ymax / 2 (i.e., half way up the plot):
ymax <- 20; shift <- ymax / 2
ybreaks <- -5:5 + shift
ylabels <- 2^(-5:5)

ggplot() +
  theme_bw() +
  geom_histogram(data = dat.df01, aes(x = edate, weight = ncas), 
     binwidth = 1, fill = "#008080", alpha = 0.45) +
  geom_line(data = dat.df01, aes(x = edate, y = ledr.500 + shift), 
     col = "black", linewidth = 1) +
  geom_line(data = dat.df01, aes(x = edate, y = ledr.025 + shift), 
     col = "black", linetype = "dashed", linewidth = 0.5) +
  geom_line(data = dat.df01, aes(x = edate, y = ledr.975 + shift), 
      col = "black", linetype = "dashed", linewidth = 0.5) +
  scale_x_date(breaks = date_breaks("7 days"), name = "Date") +
  scale_y_continuous(limits = c(0,ymax), 
     name= "Estimated dissemination ratio (EDR)",
     sec.axis = sec_axis(trans = ~ ., name = "Number of events"),
     breaks = ybreaks, labels = ylabels) + 
  geom_hline(yintercept = shift, col = "red", linetype = "dashed") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1), 
     panel.grid.minor = element_blank())
 
## End(Not run)

---

Empirical Bayes estimates of observed event counts

Description

Computes empirical Bayes estimates of observed event counts using the method of moments.

Usage

epi.empbayes(obs, pop)

Arguments

obs	a vector representing the observed event counts in each unit of interest.
pop	a vector representing the population count in each unit of interest.

Details

The gamma distribution is parameterised in terms of shape ($\alpha$) and scale ($\nu$) parameters. The mean of a given gamma distribution equals $\nu / \alpha$. The variance equals $\nu / \alpha^{2}$. The empirical Bayes estimate of event risk in each unit of interest equals $(obs + \nu) / (pop + \alpha)$.

This technique performs poorly when your data contains large numbers of zero event counts. In this situation a Bayesian approach for estimating $\alpha$ and $\nu$ would be advised.

Value

A data frame with four elements: gamma the mean event risk across all units, phi the variance of event risk across all units, alpha the estimated shape parameter of the gamma distribution, and nu the estimated scale parameter of the gamma distribution.

References

Bailey TC, Gatrell AC (1995). Interactive Spatial Data Analysis. Longman Scientific & Technical. London, pp. 303 - 308.

Langford IH (1994). Using empirical Bayes estimates in the geographical analysis of disease risk. Area 26: 142 - 149.

Meza J (2003). Empirical Bayes estimation smoothing of relative risks in disease mapping. Journal of Statistical Planning and Inference 112: 43 - 62.

Examples

## EXAMPLE 1:
data(epi.SClip)
obs <- epi.SClip$cases; pop <- epi.SClip$population

est <- epi.empbayes(obs, pop)
crude.p <- ((obs) / (pop)) * 100000
crude.r <- rank(crude.p)
ebay.p <- ((obs + est[4]) / (pop + est[3])) * 100000

dat.df01 <- data.frame(rank = c(crude.r, crude.r), 
   Method = c(rep("Crude", times = length(crude.r)), 
      rep("Empirical Bayes", times = length(crude.r))),
   est = c(crude.p, ebay.p)) 

## Scatter plot showing the crude and empirical Bayes adjusted lip cancer 
## incidence rates as a function of district rank for the crude lip 
## cancer incidence rates: 
                          
## Not run: 
library(ggplot2)

ggplot(dat = dat.df01, aes(x = rank, y = est, colour = Method)) +
  theme_bw() +
  geom_point() +
  scale_x_continuous(name = "District rank", 
     breaks = seq(from = 0, to = 60, by = 10), 
     labels = seq(from = 0, to = 60, by = 10), 
     limits = c(0,60)) +
  scale_y_continuous(limits = c(0,30), name = "Lip cancer incidence rates 
     (cases per 100,000 person years)") 

## End(Not run)

---

Rates of use of epidural anaesthesia in trials of caregiver support

Description

This data set provides results of six trials investigating rates of use of epidural anaesthesia during childbirth. Each trial is made up of a group where a caregiver (midwife, nurse) provided support intervention and a group where standard care was provided. The objective was to determine if there were higher rates of epidural use when a caregiver was present at birth.

Usage

data(epi.epidural)

Format

A data frame with 6 observations on the following 5 variables.

trial

the name and year of the trial.

ev.trt

number of births in the caregiver group where an epidural was used.

n.trt

number of births in the caregiver group.

ev.ctrl

number of births in the standard care group where an epidural was used.

n.ctrl

number of births in the standard care group.

References

Deeks JJ, Altman DG, Bradburn MJ (2001). Statistical methods for examining heterogeneity and combining results from several studies in meta-analysis. In: Egger M, Davey Smith G, Altman D (eds). Systematic Review in Health Care Meta-Analysis in Context. British Medical Journal, London, pp. 291 - 299.

---

Estimate the characteristics of diagnostic tests applied at the herd (group) level

Description

When tests are applied to individuals within a group we may wish to designate the group as being either diseased or non-diseased on the basis of the individual test results. This function estimates sensitivity and specificity of this testing regime at the group (or herd) level.

Usage

epi.herdtest(se, sp, P, N, n, k)

Arguments

se	a vector of length one defining the sensitivity of the individual test used.
sp	a vector of length one defining the specificity of the individual test used.
P	scalar, defining the estimated true prevalence.
N	scalar, defining the herd size.
n	scalar, defining the number of individuals to be tested per group (or herd).
k	scalar, defining the critical number of individuals testing positive that will denote the group as test positive.

Value

A list with one scalar and two data frames.

Scalar sfraction reports the sampling fraction (i.e., n / N). The binomial distribution is recommended if sfraction is less than 0.2.

Data frame dbinom lists APpos the probability of obtaining a positive test, APneg the probability of obtaining a negative test, HSe the estimated group (herd) sensitivity, and HSp the estimated group (herd) specificity calculated using the binomial distribution.

Data frame dhyper lists APpos the probability of obtaining a positive test, APneg the probability of obtaining a negative test, HSe the estimated group (herd) sensitivity, and HSp the estimated group (herd) specificity calculated using the hypergeometric.

Author(s)

Ron Thornton, MAF New Zealand, PO Box 2526 Wellington, New Zealand.

References

Dohoo I, Martin W, Stryhn H (2003). Veterinary Epidemiologic Research. AVC Inc, Charlottetown, Prince Edward Island, Canada, pp. 113 - 115.

Examples

## EXAMPLE 1:
## We want to estimate the herd-level sensitivity and specificity of 
## a testing regime using an individual animal test of sensitivity 0.391 
## and specificity 0.964. The estimated true prevalence of disease is 0.12.
## Assume that 60 individuals will be tested per herd and we have
## specified that two or more positive test results identify the herd
## as positive.

epi.herdtest(se = 0.391, sp = 0.964, P = 0.12, N = 1E06, n = 60, k = 2)

## This testing regime gives a herd sensitivity of 0.99 and a herd
## specificity of 0.36 (using the binomial distribution). With a herd 
## sensitivity of 0.95 we can be confident that we will declare a herd 
## as disease positive if it truly is disease positive. With a herd specficity 
## of only 0.36, we will declare 0.64 of disease negative herds as infected, 
## so false positives are a problem.

---

Laryngeal and lung cancer cases in Lancashire 1974 - 1983

Description

Between 1972 and 1980 an industrial waste incinerator operated at a site about 2 kilometres southwest of the town of Coppull in Lancashire, England. Addressing community concerns that there were greater than expected numbers of laryngeal cancer cases in close proximity to the incinerator Diggle et al. (1990) conducted a study investigating risks for laryngeal cancer, using recorded cases of lung cancer as controls. The study area is 20 km x 20 km in size and includes location of residence of patients diagnosed with each cancer type from 1974 to 1983. The site of the incinerator was at easting 354500 and northing 413600.

Usage

data(epi.incin)

Format

A data frame with 974 observations on the following 3 variables.

xcoord

easting coordinate (in metres) of each residence.

ycoord

northin coordinate (in metres) of each residence.

status

disease status: 0 = lung cancer, 1 = laryngeal cancer.

Source

Bailey TC and Gatrell AC (1995). Interactive Spatial Data Analysis. Longman Scientific & Technical. London.

References

Diggle P, Gatrell A, and Lovett A (1990). Modelling the prevalence of cancer of the larynx in Lancashire: A new method for spatial epidemiology. In: Thomas R (Editor), Spatial Epidemiology. Pion Limited, London, pp. 35 - 47.

Diggle P (1990). A point process modelling approach to raised incidence of a rare phenomenon in the viscinity of a prespecified point. Journal of the Royal Statistical Society A 153: 349 - 362.

Diggle P, Rowlingson B (1994). A conditional approach to point process modelling of elevated risk. Journal of the Royal Statistical Society A 157: 433 - 440.

---

Indirectly adjusted incidence risk estimates

Description

Compute indirectly adjusted incidence risks and standardised mortality (incidence) ratios.

Usage

epi.indirectadj(obs, pop, std, units, conf.level = 0.95)

Arguments

obs	a one column matrix representing the number of observed number of events in each strata. The dimensions of obs must be named (see the examples, below).
pop	a matrix representing population size. Rows represent strata (e.g., region); columns represent the levels of the explanatory variable to be adjusted for (e.g., age class, gender). The sum of each row will equal the total population size within each stratum. If there are no covariates pop will be a one column matrix. The dimensions of the pop matrix must be named (see the examples, below).
std	a one row matrix specifying the standard incidence risks to be applied to each level of the covariate to be adjusted for. The length of std should be one plus the number of covariates to be adjusted for (the additional value represents the incidence risk in the entire population). If there are no explanatory variables to adjust-for std is a single number representing the incidence risk in the entire population.
units	multiplier for the incidence risk estimates.
conf.level	magnitude of the returned confidence interval. Must be a single number between 0 and 1.

Details

Indirect standardisation can be performed whenever the stratum-specific incidence risk estimates are either unknown or unreliable. If the stratum-specific incidence risk estimates are known, direct standardisation is preferred.

Confidence intervals for the standardised mortality ratio estimates are based on the Poisson distribution (see Breslow and Day 1987, p 69 - 71 for details).

Value

A list containing the following:

crude.strata	the crude incidence risk estimates for each stratum.
adj.strata	the indirectly adjusted incidence risk estimates for each stratum.
smr	the standardised mortality (incidence) ratios for each stratum.

Author(s)

Thanks to Dr. Telmo Nunes (UISEE/DETSA, Faculdade de Medicina Veterinaria - UTL, Rua Prof. Cid dos Santos, 1300-477 Lisboa Portugal) for details and code for the confidence interval calculations.

References

Breslow NE, Day NE (1987). Statistical Methods in Cancer Reasearch: Volume II - The Design and Analysis of Cohort Studies. Lyon: International Agency for Cancer Research.

Dohoo I, Martin W, Stryhn H (2009). Veterinary Epidemiologic Research. AVC Inc, Charlottetown, Prince Edward Island, Canada, pp. 85 - 89.

Lash TL, VanderWeele TJ, Haneuse S, Rothman KJ (2021). Modern Epidemiology. Lippincott - Raven Philadelphia, USA, pp. 75.

Sahai H, Khurshid A (1993). Confidence intervals for the mean of a Poisson distribution: A review. Biometrical Journal 35: 857 - 867.

Sahai H, Khurshid A (1996). Statistics in Epidemiology. Methods, Techniques and Applications. CRC Press, Baton Roca.

See Also

epi.directadj

Examples

## EXAMPLE 1 (without covariates):
## Adapted from Dohoo, Martin and Stryhn (2009). In this example the frequency
## of tuberculosis is expressed as incidence risk (i.e., the number of 
## tuberculosis positive herds divided by the size of the herd population at 
## risk). In their text Dohoo et al. present the data as incidence rate (the
## number of tuberculosis positive herds per herd-year at risk).

## Data have been collected on the incidence of tuberculosis in two
## areas ("A" and "B"). Provided are the counts of (new) incident cases and
## counts of the herd population at risk. The standard incidence risk for
## the total population is 0.060 (6 cases per 100 herds at risk):

obs.m01 <- matrix(data = c(58,130), nrow = 2, byrow = TRUE,
   dimnames = list(c("A", "B"), ""))
pop.m01 <- matrix(data = c(1000,2000), nrow = 2, byrow = TRUE,
   dimnames = list(c("A", "B"), ""))
std.m01 <- 0.060

epi.indirectadj(obs = obs.m01, pop = pop.m01, std = std.m01, units = 100,
   conf.level = 0.95)


## EXAMPLE 2 (with covariates):
## We now have, for each area, data stratified by herd type (dairy, beef).
## The standard incidence risks for beef herds, dairy herds, and the total
## population are 0.025, 0.085, and 0.060 cases per herd, respectively:

obs.m02 <- matrix(data = c(58,130), nrow = 2, byrow = TRUE,
   dimnames = list(c("A", "B"), ""))
pop.m02 <- matrix(data = c(550,450,500,1500), nrow = 2, byrow = TRUE,
   dimnames = list(c("A", "B"), c("Beef", "Dairy")))
std.m02 <- matrix(data = c(0.025,0.085,0.060), nrow = 1, byrow = TRUE,
   dimnames = list("", c("Beef", "Dairy", "Total")))

epi.indirectadj(obs = obs.m02, pop = pop.m02, std = std.m02, units = 100,
   conf.level = 0.95)

## > $crude.strata
## >   est    lower    upper
## > A 5.8 4.404183 7.497845
## > B 6.5 5.430733 7.718222

## > $adj.strata
## >        est    lower    upper
## > A 6.692308 5.076923 8.423077
## > B 5.571429 4.628571 6.557143

## > $smr.strata
## >   obs exp       est     lower    upper
## > A  58  52 1.1153846 0.8461538 1.403846
## > B 130 140 0.9285714 0.7714286 1.092857

## The crude incidence risk of tuberculosis in area A was 5.8
## (95% CI 4.0 to 7.5) cases per 100 herds at risk. The crude incidence
## risk of tuberculosis in area B was 6.5 (95% CI 5.4 to 7.7) cases
## per 100 herds at risk.

## The indirectly adjusted incidence risk of tuberculosis in area A was 6.7
## (95% CI 5.1 to 8.4) cases per 100 herds at risk. The indirectly
## adjusted incidence risk of tuberculosis in area B was 5.6
## (95% CI 4.6 to 6.6) cases per 100 herds at risk.

---

Event instantaneous hazard based on Kaplan-Meier survival estimates

Description

Compute event instantaneous hazard on the basis of a Kaplan-Meier survival function.

Usage

epi.insthaz(survfit.obj, conf.level = 0.95)

Arguments

survfit.obj	a survfit object, computed using the survival package.
conf.level	magnitude of the returned confidence interval. Must be a single number between 0 and 1.

Details

Computes the instantaneous hazard of the event of interest, equivalent to the proportion of the group at risk failing per unit time.

Value

A data frame with the following variables: strata the strata identifier, time the observed event times, n.risk the number of individuals at risk at the start of the event time, n.event the number of individuals that experienced the event of interest at the event time, n.censor the number of individuals censored at the event time, sest the observed Kaplan-Meier survival estimate at the event time, slow the lower bound of the confidence interval for the observed Kaplan-Meier survival estimate at the event time, supp the upper bound of the confidence interval for the observed Kaplan-Meier survival estimate at the event time, hest the observed instantaneous hazard at the event time, hlow the lower bound of the confidence interval for the observed instantaneous hazard at the event time, and hupp the upper bound of the confidence interval for the observed instantaneous hazard at the event time.

References

Venables W, Ripley B (2002). Modern Applied Statistics with S, fourth edition. Springer, New York, pp. 353 - 385.

Singer J, Willett J (2003). Applied Longitudinal Data Analysis Modeling Change and Event Occurrence. Oxford University Press, London, pp. 348.

Examples

## EXAMPLE 1:
library(survival)
lung.df01 <- lung

lung.df01$status <- ifelse(lung.df01$status == 1, 0, lung.df01$status)
lung.df01$status <- ifelse(lung.df01$status == 2, 1, lung.df01$status)
lung.df01$sex <- factor(lung.df01$sex, levels = c(1,2), 
   labels = c("Male","Female"))

lung.km01 <- survfit(Surv(time = time, event = status) ~ 1, data = lung.df01)
lung.haz01 <- epi.insthaz(survfit.obj = lung.km01, conf.level = 0.95)

lung.shaz01 <- data.frame(
  time = lowess(lung.haz01$time, lung.haz01$hlow, f = 0.20)$x,
  shest =  lowess(lung.haz01$time, lung.haz01$hest, f = 0.20)$y,
  shlow =  lowess(lung.haz01$time, lung.haz01$hlow, f = 0.20)$y,
  shupp =  lowess(lung.haz01$time, lung.haz01$hupp, f = 0.20)$y)

## What was the maximum follow-up time? Use this to guide limits for the
## horizontal axis:
summary(lung.haz01$time)

## Maximum follow up time 883 days, so set horizontal axis limit to a 
## bit less, say 800 days. This will avoid instability in the smoothed results
## at the extremes of the data. 

plot(x = lung.haz01$time, y = lung.haz01$hest, xlab = "Time (days)", 
   ylab = "Daily probability of event", type = "s", 
   col = "grey", xlim = c(0,800), ylim = c(0, 0.05))
lines(x = lung.shaz01$time, y = lung.shaz01$shest, 
      lty = 1, lwd = 2, col = "black")
lines(x = lung.shaz01$time, y = lung.shaz01$shlow, 
      lty = 2, lwd = 1, col = "black")
lines(x = lung.shaz01$time, y = lung.shaz01$shupp, 
      lty = 2, lwd = 1, col = "black")

## Not run:  
library(ggplot2)

ggplot() +
  theme_bw() +
  geom_step(data = lung.haz01, aes(x = time, y = hest), colour = "grey") + 
  geom_line(data = lung.shaz01, aes(x = time, y = shest), colour = "black", 
     linewidth = 0.75, linetype = "solid") +
  geom_line(data = lung.shaz01, aes(x = time, y = shlow), colour = "black", 
     linewidth = 0.50, linetype = "dashed") +
  geom_line(data = lung.shaz01, aes(x = time, y = shupp), colour = "black", 
     linewidth = 0.50, linetype = "dashed") +
  scale_x_continuous(limits = c(0,800), name = "Time (days)") +
  scale_y_continuous(limits = c(0,0.10), name = "Daily probability of event") 

## End(Not run)


## EXAMPLE 2:
## Stratify the analyses by sex:
lung.km02 <- survfit(Surv(time = time, event = status) ~ sex, data = lung.df01)
lung.haz02 <- epi.insthaz(survfit.obj = lung.km02, conf.level = 0.95)

## Split the data by sex:
lung.split02 <- split(x = lung.haz02, f = lung.haz02$strata)

## Loess smooth teh data for each group using lapply:
lung.loess02 <- lapply(X = lung.split02, FUN = function(dat, span = 0.4) {
   hest <- loess(hest ~ time, data = dat, span = span)$fit
   hlow <- loess(hlow ~ time, data = dat, span = span)$fit
   hupp <- loess(hupp ~ time, data = dat, span = span)$fit
  data.frame(strata = dat$strata, time = dat$time, hest = hest, 
   hlow = hlow, hupp = hupp)
})

## Combine the loess smoothed results into a single data frame:
lung.shaz02 <- do.call(rbind, lung.loess02)
row.names(lung.shaz02) <- NULL

## Not run:  
library(ggplot2)

## Use the same horizontal axis limits calculated above:
ggplot() +
  theme_bw() +
  geom_step(data = lung.haz02, aes(x = time, y = hest), colour = "grey") + 
  facet_grid(strata ~ .) +
  geom_ribbon(data = lung.shaz02, aes(x = time, ymin = hlow, ymax = hupp), 
     alpha = 0.25) +
  scale_x_continuous(limits = c(0,800), name = "Time (days)") +
  scale_y_continuous(limits = c(0,0.10), name = "Daily probability of event")

## End(Not run)

---