#' Perform comparative risk assesment
#'
#' This function accepts up to eight parameters for use with the
#' parameteric physical activity model, or two vectors of quantiles
#' for the non-paraetric model. The function returns the population
#' attributable fraction.
#'
#' @param meanlog.baseline Mean of exposure in baseline (log-scale, parameteric model)
#' @param sdlog.baseline Standard deviation of exposure in baseline (log-scale, parameteric model)
#' @param p0.baseline Proportion of inactive travelers in baseline
#' @param meanlog.scenario Mean of exposure in scenario (log-scale, parameteric model)
#' @param sdlog.scenario Standard deviation of exposure in scenario (log-scale, parameteric model)
#' @param p0.scenario Proportion of inactive travelers in scenario
#' @param meanlog.leisure Mean of exposure in leisure physical activity (log-scale, parameteric model)
#' @param sdlog.leisure Standard deviation of leisure physical activity in baseline (log-scale, parameteric model)
#' @param n Number of quantiles (parametric)
#' @param B Sample size (parametric model)
#' @param P Quantiles of exposure in baseline (non-parametric model)
#' @param Q Quantiles of exposure in scenario (non-parametric model)
#' @param R Relative risk function in terms of total exposure (physical activity)
#' @param type A character string. Either "parametric" or "non-parametric"
#'
#' @return The population attributable fraction
CRA.function <- function(meanlog.baseline = log(5),
sdlog.baseline = 1,
p0.baseline = 0.5,
meanlog.scenario = log(5),
sdlog.scenario = 1,
p0.scenario = 0.5,
meanlog.leisure = log(10),
sdlog.leisure = 1,
n = 1e4,
B = 1e6,
P,
Q,
R = function(x) exp(-0.03*sqrt(x)), # shape parameter for colon cancer from Excel workbook
type = "parametric"){
if(type == "parametric"){
LA <- rlnorm(n = B, meanlog = meanlog.leisure, sdlog = sdlog.leisure)
runif <- runif(n = B)
TA.sample.baseline <- rlnorm(n = B, meanlog = meanlog.baseline, sdlog = sdlog.baseline)
TA.baseline <- ifelse(runif < p0.baseline, 0, TA.sample.baseline)
P <- quantile(TA.baseline + LA, probs = (1:(n-1))/n)
TA.sample.scenario <- rlnorm(n = B, meanlog = meanlog.scenario, sdlog = sdlog.scenario)
TA.scenario <- ifelse(runif < p0.scenario, 0, TA.sample.scenario)
Q <- quantile(TA.scenario + LA, probs = (1:(n-1))/n)
}else if( type == "non-parametric" ){
# Nothing to be processed. P and Q are being entered directly.
}else{
stop("Error with variable type. Must be parametric or non-parametric.")
}
return((sum(R(Q)) - sum(R(P)))/sum(R(P)))
}
#' Perform comparative risk assesment
#'
#' Foo
#'
#' @param leisue A character string, either "low", "medium" or "high"
#' indicating the amonut of leisure activity on the population
#' @param n Number of quantiles (parametric)
#'
#' @return A data frame with population attributable fractions
CRA.DF <- function(leisure = "medium"){
if(leisure == "low"){
leisure.value <- log(1)} # 3 min/day walking
else if(leisure == "medium"){
leisure.value <- log(10) # 30 min/day walking
}else if( leisure == "high"){
leisure.value <- log(100) # 24 hrs/day walking
}else{
message("Problem with leisure argument.")
}
CRA.interior <- function(x, y) {
cra.value <- CRA.function(meanlog.baseline = log(x$mean/sqrt(1+x$sd^2/x$mean^2)),
meanlog.scenario = log(y$mean/sqrt(1+y$sd^2/y$mean^2)),
sdlog.baseline = log(1+x$sd^2/x$mean^2),
sdlog.scenario = log(1+y$sd^2/y$mean^2),
p0.baseline = x$p0,
p0.scenario = y$p0,
meanlog.leisure = leisure.value,
type = "parametric"
)
return(cra.value)
}
return(CRA.interior)
}
#' Perform comparative risk assesment
#'
#' Foo
#'
#' @param activityFile.baseline Name of activity file for baseline
#' @param activityFile.scenario Name of activity file for scenario
#'
#' @return A data frame with population attributable fractions
#' @export
CRA <- function(activityFile.baseline = system.file("activity.baseline.csv", package = "ITHIM"), activityFile.scenario = system.file("activity.scenario.csv", package = "ITHIM")){
activityDF.baseline <- read.csv(file = activityFile.baseline, row.names = 1)
activityDF.scenario <- read.csv(file = activityFile.scenario, row.names = 1)
activityList.baseline <- setNames(split(activityDF.baseline,seq(nrow(activityDF.baseline))),rownames(activityDF.baseline))
activityList.scenario <- setNames(split(activityDF.scenario,seq(nrow(activityDF.scenario))),rownames(activityDF.scenario))
return(data.frame(low = mapply( CRA.DF(leisure = "low"), activityList.baseline, activityList.scenario), medium = mapply( CRA.DF(leisure = "medium"), activityList.baseline, activityList.scenario), high = mapply( CRA.DF(leisure = "high"), activityList.baseline, activityList.scenario)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.