#install.packages("AF")
# Simulate a cross-sectional sample
library(AF)
## Warning: package 'AF' was built under R version 4.0.4
## Loading required package: survival
## Loading required package: drgee
## Warning: package 'drgee' was built under R version 4.0.4
## Loading required package: stdReg
## Warning: package 'stdReg' was built under R version 4.0.4
## Loading required package: data.table
## Loading required package: ivtools
## Warning: package 'ivtools' was built under R version 4.0.4
expit <- function(x) 1 / (1 + exp( - x))
n <- 1000
Z <- rnorm(n = n)
X <- rbinom(n = n, size = 1, prob = expit(Z))
Y <- rbinom(n = n, size = 1, prob = expit(Z + X))

# Example 1: non clustered data from a cross-sectional sampling design
data <- data.frame(Y, X, Z)
head(data)
##   Y X           Z
## 1 1 0 -0.53660410
## 2 1 0 -1.11277143
## 3 1 0  0.21553832
## 4 1 0  1.40547736
## 5 0 0 -0.05393273
## 6 1 0  0.16162830
# Fit a glm object
fit <- glm(formula = Y ~ X + Z + X * Z, family = binomial, data = data)

# Estimate the attributable fraction from the fitted logistic regression
AFglm_est <- AFglm(object = fit, data = data, exposure = "X")
summary(AFglm_est)
## Call:  
## AFglm(object = fit, data = data, exposure = "X")
## 
## Estimated attributable fraction (AF) and untransformed 95% Wald CI: 
## 
##         AF  Std.Error  z value     Pr(>|z|) Lower limit Upper limit
##  0.1541473 0.02566699 6.005663 1.905515e-09   0.1038409   0.2044536
## 
## Exposure : X 
## Outcome  : Y 
## 
##  Observations Cases
##          1000   603
## 
## Method for confounder adjustment:  Logistic regression 
## Call:  
## glm(formula = Y ~ X + Z + X * Z, family = binomial, data = data)
# Example 2: clustered data from a cross-sectional sampling design
# Duplicate observations in order to create clustered data
id <- rep(1:n, 2)
data <- data.frame(id = id, Y = c(Y, Y), X = c(X, X), Z = c(Z, Z))

# Fit a glm object
fit <- glm(formula = Y ~ X + Z + X * Z, family = binomial, data = data)

# Estimate the attributable fraction from the fitted logistic regression
AFglm_clust <- AFglm(object = fit, data = data,
                     exposure = "X", clusterid = "id")
summary(AFglm_clust)
## Call:  
## AFglm(object = fit, data = data, exposure = "X", clusterid = "id")
## 
## Estimated attributable fraction (AF) and untransformed 95% Wald CI: 
## 
##         AF  Robust SE  z value     Pr(>|z|) Lower limit Upper limit
##  0.1541473 0.02566699 6.005663 1.905515e-09   0.1038409   0.2044536
## 
## Exposure : X 
## Outcome  : Y 
## 
##  Observations Cases Clusters
##          2000  1206     1000
## 
## Method for confounder adjustment:  Logistic regression 
## Call:  
## glm(formula = Y ~ X + Z + X * Z, family = binomial, data = data)
# Example 3: non matched case-control
# Simulate a sample from a non matched case-control sampling design
# Make the outcome a rare event by setting the intercept to -6

expit <- function(x) 1 / (1 + exp( - x))
NN <- 1000000
n <- 500
intercept <- -6
Z <- rnorm(n = NN)
X <- rbinom(n = NN, size = 1, prob = expit(Z))
Y <- rbinom(n = NN, size = 1, prob = expit(intercept + X + Z))
population <- data.frame(Z, X, Y)
Case <- which(population$Y == 1)
Control <- which(population$Y == 0)
# Sample cases and controls from the population
case <- sample(Case, n)
control <- sample(Control, n)
data <- population[c(case, control), ]

# Fit a glm object
fit <- glm(formula = Y ~ X + Z + X * Z, family = binomial, data = data)

# Estimate the attributable fraction from the fitted logistic regression
AFglm_est_cc <- AFglm(object = fit, data = data, exposure = "X", case.control = TRUE)
summary(AFglm_est_cc)
## Call:  
## AFglm(object = fit, data = data, exposure = "X", case.control = TRUE)
## 
## Estimated attributable fraction (AF) and untransformed 95% Wald CI: 
## 
##         AF  Std.Error  z value     Pr(>|z|) Lower limit Upper limit
##  0.6452836 0.06399425 10.08346 6.538119e-24   0.5198572     0.77071
## 
## Exposure : X 
## Outcome  : Y 
## 
##  Observations Cases
##          1000   500
## 
## Method for confounder adjustment:  Logistic regression 
## Call:  
## glm(formula = Y ~ X + Z + X * Z, family = binomial, data = data)
################################################
#REF:https://rdrr.io/cran/AF/man/AFglm.html
#https://pubmed.ncbi.nlm.nih.gov/26992709/