Description Usage Arguments Value Note Author(s) See Also Examples
View source: R/pif_confidence.R
Function that estimates confidence intervals for the Potential
Impact Fraction pif
from a cross-sectional sample of
the exposure X
with a known Relative Risk function rr
with
parameter theta
, where the Potential Impact Fraction is given
by:
PIF = (mean(rr(X; theta)) - mean(rr(cft(X);theta)))/mean(rr(X; theta)) .
1 2 3 4 5 6 7 | pif.confidence(X, thetahat, rr, thetavar, cft = NA, method = "empirical",
confidence_method = "bootstrap", confidence = 95, nsim = 1000,
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_cft = TRUE, check_rr = TRUE, check_xvar = TRUE,
check_integrals = TRUE, check_thetas = TRUE, is_paf = FALSE)
|
X |
Random sample ( |
thetahat |
Asymptotically consistent or Fisher consistent estimator
( |
rr |
|
thetavar |
Estimator of variance **Optional** |
cft |
Function |
method |
Either |
confidence_method |
Either |
confidence |
Confidence level % (default |
nsim |
Number of simulations for estimation of variance. |
weights |
Normalized survey |
Xvar |
Variance of exposure levels (for |
deriv.method.args |
|
deriv.method |
|
adjust |
Adjust bandwith parameter (for |
n |
Number of equally spaced points at which the density (for
|
ktype |
|
bw |
Smoothing bandwith parameter (for
|
check_exposure |
|
check_cft |
|
check_rr |
|
check_xvar |
|
check_integrals |
|
check_thetas |
|
is_paf |
Boolean forcing evaluation of |
pifvec Vector with lower ("Lower_CI"
), and upper
("Upper_CI"
) confidence bounds for the pif
as well as
point estimate "Point_Estimate"
and estimated variance
of log(pif)
(if confidence_method
is "loglinear"
).
For more information on kernels see density
.
Do not use the $
operator when using "approximate"
method
.
Rodrigo Zepeda-Tello rzepeda17@gmail.com
Dalia Camacho-Garc<c3><ad>a-Forment<c3><ad> daliaf172@gmail.com
See paf.confidence
for confidence interval estimation of
paf
, and pif
for only point estimate.
Sensitivity analysis plots can be done with paf.plot
, and
paf.sensitivity
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 | #Example 1: Exponential Relative Risk
#--------------------------------------------
set.seed(18427)
X <- data.frame(Exposure = rnorm(100,3,1))
thetahat <- 0.12
thetavar <- 0.02
rr <- function(X, theta){exp(theta*X)}
#Counterfactual of halving exposure
cft <- function(X){ 0.5*X }
#Using bootstrap method
pif.confidence(X, thetahat, rr, thetavar, cft)
## Not run:
#Same example with loglinear method
pif.confidence(X, thetahat, rr, thetavar, cft, confidence_method = "loglinear")
#Same example with linear method (usually the widest and least precise)
pif.confidence(X, thetahat, rr, thetavar, cft, confidence_method = "linear")
#Example 2: Linear Relative Risk
#--------------------------------------------
set.seed(18427)
X <- data.frame(Exposure = rbeta(100,3,1))
thetahat <- 0.17
thetavar <- 0.01
rr <- function(X, theta){theta*X + 1}
cft <- function(X){ 0.5*X }
weights <- runif(100)
normalized_weights <- weights/sum(weights)
pif.confidence(X, thetahat, rr, thetavar, cft, weights = normalized_weights)
#Same example with more complex counterfactual that reduces
#only the values > 0.75 are halved
cft <- function(X){
#Indentify the ones with "a lot" of exposure:
where_excess_exposure <- which(X[,"Exposure"] > 0.75)
#Halve their exposure
X[where_excess_exposure, "Exposure"] <-
X[where_excess_exposure, "Exposure"]/2
return(X)
}
pif.confidence(X, thetahat, rr, thetavar, cft, weights = normalized_weights)
#Example 3: Multivariate Linear Relative Risk
#--------------------------------------------
set.seed(18427)
X1 <- rnorm(100,4,1)
X2 <- rnorm(100,2,0.4)
thetahat <- c(0.12, 0.03)
thetavar <- diag(c(0.01, 0.02))
#But the approximate method crashes due to operator
Xmean <- data.frame(Exposure = mean(X1),
Covariate = mean(X2))
Xvar <- var(cbind(X1, X2))
#When creating relative risks avoid using the $ operator
#as it doesn't work under approximate method of PAF
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"])
}
#Creating a counterfactual.
cft <- function(X){
Y <- X
Y[,"Exposure"] <- 0.5*X[,"Exposure"]
Y[,"Covariate"] <- 1.1*X[,"Covariate"] + 1
return(Y)
}
pif.confidence(Xmean, thetahat, rr_better, thetavar, cft,
method = "approximate", Xvar = Xvar)
## End(Not run)
## Not run:
#Error: $ operator in rr definitions don't work in approximate
pif.confidence(Xmean, thetahat, rr_not, thetavar, cft, method = "approximate", Xvar = Xvar)
## End(Not run)
## 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)
thetavar <- diag(c(0.1, 0.2, 0.3))
#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)
}
#Counterfactual of reducing all obesity to normality
cft <- function(X){
X[which(X[,"Exposure"] == "Obese"),] <- "Normal"
return(X)
}
pif.confidence(X, thetahat, rr, thetavar, cft, check_rr = FALSE)
#Example 5: Categorical Relative Risk & continuous exposure
#----------------------------------------------------------
set.seed(18427)
BMI <- data.frame(Exposure = rlnorm(100, 3.1, sdlog = 0.1))
#Theoretical minimum risk exposure is 20kg/m^2 in borderline "Normal" category
BMI_adjusted <- BMI - 20
thetahat <- c(Malnourished = 2.2, Normal = 1, Overweight = 1.8,
Obese = 2.5)
thetavar <- diag(c(0.1, 0.2, 0.2, 0.1))
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)
}
#Counterfactual of everyone in normal range
cft <- function(bmi){
bmi <- data.frame(rep(2.5, nrow(bmi)), ncol = 1)
colnames(bmi) <- c("Exposure")
return(bmi)
}
pif.confidence(BMI_adjusted, thetahat, rr, thetavar, cft,
check_exposure = FALSE, method = "empirical")
#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)
thetavar <- matrix(c(0.04, 0.02, 0.02, 0.03), ncol = 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)
}
#Counterfactual of reducing the exposure in half of the individuals
cft <- function(X){
#Find out which ones are exposed
Xexp <- which(X[,"Exposure"] == "Exposed")
#Use only half of the exposed randomly
reduc <- sample(Xexp, length(Xexp)/2)
#Unexpose those individuals
X[reduc, "Exposure"] <- "Unexposed"
return(X)
}
pif.confidence(X, theta, rr, thetavar, cft)
#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)
#Create variance of theta
almostvar <- matrix(runif(6^2), ncol = 6)
thetavar <- t(almostvar) %*% almostvar
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)
}
#Counterfactual of reducing BMI
cft <- function(X){
excess_bmi <- which(X[,"BMI"] > 25)
X[excess_bmi,"BMI"] <- 25
return(X)
}
pif.confidence(X, thetahat, rr, thetavar, cft)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.