Difference between revisions of "Attributable risk"

Revision as of 07:02, 25 April 2016

Attributable risk is a fraction of total risk that can be attributed to a particular cause. There are a few different ways to calculate it. Population attributable fraction of an exposure agent is the fraction of disease that would disappear if the exposure to that agent would disappear in a population. Etiologic fraction is the fraction of cases that have occurred earlier than they would have occurred (if at all) without exposure. Etiologic fracion cannot typically be calculated based on risk ratio (RR) alone, but it requires knowledge about biological mechanisms.

Question

How to calculate attributable risk? What different approaches are there, and what are their differences in interpretation and use?

Answer

Risk ratio (RR)

risk among the exposed divided by the risk among the non-exposed

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): RR = \frac{R_1}{R_0}.

Attributable fraction

(aka excess fraction) the fraction of cases among the exposed that would not have occurred if the exposure would not have taken place:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): AF = \frac{RR - 1}{RR}

Population attributable fraction

the fraction of cases among the total population that would not have occurred if the exposure would not have taken place. The most useful formulas are

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = 1 - \frac{1}{\sum_{i=0}^k p_i (RR_i)}

for use with several population subgroups (typically with different exposure levels). Not valid when confounding exists. Subscript i refers to the ith subgroup. p_i = proportion of total population in ith subgroup.

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = 1- \sum_{i=0}^k \frac{p_{ci}}{RR_i} = \sum_i p_{ci} \frac{p_{ei}(RR_i - 1)}{p_{ei}(RR_i - 1) + 1}

which produces valid estimates when confounding exists but with a problem that parameters are often not known. p_ci is the proportion of cases falling in subgroup i (so that Σ_ip_ci = 1), p_ei is the proportion of exposed people within subgroup i (and 1-p_i is the fraction of unexposed)

Etiologic fraction

Fraction of cases among the exposed that would have occurred later (if at all) if the exposure would not have taken place. It cannot be calculated without understanding of the biological mechanism, but there are equations for several different models:

With this code, you can compare attributable fraction and lower (assuming exponential survival distribution) and upper bounds of etiological fraction.

library(OpasnetUtils)
AF <- function(x) {return(data.frame(RR = x, AF = (x-1)/x, EF_exp_lower = (x-1)/x^(x/(x-1)), EF_upper = 1))}

oprint(AF(RR))

This code creates a simulated population of 180 individuals and calculates their survival and attributable and etiologic fractions in different mechanistic settings.

#This is code 6211/ on page [[Attributable risk]]
library(OpasnetUtils)
library(reshape2)
library(ggplot2)

lifetime <- data.frame(
  Id = 1:180,
  BAU = 10 - 0:179 * 8 / 180
)

scenarios <- c(
	"BAU", 
	"Orderly", 
	"Relative", 
	"Maxloss", 
	"Competcause", 
	"EF0.25", 
	"BAU_expon", 
	"Rel_expon",
	"Rel_exMaxloss"
)

bau <- scenarios[bau]
scenario <- scenarios[scenario]
scenario <- union(bau, scenario)

scenario <- c("Id", scenario)

lifetime$Orderly <- lifetime$BAU -1
#RR <- sum(lifetime$BAU) / sum(lifetime$Orderly)
lifetime$Relative <- lifetime$BAU / 1.2
lifetime$Maxloss <- lifetime$Orderly[c(156:180, 1:155)]
#lifetime$Submaxloss <- lifetime$Orderly[c((1:30)*6, rep(0:29, each = 5)*6+(1:5))]
temp <- numeric()
for(i in 0:29) {
  temp <- c(temp, i * 5 + 1:5, 151 + i)
}
lifetime$Competcause <- lifetime$Orderly[temp]
a <- 0.3
lifetime$EF0.25 <- ifelse(lifetime$Id/4 == round(lifetime$Id/4), lifetime$BAU * a, lifetime$BAU)
lifetime$BAU_expon <- qexp((180:1)/181, 1/6)
lifetime$Rel_expon <- qexp((180:1)/181, 1.2/6)
lifetime$Rel_exMaxloss <- lifetime$Rel_expon[c(168:180, 1:167)]
lifetime <- lifetime[scenario]

objects.latest("Op_en6211", code_name = "EF")

