## 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.
Run the code above in your browser using DataLab