## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----setup-------------------------------------------------------------------- library(deltapif) ## ----------------------------------------------------------------------------- data(dementiarisk) ## ----echo = FALSE------------------------------------------------------------- dementiarisk[,c("risk_factor","logrr","sdlog", "total", "hispanic", "asian", "black", "white")] ## ----------------------------------------------------------------------------- paf_hispanic <- paf(p = 0.069, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "Hispanic") paf_hispanic ## ----------------------------------------------------------------------------- paf_asian <- paf(p = 0.049, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "Asian") paf_black <- paf(p = 0.117, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "Black") paf_white <- paf(p = 0.084, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "White") ## ----------------------------------------------------------------------------- weights <- c("white" = 0.5784, "hispanic" = 0.1873, "black" = 0.1205, "asian" = 0.0592) #Normalized weights to sum to 1 weights <- weights / sum(weights) ## ----------------------------------------------------------------------------- paf_population <- paf_total(paf_white, paf_hispanic, paf_black, paf_asian, weights = weights, var_weights = 0, label = "All") paf_population ## ----------------------------------------------------------------------------- as.data.frame(paf_population) ## ----------------------------------------------------------------------------- attributable_cases(426.5, paf_population, variance = 2647) ## ----------------------------------------------------------------------------- 2.20 - 1.59 ## ----------------------------------------------------------------------------- 1.59 - 1.15 ## ----------------------------------------------------------------------------- log(2.20) - log(1.59) ## ----------------------------------------------------------------------------- log(1.59) - log(1.15) ## ----------------------------------------------------------------------------- pif_hispanic <- pif(p = 0.069, p_cft = 0.069*0.85, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "Hispanic") pif_asian <- pif(p = 0.049, p_cft = 0.049*0.85, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "Asian") pif_black <- pif(p = 0.117, p_cft = 0.117*0.85, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "Black") pif_white <- pif(p = 0.084, p_cft = 0.084*0.85, beta = 0.4637340, var_beta = 0.0273858, var_p = 0, rr_link = exp, label = "White") ## ----------------------------------------------------------------------------- pif_population <- pif_total(pif_white, pif_hispanic, pif_black, pif_asian, weights = weights, var_weights = 0, label = "All") pif_population ## ----------------------------------------------------------------------------- df <- as.data.frame(pif_population, pif_white, pif_black, pif_hispanic, pif_asian) df ## ----------------------------------------------------------------------------- averted_cases(426.5, pif_population, variance = 2647) ## ----------------------------------------------------------------------------- data(alcohol) ## ----echo = FALSE------------------------------------------------------------- alcohol[1:12,c("sex", "alcohol_g", "median_intake", "age_18_plus")] ## ----------------------------------------------------------------------------- #Get alcohol intake for men divided by 10 cause RR is per 10g/day alcohol_men_intake <- c(0, 0.8, 1.6, 3.2, 7.1, 12.6, 17.4, 24.5, 34.7, 44.2, 54.4, 85.2) / 10 #Divided by 10 cause risk is per 10 g/day relative_risk <- exp(log(1.10)*alcohol_men_intake) #Divide by 100 to get the prevalence in decimals prevalence <- c(0.472, 0.014, 0.022, 0.081, 0.081, 0.069, 0.046, 0.073, 0.042, 0.032, 0.02, 0.048) #Calculate the paf paf(prevalence, beta = relative_risk, var_p = 0, var_b = 0) ## ----------------------------------------------------------------------------- #Calculate covariance of the betas given by variance*intake_i*intake_j beta_sd <- (1.12 - 1.07) / (2*qnorm(0.975)) beta_var <- matrix(NA, ncol = length(prevalence), nrow = length(prevalence)) for (k in 1:length(prevalence)){ for (j in 1:length(prevalence)){ beta_var[k,j] <- beta_sd^2*(alcohol_men_intake[j])*(alcohol_men_intake[k]) } } #Calculate the paf paf_males <- paf(prevalence, beta = relative_risk, var_p = 0, var_b = beta_var, label = "Males") paf_males ## ----------------------------------------------------------------------------- #Get alcohol intake for women divided by 10 cause RR is per 10g/day alcohol_female_intake <- c(0, 0.8, 1.6, 3.9, 7.1, 12.6, 17.4, 24.5, 34.7, 44.2, 54.4, 78.1) / 10 #Divided by 10 cause risk is per 10 g/day relative_risk <- exp(log(1.10)*alcohol_female_intake) #Divide by 100 to get the prevalence in decimals prevalence <- c(0.603, 0.02, 0.054, 0.091, 0.081, 0.061, 0.026, 0.034, 0.015, 0.007, 0.003, 0.005) #Calculate covariance of the betas given by variance*intake_i*intake_j beta_sd <- (1.12 - 1.07) / (2*qnorm(0.975)) beta_var <- matrix(NA, ncol = length(prevalence), nrow = length(prevalence)) for (k in 1:length(prevalence)){ for (j in 1:length(prevalence)){ beta_var[k,j] <- beta_sd^2*(alcohol_female_intake[j])*(alcohol_female_intake[k]) } } #Calculate the paf paf_females <- paf(prevalence, beta = relative_risk, var_p = 0, var_b = beta_var, label = "Females") paf_females ## ----------------------------------------------------------------------------- paf_total(paf_males, paf_females, weights = c(0.5, 0.5)) ## ----------------------------------------------------------------------------- data(dementiarisk) ## ----echo = FALSE------------------------------------------------------------- dementiarisk[which(dementiarisk$risk_factor %in% c("Excessive alcohol","Smoking","Physical inactivity")), c("risk_factor", "logrr", "sdlog", "total")] ## ----------------------------------------------------------------------------- paf_alcohol <- paf(p = 0.036, beta = 0.1655144, var_beta = 0.05402095^2, var_p = 0, label = "Alcohol") paf_smoking <- paf(p = 0.085, beta = 0.4637340, var_beta = 0.16548657^2, var_p = 0, label = "Smoking") paf_physical <- paf(p = 0.628, beta = 0.3293037, var_beta = 0.09296182^2, var_p = 0, label = "Physical Activity") ## ----------------------------------------------------------------------------- paf_ensemble(paf_alcohol, paf_physical, paf_smoking, weights = c(0.038, 0.311, 0.059), label = "Behavioural factors") ## ----------------------------------------------------------------------------- #Variance covariance matrix weight_variance <- matrix( c( 1.00, 0.21, -0.02, 0.21, 1.00, 0.09, -0.02, 0.09, 1.00), byrow = TRUE, ncol = 3 ) colnames(weight_variance) <- c("Alcohol", "Physical inactivity", "Smoking") rownames(weight_variance) <- c("Alcohol", "Physical inactivity", "Smoking") paf_ensemble(paf_alcohol, paf_physical, paf_smoking, weights = c(0.038, 0.311, 0.059), var_weights = weight_variance, label = "Behavioural factors") ## ----------------------------------------------------------------------------- #Calculate each of the fractions paf_alcohol <- paf(p = 0.036, beta = 0.1655144, var_beta = 0.05402095^2, var_p = 0, label = "Alcohol") paf_smoking <- paf(p = 0.085, beta = 0.4637340, var_beta = 0.16548657^2, var_p = 0, label = "Smoking") paf_physical <- paf(p = 0.628, beta = 0.3293037, var_beta = 0.09296182^2, var_p = 0, label = "Physical Activity") #Get the communality weights cweights <- c(0.038, 0.311, 0.059) #Calculate the adjusted versions weighted_adjusted_paf(paf_alcohol, paf_smoking, paf_physical, weights = cweights) ## ----------------------------------------------------------------------------- paf(p = 0.05, beta = 1.94, var_beta = 0.591^2, var_p = 0) ## ----------------------------------------------------------------------------- paf(p = 0.05, beta = 1.94, var_beta = 0.591^2, var_p = 0, link = "logit") ## ----------------------------------------------------------------------------- pif_smoking <- pif( p = 0.085, p_cft = 0.085 * (1 - 0.15), # 15% reduction beta = log(1.59), var_beta = ((log(2.20) - log(1.15)) / (2 * 1.96))^2, var_p = 0 ) pif_smoking ## ----------------------------------------------------------------------------- averted_cases(426.5, pif_smoking, variance = 2647.005) ## ----------------------------------------------------------------------------- averted_cases(426.5, pif_smoking, variance = 100000) ## ----------------------------------------------------------------------------- averted_cases(426.5, pif_smoking, variance = 100000, link = "log") ## ----echo = FALSE------------------------------------------------------------- upper_rr <- 0.7965867 lower_rr <- 0.7579948 rr_value <- 0.7735267 ## ----eval = FALSE------------------------------------------------------------- # library(EValue) # # upper_rr <- HR(0.72, rare = FALSE) |> toRR() |> as.numeric() # lower_rr <- HR(0.67, rare = FALSE) |> toRR() |> as.numeric() # rr_value <- HR(0.69, rare = FALSE) |> toRR() |> as.numeric() ## ----------------------------------------------------------------------------- pif( p = 0.438, beta = log(rr_value), p_cft = 0.876, var_beta = ((log(upper_rr) - log(lower_rr)) / (2 * 1.96))^2, quiet = TRUE )