metrices <- EvalOutput(metrices)

cat("Relative risks.\n")
oprint(RR@output)

cat("Different etiologic and attributable fractions.\n")
oprint(unkeep(metrices, sources = TRUE))

BS <- 24
ggplot(lif@output, aes(x = Id, y = lifResult, colour = Scenario))+geom_point()+
  coord_cartesian(ylim=c(0,10)) + theme_gray(base_size = BS) + labs(title = "Life expectancies of 180 individuals", y = "Age at death", x = "Individual")
ggplot(fr@output, aes(x = Time, y = frResult, colour = Scenario, group = Scenario))+geom_line() + 
  theme_gray(base_size = BS) + labs(title = "fraction of people dying at different time groups")
ggplot(surv@output, aes(x = Time, y = survResult, colour = Scenario, group = Scenario))+geom_line() + 
  theme_gray(base_size = BS) + labs(title = "Survival curves in different scenarios")
ggplot(EF_eq9@output, aes(x = Time, y = EF_eq9Result, colour = Scenario, group = Scenario))+geom_line()+
  theme_gray(base_size = BS) + labs(title = "Development of etiologic fraction in time")

Rationale

Definitions of terms

There are several different kinds of proportions that sound alike but are not. Therefore, we explain the specific meaning of several terms.

Total population (N)

The number of people in the total population considered, including cases, non-cases, exposed and non-exposed.

Classifications

There are three classifications, and every person in the total population belongs to exactly one group in each classification.

• Disease (D): classes case (C) and non-case (nc)
• Exposure (E): classes exposed (1) and non-exposed (0)
• Population subgroup (S): classes i = 1, 2, ..., k (typically based on different exposure levels)

Attributable fraction (AF)

The proportion of cases caused by exposure among all cases (in the subgroup)

Proportion exposed (p_e, p_ei)

proportion of exposed among the total population or within subgroup i: p_e = N(E=1)/N, p_ei = N(E=1,S=i)/N(S=i)

Proportion of population (p_i)

proportion of population in subgroups i among the total population: N(S=i)/N

Proportion of cases (p_ci)

proportion of cases in subgroups i among the total cases: N(D=c,S=i)/N(D=c)

Attributable fraction

Rockhill et al.^[1] give an extensive description about different ways to calculate attributable fraction (AF) and population attributable fraction (PAF) and assumptions needed in each approach. Modern Epidemiology ^[2] is the authoritative source of epidemiology. They first define attributable fraction AF for a cohort of people (pages 295-297). It is the fraction of cases among the exposed that would not have occurred if the exposure would not have taken place.R↻

Impact of confounders

