paf: Population Attributable Fraction

In pifpaf: Potential Impact Fraction and Population Attributable Fraction for Cross-Sectional Data

Description

Function for estimating the Population Attributable Fraction paf from a cross-sectional sample of the exposure X with a known Relative Risk function rr with meta-analytical parameter theta, where the Population Attributable Fraction is given by:

PAF = mean(rr(X; theta) - 1)/mean(rr(X; theta)).

Usage

paf(X, thetahat, rr, method = "empirical",
  weights = rep(1/nrow(as.matrix(X)), nrow(as.matrix(X))), Xvar = var(X),
  deriv.method.args = list(), deriv.method = "Richardson", adjust = 1,
  n = 512, ktype = "gaussian", bw = "SJ", check_exposure = TRUE,
  check_rr = TRUE, check_integrals = TRUE)

Arguments

X	Random sample (data.frame) which includes exposure and covariates or sample mean if "approximate" method is selected.
thetahat	Asymptotically consistent or Fisher consistent estimator (vector) of theta for the Relative Risk function rr.
rr	function for Relative Risk which uses parameter theta. The order of the parameters should be rr(X, theta). **Optional**
method	Either "empirical" (default), "kernel" or "approximate". For details on estimation methods see pif.
weights	Normalized survey weights for the sample X.
Xvar	Variance of exposure levels (for "approximate" method).
deriv.method.args	method.args for hessian (for "approximate" method).
deriv.method	method for hessian. Don't change this unless you know what you are doing (for "approximate" method).
adjust	Adjust bandwith parameter (for "kernel" method) from density.
n	Number of equally spaced points at which the density (for "kernel" method) is to be estimated (see density).
ktype	kernel type: "gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine" (for "kernel" method). Additional information on kernels in density.
bw	Smoothing bandwith parameter (for "kernel" method) from density. Default "SJ".
check_exposure	boolean Check that exposure X is positive and numeric.
check_rr	boolean Check that Relative Risk function rr equals 1 when evaluated at 0.
check_integrals	boolean Check that counterfactual cft and relative risk's rr expected values are well defined for this scenario.

Details

The Relative Risk function rr and counterfactual cft should consider X to be a data.frame. Each row of X is an "individual" and each column a variable (exposure or covariate) of said individual.

Value

paf Estimate of Population Attributable Fraction.

Note

"approximate" method should be the last choice. In practice "empirical" should be preferred as convergence is faster in simulations than "kernel". In addition, the scope of "kernel" is limited as it does not work with multivariate exposures X.

paf is a wrapper for pif with counterfactual of theoretical minimum risk exposure (rr = 1).

For more information on kernels see density.

Do not use the $ operator when using "approximate" method.

Author(s)

Rodrigo Zepeda-Tello [e-mail]

Dalia Camacho-Garc<c3><ad>a-Forment<c3><ad> [e-mail]

References

Vander Hoorn, S., Ezzati, M., Rodgers, A., Lopez, A. D., & Murray, C. J. (2004). Estimating attributable burden of disease from exposure and hazard data. Comparative quantification of health risks: global and regional burden of disease attributable to selected major risk factors. Geneva: World Health Organization, 2129-40.

See Also

paf.confidence for confidence interval estimation, pif for Potential Impact Fraction estimation.

See paf.exponential and paf.linear for fractions with ready-to-use exponential and linear Relative Risks respectively.

Sensitivity analysis plots can be done with paf.plot, and paf.sensitivity.

Examples

#Example 1: Exponential Relative Risk
#--------------------------------------------
set.seed(18427)
X        <- data.frame(Exposure = rnorm(100,3,1))
thetahat <- 0.12
rr       <- function(X, theta){exp(theta*X)}

#Using the empirical method
paf(X, thetahat, rr)

#Same example with kernel method
paf(X, thetahat, rr, method = "kernel")

#Same example with approximate method
Xmean <- data.frame(Exposure = mean(X[,"Exposure"]))
Xvar  <- var(X[,"Exposure"])
paf(Xmean, thetahat, rr, method = "approximate", Xvar = Xvar)

#Additional options for approximate:
paf(Xmean, thetahat, rr, method = "approximate", Xvar = Xvar, 
   deriv.method = "Richardson",  deriv.method.args = list(eps=1e-3, d=0.1))

#Example 2: Linear Relative Risk with weighted sample
#--------------------------------------------
set.seed(18427)
X                   <- data.frame(Exposure = rbeta(100,3,1))
weights             <- runif(100)
normalized_weights  <- weights/sum(weights)
thetahat            <- 0.12
rr                  <- function(X, theta){theta*X^2 + 1}
paf(X, thetahat, rr, weights = normalized_weights)

   
#Additional options for kernel:
paf(X, thetahat, rr, weights = normalized_weights, 
     method = "kernel", ktype = "cosine", bw = "nrd0")


#Example 3: Multivariate Linear Relative Risk
#--------------------------------------------
set.seed(18427)
X1       <- rnorm(100,4,1)
X2       <- rnorm(100,2,0.4)
X        <- data.frame(Exposure = X1, Covariate = X2)
thetahat <- c(0.12, 0.03)