The problem with the two PAF equations (see [[#Answer|]]) is that the former has easier-to-collect input, but it is not valid if there is confounding. It is still often mistakenly used. The latter equation would produce an unbiased estimate, but the data needed is harder to collect. Darrow and Steenland^[5] have studied the impact of confounding on the bias in attributable fraction. This is their summary:

Population attributable fraction

The population attributable fraction PAF is that fraction among the whole cohort:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = \frac{N_1 (R_1 - R_0)}{N_1 R_1 + N_0 R_0} = \frac{N_1 (R_1 - R_0)/R_0}{N_1 R_1/R_0 + N_0 R_0/R_0} = \frac{N_1 (RR - 1)}{N_1 RR + N_0}

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): = \frac{ \frac{N_1 (RR - 1)}{N_1 + N_0} }{ \frac{N_1 RR + N_0}{N_1 + N_0}} = \frac{ p (RR - 1) }{ \frac{N_1 RR - N_1 + (N_1 + N_0)}{N_1 + N_0}} = \frac{p (RR - 1)}{p RR - p + 1} = \frac{p (RR - 1)}{p (RR - 1) + 1},

where

• N₁ and N₀ are the numbers of exposed and unexposed people, respectively,
• R₁ and R₀ are the risks of disease in the exposed and unexposed group, respectively, and RR = R₁ / R₀,
• p is the fraction of exposed people among the whole cohort.

Note that there is a typo in the Modern Epidemiology book: the denominator should be p(RR-1)+1, not p(RR-1)-1.

Population attributable fraction can be calculated as a weighted average based on subgroup data:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = \Sigma_i p_{ci} PAF_{i},

where

• p_ci is the proportion of cases falling in stratum (subgroup) i,
• PAF_i is the population attributable fraction calculated for the subgroup.

Specifically, we can divide the cohort into subgroups based on exposure (in the simplest case exposed and unexposed), so we get

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = p_c \frac{1(RR - 1)}{1(RR - 1) + 1} + (1 - p_c) \frac{0(RR - 1)}{0(RR - 1) +1} = p_c \frac{RR - 1}{RR},

where p_c is the proportion of cases in the exposed group among all cases; this is the same as exposure prevalence among cases.

WHO approach

PAF is ^[9]

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = \frac{\sum_{i=0}^k P_i RR_i - \Sigma_{i=0}^k P'_i RR_i}{\Sigma_{i=0}^k P_i RR_i}

where i is a certain exposure level, P is the fraction of population in that exposure level, RR is the relative risk at that exposure level, and P' is the fraction of population in a counterfactual ideal situation (where the exposure is typically lower).

Based on this, we can limit our examination to a situation where there are only two population groups, one exposed to background level (with relative risk 1) and the other exposed to a higher level (with relative risk RR). In the counterfactual situation nobody is exposed. Thus, we get

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = \frac{(P RR + (1-P)*1) - (0*RR + 1*1)}{P RR + (1-P)*1}

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = \frac{P RR - P}{P RR + 1 - P}

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): PAF = \frac{P(RR - 1)}{P(RR -1) + 1}

This equation is used in e.g. Health impact assessment.

Constant background assumption

⇤# : Is this section necessary? --Jouni (talk) 20:53, 7 April 2016 (UTC)

p_ci can be calculated for each subgroup with the following equation if the background risk of disease is equal in all subgroups (and thus cancels out):

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): p_{ci} = \frac{N_i \Pi_j RR_{i,j}}{\Sigma_i N_i \Pi_j RR_{i,j}},

where

• N_i is the number of people in each subgroup i,
• RR_i,j is the risk ratio in subgroup i due to pollutant j (accounting for the estimated exposure in the subgroup). Note that this assumes that multiplicative assumption holds between different pollutant effects.

This page does not contain R code. Instead, it is written as part of the model in Health impact assessment.

p_c can be calculated by first calculating number of cases in each subgroup:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): cases_i = N_i * background * \Pi_j e^{ln(ERF_{j}) exposure_{i,j}},

where

• cases_i is the number of cases in subgroup i,
• N_i is the number of people in subgroup i,
• background is the background risk of the disease in the unexposed; we assume that it is the same in all subgroups,
• ERF_j is the risk ratio for unit exposure for each pollutant j (if the exposure response function ERF assumes another form than relative risk, i.e. exponential, then another equations must be used),
• exposure_i,j is the amount of exposure in a subgroup i to pollutant j.

Therefore,

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): p_{ci} = \frac{cases_i}{\Sigma_i cases_i} = \frac{N_i * background * \Pi e^{ln(ERF_{j}) exposure_{i,j}}}{background \Sigma N_i \Pi e^{ln(ERF_{j}) exposure_{i,j}}}

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): p_{ci} = \frac{N_i \Pi_j RR_{i,j}}{\Sigma_i N_i \Pi_j RR_{i,j}},

where RR_i,j = exp(ln(ERF_j) exposure_i,j).

In addition, if only fraction p of the population is exposed, for the total population we get

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): RR = \frac{p * N * background * RR_{exposed} + (1-p) * N * background * RR_{unexposed}}{N * background * RR_{unexposed}}

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): = \frac{p e^{ln(ERF)exposure} + (1-p)1}{1} = p e^{ln(ERF)exposure} -p + 1

Etiologic fraction

Etiologic fraction is defined as the fraction of cases that is advanced in time because of exposure.^[10]R↻ It can also be called probability of causation, which has importance in court. It can also be used to calculate premature cases, but that word is ambiguous and sometimes attributable fraction is used instead.R↻ Therefore, it is important to explicitly explain what is meant by the work premature.

The exact value of etiologic fraction cannot be estimated directly from risk ratio (RR) because some knowledge is needed about biological mechanisms (more precisely: timing of disease). In any case, the etiologic fraction always lies between f and 1, when f is

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): \frac{RR - 1}{RR^{RR/(RR-1)}}.