#When creating relative risks avoid using the $ operator
#as it doesn't work under approximate method
rr_not    <- function(X, theta){
               exp(theta[1]*X$Exposure + theta[2]*X$Covariate)
             }
rr_better <- function(X, theta){
               exp(theta[1]*X[,"Exposure"] + theta[2]*X[,"Covariate"])
             }
             
#For the empirical method it makes no difference:              
paf(X, thetahat, rr_better) 
paf(X, thetahat, rr_not) 


#But the approximate method crashes due to operator
Xmean <- data.frame(Exposure = mean(X[,"Exposure"]), 
                    Covariate = mean(X[,"Covariate"]))
Xvar  <- var(X)

paf(Xmean, thetahat, rr_better, method = "approximate", Xvar = Xvar)
## Not run: 
#Error: $ operator in rr definitions don't work in approximate
paf(Xmean, thetahat, rr_not, method = "approximate", Xvar = Xvar)

## End(Not run)


## Not run: 
#Error: Multivariate cases cannot be evaluated with kernel method
paf(X, thetahat, rr, method = "kernel") 

## End(Not run)

#Example 4: Categorical Relative Risk & Exposure
#--------------------------------------------
set.seed(18427)
mysample  <- sample(c("Normal","Overweight","Obese"), 100, 
                   replace = TRUE, prob = c(0.4, 0.1, 0.5))
X        <- data.frame(Exposure = mysample)

thetahat <- c(1, 1.2, 1.5)

#Categorical relative risk function
rr <- function(X, theta){

   #Create return vector with default risk of 1
   r_risk <- rep(1, nrow(X))
   
   #Assign categorical relative risk
   r_risk[which(X[,"Exposure"] == "Normal")]      <- thetahat[1]
   r_risk[which(X[,"Exposure"] == "Overweight")]  <- thetahat[2]
   r_risk[which(X[,"Exposure"] == "Obese")]       <- thetahat[3]
   
   return(r_risk)
}

paf(X, thetahat, rr, check_rr = FALSE)

#Example 5: Continuous Exposure and Categorical Relative Risk
#------------------------------------------------------------------
set.seed(18427)

#Assume we have BMI from a sample
BMI          <- data.frame(Exposure = rlnorm(100, 3.1, sdlog = 0.1))

#Theoretical minimum risk exposure is at 20kg/m^2 in borderline "Normal" category
BMI_adjusted <- BMI - 20

thetahat <- c(Malnourished = 2.2, Normal = 1, Overweight = 1.8, 
              Obese = 2.5)
rr       <- function(X, theta){
     
     #Create return vector with default risk of 1
     r_risk <- rep(1, nrow(X))
   
     #Assign categorical relative risk
     r_risk[which(X[,"Exposure"] < 0)]             <- theta[1] #Malnourished
     r_risk[intersect(which(X[,"Exposure"] >= 0), 
                      which(X[,"Exposure"] < 5))]  <- theta[2] #Normal
     r_risk[intersect(which(X[,"Exposure"] >= 5), 
                      which(X[,"Exposure"] < 10))] <- theta[3] #Overweight
     r_risk[which(X[,"Exposure"] >= 10)]           <- theta[4] #Obese
   
   return(r_risk)
}

paf(BMI_adjusted, thetahat, rr, check_exposure = FALSE)

#Example 6: Bivariate exposure and rr ("classical PAF")
#------------------------------------------------------------------
set.seed(18427)
mysample  <- sample(c("Exposed","Unexposed"), 1000, 
                replace = TRUE, prob = c(0.1, 0.9))
X         <- data.frame(Exposure = mysample)
theta     <- c("Exposed" = 2.5, "Unexposed" = 1.2)  
rr        <- function(X, theta){
   
   #Create relative risk function
   r_risk <- rep(1, nrow(X))
   
   #Assign values of relative risk
   r_risk[which(X[,"Exposure"] == "Unexposed")] <- theta["Unexposed"]
   r_risk[which(X[,"Exposure"] == "Exposed")]   <- theta["Exposed"]
   
   return(r_risk)
}    

paf(X, theta, rr)

#Example 7: Continuous exposure, several covariates
#------------------------------------------------------------------
X <- data.frame(Exposure = rbeta(100, 2, 3),
                Age      = runif(100, 20, 100),
                Sex      = sample(c("M","F"), 100, replace = TRUE),
                BMI      = rlnorm(100, 3.2, 0.2))
thetahat <- c(-0.1, 0.05, 0.2, -0.4, 0.3, 0.1)

rr <- function(X, theta){

     #Create risk vector
     Risk    <- rep(1, nrow(X))
     
     #Identify subpopulations
     males   <- which(X[,"Sex"] == "M")
     females <- which(X[,"Sex"] == "F")
     
     #Calculate population specific rr
     Risk[males] <- theta[1]*X[males,"Exposure"] + 
                                      theta[2]*X[males,"Age"]^2 + 
                                      theta[3]*X[males,"BMI"]/2 
                                     
     Risk[females] <- theta[4]*X[females,"Exposure"] + 
                                      theta[5]*X[females,"Age"]^2 + 
                                      theta[6]*X[females,"BMI"]/2 
                                     
    return(Risk)
}

paf(X, thetahat, rr)

pifpaf documentation built on May 1, 2019, 9:11 p.m.