The code below calculates the attributable fraction and lower and upper bounds of the etiological fraction for user-defined RRs.

Estimating etiologic fraction

Etiologic fraction (EF) tells which fraction of the population dies earlier because of the exposure considered, compared with the situation that they would not be exposed. Often that is calculated is the attributable fraction AF:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): AF = \frac{RR - 1}{RR},

where RR is the relative risk: risk of exposed people divided by the risk of non-exposed. Robins and Greenland^[10] studied the estimability of etiologic fraction. They concluded that observations are not enough to conclude about the precise value of EF, because irrespective of observation, the same amount of observed life years lost may be due to many people losing a short time each, or due to a few losing a long time each. The upper limit in theory is always 1, and the lower bound they estimated by this equation (equation 9 in the article):

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): \frac{\int_G [f_1(u) - f_0(u)]du}{1 - S_1(t)},

where 1 means the exposed group, 0 means the non-exposed group, f is the proportion of population dying particular time points, S is the survival function, t is the length of the observation time, u the observation time and G is the set of all u < t such that f₁(u) > f₀(u).

From this, they derive the following equation (number 11 in the article):

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): \frac{RR-1}{RR^{RR/(RR-1)}}

It is not clear how they got from the first equation to the second equation, especially because the second gives values that are typically less than half of that if the first one. Both are used in the code below. In addition, the true etiologic fraction is calculated for this simulated population, because in the simulation we assume that we know exactly what happens to each individual in each scenario and how much their lengths of lives change.

Calculations

⇤# : UPDATE AF TO REFLECT THE CURRENT IMPLEMENTATION OF ERF Exposure-response function --Jouni (talk) 05:20, 13 June 2015 (UTC)

# This is code Op_en6211/AF on page [[Attributable risk]]
# Parameters: none

library(OpasnetUtils)

# AF = attributable fraction
# EF = etiologic fraction
# PAF = population attributable fraction using 
EF <- Ovariable("EF", 
	dependencies = data.frame(Name = c(
		"RR" # Risk ratio
	)),
	
	formula = function(...) {

		R <- unkeep(RR, sources = TRUE, prevresults = TRUE)
		EF <- (RR - 1) / R^(R/(R-1))
		EF <- EF * Ovariable("temp", data = data.frame(
			EFestimate = c("Low", "High"),
			Result = 1
		))
		result(EF)[EF$EFestimate == "High"] <- 1

		return(EF)
	}
)

AF <- Ovariable("AF", 
	dependencies = data.frame(Name = c(
		"RR" # Risk ratio
	)),
	
	formula = function(...) {

		AF <- (RR - 1) / unkeep(RR, sources = TRUE, prevresults = TRUE)

		return(AF)
	}
)

PAF <- Ovariable("PAF", 
	dependencies = data.frame(Name = c(
		"RR", # Risk ratio
		"pci", # proportion of cases falling subgroup i among all cases
		"pei" # proportion of exposed people within subgroup i
	)),
	
	formula = function(...) {

		peirri <- pei * (RR - 1)
		peirri <- unkeep(peirri, sources = TRUE, prevresults = TRUE)

		PAF <- pci * peirri / (peirri + 1) # The population subgroup could be summed up.

		return(PAF)
	}
)

objects.store(EF, AF, PAF)
cat("Ovariables EF, AF, PAF stored.\n")

A previous version of code looked at RRs of all exposure agents and summed PAFs up.

Some interesting model runs:

• Population with exponentially distributed lifetimes
• Life loss to a fraction of people across the whole population
• [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=MR9wpPlPiKmsr7mD Give all life loss to the people living the longest9

#This is code Op_en6211/EF on page [[Attributable risk]]

library(OpasnetUtils)

lif <- Ovariable(
  "lif",
  dependencies = data.frame(Name = "lifetime"),
  formula = function(...) {
    out <- melt(
       lifetime, 
       id.vars = "Id", 
       value.name = "Result", 
       variable.name = "Scenario"
    )
    out <- Ovariable(
      output = out, 
      marginal = c(TRUE, TRUE, FALSE)
    )
    return(out)
  }
)

le <- Ovariable(
  "le",
  dependencies = data.frame(Name = "lif"),
  formula = function(...) {
    le <- oapply(lif, INDEX = "Scenario", FUN = sum) / 
    oapply(lif, INDEX = "Scenario", FUN = length)
    return(le)
  }
)

RR <- Ovariable(
  "RR",
  dependencies = data.frame(Name = "le"),
  formula = function(...) {
    RR <- le[le$Scenario == bau , ]
    RR <- unkeep(RR, cols = c("Scenario", "lifResult"))
    RR <- RR / le
    RR <- unkeep(RR, prevresults = TRUE)
    return(RR)
  }
)

fr <- Ovariable("fr", 
  dependencies = data.frame(
    Name = "lif"
  ),
  formula = function(...) {
    out <- lif
    temp2 <- cut(result(out), breaks = 12)
    out$Time <- temp2
    out <- out * 0 + 1/oapply(out, cols = c("Id", "Time"), FUN = length)
    temp <- Ovariable(
      "temp", 
      data = data.frame(
        Time = levels(temp2),
        Result = 0
      )
    )
    out <- combine(EvalOutput(temp), out)
    out <- oapply(out, cols = "Id", FUN = sum) # Automatic fillna is OK.
    return(out)
  }
)

surv <- Ovariable(
  "surv",
  dependencies = data.frame(Name = "fr"),
  formula = function(...) {
    out <- fr[order(fr$Time) , ]
    temp <- data.frame()
    for(i in unique(out$Scenario)) {
      temp2 <- out[out$Scenario == i , ]
      result(temp2) <- 1 - cumsum(result(temp2))
      temp <- rbind(temp, temp2@output)
    }
    out@output <- temp
    return(out)
  }
)

EF_eq9 <- Ovariable(
  "EF_eq9",
  dependencies = data.frame(Name = c("fr", "surv")),
  formula = function(...) {
    BAU <- fr[fr$Scenario == bau , ]
    BAU <- unkeep(BAU, prevresults = TRUE, sources = TRUE, cols = "Scenario")
    
    out <- fr
    result(out) <- pmax(0, result(out - BAU))
    
    out <- out[order(out$Time) , ]
    temp <- data.frame()
    for(i in unique(out$Scenario)) {
      temp2 <- out[out$Scenario == i , ]
      result(temp2) <- cumsum(result(temp2))
      temp <- rbind(temp, temp2@output)
    }
    out@output <- temp
    out <- out / (1 - surv)
    
    return(out)
  }
)

EF_true <- Ovariable(
  "EF_true", 
  dependencies = data.frame(Name = "lif"),
  formula = function(...) {
    BAU <- lif[lif$Scenario == bau , ]
    BAU <- unkeep(BAU, cols = "Scenario", prevresults = TRUE, sources = TRUE)
    out <- lif < BAU
    out <- oapply(out, cols = "Id", FUN = sum) / 
      oapply(out, cols = "Id", FUN = length)
    
    return(out)
  }
)

metrices <- Ovariable(
  "metrices", 
  dependencies = data.frame(Name = c("RR", "lif", "EF_true", "EF_eq9")),
  formula = function(...) {
    out <- (RR - 1) / RR
    out$Metric <- "Attributable fraction"
    temp <- (RR - 1)/(RR^(RR/(RR-1)))
    temp$Metric <- "EF_low from eq 11"
    out <- combine(out, temp)
#    result(temp) <- 1
#    temp$Metric <- "EF_up theoretical"
#    out <- combine(out, temp)
    temp <- unkeep(EF_true, sources = TRUE, prevresults = TRUE)
    temp$Metric <- "EF_true"
    out <- combine(out, temp)
    temp <- unkeep(EF_eq9[EF_eq9$Time == levels(EF_eq9$Time)[length(levels(EF_eq9$Time))] , ],
           cols = "Time", sources = TRUE, prevresults = TRUE
    )
    temp$Metric <- "EF_low from Eq 9"
    out <- combine(out, temp)
    

    return(out)
  }
)

objects.store(lif, le, RR, fr, surv, EF_eq9, EF_true, metrices)
cat("Ovariables lif, le, RR, fr, surv, EF_eq9, EF_true, metrices stored.\n")

See also

• Attributable risk in Wikipedia.

References