R/prostata.R
In mclements/prostata: Prostate Cancer Microsimulation for Sweden
Defines functions
nn.fhcrc
nn
lines.fhcrc
plot.fhcrc
predict.fhcrc
print.fhcrc
ICER.fhcrc
print.summary.fhcrc
summary.fhcrc
callFhcrc
.xtabsDFCheck
Documented in
callFhcrc
ICER.fhcrc
lines.fhcrc
nn
plot.fhcrc
predict.fhcrc
print.fhcrc
print.summary.fhcrc
summary.fhcrc
#' prostata
#'
#' @section Introduction:
#'
#' The packages provides implementations of discrete event simulation in both R
#' and C++.
#'
#' @section Background:
#'
#' Prostata is a natural history model of prostate cancer model which extends
#' the model developed by Ruth Etzioni and colleagues at the
#' [[http://www.fredhutch.org][Fred Hutchinson Cancer Research Center
#' (FHCRC)]]. This is a well validated prostate cancer natural history model
#' developed within the
#' [[http://cisnet.cancer.gov/prostate/profiles.html][Cancer Intervention and
#' Surveillance Modeling Network (CISNET)]]. We here provide our extensions of
#' the model with extended states and finer calibration. The model is provided
#' as an R package, depending on our
#' [[https://github.com/mclements/microsimulation][microsimulation]] frame work.
#' We specifically developed this package for modelling the cost-effectiveness
#' of prostate cancer screening, where many (e.g. 10^7) men are followed from
#' birth to death.
#'
#' @section Running the simulation:
#'
#' There are a number of available testing scenarios. They determine testing
#' frequencies and re-testing intervals over calendar period and ages.
#'
#' @docType package
#' @name prostata-package
#' @aliases prostata
#' @author Mark Clements \email{mark.clements@ki.se}
#' @references \url{https://github.com/mclements/prostata}
#' @references \url{https://github.com/mclements/microsimulation}
#' @seealso \code{\link{Rcpp}}
#' @useDynLib prostata, .registration=TRUE
#' @import microsimulation
#' @importFrom Rcpp evalCpp compileAttributes
#' @importFrom utils packageName
NULL
if(.Platform$OS.type == "unix") {
pkg <- packageName()
## http://r.789695.n4.nabble.com/How-to-construct-a-valid-seed-for-l-Ecuyer-s-method-with-given-Random-seed-td4656340.html
#' @param seed Random seed
#' @keywords internal
#' @importFrom microsimulation set.user.Random.seed
set.user.Random.seed <- function (seed)
microsimulation::set.user.Random.seed(seed,pkg)
#' @keywords internal
#' @importFrom microsimulation next.user.Random.substream
next.user.Random.substream <- function ()
microsimulation::next.user.Random.substream(pkg)
#' @keywords internal
#' @importFrom microsimulation user.Random.seed
user.Random.seed <- function()
microsimulation::user.Random.seed(pkg)
}
##
#' @title Initial values for the FHCRC model
#' @description Initial values for the FHCRC model
#' @format A list
#' \describe{
#'}
#' @details Currently, see the R documentation
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname FhcrcParameters
#' @export
"FhcrcParameters"
FhcrcParameters <- list(
Andreas = TRUE, # version for Andreas's CEA paper
revised_natural_history=TRUE,
## panel=FALSE,
grade.clinical.rate.high=0.3042454700,
startFullUptake = 1932.0,
fullUptakePortion = 0.9,
yearlyUptakeIncrease = 0.03,
startUptakeMixture = 1945.0,
endUptakeMixture = 1960.0,
uptakeStartAge=35.0,
screeningIntroduced = 1995.0,
shapeA = 3.8,
scaleA = 15.0,
shapeT = 2.16,
scaleT = 11.7,
tau2 = 0.0829, # log PSA measurement variance - normal
susceptible = 1.0, # portion susceptible
g0=0.0005, # onset parameter
g3p = exp(-6.73546100), # T3+ parameter
gm = exp(-7.46200003), # metastatic parameter
gc=0.0015, # clinical diagnosis parameter
thetac=19.1334, # clinical diagnosis parameter after metastatic
mubeta0=-1.609, # mean of beta0, where beta0 is the log PSA intercept at age 35 years
sebeta0=0.2384, # SE of beta0
mubeta1=0.04463, # mean of beta1, where beta1 is the log PSA slope prior to cancer onset
sebeta1=0.0430, # SE of beta1
mubeta2=c(0.0397,0.1678, 0.0), # base::grade: the gleason specific change in slope of log PSA after cancer onset
sebeta2=c(0.0913,0.3968, 0.0), # base::grade: variance of beta2
rev_mubeta2=c(0.051, 0.129, 0.1678), # ext::grade: same as above for extended gleason grade (6-, 7, 8+)
rev_sebeta2=c(0.064, 0.087, 0.3968), # ext::grade
alpha7 = log(0.2), # log of the proportion gleason 7 at age 35
beta7 = 0.06840404, # slope of log proportion of gleason 7
alpha8 = log(0.002), # log of the proportion of gleason 8+ at age 35
beta8 = 0.19894431, # slope of log proportion of gleason 8+
RR_T3plus=2.0, # prostate cancer mortality rate ratio comparing T3+ with T1-T2, add lit reference
## mubeta2.scale=1.0, # cf. 2.1
## beta.rho=0.62,
c_txlt_interaction = 1.0, # treatment lead time interaction, default to no effect
c_baseline_specific = 1.0, # scale factor on survival, default to no effect
c_benefit_value0 = 10, # (value -> reduction): 0.04 -> 10%; 0.18 -> 20%; 10 -> 28%, generalised stage-shift parameter
sxbenefit = 1.0, # hazard rate ratio for screening benefit, defaults to no effect
c_benefit_type = 0, # 0=stage-shift (=> c_benefit_value0=10), 1=lead-time based (=> c_benefit_value1=0.1)
c_benefit_value1 = 0.1, # approximate increase in cure for each year of lead-time, used for the lead-time based screening effect
RP_mortHR = 0.56, # mortality hazard ratio for surgery over watchful-waiting from SPCG-4 Bill-Axelson 2014
screeningParticipation = 0.75, # probability of actually having the first PSA test
rescreeningParticipation = 0.95, # probability of actually having the re-screening PSA tests
biopsyCompliance = 0.856, # Formal biopsy compliance, Schroder ERSPC 2014
biopsySensitivityTimeProportionT1T2 = 0.5314886, # time portion when T1-T2 cancers are sensitivity to biopsies (expit from calibration). The remaining part, starting at onset, is not detectable.
studyParticipation = 50.0/260.0, # observed fraction of population who participated in STHLM3 study
nLifeHistories = 10L, screen = 0L, ## integers
psaThreshold = 3.0,
psaThresholdBiopsyFollowUp = 4.0, # revised PSA threshold for negative follow-up screen
biomarker_model = 1, # biomarker_model = 0 random, biomarker_model = 1 psa/risk based correction of FP
PSA_FP_threshold_nCa=4.21, # reduce FP in no cancers with PSA threshold
PSA_FP_threshold_GG6=3.76, # reduce FP in GG 6 with PSA threshold
PSA_FP_threshold_GG7plus=3, # reduce FP in GG >= 7 with PSA threshold
panelReflexThreshold = 1.0,
## Natural history calibration
rTPF=1.0,
rFPF=0.6,
c_low_grade_slope=-0.006,
stockholmTreatment = TRUE,
discountRate.effectiveness = 0.03,
discountRate.costs = 0.03,
full_report = 1.0,
formal_costs = 0.0,
formal_compliance = 0.0,
start_screening = 50.0, # start of organised screening
stop_screening = 70.0, # end of organised screening
screening_interval = 2.0, # screening interval for regular_screening and introduced_screening*
screening_interval1 = 2.0, # screening interval for grs_stratified_age for first age group
screening_interval2 = 4.0, # screening interval for grs_stratified_age for second age group
screening_interval_split = 60.0, # age to split screening interval for grs_stratified_age
introduction_year = 2015.0, # year to start organised screening for introduced_screening* & cancel the opportunistic testing under stopped_screening
risk_psa_threshold=1, # PSA threshold for risk-stratified screening
risk_lower_interval=8, # re-screening interval for lower risk
risk_upper_interval=4, # re-screening interval for higher risk
mu0=c(0.00219, 0.000304, 5.2e-05, 0.000139, 0.000141, 3.6e-05, 7.3e-05,
0.000129, 3.8e-05, 0.000137, 6e-05, 8.1e-05, 6.1e-05, 0.00012,
0.000117, 0.000183, 0.000185, 0.000397, 0.000394, 0.000585, 0.000448,
0.000696, 0.000611, 0.000708, 0.000659, 0.000643, 0.000654, 0.000651,
0.000687, 0.000637, 0.00063, 0.000892, 0.000543, 0.00058, 0.00077,
0.000702, 0.000768, 0.000664, 0.000787, 0.00081, 0.000991, 9e-04,
0.000933, 0.001229, 0.001633, 0.001396, 0.001673, 0.001926, 0.002217,
0.002562, 0.002648, 0.002949, 0.002729, 0.003415, 0.003694, 0.004491,
0.00506, 0.004568, 0.006163, 0.006988, 0.006744, 0.00765, 0.007914,
0.009153, 0.010231, 0.011971, 0.013092, 0.013839, 0.015995, 0.017693,
0.018548, 0.020708, 0.022404, 0.02572, 0.028039, 0.031564, 0.038182,
0.042057, 0.047361, 0.05315, 0.058238, 0.062619, 0.074934, 0.089776,
0.099887, 0.112347, 0.125351, 0.143077, 0.153189, 0.179702, 0.198436,
0.240339, 0.256215, 0.275103, 0.314157, 0.345252, 0.359275, 0.41768,
0.430279, 0.463636, 0.491275, 0.549738, 0.354545, 0.553846, 0.461538,
0.782609),
hr_locoregional = transform(expand.grid(age=c(50,60,70), ext_grade=0:2,
psa10=0:1),
hr = c(0.2640235, 0.5107137,
1.0917129, 1.9397446,
1.6671078, 1.9372717,
1.7998275, 1.0379953,
1.7205986, 0.8606931,
0.8966157, 1.1941040,
3.6151038, 2.9625742 ,
2.1622938, 1.0185837,
0.9337168, 0.8849155)),
hr_metastatic = data.frame(age = c(50, 60, 70),
hr = c(0.893, 0.779, 0.698)),
biopsy_sensitivity = data.frame(Year = c(1987, 1988, 1989, 1990, 1991, 1992,
1993, 1994, 1995, 1996, 1997, 1998,
1999, 2000),
Sensitivity = c(0.502, 0.574, 0.642, 0.702,
0.749, 0.782, 0.803, 0.814,
0.821, 0.830, 0.844, 0.866,
0.894, 0.925)),
rand_biopsy_sensitivityG6 = 1.0,
neg_biopsy_to_psa=structure(list(age = c(50, 60, 70, 80),
meanlog = c(-0.314249298970759, -0.36332284250614,
-0.333753007363153, -0.175063351183678),
sdlog = c(0.996188493950768, 0.974965602213158,
1.0490009138328, 1.23085968935057)),
.Names = c("age", "meanlog", "sdlog"),
row.names = c(NA, -4L),
class = "data.frame"),
neg_biopsy_to_biopsy = structure(list(age = c(50, 60, 70, 80),
meanlog = c(3.05337156798731, 2.65573907207411,
1.78297474442106, 1.11581637027693),
sdlog = c(2.41326118058954, 2.3225660366395,
1.86615343381663, 2.04687048711335)),
.Names = c("age", "meanlog", "sdlog"),
row.names = c(NA, -4L), class = "data.frame"),
cure_m_CM_to_RP = structure(list(age = c(50, 60, 70, 80),
pnever = c(0.485026452157334, 0.605786246502961,
0.922871191840356, 0.999745963308522),
meanlog = c(1.14257205575071, 1.04021691430639,
0.519760988362341, 4.35555555555556),
sdlog = c(0.758280639791713, 0.875661155361679,
0.782118228501619, 0.255901426097166)),
.Names = c("age", "pnever", "meanlog", "sdlog"),
row.names = c(NA, -4L), class = "data.frame"),
cure_m_CM_to_RT =structure(list(age = c(50, 60, 70, 80),
pnever = c(0.60082830813107,
2.71246734638524e-06,
0.493326620739727,
0.739523343815296),
meanlog = c(3.14746283053335,
4.24925205843398,
2.0235406338416,
0.705091563838411),
sdlog = c(2.0453909799028,
2.26870437384556,
2.1703325154377,
1.81819516420271)),
.Names = c("age", "pnever", "meanlog", "sdlog"),
row.names = c(NA, -4L), class = "data.frame"),
currency_rate = 0.750/9.077, # PPP EU19@2017/Swe@2016 https://data.oecd.org/conversion/purchasing-power-parities-ppp.htm
## Should we add cost and dis-utilities for using test characteristics based on prostate volume 7% from supplement
## Swedish governmental report on organised PSA testing (p.22):
## https://www.socialstyrelsen.se/SiteCollectionDocuments/2018-2-13-halsoekonomisk-analys.pdf
## Based on the Swedish south region 2017:
## http://sodrasjukvardsregionen.se/avtal-priser/regionala-priser-och-ersattningar-foregaende-ar/
## S3M cost from Karolinska University Laboratory (KUL):
## https://www.karolinska.se/KUL/Alla-anvisningar/Anvisning/10245
cost_parameters = c("Invitation" = 7 # Invitation letter
+ 7, # Results letter
"Formal PSA" = 349 # test sampling, primary care
+ 45 # PSA analysis
+ 0 * 1539, # No GP primary care
"Formal panel" = 349 # test sampling, primary care
+ 45 # PSA analysis not included in panel price
+ 2300 # From KUL price list
+ 0 * 1539, # No GP for formal
"Opportunistic PSA" = 349 # test sampling, primary care
+ 45 # PSA analysis
+ 0.2 * 1539, # GP primary care
"Opportunistic panel" = 349 # test sampling, primary care
+ 45 # PSA analysis not included in panel price
+ 2300 # From KUL price list
+ 0.2 * 1539, # GP primary care
"Biopsy" = 4733, # Biopsy cost
"Assessment" = 1794, # Urology assessment
"Prostatectomy" = 80000 # Surgery
+ 100000 * 0.25 # Radiation therapy
+ 1794 * 2 # Urology visit
+ 1205 * 2, # Nurse visit
"Radiation therapy" = 100000 # Radiation therapy
+ 1794 * 2 # Urology visit
+ 1205 * 2, # Nurse visit
"Active surveillance - yearly" = 1794 # Urology visit
+ 45 * 2 # PSA analysis
+ 4733 * 0.5, # Biopsy
"Active surveillance - single MR" = 3090, # Used once for active surveillance
"Post-Tx follow-up - yearly" = 349 # PSA test sampling
+ 45 # PSA analysis,
+ 474, # Telefollow-up by urologist
"Cancer death" = 100160 * 3 # Care for spread disease
+ 68000 * 3, # Drugs for spread disease
"Polygenic risk stratification" = 250*12, # Callender et al (2021) with exchange rate of approximately 12
"Opportunistic DRE" = 349 # DRE procedure in primary care
+ 0.2 * 1539), # part of primary care visit
## Swedish governmental report on organised PSA testing (p.23):
## https://www.socialstyrelsen.se/SiteCollectionDocuments/2018-2-13-halsoekonomisk-analys.pdf
production = data.frame(ages = c(0, 54, 64, 74),
values=apply(
rbind(0.878, 0.756, 0.158, 0), # employment portion, full time
1,
function(empl) empl
* 30.7/40 # average working hours
* 32800*12 # average salary
* 1.485)), # including non-optional social fees
lost_production_years= c("Formal PSA"=2/24/365.25,
"Formal panel"=2/24/365.25,
"Opportunistic PSA"=2/24/365.25,
"Opportunistic panel"=2/24/365.25,
"Assessment"=2/24/365.25, # Urology assessment
"Biopsy"=2/24/365.25,
"Prostatectomy"=6/52,
"Radiation therapy"=8/52,
"Active surveillance - yearly"=
1 * 2/24/365.25 # Urology assessment
+ 2 * 2/24/365.25 # PSA tests
+ 0.5 * 2/24/365.25, # Biopsy
"Metastatic cancer"=6/12,
"Terminal illness" = 6/12), # NOTE: potentially add follow-up Tx production losses
## Heijnsdijk 2012
utility_estimates = c("Invitation" = 1,
"Formal PSA" = 0.99,
"Formal panel" = 0.99,
"Opportunistic PSA" = 0.99,
"Opportunistic panel" = 0.99,
"Biopsy" = 0.90,
"Cancer diagnosis" = 0.80,
"Prostatectomy part 1" = 0.67,
"Prostatectomy part 2" = 0.77,
"Radiation therapy part 1" = 0.73,
"Radiation therapy part 2" = 0.78,
"Active surveillance" = 0.97,
"Postrecovery period" = 0.95,
"Palliative therapy" = 0.60,
"Terminal illness" = 0.40,
"Death" = 0.00),
## Utility duration is given in years.
utility_duration = c("Invitation" = 0.0,
"Formal PSA" = 1/52,
"Formal panel" = 1/52,
"Opportunistic PSA" = 1/52,
"Opportunistic panel" = 1/52,
"Biopsy" = 3/52,
"Cancer diagnosis" = 1/12,
"Prostatectomy part 1" = 2/12,
"Prostatectomy part 2" = 10/12,
"Radiation therapy part 1" = 2/12,
"Radiation therapy part 2" = 10/12,
"Active surveillance" = 7,
"Postrecovery period" = 9,
"Palliative therapy" = 30/12,
"Terminal illness" = 6/12),
utility_truncate = TRUE, # should utilities be truncated at zero?
utility_scales = c("UtilityAdditive"=0,"UtilityMultiplicative"=1,"UtilityMinimum"=2), # encoding for the utility scales
utility_scale = as.double(1), # default scale = UtilityMultiplicative
includeEventHistories = TRUE,
includePSArecords = FALSE,
includeBxrecords = FALSE,
includeDiagnoses = FALSE,
MRI_screen = FALSE, # defines whether MRI pathway is used for screen-positive patients
MRI_clinical = FALSE, # defines whether MRI pathway is used for clinically detected cancers
MRInegSBx=FALSE, # No SBx for MRI- (by default)
rescreenDoubleNeg = FALSE, # defines whether MRI-/Bx- men should move to rescreening (cf. repeat PSA and Bx within ~12 months)
indiv_reports = FALSE, # should the cost and standard reports include individual values?
startReportAge = 0.0, # Age to start reporting (currently only means_utilities and means_costs)
weibull_onset = FALSE, # flag to use Weibull-distributed onset
weibull_onset_shape = 2, # shape; values (0,Inf)
weibull_onset_scale= 40, # scale; values (0,Inf)
frailty = FALSE, # assume a frailty distribution on the onset distribution?
grs_variance = 0.68, # Callender et al (2019)
other_variance = 1.14, # total variance = 1.82 from Kicinski et al (2011; https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0027130)
grs_risk_threshold = 0.04 # ten-year risk threshold for GRS-based risk stratified screening
)
IHE <- list(prtx=data.frame(Age=50.0,DxY=1973.0,G=1:2,CM=0.6,RP=0.26,RT=0.14)) ## assumed constant across ages and periods
ParameterNV <- FhcrcParameters[sapply(FhcrcParameters,class)=="numeric" & sapply(FhcrcParameters,length)==1]
## ParameterIV <- FhcrcParameters[sapply(FhcrcParameters,class)=="integer" & sapply(FhcrcParameters,length)==1]
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 15 rows and 3 variables:
#' \describe{
#' \item{\code{psa}}{double COLUMN_DESCRIPTION}
#' \item{\code{age}}{double COLUMN_DESCRIPTION}
#' \item{\code{compliance}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"swedenOpportunisticBiopsyCompliance"
swedenOpportunisticBiopsyCompliance <- data.frame(
psa = c(3, 5, 10, 3, 5, 10, 3, 5, 10, 3, 5, 10, 3, 5, 10),
age = c(40, 40, 40, 50, 50, 50, 60, 60, 60, 70, 70, 70, 80, 80, 80),
compliance = c(0.3764045, 0.5680751, 0.7727273, 0.3110770, 0.5726548, 0.7537372, 0.2385155, 0.4814588, 0.6929770, 0.1754264, 0.3685056, 0.5602030, 0.1629213, 0.2697368, 0.5010052))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 15 rows and 3 variables:
#' \describe{
#' \item{\code{psa}}{double COLUMN_DESCRIPTION}
#' \item{\code{age}}{double COLUMN_DESCRIPTION}
#' \item{\code{compliance}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"swedenFormalBiopsyCompliance"
swedenFormalBiopsyCompliance <- cbind(expand.grid(psa=c(3,5,10),age=seq(40,80,10)),
compliance=FhcrcParameters$biopsyCompliance)
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 24 rows and 6 variables:
#' \describe{
#' \item{\code{DxY}}{double COLUMN_DESCRIPTION}
#' \item{\code{Age}}{double COLUMN_DESCRIPTION}
#' \item{\code{G}}{integer COLUMN_DESCRIPTION}
#' \item{\code{CM}}{double COLUMN_DESCRIPTION}
#' \item{\code{RP}}{double COLUMN_DESCRIPTION}
#' \item{\code{RT}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"stockholmTreatment"
stockholmTreatment <-
data.frame(DxY=2008,
Age=c(50,50,50,55,55,55,60,60,60,65,65,65,70,70,70,75,75,75,80,80,80,85,85,85),
G=as.integer(c(6,7,8,6,7,8,6,7,8,6,7,8,6,7,8,6,7,8,6,7,8,6,7,8)-6),
CM=c(0.231023,0.044872,0.25,0.333333,0.09542,0.328358,0.409439,0.12825,0.348101,0.479167,0.178182,0.401639,0.689013,0.359143,0.56087,0.876543,0.744444,0.809524,1,0.970711,0.952096,0.9375,1,1),
RP=c(0.700623,0.815839,0.5,0.553592,0.740111,0.326226,0.49981,0.6866,0.291437,0.409879,0.552483,0.279235,0.210949,0.318374,0.179363,0.041152,0.058652,0.049689,0,0.009763,0.023952,0.0625,0,0),
RT=c(0.068354,0.13929,0.25,0.113074,0.164469,0.345416,0.090751,0.185151,0.360462,0.110954,0.269335,0.319126,0.100038,0.322482,0.259767,0.082305,0.196903,0.140787,0,0.019526,0.023952,0,0,0))
secularTrendTreatment2008OR <-
data.frame(year=as.numeric(1988:2009),
OR=c(0.194409,0.155874,0.130704,0.154104,0.148039,0.237752,0.299491,0.330934,0.407047,0.377968,0.443079,0.551047,0.686101,0.700657,0.846134,0.93956,1.03026,1.191726,1.067736,1.012136,1,0.911471))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 52 rows and 5 variables:
#' \describe{
#' \item{\code{age5}}{double COLUMN_DESCRIPTION}
#' \item{\code{total_cat}}{double COLUMN_DESCRIPTION}
#' \item{\code{shape}}{double COLUMN_DESCRIPTION}
#' \item{\code{cure}}{double COLUMN_DESCRIPTION}
#' \item{\code{scale}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"rescreening"
rescreening <- data.frame(age5 = c(30, 30, 30, 30, 35, 35, 35, 35, 40, 40,
40, 40, 45, 45, 45, 45, 50, 50, 50, 50, 55, 55, 55, 55, 60, 60,
60, 60, 65, 65, 65, 65, 70, 70, 70, 70, 75, 75, 75, 75, 80, 80,
80, 80, 85, 85, 85, 85, 90, 90, 90, 90), total_cat = c(0, 1,
3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0,
1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10,
0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10),
shape = c(1.04293463521838,
0.825156738737092, 0.899003822110993, 0.589358680965725, 1.05948319345651,
0.895913335495644, 0.84619233996788, 0.636706369588122, 1.16028493785731,
0.97919993225436, 0.68848363112076, 0.667566641635962, 1.25727757690268,
1.11021543939268, 0.717059585208267, 0.724574127291869, 1.3247305288869,
1.18510576695807, 0.858171614262325, 0.749998577524088, 1.3340371479645,
1.20792880820659, 0.920717960779785, 0.852052206426261, 1.32268361303857,
1.24450960892428, 0.988946105418831, 0.919376727690091, 1.31156083562523,
1.26972967056357, 1.03454490825271, 0.954911145157684, 1.27126976718353,
1.26045834296745, 1.06684366805728, 0.968088650561122, 1.18494141305856,
1.20593215449012, 1.05998296831204, 0.989315997334925, 1.12395322065035,
1.13733607943079, 1.04483735235382, 0.988385269358649, 1.1220716456468,
1.09128664878198, 1.01326552242001, 0.948730052972523, 1.08315738475919,
1.09407460180658, 1.05197207574662, 0.978397483932468),
cure = c(0.585284664126963,
0.545117838495959, 0.402804166722055, 0.14283955595853, 0.398430641565923,
0.327274451509188, 0.148468826128857, 0.00333583482942569, 0.257384245965907,
0.208165119313289, 0.157385208407254, 0.0696401481749125, 0.166635358789789,
0.133819199845029, 0.0610524867384188, 0.0535645284203467, 0.134003946557746,
0.114600739520109, 0.0537234127141967, 0.0308512461181499, 0.104553332825899,
0.0862550079726627, 0.0399447271614142, 0.0259518317972345, 0.10432761596883,
0.0823122881511086, 0.0355479853523502, 0.0190476453614591, 0.103738700348739,
0.0816458941897609, 0.0412696810847696, 0.0319927106467962, 0.108815314898247,
0.0867189417281315, 0.0472560235806836, 0.0300347697835532, 0.124252895869383,
0.111313062025094, 0.0675699701502124, 0.0373036736422192, 0.15590486866906,
0.153015363026187, 0.0974421665330981, 0.0563274705516883, 0.228204960789236,
0.2051853623161, 0.169506663277659, 0.0887344403643672, 0.334252927026757,
0.295865002768107, 0.24644310712191, 0.151125704749819),
scale = c(3.33440594470802,
3.56052439703143, 0.247924244305486, 0.150095610340916, 4.05027295845979,
3.78800058764034, 0.35102275649495, 0.486681740716857, 3.62361340129802,
3.48003996256557, 0.661105185961319, 0.247063080430436, 2.80272048406229,
2.6202749242671, 0.764421647692264, 0.268131101357167, 2.29273349177322,
2.11664063949949, 0.749125964202465, 0.417294230771058, 2.06795900869592,
1.80352742984593, 0.78399749774511, 0.451625737966194, 1.82650979495336,
1.62873839529095, 0.806784448546511, 0.478644431229979, 1.61927028533258,
1.46189762327418, 0.810126410765397, 0.500531738601558, 1.53243134315106,
1.42674067036488, 0.869134221688119, 0.569591312385404, 1.44388696689722,
1.41174931932889, 0.938154963786188, 0.631891327069348, 1.47052388447475,
1.4543783566239, 1.01084363142901, 0.644940081254986, 1.29365253055963,
1.4079170674224, 1.03720449704243, 0.643101192871478, 1.08541958643012,
1.29524623033074, 1.02176143186057, 0.673175882772333))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 136 rows and 2 variables:
#' \describe{
#' \item{\code{cohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{pop}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"pop1"
pop1 <- data.frame(cohort=as.double(2035:1900),
pop=c(rep(17239,32), 16854, 16085, 15504, 15604, 16381, 16705,
16762, 16853, 15487, 14623, 14066, 13568, 13361, 13161, 13234,
13088, 12472, 12142, 12062, 12078, 11426, 12027, 11963, 12435,
12955, 13013, 13125, 13065, 12249, 11103, 9637, 9009, 8828,
8350, 7677, 7444, 7175, 6582, 6573, 6691, 6651, 6641, 6268,
6691, 6511, 6857, 7304, 7308, 7859, 7277, 8323, 8561, 7173,
6942, 7128, 6819, 5037, 6798, rep(6567,46)))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 4 rows and 4 variables:
#' \describe{
#' \item{\code{cohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{pop}}{double COLUMN_DESCRIPTION}
#' \item{\code{screened}}{double Number screened at baseline}
#' \item{\code{pscreened}}{double Probability of being screened at baseline}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname cap
#' @export
"cap_study"
cap_study <- data.frame(cohort=c(1955,1950,1945,1940),
pop=c(55229, 55077,44057,35023),
screened=c(17869,19448,15783,11335))
cap_study <- transform(cap_study, pscreened=screened/pop)
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 4 rows and 4 variables:
#' \describe{
#' \item{\code{cohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{pop}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname cap
#' @export
"cap_control"
cap_control <- data.frame(cohort=c(1955,1950,1945,1940),
pop=c(63423,63285,51507,41224))
#' @title background_utilities
#' @description Swedish background utilities from Burstrom and Rehnberg (2006)
#' @format A data frame with 14 rows and 3 variables:
#' \describe{
#' \item{\code{lower}}{double lower age}
#' \item{\code{lower}}{double upper age}
#' \item{\code{utility}}{double background utility value}
#'}
#' @details http://libris.kb.se/bib/10708053?vw=short
#' @examples
#' \dontrun{
#' if(interactive()){
#' background_utilities
#' }
#' }
#' @rdname background_utilities
#' @export
"background_utilities"
background_utilities <-
data.frame(lower=c(0, 18, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80),
upper=c(18, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 1.0e55),
utility=c(1, 0.89, 0.89, 0.88, 0.87, 0.84, 0.84, 0.83, 0.83, 0.82, 0.83, 0.81,
0.79, 0.74))
#' @title Population for the STHLM3-MRI MRI arm
#' @description Birth cohorts for the study, assuming 2019.5 is the middle of the study
#' @format A data frame with 30 rows and 2 variables:
#' \describe{
#' \item{\code{cohort}}{birt cohort}
#' \item{\code{pop}}{population count}
#'}
#' @examples
#' \dontrun{
#' if(interactive()){
#' background_utilities
#' }
#' }
#' @rdname STHLM3
#' @export
"sthlm3_mri_arm"
sthlm3_mri_arm <- data.frame(cohort = 2019 - 45:74,
pop = c(1, 4, 6, 9, 10, 11, 21, 22, 20, 28, 23,
34, 44, 37, 48, 58, 59, 60, 71, 60, 71, 79, 76, 71, 79, 91, 95,
77, 83, 54))
#' @title List tables for the simulation
#' @description DATASET_DESCRIPTION
#' \strong{all_cause_mortality}
#' @format A data frame with 12120 rows and 4 variables:
#' \describe{
#' \item{\code{BirthCohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Age}}{integer COLUMN_DESCRIPTION}
#' \item{\code{AnnualMortality}}{double COLUMN_DESCRIPTION}
#' \item{\code{Survival}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{biopsy_frequency}
#' @format A data frame with 3 rows and 7 variables:
#' \describe{
#' \item{\code{PSA.beg}}{integer COLUMN_DESCRIPTION}
#' \item{\code{PSA.end}}{integer COLUMN_DESCRIPTION}
#' \item{\code{X55.59}}{double COLUMN_DESCRIPTION}
#' \item{\code{X60.64}}{double COLUMN_DESCRIPTION}
#' \item{\code{X65.69}}{double COLUMN_DESCRIPTION}
#' \item{\code{X70.74}}{double COLUMN_DESCRIPTION}
#' \item{\code{X75.79}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{biopsy_sensitivity}
#' @format A data frame with 14 rows and 2 variables:
#' \describe{
#' \item{\code{Year}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Sensitivity}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{neg_biopsy_to_psa}
#' @format A data frame with 4 rows and 3 variables. Describing the
#' time to the next PSA test following a negative biopsy. Modelled
#' as a competing risk with a biopsy. Informed by the Stockholm
#' PSA and Biopsy Register (SPBR).:
#' \describe{
#' \item{\code{age}}{integer with age groups}
#' \item{\code{meanlog}}{double for the mean of a log-normal time to
#' a PSA test following a negative biopsy.}
#' \item{\code{sdlog}}{double for the standard deviation of a
#' log-normal.}
#' }
#' \strong{neg_biopsy_to_biopsy}
#' @format A data frame with 4 rows and 3 variables. Describing the
#' time to the next biopsy following a negative biopsy. Modelled
#' as a competing risk with a PSA test. Informed by the Stockholm
#' PSA and Biopsy Register (SPBR).:
#' \describe{
#' \item{\code{age}}{integer with age groups}
#' \item{\code{meanlog}}{double describing the mean of a log-normal}
#' \item{\code{sdlog}}{double describing the standard deviation of a
#' log-normal}
#' }
#' \strong{cure_m_CM_to_RP}
#' @format A data frame with 4 rows and 4 variables. Cure model giving
#' the probability of not having a radical prostatectomy following
#' assignment to conservative management (active surveillance &
#' watchfull waiting). Log-normal distribution of the time to
#' radical prostatectomy for those that do have that. A future
#' extension would be to model active surveillance and watchfull
#' waiting separately:
#' \describe{
#' \item{\code{age}}{double with age groups}
#' \item{\code{pnever}}{double probability of not having a radical
#' prostatectomy}
#' \item{\code{meanlog}}{double mean of log-normal time to radical
#' prostatectomy}
#' \item{\code{sdlog}}{double standard deviation of log-normal time to
#' radical prostatectomy}
#' }
#' \strong{cure_m_CM_to_RT}
#' @format A data frame with 4 rows and 4 variables. Cure model giving
#' the probability of not having a radiation therapy following
#' assignment to conservative management (active surveillance &
#' watchfull waiting). Log-normal distribution of the time to
#' radiation therapy for those that do have that. A future
#' extension would be to model active surveillance and watchfull
#' waiting separately:
#' \describe{
#' \item{\code{age}}{double with age groups}
#' \item{\code{pnever}}{double probability of not having a radical
#' prostatectomy}
#' \item{\code{meanlog}}{double mean of log-normal time to
#' radiation therapy}
#' \item{\code{sdlog}}{double standard deviation of log-normal time to
#' radiation therapy}
#' }
#' @format A data frame with 4 rows and 3 variables. Describing the
#' time to the next biopsy following a negative biopsy. Modelled
#' as a competing risk with a PSA test. Informed by the Stockholm
#' PSA and Biopsy Register (SPBR).:
#' \describe{
#' \item{\code{age}}{integer with age groups}
#' \item{\code{meanlog}}{double describing the mean of a log-normal}
#' \item{\code{sdlog}}{double describing the standard deviation of a
#' log-normal}
#' }
#' \strong{dre}
#' @format A data frame with 4 rows and 4 variables:
#' \describe{
#' \item{\code{psa.low}}{integer COLUMN_DESCRIPTION}
#' \item{\code{psa.high}}{double COLUMN_DESCRIPTION}
#' \item{\code{sensitivity}}{double COLUMN_DESCRIPTION}
#' \item{\code{specificity}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{prob_grade7}
#' @format A data frame with 51 rows and 2 variables:
#' \describe{
#' \item{\code{slope}}{double COLUMN_DESCRIPTION}
#' \item{\code{Pr.grade.7.}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{pradt}
#' @format A data frame with 8208 rows and 6 variables:
#' \describe{
#' \item{\code{Tx}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Age}}{integer COLUMN_DESCRIPTION}
#' \item{\code{DxY}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{NoADT}}{double COLUMN_DESCRIPTION}
#' \item{\code{ADT}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{prtx}
#' @format A data frame with 2664 rows and 6 variables:
#' \describe{
#' \item{\code{Age}}{integer COLUMN_DESCRIPTION}
#' \item{\code{DxY}}{integer COLUMN_DESCRIPTION}
#' \item{\code{G}}{integer COLUMN_DESCRIPTION}
#' \item{\code{CM}}{double COLUMN_DESCRIPTION}
#' \item{\code{RP}}{double COLUMN_DESCRIPTION}
#' \item{\code{RT}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{survival_dist}
#' @format A data frame with 42 rows and 3 variables:
#' \describe{
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Time}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Survival}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{survival_local}
#' @format A data frame with 294 rows and 5 variables:
#' \describe{
#' \item{\code{AgeLow}}{integer COLUMN_DESCRIPTION}
#' \item{\code{AgeHigh}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Time}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Survival}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{prob_earlystage}
#' @format A data frame with 84 rows and 6 variables:
#' \describe{
#' \item{\code{BeginAge}}{integer COLUMN_DESCRIPTION}
#' \item{\code{EndAge}}{integer COLUMN_DESCRIPTION}
#' \item{\code{BeginPSA}}{integer COLUMN_DESCRIPTION}
#' \item{\code{EndPSA}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Prob}}{double COLUMN_DESCRIPTION}
#'}
"fhcrcData"
fhcrcData$biopsyOpportunisticComplianceTable <- swedenOpportunisticBiopsyCompliance
fhcrcData$biopsyFormalComplianceTable <- swedenFormalBiopsyCompliance
fhcrcData$secularTrendTreatment2008OR <- secularTrendTreatment2008OR
## https://www.socialstyrelsen.se/Lists/Artikelkatalog/Attachments/20008/2015-12-26.pdf
fhcrcData$background_utilities <- background_utilities
## fhcrcData$uk_screen_uptake <-
## data.frame(age = as.double(44:76),
## H = c(0, 0.00279457395690957, 0.0111782958276383,
## 0.0251511656121861, 0.0447131833105531, 0.0698643489227393, 0.098027986598169,
## 0.126627420486267, 0.155662650587033, 0.185133676900467, 0.215040499426569,
## 0.247663838926632, 0.285284416161949, 0.32790223113252, 0.375517283838344,
## 0.428129574279422, 0.483488769981715, 0.539344538471183, 0.595696879747827,
## 0.652545793811647, 0.709891280662642, 0.769265269782334, 0.832199690652243,
## 0.89869454327237, 0.968749827642714, 1.04236554376328, 1.11776147575895,
## 1.19315740775462, 1.26855333975029, 1.34394927174596, 1.41934520374163,
## 1.42934520374163, 1.43934520374163))
#' @title UK rescreening rates
#' @description DATASET_DESCRIPTION
#' @format A data frame with 48 rows and 5 variables:
#' \describe{
#' \item{\code{age5}}{double COLUMN_DESCRIPTION}
#' \item{\code{total_cat}}{double COLUMN_DESCRIPTION}
#' \item{\code{shape}}{double COLUMN_DESCRIPTION}
#' \item{\code{cure}}{double COLUMN_DESCRIPTION}
#' \item{\code{scale}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname uk_rescreening
#' @export
"uk_rescreening"
uk_rescreening <- data.frame(age5 = c(30, 30, 30, 30, 30, 30, 50, 50, 50, 50,
50, 50, 55, 55, 55, 55, 55, 55, 60, 60, 60, 60, 60, 60, 65, 65,
65, 65, 65, 65, 80, 80, 80, 80, 80, 80, 85, 85, 85, 85, 85, 85,
90, 90, 90, 90, 90, 90),
total_cat = c(0, 3, 4, 6, 10, 20, 0,
3, 4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3,
4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3, 4,
6, 10, 20),
shape = c(1.18510576695807, 0.858171614262325, 0.858171614262325,
0.858171614262325, 0.749998577524088, 0.749998577524088, 1.18510576695807,
0.858171614262325, 0.858171614262325, 0.858171614262325, 0.749998577524088,
0.749998577524088, 1.20792880820659, 0.920717960779785, 0.920717960779785,
0.920717960779785, 0.852052206426261, 0.852052206426261, 1.24450960892428,
0.988946105418831, 0.988946105418831, 0.988946105418831, 0.919376727690091,
0.919376727690091, 1.26972967056357, 1.03454490825271, 1.03454490825271,
1.03454490825271, 0.954911145157684, 0.954911145157684, 1.13733607943079,
1.04483735235382, 1.04483735235382, 1.04483735235382, 0.988385269358649,
0.988385269358649, 1.09128664878198, 1.01326552242001, 1.01326552242001,
1.01326552242001, 0.948730052972523, 0.948730052972523, 1.09407460180658,
1.05197207574662, 1.05197207574662, 1.05197207574662, 0.978397483932468,
0.978397483932468),
cure = c(0.114600739520109, 0.0537234127141967,
0.0537234127141967, 0.0537234127141967, 0.0308512461181499, 0.0308512461181499,
0.114600739520109, 0.0537234127141967, 0.0537234127141967, 0.0537234127141967,
0.0308512461181499, 0.0308512461181499, 0.0862550079726627, 0.0399447271614142,
0.0399447271614142, 0.0399447271614142, 0.0259518317972345, 0.0259518317972345,
0.0823122881511086, 0.0355479853523502, 0.0355479853523502, 0.0355479853523502,
0.0190476453614591, 0.0190476453614591, 0.0816458941897609, 0.0412696810847696,
0.0412696810847696, 0.0412696810847696, 0.0319927106467962, 0.0319927106467962,
0.153015363026187, 0.0974421665330981, 0.0974421665330981, 0.0974421665330981,
0.0563274705516883, 0.0563274705516883, 0.2051853623161, 0.169506663277659,
0.169506663277659, 0.169506663277659, 0.0887344403643672, 0.0887344403643672,
0.295865002768107, 0.24644310712191, 0.24644310712191, 0.24644310712191,
0.151125704749819, 0.151125704749819),
scale = c(6.58183283712968,
3.56486060274202, 0.96113775302889, 0.784828252112451, 0.937020402040802,
1.59551751310995, 6.58183283712968, 3.56486060274202, 0.96113775302889,
0.784828252112451, 0.937020402040802, 1.59551751310995, 5.47461832642227,
4.85185826844681, 1.19197420365362, 0.960420044179229, 1.12429005910723,
1.6320163836042, 4.5256676489079, 4.92695557504941, 1.37335505106699,
1.08828035047345, 1.16345156553308, 1.88292622499511, 3.88781952822644,
4.07180517697073, 1.814341631623, 1.10397958554216, 1.27491455185655,
2.0730535296673, 1.4543783566239, 1.01084363142901, 1.01084363142901,
1.01084363142901, 0.644940081254986, 0.644940081254986, 1.4079170674224,
1.03720449704243, 1.03720449704243, 1.03720449704243, 0.643101192871478,
0.643101192871478, 1.29524623033074, 1.02176143186057, 1.02176143186057,
1.02176143186057, 0.673175882772333, 0.673175882772333))
#' @title ageStandards
#' @description Age standardisation weights
#' @format A data frame with 18 rows and 3 variables:
#' \describe{
#' \item{\code{Age}}{integer Left-hand age for five-year age groups to 85+}
#' \item{\code{Sweden2000}}{double Sweden 2000 standard weights}
#' \item{\code{World}}{double World standard weights}
#'}
#' @details Age groups
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname ageStandards
#' @export
"ageStandards"
ageStandards <- data.frame(Age = cut(seq(0, 85, 5),
breaks = c(seq(0,85,5), Inf),
right = FALSE),
Sweden2000 = c(5.3, 6.9, 6.4, 5.7, 5.9, 6.7, 7.2, 6.9,
6.6, 6.6, 7.6, 6.3, 4.9, 4.3, 4.1, 3.9,
2.6, 2.3),
World = c(12.0, 10.0, 9.0, 9.0, 8.0, 8.0, 6.0, 6.0, 6.0,
6.0, 5.0, 4.0, 4.0, 3.0, 2.0, 1.0, 0.5, 0.5))
#' @title dre
#' @description Sensitivity and specificity for digital rectal examination (DRE)
#' @format A data frame with 6 rows and 3 variables:
#' \describe{
#' \item{\code{psa_low}}{double Left-hand interval for PSA}
#' \item{\code{se}}{double sensitivity}
#' \item{\code{sp}}{double specificity}
#'}
#' @details See https://doi.org/10.1093/jnci/90.23.1817, table 2.
#' @rdname dre
#' @export
dre = data.frame(psa_low=c(0:4,10),
psa_high=c(1:4,10,1e6),
dre=c(1702,3305,1314,734,1095,217),
appa=c(19,119,97,89,254,128), # predicted!!
bx=c(109,296,128,106,249,82),
pc=c(4,29,14,35,115,67))
dre = transform(dre, sensitivity=pc/appa, specificity=(dre-appa-(bx-pc))/(dre-appa))
.xtabsDFCheck <- function(df, dim_cols) {
stopifnot(is.data.frame(df))
stopifnot(all(dim_cols %in% names(df)))
stopifnot(prod(sapply(dim_cols, function(name) length(unique(df[[name]])))) == nrow(df))
}
## fix pradt: repeated values in 2004
fhcrcData$pradt <- unique(fhcrcData$pradt)
#' @title Define and run simulation
#' @description Function to specify and run the microsimulation. A large number
#' of simulated men, \code{n}, will cause the simulation to take longer time
#' but will reduce the stochastic variation particularly for rare events.
#' @param n Integer number of men to simulate. Default: 10
#' @param screen String with one of the following simulated screening scenarios:
#' \describe{
#' \item{\code{noScreening}}{no screening test, only diagnosis from symptoms}
#' \item{\code{randomScreen50to70}}{TBA}
#' \item{\code{twoYearlyScreen50to70}}{two-yearly screening from age 50 to 70}
#' \item{\code{fourYearlyScreen50to70}}{four-yearly screening from age 50 to 70}
#' \item{\code{screen50}}{one screen at age 50}
#' \item{\code{screen60}}{one screen at age 60}
#' \item{\code{screen70}}{one screen at age 70}
#' \item{\code{screenUptake}}{current testing pattern in Sweden}
#' \item{\code{stockholm3_goteborg}}{TBA}
#' \item{\code{stockholm3_risk_stratified}}{TBA}
#' \item{\code{goteborg}}{risk stratified re-screening 2+4 from age 50 to 70}
#' \item{\code{risk_stratified}}{risk stratified re-screening 4+8 from age 50 to 70}
#' \item{\code{mixed_screening}}{risk stratified re-screening 2+4 from age 50 to 70 & opportunistic testing for other ages}
#' \item{\code{regular_screen}}{TBA}
#' \item{\code{single_screen}}{TBA}
#' \item{\code{introduced_screening_only}}{TBA}
#' \item{\code{introduced_screening_preference}}{TBA}
#' \item{\code{introduced_screening}}{TBA}
#' \item{\code{stopped_screening}}{TBA}
#' } . Default: 'noScreening'
#'
#' @param nLifeHistories Integer with number of men for all events should be
#' recorded, Default: 10
#' @param seed Integer random number seed, Default: 12345
#' @param panel Boolean to use the Stockholm3 biomarker panel test
#' characteristics, otherwise PSA is used, Default: FALSE
#' @param flatPop Boolean to set all birth cohorts of equal size, Default: FALSE
#' @param pop Data.frame with two integer columns \code{cohort} with year of the
#' birth cohorts and \code{pop} with the size of the corresponding birth
#' cohorts, Default: pop1
#' @param tables List of data.frames to update the tables in fhcrcData, Default:
#' IHE (NB: fhcrcData$ptrx gets changed below to stockholmTreatment if parameters$stockholmTreatment = TRUE -- which is the default! This implies that IHE is the default when parameters$stockholmTreatment = FALSE)
#' @param debug Boolean to print debugging, Default: FALSE
#' @param parms List to update FhcrcParameters, Default: NULL
#' @param mc.cores Integer with the number of cores to use for the computation,
#' Default: 1
#' @param cl a cluster object, created by the \code{parallel} package or by the \code{snow}
#' package. If NULL, use \code{mclapply} (i.e. does not use the registered
#' default cluster).
#' @param print.timing Boolean should the required time be printed after the
#' simulation run, Default: TRUE
#' @param ... TBA
#' @return A fhcrc object
#' @details TBA
#' @examples
#' \dontrun{
#' if(interactive()){
#' library(prostata)
#' sim1 <- callFhcrc(n=1e6, screen="screenUptake", mc.cores=3)
#' }
#' }
#' @seealso
#' \code{\link[parallel]{mclapply}}
#' @rdname callFhcrc
#' @export
#' @importFrom parallel mclapply
callFhcrc <- function(n=10, screen= "noScreening", nLifeHistories=10,
seed=12345, panel=FALSE, flatPop = FALSE, pop = pop1,
tables = IHE, debug=FALSE, parms = NULL, mc.cores = 1,
cl = NULL,
print.timing = TRUE,...) {
## save the random number state for resetting later
state <- RNGstate(); on.exit(state$reset())
## yes, we use the user-defined RNG
RNGkind("user")
set.user.Random.seed(seed)
call <- match.call()
## birth cohorts that should give approximately the number of men alive in Stockholm in 2012
## check the input arguments
screenT <- c("noScreening", "randomScreen50to70",
"twoYearlyScreen50to70", "fourYearlyScreen50to70",
"screen50", "screen60", "screen70", "screenUptake",
"stockholm3_goteborg", "stockholm3_risk_stratified",
"goteborg", "risk_stratified", "mixed_screening",
"regular_screen", "single_screen",
"introduced_screening_only", "introduced_screening_preference",
"introduced_screening", "stopped_screening",
"cap_control", "cap_study", "sthlm3_mri_arm", "grs_stratified", "grs_stratified_age",
"germany_2018")
screen <- match.arg(screen, screenT)
stopifnot(is.na(n) || is.integer(as.integer(n)))
stopifnot(is.integer(as.integer(nLifeHistories)))
## these enum strings should be moved to C++
stateT <- c("Healthy","Localised","Metastatic")
ext_stateT <- c("Healthy","T1_T2","T3plus","Metastatic")
gradeT <- c("Gleason_le_6","Gleason_7","Gleason_ge_8","Healthy")
eventT <- c("toLocalised","toMetastatic","toClinicalDiagnosis", "toCancerDeath",
"toOtherDeath","toScreen","toBiopsyFollowUpScreen",
"toScreenInitiatedBiopsy","toClinicalDiagnosticBiopsy",
"toScreenDiagnosis", "toOverDiagnosis", "toOrganised","toTreatment",
"toCM","toRP", "toRT","toADT","toUtilityChange","toUtilityRemove",
"toSTHLM3", "toOpportunistic","toT3plus", "toCancelScreens",
"toYearlyActiveSurveillance", "toYearlyPostTxFollowUp", "toMRI",
"toPalliative", "toTerminal", "toDRE")
diagnosisT <- c("NotDiagnosed","ClinicalDiagnosis","ScreenDiagnosis")
treatmentT <- c("no_treatment","CM","RP","RT")
psaT <- c("PSA<3","PSA>=3") # not sure where to put this...
screenIndex <- which(screen == screenT) - 1
timingfunction <- if (print.timing) function(x) print(system.time(x)) else function(x) x
## NB: sample() calls the random number generator (!)
if (is.vector(pop)) {
flatPop <- TRUE
pop <- data.frame(cohort=pop,pop=1)
}
if (is.na(n)) {
cohort <- pop$cohort[rep.int(1:nrow(pop),times=pop$pop)]
n <- length(cohort)
} else {
if (flatPop) {
cohort <- rep(pop$cohort,each=ceiling(n/nrow(pop))) #Need ceiling so int n=!0
n <- ceiling(n/nrow(pop)) * nrow(pop) #To get the chunks right
} else
cohort <- sample(pop$cohort,n,prob=pop$pop/sum(pop$pop),replace=TRUE)
}
cohort <- sort(cohort)
## now separate the data into chunks and set the initial random numbers
if (!is.null(cl)) mc.cores <- length(cl) # HACK
currentSeed <- user.Random.seed()
if (mc.cores==1) {
chunks <- list(cohort)
ns <- 0
initialSeeds <- list(currentSeed)
} else {
ns <- floor((0:mc.cores)/mc.cores*n)
chunks <- lapply(1:mc.cores, function(i) cohort[(ns[i]+1):ns[i+1]])
initialSeeds <- c(list(currentSeed), lapply(ns[-c(1,mc.cores+1)], function(i) advance.substream(currentSeed, i)))
ns <- ns[-length(ns)]
}
## Minor changes to fhcrcData
fhcrcData$biopsyOpportunisticComplianceTable <- swedenOpportunisticBiopsyCompliance
fhcrcData$biopsyFormalComplianceTable <- swedenFormalBiopsyCompliance
fhcrcData$secularTrendTreatment2008OR <- secularTrendTreatment2008OR
fhcrcData$rescreening <- rescreening
fhcrcData$newdre = with(dre, data.frame(psa=psa_low, sensitivity, specificity))
if (!is.null(tables))
for (name in names(tables))
fhcrcData[[name]] <- tables[[name]]
updateParameters <- parms
updateParameters$nLifeHistories <- as.integer(nLifeHistories)
updateParameters$screen <- as.integer(screenIndex)
parameter <- FhcrcParameters
for (name in names(updateParameters)){
parameter[[name]] <- updateParameters[[name]]
}
parameter$g0 <- parameter$g0 / parameter$susceptible
parameter$cost_parameters <- parameter$currency_rate * parameter$cost_parameters
parameter$production <- data.frame(ages = parameter$production$ages,
values = parameter$currency_rate * parameter$production$values)
if (parameter$stockholmTreatment)
fhcrcData$prtx <- stockholmTreatment
if (any(.intersection <- (names(parameter) %in% names(fhcrcData))))
warning("The following tables are being over-written by values from parms: ",
paste0(names(parameter)[.intersection], collapse=", "))
temp <- parameter
parameter <- fhcrcData
for (name in names(temp))
parameter[[name]] <- temp[[name]]
## table checks
parameter$rescreening$total <- parameter$rescreening$total_cat
parameter$prtx$Age <- as.double(parameter$prtx$Age)
parameter$prtx$DxY <- as.double(parameter$prtx$DxY)
parameter$prtx$G <- parameter$prtx$G - 1L
parameter$pradt$Grade <- parameter$pradt$Grade - 1L
parameter$biopsy_sensitivity$Year <- as.double(parameter$biopsy_sensitivity$Year)
parameter$neg_biopsy_to_psa$age <- as.double(parameter$neg_biopsy_to_psa$age)
parameter$neg_biopsy_to_biopsy$age <- as.double(parameter$neg_biopsy_to_biopsy$age)
parameter$cure_m_CM_to_RP$age <- as.double(parameter$cure_m_CM_to_RP$age)
parameter$cure_m_CM_to_RT$age <- as.double(parameter$cure_m_CM_to_RT$age)
parameter$pradt$Age <- as.double(parameter$pradt$Age)
parameter$pradt$DxY <- as.double(parameter$pradt$DxY)
## parameter$biopsyComplianceTable <-
## data.frame(expand.grid(psa=c(4,7,10),age=seq(55,75,by=5)),
## compliance=unlist(parameter$biopsy_frequency[,-(1:2),]))
parameter$survival_local <-
with(parameter$survival_local,
data.frame(Age=as.double(AgeLow),Grade=Grade,Time=as.double(Time),
Survival=Survival))
parameter$survival_dist <-
with(parameter$survival_dist,
data.frame(Grade=Grade,Time=as.double(Time),
Survival=Survival))
.xtabsDFCheck(parameter$prtx,c("Age","DxY","G"))
.xtabsDFCheck(parameter$pradt,c("Tx","Age","DxY","Grade"))
.xtabsDFCheck(parameter$hr_locoregional,c("age","ext_grade","psa10")) # should this be in fhcrcData?
.xtabsDFCheck(parameter$biopsyOpportunisticComplianceTable, c("psa","age"))
.xtabsDFCheck(parameter$rescreening,c("age5","total_cat"))
.xtabsDFCheck(parameter$survival_dist,c("Grade","Time"))
.xtabsDFCheck(parameter$survival_local,c("Age","Grade","Time"))
## .xtabsDFCheck(parameter$dre,c("psa_low","sensitivity","specificity"))
## parameter <- modifyList(fhcrcData, parameter)
pind <- sapply(parameter,class)=="numeric" & sapply(parameter,length)==1
bInd <- sapply(parameter,class)=="logical" & sapply(parameter,length)==1
## check some parameters for sanity
if (panel && parameter$rTPF>1) stop("Panel: rTPF>1 (not currently implemented)")
if (panel && parameter$rFPF>1) stop("Panel: rFPF>1 (not currently implemented)")
if (parameter$Andreas && parameter$MRI_screen)
stop("For 'MRI_screen=TRUE', use 'parms=c(prostata:::ShuangParameters, MRI_screen=TRUE)'")
## if (panel && parameter$MRI_screen)
## stop("Scenarios for 'panel=TRUE' and 'MRI_screen=TRUE' have not been defined")
if (parameter$MRI_clinical && !parameter$MRI_screen)
stop("Scenarios for 'MRI_clinical=TRUE' and 'MRI_screen=FALSE' have not been defined")
## now run the chunks separately
step <- function(i) {
chunk <- chunks[[i]]
set.user.Random.seed(initialSeeds[[i]])
.Call("callFhcrc",
parms=list(n=as.integer(length(chunk)),
firstId=ns[i],
panel=panel, # bool
debug=debug, # bool
cohort=as.double(chunk),
parameter=unlist(parameter[pind]),
bparameter=unlist(parameter[bInd]),
otherParameters=parameter[!pind & !bInd]),
PACKAGE="prostata")
}
if (is.null(cl)) {
timingfunction(out <- parallel::mclapply(1:mc.cores, step, mc.cores=mc.cores))
} else {
clusterEvalQ(cl, {library(prostata); RNGkind("user")})
clusterExport(cl, c("chunks", "initialSeeds", "ns", "panel", "debug", "pind",
"bInd", "parameter"), envir=environment())
timingfunction(out <- parallel::parLapply(cl, 1:length(cl), step))
}
## Apologies: we now need to massage the chunks from C++
## reader <- function(obj) {
## out <- cbind(data.frame(state=enum(obj$state[[1]],stateT),
## dx=enum(obj$state[[2]],diagnosisT),
## psa=enum(obj$state[[3]],psaT),
## cohort=obj$state[[4]]),
## data.frame(obj[-1]))
## out$year <- out$cohort + out$age
## out
## }
ext_state2state <- function(obj) # collapses T-stages to localised
`levels<-`(factor(obj),list(Healthy="Healthy",Localised=list("T1_T2","T3plus"),Metastatic="Metastatic"))
cbindList <- function(obj) # recursive (not used)
if (is.list(obj)) do.call("cbind",lapply(obj,cbindList)) else data.frame(obj)
rbindList <- function(obj) # recursive (not used)
if (is.list(obj)) do.call("rbind",lapply(obj,rbindList)) else data.frame(obj)
rbindExtract <- function(obj,name)
do.call("rbind",lapply(obj, function(obji) data.frame(obji[[name]])))
reader <- function(obj) {
cbind(data.frame(state=ext_state2state(enum(obj[[1]],ext_stateT)),
ext_state=enum(obj[[1]],ext_stateT),
grade=enum(obj[[2]],gradeT),
dx=enum(obj[[3]],diagnosisT),
psa=enum(obj[[4]],psaT),
cohort=obj[[5]]),
data.frame(obj[,-(1:5)]))
}
## grab all of the pt, prev, ut, events from summary
## pt <- lapply(out, function(obj) obj$summary$pt)
if (length(out[[1]]$summary) > 0) {
summary <- lapply(names(out[[1]]$summary),
function(name) do.call(rbind,
lapply(out, function(obj) reader(obj$summary[[name]]))))
names(summary) <- names(out[[1]]$summary)
states <- c("state","ext_state","grade","dx","psa","cohort")
names(summary$prev) <- c(states,"age","count")
names(summary$pt) <- c(states,"age","pt")
names(summary$ut) <- c(states,"age","ut")
names(summary$events) <- c(states,"event","age","n")
if(FALSE) age <- NULL # To pass false-positive check note
summary <- lapply(summary,function(obj) within(obj,year <- cohort+age))
enum(summary$events$event) <- eventT
}
else if (length(out[[1]]$shortSummary) > 0) {
## use shortSummary
summary <- lapply(names(out[[1]]$shortSummary),
function(name) do.call(rbind,
lapply(out, function(obj) obj$shortSummary[[name]])))
names(summary) <- names(out[[1]]$shortSummary)
states <- "state"
names(summary$prev) <- c(states,"age","count")
names(summary$pt) <- c(states,"age","pt")
names(summary$ut) <- c(states,"age","ut")
names(summary$events) <- c(states,"event","age","n")
if(FALSE) age <- NULL # To pass false-positive check note
if ("cohort" %in% names(summary))
summary <- lapply(summary,function(obj) within(obj,year <- cohort+age))
enum(summary$events$event) <- eventT
}
else {
summary <- list()
}
## lifeHistories <- do.call("rbind",lapply(out,function(obj) data.frame(obj$lifeHistories)))
## psarecord <- do.call("rbind",lapply(out,function(obj) data.frame(obj$psarecord)))
## diagnoses <- do.call("rbind",lapply(out,function(obj) data.frame(obj$diagnoses)))
## falsePositives <- do.call("rbind",lapply(out,function(obj) data.frame(obj$falsePositives)))
## parameters <- do.call("rbind",lapply(out,function(obj) data.frame(obj$parameters)))
lifeHistories <- rbindExtract(out,"lifeHistories")
psarecord <- rbindExtract(out,"psarecord")
bxrecord <- rbindExtract(out,"bxrecord")
diagnoses <- rbindExtract(out,"diagnoses")
falsePositives <- rbindExtract(out,"falsePositives")
parameters <- rbindExtract(out,"parameters")
indiv_costs <- do.call(c, lapply(out, "[[", "indiv_costs"))
indiv_utilities <- do.call(c, lapply(out, "[[", "indiv_utilities"))
combine.Means <- function(df) {
n <- sum(df$n)
sum <- sum(df$sum)
sumsq <- sum(df$sumsq)
mean <- sum/n
var <- n/(n-1)*(sumsq/n-mean*mean)
sd <- sqrt(var)
se <- sd/sqrt(n)
data.frame(n,mean,var,sd,se,sum,sumsq)
}
mean_utilities <- combine.Means(df=rbindExtract(out, "mean_utilities"))
mean_costs <- combine.Means(df=rbindExtract(out, "mean_costs"))
appendMeans <- function(x) c(x,
mean.sum = x[["sum"]] / x[["n"]],
mean.sumsq = x[["sumsq"]] / x[["n"]])
natural.history.summary <- data.frame(tmc_minus_t0 = appendMeans(sapply(rbindExtract(out,"tmc_minus_t0"), sum)))
## Identifying elements without name which also need to be rbind:ed
societal.costs <- do.call("rbind",lapply(out,function(obj) data.frame(obj$costs))) #split in sociatal and healthcare perspective
## names(costs) <- c("type","item","cohort","age","costs")
names(societal.costs) <- c("type","item","age","costs")
societal.costs$type <- factor(ifelse(societal.costs$type,
"Productivity loss",
"Health sector cost")) # societal perspective
healthsector.costs <- societal.costs[societal.costs["type"] == "Health sector cost", c("item", "age", "costs")] # healthcare perspective
names(lifeHistories) <- c("id", "ext_state", "ext_grade", "dx", "event", "begin", "end", "year", "psa", "utility")
enum(lifeHistories$ext_state) <- ext_stateT
lifeHistories$state <- ext_state2state(lifeHistories$ext_state)
lifeHistories <- lifeHistories[c(names(lifeHistories)[1], "state", names(lifeHistories)[-1])] # shift col order
enum(lifeHistories$dx) <- diagnosisT
enum(lifeHistories$event) <- eventT
enum(diagnoses$ext_state) <- ext_stateT
diagnoses$state <- ext_state2state(diagnoses$ext_state)
enum(diagnoses$ext_grade) <- gradeT
enum(diagnoses$dx) <- diagnosisT
enum(diagnoses$tx) <- treatmentT
enum <- list(stateT = stateT, ext_stateT = ext_stateT, eventT = eventT, screenT = screenT,
diagnosisT = diagnosisT, psaT = psaT)
out <- list(n=n,screen=screen,enum=enum,lifeHistories=lifeHistories,
parameters=parameters, summary=summary,
healthsector.costs=healthsector.costs, societal.costs=societal.costs,
psarecord=psarecord, diagnoses=diagnoses, bxrecord=bxrecord,
cohort=data.frame(table(cohort)),simulation.parameters=parameter,
falsePositives=falsePositives, panel=panel, call = call,
natural.history.summary=natural.history.summary,
indiv_costs=indiv_costs,
indiv_utilities=indiv_utilities,
mean_utilities=mean_utilities,
mean_costs=mean_costs)
class(out) <- "fhcrc"
out
}
## R --slave -e "options(width=200); require(microsimulation); callFhcrc(100,nLifeHistories=1e5,screen=\"screen50\")[[\"parameters\"]]"
#' @title Summarise simulation results
#' @description FUNCTION_DESCRIPTION
#' @param object PARAM_DESCRIPTION
#' @param from Age from which to report
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname summary.fhcrc
#' @export
summary.fhcrc <- function(object, from=0, ...) {
if (from<0) {
warning("from argument should be non-negative - changed to 0")
from <- 0
}
if (from != round(from)) {
warning("from argument should be integer - changed to round(from)")
from <- round(from)
}
if (from>0) {
adjust.costs <- (1+object$simulation.parameters$discountRate.costs)^from
adjust.effectiveness <- (1+object$simulation.parameters$discountRate.effectiveness)^from
object$n <- sum(subset(object$summary$prev,age==from)$count)
object$summary$ut <- subset(object$summary$ut, age>=from)
object$summary$pt <- subset(object$summary$pt, age>=from)
object$summary$prev <- subset(object$summary$prev, age>=from)
object$summary$events <- subset(object$summary$events, age>=from)
object$societal.costs <- subset(object$societal.costs, age>=from)
object$healthsector.costs <- subset(object$healthsector.costs, age>=from)
} else {
adjust.costs <- adjust.effectiveness <- 1
}
newobj <- object[c("n","screen")]
with(object,
structure(.Data=c(newobj,
with(object, list(
from=from,
discountRate.costs = simulation.parameters$discountRate.costs,
discountRate.effectiveness = simulation.parameters$discountRate.effectiveness,
screening.tests = with(summary$events,
sum(n[event == "toScreen"])) / n,
biopsies = with(summary$events,
sum(n[event %in% c("toScreenInitiatedBiopsy",
"toClinicalDiagnosticBiopsy")])) / n,
biopsies.screen.initiated = with(summary$events,
sum(n[event %in% c("toScreenInitiatedBiopsy")])) / n,
biopsies.clinical.diagnostic = with(summary$events,
sum(n[event %in% c("toClinicalDiagnosticBiopsy")])) / n,
negative.biopsies = with(summary$events,
sum(n[event %in% c("toScreenInitiatedBiopsy",
"toClinicalDiagnosticBiopsy")])
- sum(n[event %in% c("toScreenDiagnosis",
"toClinicalDiagnosis")])) / n,
screen.diagnosis = with(summary$events,
sum(n[event == "toScreenDiagnosis"])) / n,
clinical.diagnosis = with(summary$events,
sum(n[event == "toClinicalDiagnosis"])) / n,
over.diagnosis = with(summary$events,
sum(n[event == "toOverDiagnosis"])) / n,
diagnosis = with(summary$events,
sum(n[event %in% c("toClinicalDiagnosis",
"toScreenDiagnosis")])) / n,
cancer.deaths = with(summary$events,
sum(n[event == "toCancerDeath"])) / n,
LE = sum(summary$pt$pt) / n,
QALE = sum(summary$ut$ut) / n * adjust.effectiveness,
healthsector.costs = sum(healthsector.costs$costs) / n * adjust.costs,
productivity.loss = tapply(societal.costs$costs,
societal.costs$type=="Productivity loss",
sum)[["TRUE"]] / n * adjust.costs,
societal.costs = sum(societal.costs$costs) / n * adjust.costs,
healthsector.cost.per.qaly = sum(healthsector.costs$costs) /
sum(summary$ut$ut),
productivity.loss.per.qaly = tapply(societal.costs$costs,
societal.costs$type=="Productivity loss",
sum)[["TRUE"]] /
sum(summary$ut$ut),
societal.cost.per.qaly = sum(societal.costs$costs) /
sum(summary$ut$ut)))),
class="summary.fhcrc"))
}
#' @title Minus operator method for \code{summary.fhcrc} objects
#' @description Calculates the difference between two \code{summary.fhcrc}
#' objects.
#' @param object1 The positive reference \code{summary.fhcrc} object.
#' @param object2 The negative \code{summary.fhcrc} object which is subtracted
#' from \code{object1}.
#' @return A \code{summary.fhcrc} object representing the difference between the
#' two compared \code{summary.fhcrc}.
#' @details Calculates the difference between two \code{summary.fhcrc} objects,
#' e.g. the difference in life expectancy, quality of life etc. The
#' exceptions are \code{discountRate.costs} and
#' \code{discountRate.effectiveness} which are required two be the same for
#' the two scenarios and the original value is returned. If \code{n} has the
#' same value the original value is returned if they differ a list with both
#' is returned.
#' @examples
#' \dontrun{
#' if(interactive()){
#' scenario1 <- summary(callFhcrc(screen = "screenUptake"))
#' scenario2 <-summary(callFhcrc(screen = "noScreening"))
#' scenario1 - scenario2
#' }
#' }
#' @rdname summary.fhcrc
#' @export
'-.summary.fhcrc' <- function(object1, object2) {
if(object1$from != object2$from)
stop("The two scenarios have different 'from' arguments.")
if(object1$discountRate.costs != object2$discountRate.costs)
stop("The two scenarios have different discount rates for the costs.")
if(object1$discountRate.effectiveness != object2$discountRate.effectiveness)
stop("The two scenarios have different discount rates for the effectiveness.")
if(object1$n != object2$n)
warning("The two scenarios have different number of men.")
structure(.Data= c(list(n = if(object1$n == object2$n) {object1$n} else {list(object1$n, object2$n)},
screen = sprintf("%s-%s", object1$screen, object2$screen),
discountRate.costs = object1$discountRate.costs,
discountRate.effectiveness = object1$discountRate.effectiveness),
## exempt special elements calc difference on all others
sapply(names(object1)[!names(object1) %in% c("n", "screen",
"discountRate.costs",
"discountRate.effectiveness")],
function(x) object1[[x]] - object2[[x]],
simplify = FALSE, USE.NAMES = TRUE)),
class="summary.fhcrc")
}
#' @title Multiplaction with scalar operator method for \code{summary.fhcrc} objects
#' @description Calculates the results per a number of persons of a
#' \code{summary.fhcrc} object.
#' @param obj The \code{summary.fhcrc} object.
#' @param perpersons a scalar to multiply the summary results with.
#' @return A \code{summary.fhcrc} object scaled for a number of persons.
#' @details N.b. this is unfortunatly not a commutative
#' operator. Therefor the first parameter must be the
#' \code{summary.fhcrc} object and the second parameter a nummeric
#' scalar.
#' @examples
#' \dontrun{
#' if(interactive()){
#' ## Example 1
#' scenario1 <- summary(callFhcrc(screen = "screenUptake"))
#' scenario1 * 1000 # the summery per thousand persons
#'
#' ## Example 2
#' scenario2 <-summary(callFhcrc(screen = "noScreening"))
#' (scenario1 - scenario2) * 1000 # the difference per thousand persons
#' }
#' }
#' @rdname summary.fhcrc
#' @export
'*.summary.fhcrc' <- function(obj, perpersons) {
if(!is.numeric(perpersons))
stop("The 'perpersons' variable needs to be a numeric.")
structure(.Data= c(list(n = ifelse(is.list(obj$n),
list(lapply(obj$n, function(x) x * perpersons)),
obj$n * perpersons),
screen = sprintf("(%s)*%d", obj$screen, perpersons),
discountRate.costs = obj$discountRate.costs,
discountRate.effectiveness = obj$discountRate.effectiveness),
## exempt special elements scale all others
lapply(obj[!names(obj) %in% c("n", "screen",
"discountRate.costs",
"discountRate.effectiveness")],
function(x) x * perpersons)),
class="summary.fhcrc")
}
#' @title Divide with scalar operator method for \code{summary.fhcrc} objects
#' @description Calculates the results divided by a denominator \code{summary.fhcrc} object.
#' @param obj The \code{summary.fhcrc} object.
#' @param denominator a scalar to divide the summary results with.
#' @return A \code{summary.fhcrc} object scaled by a denominator.
#' @details N.b. this is unfortunatly not a commutative
#' operator. Therefor the first parameter must be the
#' \code{summary.fhcrc} object and the second parameter a nummeric
#' scalar.
#' @examples
#' \dontrun{
#' if(interactive()){
#' }
#' }
#' @rdname summary.fhcrc
#' @export
'/.summary.fhcrc' <- function(obj, denominator) {
if(!is.numeric(denominator))
stop("The 'denominator' variable needs to be a numeric.")
structure(.Data= c(list(n = ifelse(is.list(obj$n),
list(lapply(obj$n, function(x) x / denominator)),
obj$n * denominator),
screen = sprintf("(%s)/%d", obj$screen, denominator),
discountRate.costs = obj$discountRate.costs,
discountRate.effectiveness = obj$discountRate.effectiveness),
## exempt special elements scale all others
lapply(obj[!names(obj) %in% c("n", "screen",
"discountRate.costs",
"discountRate.effectiveness")],
function(x) x / denominator)),
class="summary.fhcrc")
}
#' @title Print a selection of the summarised simulation
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname print.summary.fhcrc
#' @export
print.summary.fhcrc <- function(x, ...) {
obj <- x
cat(sprintf(
"Screening scenario: %s
Life expectancy: %f
Discounted QALE: %f
Discounted health sector costs: %f
Discounted societal costs: %f
Discounted rate (effect.): %f
Discounted rate (costs): %f
", obj$screen, obj$LE, obj$QALE,
obj$healthsector.costs, obj$societal.costs,
obj$discountRate.effectiveness,
obj$discountRate.costs))
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param object1 PARAM_DESCRIPTION
#' @param object2 PARAM_DESCRIPTION
#' @param perspective PARAM_DESCRIPTION, Default: c("societal.costs", "healthsector.costs")
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname ICER.fhcrc
#' @export
ICER.fhcrc <- function(object1, object2,
perspective = c("societal.costs",
"healthsector.costs"), from=0, ...) {
perspective <- match.arg(perspective)
p1 <- object1$simulation.parameters
p2 <- object2$simulation.parameters
stopifnot(p1$discountRate.costs == p2$discountRate.costs)
stopifnot(p1$discountRate.effectiveness == p2$discountRate.effectiveness)
summary1 <- summary(object1, from=from, ...)
summary2 <- summary(object2, from=from, ...)
out <- list(ICER.QALE=(summary1[[perspective]] - summary2[[perspective]]) /
(summary1$QALE - summary2$QALE),
delta.QALE=summary1$QALE - summary2$QALE,
delta.costs=summary1[[perspective]] - summary2[[perspective]],
from=from, perspective=perspective)
if (p1$discountRate.costs == 0 && p2$discountRate.costs == 0 &&
p1$discountRate.effectiveness == 0 && p2$discountRate.effectiveness == 0) {
out <- c(out,
list(ICER.LE = (summary1[[perspective]] - summary2[[perspective]]) /
(summary1$LE - summary2$LE),
delta.LE = summary1$LE - summary2$LE))
}
out
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname print.fhcrc
#' @export
print.fhcrc <- function(x, ...)
cat(sprintf("FHCRC prostate cancer model with %i individual(s) under scenario '%s'.\n",
x$n, x$screen),
...)
## TODO: better solve issue below with testing.rate for noScreening scenario
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param object PARAM_DESCRIPTION
#' @param scenarios PARAM_DESCRIPTION, Default: NULL
#' @param type PARAM_DESCRIPTION, Default: 'incidence.rate'
#' @param group PARAM_DESCRIPTION, Default: 'age'
#' @param age.breaks PARAM_DESCRIPTION, Default: NULL
#' @param year.breaks PARAM_DESCRIPTION, Default: NULL
#' @param cohort.breaks PARAM_DESCRIPTION, Default: NULL
#' @param age.weights PARAM_DESCRIPTION, Default: NULL
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @importFrom stats predict
#' @rdname predict.fhcrc
#' @export
predict.fhcrc <- function(object, scenarios=NULL, type = "incidence.rate",
group = "age", group.pt = group, age.breaks = NULL, year.breaks = NULL,
cohort.breaks = NULL, age.weights = NULL, ...) {
if(!inherits(object,"fhcrc")) stop("Expecting object to be an fhcrc object")
if(!(is.null(scenarios) || all(sapply(scenarios,inherits,"fhcrc")) || inherits(object,"fhcrc")))
stop("Expecting scenarios is NULL, a fhcrc object or a list of fhcrc objects")
if(is.null(scenarios) && (grepl(".?rate.?ratio$|.?rr$", type, ignore.case = TRUE)))
stop("More than one simulation object is needed to calculate a rate-ratio")
## Stripping of potential rate ratio option before matching
abbr_type <- match.arg(sub(".?rr$|.?rate.?ratio$",
"", type, ignore.case = TRUE),
c("incidence.rate", "symptomatic.incidence.rate",
"screen.incidence.rate", "overdiagnosis.rate",
"testing.rate", "biopsy.rate", "metastasis.rate",
"pc.mortality.rate", "allcause.mortality.rate",
"prevalence", "ly", "life.years"))
## Allowing for several groups
group <- match.arg(group,
c("state", "ext_state", "grade", "dx", "psa", "age", "year", "cohort"),
several.ok = TRUE)
event_types <- switch(abbr_type,
incidence.rate = c("toClinicalDiagnosis", "toScreenDiagnosis"),
symptomatic.incidence.rate = "toClinicalDiagnosis",
screen.incidence.rate = "toScreenDiagnosis",
overdiagnosis.rate = "toOverDiagnosis",
testing.rate = "toScreen",
biopsy.rate = c("toClinicalDiagnosticBiopsy", "toScreenInitiatedBiopsy"),
metastasis.rate = "toMetastatic",
pc.mortality.rate = "toCancerDeath",
allcause.mortality.rate = c("toCancerDeath", "toOtherDeath"))
interval2break <- function(interval_vector) {
as.numeric(unique(unlist(regmatches(interval_vector,
gregexpr("[0-9]+|Inf", interval_vector)))))
}
if(FALSE) {age <- year <- cohort <- Freq <- event <- Freq.y <- Freq.x <-
scenario.x <- rate.x <- rate.y <- NULL} # To pass false-positive check note
## Manipulate input, make sure age.breaks exists from age standardisation
if(!is.null(age.weights)) age.breaks <- interval2break(age.weights[,1])
## Manipulate input, make sure time group exists for specified time interval
for (scale in c("age", "year", "cohort")) {
if(!is.null(eval(parse(text = paste(scale, "breaks",sep="."))))) {
if(!(scale %in% group)) group <- c(group, scale)
if(!(scale %in% group.pt)) group.pt <- c(group.pt, scale)
}
}
## fast operations by group using base-R
## http://stackoverflow.com/questions/3685492/r-speeding-up-group-by-operations
grp_apply = function(XS, INDEX, FUN, ..., simplify=T) {
FUN = match.fun(FUN)
if (!is.list(XS))
XS = list(XS)
as.data.frame(as.table(tapply(1:length(XS[[1L]]), INDEX, function(s, ...)
do.call(FUN, c(lapply(XS, `[`, s), list(...))), ..., simplify=simplify)))
}
## Fixes colnames after group operation
name_grp <- function(x, grp = group) {names(x)[grep("^Var[0-9]+$", names(x))] <- grp; x}
## After a group operation the group will become a factor, if the time scale
## is not reported as an interval we convert it back to numeric. N.b. This
## needs to be called after name_grp().
numeric_time_scale <- function(df, group) {
within(df, {
if("age" %in% group && is.null(age.breaks)) age <- as.numeric(levels(age))[age] #important factor conversion
if("year" %in% group && is.null(year.breaks)) year <- as.numeric(levels(year))[year] #important factor conversion
if("cohort" %in% group && is.null(cohort.breaks)) cohort <- as.numeric(levels(cohort))[cohort] #important factor conversion
})
}
categorise_time <- function(df, age.breaks, year.breaks, cohort.breaks) {
cut_t <- function(time, time.breaks) {
## using first column or single vector as breaks
as.factor(cut(time, breaks = data.frame(time.breaks)[,1],
right = FALSE, dig.lab=4, ordered_result = TRUE))
}
if(!is.null(age.breaks)) df <- transform(df, age = cut_t(age, age.breaks))
if(!is.null(year.breaks)) df <- transform(df, year = cut_t(year, year.breaks))
if(!is.null(cohort.breaks)) df <- transform(df, cohort = cut_t(cohort, cohort.breaks))
df
}
## For standardising over a time scale e.g. ages
standardise_time <- function(df, timestr = "age", parstr = "rate",
time.weights = age.weights) {
if(is.null(time.weights)) {
return(df)
} else {
collapsed_grps <- c(group[!group == timestr], "scenario")
## Scale with weight
weight <- function(time, par, time.weights) {
par * sapply(t(time), {function(x) time.weights[x == time.weights[,1],2]})
}
## Sum over the groups (at least scenario)
collapse <- function(df, collapsed_grps, parstr) {
within(with(df, grp_apply(eval(parse(text = parstr)),
lapply(as.list(collapsed_grps),
{function(x) eval(parse(text = x))}),
sum, na.rm = TRUE)) ,{
assign(eval(substitute(parstr)),
Freq / sum(time.weights[,2]))
rm(Freq)
})
}
## Standardise and rename
return(numeric_time_scale(name_grp(collapse(
within(df,
{assign(eval(substitute(parstr)), weight(eval(parse(text =timestr)),
eval(parse(text = parstr)),
time.weights))}),
collapsed_grps, parstr), collapsed_grps), collapsed_grps))
}
}
## Calculates rates of specific events by specified groups
calc_rate <- function(object, event_types, group){
pt <- with(categorise_time(object$summary$pt, age.breaks, year.breaks,
cohort.breaks),
numeric_time_scale(name_grp(grp_apply(pt,
lapply(as.list(group.pt), function(x) eval(parse(text = x))),
sum), grp = group.pt), group.pt))
## TODO: if subset has no dim replace with zeros, temp fix below:
if(!any(object$summary$event$event %in% event_types)) {
stop(paste("The event(s)", paste(event_types, collapse = ", "),
"was not found in the", object$screen, "scenario"))
}
events <- with(subset(
categorise_time(object$summary$events, age.breaks, year.breaks,
cohort.breaks), event %in% event_types),
numeric_time_scale(name_grp(grp_apply(n,
lapply(as.list(group),
{function(x) eval(parse(text = x))}),
sum)), group))
within(merge(pt, events, by = group.pt, all.y = TRUE),{
rate <- ifelse(is.na(Freq.y) & !is.na(Freq.x), 0, Freq.y/Freq.x) #no events but some pt -> 0
n <- Freq.y
pt <- Freq.x
rm(Freq.x,Freq.y)})
}
## Calculate prevalences by specified groups
calc_prev <- function(object, group){
within(with(categorise_time(object$summary$prev, age.breaks,
year.breaks, cohort.breaks),
numeric_time_scale(name_grp(grp_apply(count,
lapply(as.list(group),
function(x) eval(parse(text = x))), sum)), group)), {
prevalence <- Freq/object$n
rm("Freq")})
}
## Calculate life-years by specified groups
calc_ly <- function(object, group){
within(with(categorise_time(object$summary$pt, age.breaks,
year.breaks, cohort.breaks),
numeric_time_scale(name_grp(grp_apply(pt,
lapply(as.list(group),
function(x) eval(parse(text = x))), sum)), group)), {
ly <- Freq/object$n
rm("Freq")})
}
## Calculates the outcome in the passed function for all
## simulation objects in the 'scenarios' list. Then the object
## outcome (e.g. rates or prev) are for the scenarios are added as
## rows and the scenario name as a column.
predict_scenarios <- function(scenarios, calc_outcome, ...) {
do.call(rbind, lapply(scenarios,
{function(object, ...)
cbind(calc_outcome(object, ...), scenario = object$screen)}, ...))
}
## Input checks allow for scenarios to be a single fhcrc object or
## list of fhcrc objects. Now make sure scenarios is a list.
if(inherits(scenarios, "fhcrc")) {scenarios <- list(scenarios)}
## Rate-ratio if type ends with rate.ratio or RR
if(grepl(".?rate.?ratio$|.?rr$", type, ignore.case = TRUE)){
scenario_rates <- standardise_time(predict_scenarios(unique(scenarios),
calc_rate,
event_types, group),
timestr = "age",
parstr = "rate",
time.weights = age.weights)
reference_rate <- standardise_time(predict_scenarios(list(object),
calc_rate,
event_types, group),
timestr = "age",
parstr = "rate",
time.weights = age.weights)
## Alt. check if it is part of the rates
if(!is.null(age.weights)) group <- group[!group == "age"]
within(merge(scenario_rates, reference_rate, by = group),{
scenario <- scenario.x
rate.ratio <- rate.x/rate.y
rate.ratio[!is.finite(rate.ratio)] <- NaN
rm(list=ls(pattern=".x$|.y$"))})
## Prevalence if type ends with rate.ratio or RR
} else if(grepl(".?prev$|.?prevalence$", type, ignore.case = TRUE)){
standardise_time(predict_scenarios(unique(c(list(object),scenarios)), calc_prev, group),
timestr = "age", parstr = "prevalence", time.weights = age.weights)
## Life-years if type is ly or life*years
} else if(grepl("^ly$|^life.?years$", type, ignore.case = TRUE)){
standardise_time(predict_scenarios(unique(c(list(object),scenarios)), calc_ly, group),
timestr = "age", parstr = "LY", time.weights = age.weights)
## Defaults to plain rates. If reference object exist add it to
## scenario list and remove duplicates.
} else {
standardise_time(predict_scenarios(unique(c(list(object),scenarios)), calc_rate, event_types, group),
timestr = "age", parstr = "rate", time.weights = age.weights)
}
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param type PARAM_DESCRIPTION, Default: c("incidence.rate", "testing.rate", "biopsy.rate", "metastasis.rate",
#' "pc.mortality.rate", "allcause.mortality.rate")
#' @param plot.type PARAM_DESCRIPTION, Default: 'l'
#' @param add PARAM_DESCRIPTION, Default: FALSE
#' @param xlab PARAM_DESCRIPTION, Default: 'Age (years)'
#' @param ylab PARAM_DESCRIPTION, Default: NULL
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @importFrom graphics plot
#' @rdname plot.fhcrc
#' @export
plot.fhcrc <- function(x, type=c("incidence.rate", "testing.rate",
"biopsy.rate", "metastasis.rate",
"pc.mortality.rate",
"allcause.mortality.rate"),
plot.type="l", add=FALSE, xlab="Age (years)",
ylab=NULL, ...) {
type <- match.arg(type)
if (is.null(ylab)) {ylab <- switch(type,
incidence.rate="Prostate cancer incidence rates per 100,000",
testing.rate="PSA rates per 1000",
biopsy.rate="Biopsies per 1000",
metastasis.rate="Metastatic onset per 100,000",
pc.mortality.rate="Cancer mortality rates per 100,000",
allcause.mortality.rate="All cause mortality rates per 100,000")}
rates <- predict(object = x, type = type)
rates$rate = rates$rate*switch(type, testing.rate=1000, biopsy.rate=1000, incidence.rate=1e5, metastasis.rate=1e5,pc.mortality.rate=1e5,allcause.mortality.rate=1e5)
if (!add) plot(rate~age, data=rates, type=plot.type, xlab=xlab, ylab=ylab, ...) else lines(rate~age, data=rates, ...)
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname lines.fhcrc
#' @importFrom graphics lines
#' @export
lines.fhcrc <- function(x,...) {
plot(x, ..., add=TRUE)
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param object1 PARAM_DESCRIPTION
#' @param object2 PARAM_DESCRIPTION
#' @param startAge PARAM_DESCRIPTION, Default: 50
#' @param stopAge PARAM_DESCRIPTION, Default: Inf
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname nn.fhcrc
#' @export
nn <- function(object1, object2, startAge = 50, stopAge = Inf, ...){
UseMethod("nn")
}
#' @export
nn.fhcrc <- function(object1, object2, startAge = 50, stopAge = Inf, ...) {
pNNS <- function(thisScenario) {
with(thisScenario$summary$events,
sum(n[event=="toCancerDeath" & age>=startAge
& age<stopAge])) / # divided by
with(thisScenario$summary$prev,
sum(count[age==round(startAge)]))
}
pNND <- function(thisScenario) {
with(thisScenario$summary$events,
sum(n[event=="toCancerDeath" & age>=startAge & age<stopAge])) / # divided by
with(thisScenario$summary$events,
sum(n[event %in% c("toScreenDiagnosis","toClinicalDiagnosis")
& age>=startAge & age<stopAge]))
}
NNS <- 1 / (pNNS(object2) - pNNS(object1)) #number needed to screen to prevent 1 PCa death
NND <- 1 / (pNND(object2) - pNND(object1)) #number needed to detect to prevent 1 PCa death
## Include additional number needed to treat (NNT) [Gulati 2011] to show overdiagnosis?
structure(.Data = list(NNS=NNS,NND=NND), class = "nn.fhcrc")
}
mclements/prostata documentation built on Feb. 1, 2023, 1:20 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions nn.fhcrc nn lines.fhcrc plot.fhcrc predict.fhcrc print.fhcrc ICER.fhcrc print.summary.fhcrc summary.fhcrc callFhcrc .xtabsDFCheck
Documented in callFhcrc ICER.fhcrc lines.fhcrc nn plot.fhcrc predict.fhcrc print.fhcrc print.summary.fhcrc summary.fhcrc
#' prostata
#'
#' @section Introduction:
#'
#' The packages provides implementations of discrete event simulation in both R
#' and C++.
#'
#' @section Background:
#'
#' Prostata is a natural history model of prostate cancer model which extends
#' the model developed by Ruth Etzioni and colleagues at the
#' [[http://www.fredhutch.org][Fred Hutchinson Cancer Research Center
#' (FHCRC)]]. This is a well validated prostate cancer natural history model
#' developed within the
#' [[http://cisnet.cancer.gov/prostate/profiles.html][Cancer Intervention and
#' Surveillance Modeling Network (CISNET)]]. We here provide our extensions of
#' the model with extended states and finer calibration. The model is provided
#' as an R package, depending on our
#' [[https://github.com/mclements/microsimulation][microsimulation]] frame work.
#' We specifically developed this package for modelling the cost-effectiveness
#' of prostate cancer screening, where many (e.g. 10^7) men are followed from
#' birth to death.
#'
#' @section Running the simulation:
#'
#' There are a number of available testing scenarios. They determine testing
#' frequencies and re-testing intervals over calendar period and ages.
#'
#' @docType package
#' @name prostata-package
#' @aliases prostata
#' @author Mark Clements \email{mark.clements@ki.se}
#' @references \url{https://github.com/mclements/prostata}
#' @references \url{https://github.com/mclements/microsimulation}
#' @seealso \code{\link{Rcpp}}
#' @useDynLib prostata, .registration=TRUE
#' @import microsimulation
#' @importFrom Rcpp evalCpp compileAttributes
#' @importFrom utils packageName
NULL
if(.Platform$OS.type == "unix") {
pkg <- packageName()
## http://r.789695.n4.nabble.com/How-to-construct-a-valid-seed-for-l-Ecuyer-s-method-with-given-Random-seed-td4656340.html
#' @param seed Random seed
#' @keywords internal
#' @importFrom microsimulation set.user.Random.seed
set.user.Random.seed <- function (seed)
microsimulation::set.user.Random.seed(seed,pkg)
#' @keywords internal
#' @importFrom microsimulation next.user.Random.substream
next.user.Random.substream <- function ()
microsimulation::next.user.Random.substream(pkg)
#' @keywords internal
#' @importFrom microsimulation user.Random.seed
user.Random.seed <- function()
microsimulation::user.Random.seed(pkg)
}
##
#' @title Initial values for the FHCRC model
#' @description Initial values for the FHCRC model
#' @format A list
#' \describe{
#'}
#' @details Currently, see the R documentation
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname FhcrcParameters
#' @export
"FhcrcParameters"
FhcrcParameters <- list(
Andreas = TRUE, # version for Andreas's CEA paper
revised_natural_history=TRUE,
## panel=FALSE,
grade.clinical.rate.high=0.3042454700,
startFullUptake = 1932.0,
fullUptakePortion = 0.9,
yearlyUptakeIncrease = 0.03,
startUptakeMixture = 1945.0,
endUptakeMixture = 1960.0,
uptakeStartAge=35.0,
screeningIntroduced = 1995.0,
shapeA = 3.8,
scaleA = 15.0,
shapeT = 2.16,
scaleT = 11.7,
tau2 = 0.0829, # log PSA measurement variance - normal
susceptible = 1.0, # portion susceptible
g0=0.0005, # onset parameter
g3p = exp(-6.73546100), # T3+ parameter
gm = exp(-7.46200003), # metastatic parameter
gc=0.0015, # clinical diagnosis parameter
thetac=19.1334, # clinical diagnosis parameter after metastatic
mubeta0=-1.609, # mean of beta0, where beta0 is the log PSA intercept at age 35 years
sebeta0=0.2384, # SE of beta0
mubeta1=0.04463, # mean of beta1, where beta1 is the log PSA slope prior to cancer onset
sebeta1=0.0430, # SE of beta1
mubeta2=c(0.0397,0.1678, 0.0), # base::grade: the gleason specific change in slope of log PSA after cancer onset
sebeta2=c(0.0913,0.3968, 0.0), # base::grade: variance of beta2
rev_mubeta2=c(0.051, 0.129, 0.1678), # ext::grade: same as above for extended gleason grade (6-, 7, 8+)
rev_sebeta2=c(0.064, 0.087, 0.3968), # ext::grade
alpha7 = log(0.2), # log of the proportion gleason 7 at age 35
beta7 = 0.06840404, # slope of log proportion of gleason 7
alpha8 = log(0.002), # log of the proportion of gleason 8+ at age 35
beta8 = 0.19894431, # slope of log proportion of gleason 8+
RR_T3plus=2.0, # prostate cancer mortality rate ratio comparing T3+ with T1-T2, add lit reference
## mubeta2.scale=1.0, # cf. 2.1
## beta.rho=0.62,
c_txlt_interaction = 1.0, # treatment lead time interaction, default to no effect
c_baseline_specific = 1.0, # scale factor on survival, default to no effect
c_benefit_value0 = 10, # (value -> reduction): 0.04 -> 10%; 0.18 -> 20%; 10 -> 28%, generalised stage-shift parameter
sxbenefit = 1.0, # hazard rate ratio for screening benefit, defaults to no effect
c_benefit_type = 0, # 0=stage-shift (=> c_benefit_value0=10), 1=lead-time based (=> c_benefit_value1=0.1)
c_benefit_value1 = 0.1, # approximate increase in cure for each year of lead-time, used for the lead-time based screening effect
RP_mortHR = 0.56, # mortality hazard ratio for surgery over watchful-waiting from SPCG-4 Bill-Axelson 2014
screeningParticipation = 0.75, # probability of actually having the first PSA test
rescreeningParticipation = 0.95, # probability of actually having the re-screening PSA tests
biopsyCompliance = 0.856, # Formal biopsy compliance, Schroder ERSPC 2014
biopsySensitivityTimeProportionT1T2 = 0.5314886, # time portion when T1-T2 cancers are sensitivity to biopsies (expit from calibration). The remaining part, starting at onset, is not detectable.
studyParticipation = 50.0/260.0, # observed fraction of population who participated in STHLM3 study
nLifeHistories = 10L, screen = 0L, ## integers
psaThreshold = 3.0,
psaThresholdBiopsyFollowUp = 4.0, # revised PSA threshold for negative follow-up screen
biomarker_model = 1, # biomarker_model = 0 random, biomarker_model = 1 psa/risk based correction of FP
PSA_FP_threshold_nCa=4.21, # reduce FP in no cancers with PSA threshold
PSA_FP_threshold_GG6=3.76, # reduce FP in GG 6 with PSA threshold
PSA_FP_threshold_GG7plus=3, # reduce FP in GG >= 7 with PSA threshold
panelReflexThreshold = 1.0,
## Natural history calibration
rTPF=1.0,
rFPF=0.6,
c_low_grade_slope=-0.006,
stockholmTreatment = TRUE,
discountRate.effectiveness = 0.03,
discountRate.costs = 0.03,
full_report = 1.0,
formal_costs = 0.0,
formal_compliance = 0.0,
start_screening = 50.0, # start of organised screening
stop_screening = 70.0, # end of organised screening
screening_interval = 2.0, # screening interval for regular_screening and introduced_screening*
screening_interval1 = 2.0, # screening interval for grs_stratified_age for first age group
screening_interval2 = 4.0, # screening interval for grs_stratified_age for second age group
screening_interval_split = 60.0, # age to split screening interval for grs_stratified_age
introduction_year = 2015.0, # year to start organised screening for introduced_screening* & cancel the opportunistic testing under stopped_screening
risk_psa_threshold=1, # PSA threshold for risk-stratified screening
risk_lower_interval=8, # re-screening interval for lower risk
risk_upper_interval=4, # re-screening interval for higher risk
mu0=c(0.00219, 0.000304, 5.2e-05, 0.000139, 0.000141, 3.6e-05, 7.3e-05,
0.000129, 3.8e-05, 0.000137, 6e-05, 8.1e-05, 6.1e-05, 0.00012,
0.000117, 0.000183, 0.000185, 0.000397, 0.000394, 0.000585, 0.000448,
0.000696, 0.000611, 0.000708, 0.000659, 0.000643, 0.000654, 0.000651,
0.000687, 0.000637, 0.00063, 0.000892, 0.000543, 0.00058, 0.00077,
0.000702, 0.000768, 0.000664, 0.000787, 0.00081, 0.000991, 9e-04,
0.000933, 0.001229, 0.001633, 0.001396, 0.001673, 0.001926, 0.002217,
0.002562, 0.002648, 0.002949, 0.002729, 0.003415, 0.003694, 0.004491,
0.00506, 0.004568, 0.006163, 0.006988, 0.006744, 0.00765, 0.007914,
0.009153, 0.010231, 0.011971, 0.013092, 0.013839, 0.015995, 0.017693,
0.018548, 0.020708, 0.022404, 0.02572, 0.028039, 0.031564, 0.038182,
0.042057, 0.047361, 0.05315, 0.058238, 0.062619, 0.074934, 0.089776,
0.099887, 0.112347, 0.125351, 0.143077, 0.153189, 0.179702, 0.198436,
0.240339, 0.256215, 0.275103, 0.314157, 0.345252, 0.359275, 0.41768,
0.430279, 0.463636, 0.491275, 0.549738, 0.354545, 0.553846, 0.461538,
0.782609),
hr_locoregional = transform(expand.grid(age=c(50,60,70), ext_grade=0:2,
psa10=0:1),
hr = c(0.2640235, 0.5107137,
1.0917129, 1.9397446,
1.6671078, 1.9372717,
1.7998275, 1.0379953,
1.7205986, 0.8606931,
0.8966157, 1.1941040,
3.6151038, 2.9625742 ,
2.1622938, 1.0185837,
0.9337168, 0.8849155)),
hr_metastatic = data.frame(age = c(50, 60, 70),
hr = c(0.893, 0.779, 0.698)),
biopsy_sensitivity = data.frame(Year = c(1987, 1988, 1989, 1990, 1991, 1992,
1993, 1994, 1995, 1996, 1997, 1998,
1999, 2000),
Sensitivity = c(0.502, 0.574, 0.642, 0.702,
0.749, 0.782, 0.803, 0.814,
0.821, 0.830, 0.844, 0.866,
0.894, 0.925)),
rand_biopsy_sensitivityG6 = 1.0,
neg_biopsy_to_psa=structure(list(age = c(50, 60, 70, 80),
meanlog = c(-0.314249298970759, -0.36332284250614,
-0.333753007363153, -0.175063351183678),
sdlog = c(0.996188493950768, 0.974965602213158,
1.0490009138328, 1.23085968935057)),
.Names = c("age", "meanlog", "sdlog"),
row.names = c(NA, -4L),
class = "data.frame"),
neg_biopsy_to_biopsy = structure(list(age = c(50, 60, 70, 80),
meanlog = c(3.05337156798731, 2.65573907207411,
1.78297474442106, 1.11581637027693),
sdlog = c(2.41326118058954, 2.3225660366395,
1.86615343381663, 2.04687048711335)),
.Names = c("age", "meanlog", "sdlog"),
row.names = c(NA, -4L), class = "data.frame"),
cure_m_CM_to_RP = structure(list(age = c(50, 60, 70, 80),
pnever = c(0.485026452157334, 0.605786246502961,
0.922871191840356, 0.999745963308522),
meanlog = c(1.14257205575071, 1.04021691430639,
0.519760988362341, 4.35555555555556),
sdlog = c(0.758280639791713, 0.875661155361679,
0.782118228501619, 0.255901426097166)),
.Names = c("age", "pnever", "meanlog", "sdlog"),
row.names = c(NA, -4L), class = "data.frame"),
cure_m_CM_to_RT =structure(list(age = c(50, 60, 70, 80),
pnever = c(0.60082830813107,
2.71246734638524e-06,
0.493326620739727,
0.739523343815296),
meanlog = c(3.14746283053335,
4.24925205843398,
2.0235406338416,
0.705091563838411),
sdlog = c(2.0453909799028,
2.26870437384556,
2.1703325154377,
1.81819516420271)),
.Names = c("age", "pnever", "meanlog", "sdlog"),
row.names = c(NA, -4L), class = "data.frame"),
currency_rate = 0.750/9.077, # PPP EU19@2017/Swe@2016 https://data.oecd.org/conversion/purchasing-power-parities-ppp.htm
## Should we add cost and dis-utilities for using test characteristics based on prostate volume 7% from supplement
## Swedish governmental report on organised PSA testing (p.22):
## https://www.socialstyrelsen.se/SiteCollectionDocuments/2018-2-13-halsoekonomisk-analys.pdf
## Based on the Swedish south region 2017:
## http://sodrasjukvardsregionen.se/avtal-priser/regionala-priser-och-ersattningar-foregaende-ar/
## S3M cost from Karolinska University Laboratory (KUL):
## https://www.karolinska.se/KUL/Alla-anvisningar/Anvisning/10245
cost_parameters = c("Invitation" = 7 # Invitation letter
+ 7, # Results letter
"Formal PSA" = 349 # test sampling, primary care
+ 45 # PSA analysis
+ 0 * 1539, # No GP primary care
"Formal panel" = 349 # test sampling, primary care
+ 45 # PSA analysis not included in panel price
+ 2300 # From KUL price list
+ 0 * 1539, # No GP for formal
"Opportunistic PSA" = 349 # test sampling, primary care
+ 45 # PSA analysis
+ 0.2 * 1539, # GP primary care
"Opportunistic panel" = 349 # test sampling, primary care
+ 45 # PSA analysis not included in panel price
+ 2300 # From KUL price list
+ 0.2 * 1539, # GP primary care
"Biopsy" = 4733, # Biopsy cost
"Assessment" = 1794, # Urology assessment
"Prostatectomy" = 80000 # Surgery
+ 100000 * 0.25 # Radiation therapy
+ 1794 * 2 # Urology visit
+ 1205 * 2, # Nurse visit
"Radiation therapy" = 100000 # Radiation therapy
+ 1794 * 2 # Urology visit
+ 1205 * 2, # Nurse visit
"Active surveillance - yearly" = 1794 # Urology visit
+ 45 * 2 # PSA analysis
+ 4733 * 0.5, # Biopsy
"Active surveillance - single MR" = 3090, # Used once for active surveillance
"Post-Tx follow-up - yearly" = 349 # PSA test sampling
+ 45 # PSA analysis,
+ 474, # Telefollow-up by urologist
"Cancer death" = 100160 * 3 # Care for spread disease
+ 68000 * 3, # Drugs for spread disease
"Polygenic risk stratification" = 250*12, # Callender et al (2021) with exchange rate of approximately 12
"Opportunistic DRE" = 349 # DRE procedure in primary care
+ 0.2 * 1539), # part of primary care visit
## Swedish governmental report on organised PSA testing (p.23):
## https://www.socialstyrelsen.se/SiteCollectionDocuments/2018-2-13-halsoekonomisk-analys.pdf
production = data.frame(ages = c(0, 54, 64, 74),
values=apply(
rbind(0.878, 0.756, 0.158, 0), # employment portion, full time
1,
function(empl) empl
* 30.7/40 # average working hours
* 32800*12 # average salary
* 1.485)), # including non-optional social fees
lost_production_years= c("Formal PSA"=2/24/365.25,
"Formal panel"=2/24/365.25,
"Opportunistic PSA"=2/24/365.25,
"Opportunistic panel"=2/24/365.25,
"Assessment"=2/24/365.25, # Urology assessment
"Biopsy"=2/24/365.25,
"Prostatectomy"=6/52,
"Radiation therapy"=8/52,
"Active surveillance - yearly"=
1 * 2/24/365.25 # Urology assessment
+ 2 * 2/24/365.25 # PSA tests
+ 0.5 * 2/24/365.25, # Biopsy
"Metastatic cancer"=6/12,
"Terminal illness" = 6/12), # NOTE: potentially add follow-up Tx production losses
## Heijnsdijk 2012
utility_estimates = c("Invitation" = 1,
"Formal PSA" = 0.99,
"Formal panel" = 0.99,
"Opportunistic PSA" = 0.99,
"Opportunistic panel" = 0.99,
"Biopsy" = 0.90,
"Cancer diagnosis" = 0.80,
"Prostatectomy part 1" = 0.67,
"Prostatectomy part 2" = 0.77,
"Radiation therapy part 1" = 0.73,
"Radiation therapy part 2" = 0.78,
"Active surveillance" = 0.97,
"Postrecovery period" = 0.95,
"Palliative therapy" = 0.60,
"Terminal illness" = 0.40,
"Death" = 0.00),
## Utility duration is given in years.
utility_duration = c("Invitation" = 0.0,
"Formal PSA" = 1/52,
"Formal panel" = 1/52,
"Opportunistic PSA" = 1/52,
"Opportunistic panel" = 1/52,
"Biopsy" = 3/52,
"Cancer diagnosis" = 1/12,
"Prostatectomy part 1" = 2/12,
"Prostatectomy part 2" = 10/12,
"Radiation therapy part 1" = 2/12,
"Radiation therapy part 2" = 10/12,
"Active surveillance" = 7,
"Postrecovery period" = 9,
"Palliative therapy" = 30/12,
"Terminal illness" = 6/12),
utility_truncate = TRUE, # should utilities be truncated at zero?
utility_scales = c("UtilityAdditive"=0,"UtilityMultiplicative"=1,"UtilityMinimum"=2), # encoding for the utility scales
utility_scale = as.double(1), # default scale = UtilityMultiplicative
includeEventHistories = TRUE,
includePSArecords = FALSE,
includeBxrecords = FALSE,
includeDiagnoses = FALSE,
MRI_screen = FALSE, # defines whether MRI pathway is used for screen-positive patients
MRI_clinical = FALSE, # defines whether MRI pathway is used for clinically detected cancers
MRInegSBx=FALSE, # No SBx for MRI- (by default)
rescreenDoubleNeg = FALSE, # defines whether MRI-/Bx- men should move to rescreening (cf. repeat PSA and Bx within ~12 months)
indiv_reports = FALSE, # should the cost and standard reports include individual values?
startReportAge = 0.0, # Age to start reporting (currently only means_utilities and means_costs)
weibull_onset = FALSE, # flag to use Weibull-distributed onset
weibull_onset_shape = 2, # shape; values (0,Inf)
weibull_onset_scale= 40, # scale; values (0,Inf)
frailty = FALSE, # assume a frailty distribution on the onset distribution?
grs_variance = 0.68, # Callender et al (2019)
other_variance = 1.14, # total variance = 1.82 from Kicinski et al (2011; https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0027130)
grs_risk_threshold = 0.04 # ten-year risk threshold for GRS-based risk stratified screening
)
IHE <- list(prtx=data.frame(Age=50.0,DxY=1973.0,G=1:2,CM=0.6,RP=0.26,RT=0.14)) ## assumed constant across ages and periods
ParameterNV <- FhcrcParameters[sapply(FhcrcParameters,class)=="numeric" & sapply(FhcrcParameters,length)==1]
## ParameterIV <- FhcrcParameters[sapply(FhcrcParameters,class)=="integer" & sapply(FhcrcParameters,length)==1]
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 15 rows and 3 variables:
#' \describe{
#' \item{\code{psa}}{double COLUMN_DESCRIPTION}
#' \item{\code{age}}{double COLUMN_DESCRIPTION}
#' \item{\code{compliance}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"swedenOpportunisticBiopsyCompliance"
swedenOpportunisticBiopsyCompliance <- data.frame(
psa = c(3, 5, 10, 3, 5, 10, 3, 5, 10, 3, 5, 10, 3, 5, 10),
age = c(40, 40, 40, 50, 50, 50, 60, 60, 60, 70, 70, 70, 80, 80, 80),
compliance = c(0.3764045, 0.5680751, 0.7727273, 0.3110770, 0.5726548, 0.7537372, 0.2385155, 0.4814588, 0.6929770, 0.1754264, 0.3685056, 0.5602030, 0.1629213, 0.2697368, 0.5010052))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 15 rows and 3 variables:
#' \describe{
#' \item{\code{psa}}{double COLUMN_DESCRIPTION}
#' \item{\code{age}}{double COLUMN_DESCRIPTION}
#' \item{\code{compliance}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"swedenFormalBiopsyCompliance"
swedenFormalBiopsyCompliance <- cbind(expand.grid(psa=c(3,5,10),age=seq(40,80,10)),
compliance=FhcrcParameters$biopsyCompliance)
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 24 rows and 6 variables:
#' \describe{
#' \item{\code{DxY}}{double COLUMN_DESCRIPTION}
#' \item{\code{Age}}{double COLUMN_DESCRIPTION}
#' \item{\code{G}}{integer COLUMN_DESCRIPTION}
#' \item{\code{CM}}{double COLUMN_DESCRIPTION}
#' \item{\code{RP}}{double COLUMN_DESCRIPTION}
#' \item{\code{RT}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"stockholmTreatment"
stockholmTreatment <-
data.frame(DxY=2008,
Age=c(50,50,50,55,55,55,60,60,60,65,65,65,70,70,70,75,75,75,80,80,80,85,85,85),
G=as.integer(c(6,7,8,6,7,8,6,7,8,6,7,8,6,7,8,6,7,8,6,7,8,6,7,8)-6),
CM=c(0.231023,0.044872,0.25,0.333333,0.09542,0.328358,0.409439,0.12825,0.348101,0.479167,0.178182,0.401639,0.689013,0.359143,0.56087,0.876543,0.744444,0.809524,1,0.970711,0.952096,0.9375,1,1),
RP=c(0.700623,0.815839,0.5,0.553592,0.740111,0.326226,0.49981,0.6866,0.291437,0.409879,0.552483,0.279235,0.210949,0.318374,0.179363,0.041152,0.058652,0.049689,0,0.009763,0.023952,0.0625,0,0),
RT=c(0.068354,0.13929,0.25,0.113074,0.164469,0.345416,0.090751,0.185151,0.360462,0.110954,0.269335,0.319126,0.100038,0.322482,0.259767,0.082305,0.196903,0.140787,0,0.019526,0.023952,0,0,0))
secularTrendTreatment2008OR <-
data.frame(year=as.numeric(1988:2009),
OR=c(0.194409,0.155874,0.130704,0.154104,0.148039,0.237752,0.299491,0.330934,0.407047,0.377968,0.443079,0.551047,0.686101,0.700657,0.846134,0.93956,1.03026,1.191726,1.067736,1.012136,1,0.911471))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 52 rows and 5 variables:
#' \describe{
#' \item{\code{age5}}{double COLUMN_DESCRIPTION}
#' \item{\code{total_cat}}{double COLUMN_DESCRIPTION}
#' \item{\code{shape}}{double COLUMN_DESCRIPTION}
#' \item{\code{cure}}{double COLUMN_DESCRIPTION}
#' \item{\code{scale}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"rescreening"
rescreening <- data.frame(age5 = c(30, 30, 30, 30, 35, 35, 35, 35, 40, 40,
40, 40, 45, 45, 45, 45, 50, 50, 50, 50, 55, 55, 55, 55, 60, 60,
60, 60, 65, 65, 65, 65, 70, 70, 70, 70, 75, 75, 75, 75, 80, 80,
80, 80, 85, 85, 85, 85, 90, 90, 90, 90), total_cat = c(0, 1,
3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0,
1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10,
0, 1, 3, 10, 0, 1, 3, 10, 0, 1, 3, 10),
shape = c(1.04293463521838,
0.825156738737092, 0.899003822110993, 0.589358680965725, 1.05948319345651,
0.895913335495644, 0.84619233996788, 0.636706369588122, 1.16028493785731,
0.97919993225436, 0.68848363112076, 0.667566641635962, 1.25727757690268,
1.11021543939268, 0.717059585208267, 0.724574127291869, 1.3247305288869,
1.18510576695807, 0.858171614262325, 0.749998577524088, 1.3340371479645,
1.20792880820659, 0.920717960779785, 0.852052206426261, 1.32268361303857,
1.24450960892428, 0.988946105418831, 0.919376727690091, 1.31156083562523,
1.26972967056357, 1.03454490825271, 0.954911145157684, 1.27126976718353,
1.26045834296745, 1.06684366805728, 0.968088650561122, 1.18494141305856,
1.20593215449012, 1.05998296831204, 0.989315997334925, 1.12395322065035,
1.13733607943079, 1.04483735235382, 0.988385269358649, 1.1220716456468,
1.09128664878198, 1.01326552242001, 0.948730052972523, 1.08315738475919,
1.09407460180658, 1.05197207574662, 0.978397483932468),
cure = c(0.585284664126963,
0.545117838495959, 0.402804166722055, 0.14283955595853, 0.398430641565923,
0.327274451509188, 0.148468826128857, 0.00333583482942569, 0.257384245965907,
0.208165119313289, 0.157385208407254, 0.0696401481749125, 0.166635358789789,
0.133819199845029, 0.0610524867384188, 0.0535645284203467, 0.134003946557746,
0.114600739520109, 0.0537234127141967, 0.0308512461181499, 0.104553332825899,
0.0862550079726627, 0.0399447271614142, 0.0259518317972345, 0.10432761596883,
0.0823122881511086, 0.0355479853523502, 0.0190476453614591, 0.103738700348739,
0.0816458941897609, 0.0412696810847696, 0.0319927106467962, 0.108815314898247,
0.0867189417281315, 0.0472560235806836, 0.0300347697835532, 0.124252895869383,
0.111313062025094, 0.0675699701502124, 0.0373036736422192, 0.15590486866906,
0.153015363026187, 0.0974421665330981, 0.0563274705516883, 0.228204960789236,
0.2051853623161, 0.169506663277659, 0.0887344403643672, 0.334252927026757,
0.295865002768107, 0.24644310712191, 0.151125704749819),
scale = c(3.33440594470802,
3.56052439703143, 0.247924244305486, 0.150095610340916, 4.05027295845979,
3.78800058764034, 0.35102275649495, 0.486681740716857, 3.62361340129802,
3.48003996256557, 0.661105185961319, 0.247063080430436, 2.80272048406229,
2.6202749242671, 0.764421647692264, 0.268131101357167, 2.29273349177322,
2.11664063949949, 0.749125964202465, 0.417294230771058, 2.06795900869592,
1.80352742984593, 0.78399749774511, 0.451625737966194, 1.82650979495336,
1.62873839529095, 0.806784448546511, 0.478644431229979, 1.61927028533258,
1.46189762327418, 0.810126410765397, 0.500531738601558, 1.53243134315106,
1.42674067036488, 0.869134221688119, 0.569591312385404, 1.44388696689722,
1.41174931932889, 0.938154963786188, 0.631891327069348, 1.47052388447475,
1.4543783566239, 1.01084363142901, 0.644940081254986, 1.29365253055963,
1.4079170674224, 1.03720449704243, 0.643101192871478, 1.08541958643012,
1.29524623033074, 1.02176143186057, 0.673175882772333))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 136 rows and 2 variables:
#' \describe{
#' \item{\code{cohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{pop}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname Sweden
#' @export
"pop1"
pop1 <- data.frame(cohort=as.double(2035:1900),
pop=c(rep(17239,32), 16854, 16085, 15504, 15604, 16381, 16705,
16762, 16853, 15487, 14623, 14066, 13568, 13361, 13161, 13234,
13088, 12472, 12142, 12062, 12078, 11426, 12027, 11963, 12435,
12955, 13013, 13125, 13065, 12249, 11103, 9637, 9009, 8828,
8350, 7677, 7444, 7175, 6582, 6573, 6691, 6651, 6641, 6268,
6691, 6511, 6857, 7304, 7308, 7859, 7277, 8323, 8561, 7173,
6942, 7128, 6819, 5037, 6798, rep(6567,46)))
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 4 rows and 4 variables:
#' \describe{
#' \item{\code{cohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{pop}}{double COLUMN_DESCRIPTION}
#' \item{\code{screened}}{double Number screened at baseline}
#' \item{\code{pscreened}}{double Probability of being screened at baseline}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname cap
#' @export
"cap_study"
cap_study <- data.frame(cohort=c(1955,1950,1945,1940),
pop=c(55229, 55077,44057,35023),
screened=c(17869,19448,15783,11335))
cap_study <- transform(cap_study, pscreened=screened/pop)
#' @title DATASET_TITLE
#' @description DATASET_DESCRIPTION
#' @format A data frame with 4 rows and 4 variables:
#' \describe{
#' \item{\code{cohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{pop}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname cap
#' @export
"cap_control"
cap_control <- data.frame(cohort=c(1955,1950,1945,1940),
pop=c(63423,63285,51507,41224))
#' @title background_utilities
#' @description Swedish background utilities from Burstrom and Rehnberg (2006)
#' @format A data frame with 14 rows and 3 variables:
#' \describe{
#' \item{\code{lower}}{double lower age}
#' \item{\code{lower}}{double upper age}
#' \item{\code{utility}}{double background utility value}
#'}
#' @details http://libris.kb.se/bib/10708053?vw=short
#' @examples
#' \dontrun{
#' if(interactive()){
#' background_utilities
#' }
#' }
#' @rdname background_utilities
#' @export
"background_utilities"
background_utilities <-
data.frame(lower=c(0, 18, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80),
upper=c(18, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 1.0e55),
utility=c(1, 0.89, 0.89, 0.88, 0.87, 0.84, 0.84, 0.83, 0.83, 0.82, 0.83, 0.81,
0.79, 0.74))
#' @title Population for the STHLM3-MRI MRI arm
#' @description Birth cohorts for the study, assuming 2019.5 is the middle of the study
#' @format A data frame with 30 rows and 2 variables:
#' \describe{
#' \item{\code{cohort}}{birt cohort}
#' \item{\code{pop}}{population count}
#'}
#' @examples
#' \dontrun{
#' if(interactive()){
#' background_utilities
#' }
#' }
#' @rdname STHLM3
#' @export
"sthlm3_mri_arm"
sthlm3_mri_arm <- data.frame(cohort = 2019 - 45:74,
pop = c(1, 4, 6, 9, 10, 11, 21, 22, 20, 28, 23,
34, 44, 37, 48, 58, 59, 60, 71, 60, 71, 79, 76, 71, 79, 91, 95,
77, 83, 54))
#' @title List tables for the simulation
#' @description DATASET_DESCRIPTION
#' \strong{all_cause_mortality}
#' @format A data frame with 12120 rows and 4 variables:
#' \describe{
#' \item{\code{BirthCohort}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Age}}{integer COLUMN_DESCRIPTION}
#' \item{\code{AnnualMortality}}{double COLUMN_DESCRIPTION}
#' \item{\code{Survival}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{biopsy_frequency}
#' @format A data frame with 3 rows and 7 variables:
#' \describe{
#' \item{\code{PSA.beg}}{integer COLUMN_DESCRIPTION}
#' \item{\code{PSA.end}}{integer COLUMN_DESCRIPTION}
#' \item{\code{X55.59}}{double COLUMN_DESCRIPTION}
#' \item{\code{X60.64}}{double COLUMN_DESCRIPTION}
#' \item{\code{X65.69}}{double COLUMN_DESCRIPTION}
#' \item{\code{X70.74}}{double COLUMN_DESCRIPTION}
#' \item{\code{X75.79}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{biopsy_sensitivity}
#' @format A data frame with 14 rows and 2 variables:
#' \describe{
#' \item{\code{Year}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Sensitivity}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{neg_biopsy_to_psa}
#' @format A data frame with 4 rows and 3 variables. Describing the
#' time to the next PSA test following a negative biopsy. Modelled
#' as a competing risk with a biopsy. Informed by the Stockholm
#' PSA and Biopsy Register (SPBR).:
#' \describe{
#' \item{\code{age}}{integer with age groups}
#' \item{\code{meanlog}}{double for the mean of a log-normal time to
#' a PSA test following a negative biopsy.}
#' \item{\code{sdlog}}{double for the standard deviation of a
#' log-normal.}
#' }
#' \strong{neg_biopsy_to_biopsy}
#' @format A data frame with 4 rows and 3 variables. Describing the
#' time to the next biopsy following a negative biopsy. Modelled
#' as a competing risk with a PSA test. Informed by the Stockholm
#' PSA and Biopsy Register (SPBR).:
#' \describe{
#' \item{\code{age}}{integer with age groups}
#' \item{\code{meanlog}}{double describing the mean of a log-normal}
#' \item{\code{sdlog}}{double describing the standard deviation of a
#' log-normal}
#' }
#' \strong{cure_m_CM_to_RP}
#' @format A data frame with 4 rows and 4 variables. Cure model giving
#' the probability of not having a radical prostatectomy following
#' assignment to conservative management (active surveillance &
#' watchfull waiting). Log-normal distribution of the time to
#' radical prostatectomy for those that do have that. A future
#' extension would be to model active surveillance and watchfull
#' waiting separately:
#' \describe{
#' \item{\code{age}}{double with age groups}
#' \item{\code{pnever}}{double probability of not having a radical
#' prostatectomy}
#' \item{\code{meanlog}}{double mean of log-normal time to radical
#' prostatectomy}
#' \item{\code{sdlog}}{double standard deviation of log-normal time to
#' radical prostatectomy}
#' }
#' \strong{cure_m_CM_to_RT}
#' @format A data frame with 4 rows and 4 variables. Cure model giving
#' the probability of not having a radiation therapy following
#' assignment to conservative management (active surveillance &
#' watchfull waiting). Log-normal distribution of the time to
#' radiation therapy for those that do have that. A future
#' extension would be to model active surveillance and watchfull
#' waiting separately:
#' \describe{
#' \item{\code{age}}{double with age groups}
#' \item{\code{pnever}}{double probability of not having a radical
#' prostatectomy}
#' \item{\code{meanlog}}{double mean of log-normal time to
#' radiation therapy}
#' \item{\code{sdlog}}{double standard deviation of log-normal time to
#' radiation therapy}
#' }
#' @format A data frame with 4 rows and 3 variables. Describing the
#' time to the next biopsy following a negative biopsy. Modelled
#' as a competing risk with a PSA test. Informed by the Stockholm
#' PSA and Biopsy Register (SPBR).:
#' \describe{
#' \item{\code{age}}{integer with age groups}
#' \item{\code{meanlog}}{double describing the mean of a log-normal}
#' \item{\code{sdlog}}{double describing the standard deviation of a
#' log-normal}
#' }
#' \strong{dre}
#' @format A data frame with 4 rows and 4 variables:
#' \describe{
#' \item{\code{psa.low}}{integer COLUMN_DESCRIPTION}
#' \item{\code{psa.high}}{double COLUMN_DESCRIPTION}
#' \item{\code{sensitivity}}{double COLUMN_DESCRIPTION}
#' \item{\code{specificity}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{prob_grade7}
#' @format A data frame with 51 rows and 2 variables:
#' \describe{
#' \item{\code{slope}}{double COLUMN_DESCRIPTION}
#' \item{\code{Pr.grade.7.}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{pradt}
#' @format A data frame with 8208 rows and 6 variables:
#' \describe{
#' \item{\code{Tx}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Age}}{integer COLUMN_DESCRIPTION}
#' \item{\code{DxY}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{NoADT}}{double COLUMN_DESCRIPTION}
#' \item{\code{ADT}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{prtx}
#' @format A data frame with 2664 rows and 6 variables:
#' \describe{
#' \item{\code{Age}}{integer COLUMN_DESCRIPTION}
#' \item{\code{DxY}}{integer COLUMN_DESCRIPTION}
#' \item{\code{G}}{integer COLUMN_DESCRIPTION}
#' \item{\code{CM}}{double COLUMN_DESCRIPTION}
#' \item{\code{RP}}{double COLUMN_DESCRIPTION}
#' \item{\code{RT}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{survival_dist}
#' @format A data frame with 42 rows and 3 variables:
#' \describe{
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Time}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Survival}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{survival_local}
#' @format A data frame with 294 rows and 5 variables:
#' \describe{
#' \item{\code{AgeLow}}{integer COLUMN_DESCRIPTION}
#' \item{\code{AgeHigh}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Time}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Survival}}{double COLUMN_DESCRIPTION}
#'}
#' \strong{prob_earlystage}
#' @format A data frame with 84 rows and 6 variables:
#' \describe{
#' \item{\code{BeginAge}}{integer COLUMN_DESCRIPTION}
#' \item{\code{EndAge}}{integer COLUMN_DESCRIPTION}
#' \item{\code{BeginPSA}}{integer COLUMN_DESCRIPTION}
#' \item{\code{EndPSA}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Grade}}{integer COLUMN_DESCRIPTION}
#' \item{\code{Prob}}{double COLUMN_DESCRIPTION}
#'}
"fhcrcData"
fhcrcData$biopsyOpportunisticComplianceTable <- swedenOpportunisticBiopsyCompliance
fhcrcData$biopsyFormalComplianceTable <- swedenFormalBiopsyCompliance
fhcrcData$secularTrendTreatment2008OR <- secularTrendTreatment2008OR
## https://www.socialstyrelsen.se/Lists/Artikelkatalog/Attachments/20008/2015-12-26.pdf
fhcrcData$background_utilities <- background_utilities
## fhcrcData$uk_screen_uptake <-
## data.frame(age = as.double(44:76),
## H = c(0, 0.00279457395690957, 0.0111782958276383,
## 0.0251511656121861, 0.0447131833105531, 0.0698643489227393, 0.098027986598169,
## 0.126627420486267, 0.155662650587033, 0.185133676900467, 0.215040499426569,
## 0.247663838926632, 0.285284416161949, 0.32790223113252, 0.375517283838344,
## 0.428129574279422, 0.483488769981715, 0.539344538471183, 0.595696879747827,
## 0.652545793811647, 0.709891280662642, 0.769265269782334, 0.832199690652243,
## 0.89869454327237, 0.968749827642714, 1.04236554376328, 1.11776147575895,
## 1.19315740775462, 1.26855333975029, 1.34394927174596, 1.41934520374163,
## 1.42934520374163, 1.43934520374163))
#' @title UK rescreening rates
#' @description DATASET_DESCRIPTION
#' @format A data frame with 48 rows and 5 variables:
#' \describe{
#' \item{\code{age5}}{double COLUMN_DESCRIPTION}
#' \item{\code{total_cat}}{double COLUMN_DESCRIPTION}
#' \item{\code{shape}}{double COLUMN_DESCRIPTION}
#' \item{\code{cure}}{double COLUMN_DESCRIPTION}
#' \item{\code{scale}}{double COLUMN_DESCRIPTION}
#'}
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname uk_rescreening
#' @export
"uk_rescreening"
uk_rescreening <- data.frame(age5 = c(30, 30, 30, 30, 30, 30, 50, 50, 50, 50,
50, 50, 55, 55, 55, 55, 55, 55, 60, 60, 60, 60, 60, 60, 65, 65,
65, 65, 65, 65, 80, 80, 80, 80, 80, 80, 85, 85, 85, 85, 85, 85,
90, 90, 90, 90, 90, 90),
total_cat = c(0, 3, 4, 6, 10, 20, 0,
3, 4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3,
4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3, 4, 6, 10, 20, 0, 3, 4,
6, 10, 20),
shape = c(1.18510576695807, 0.858171614262325, 0.858171614262325,
0.858171614262325, 0.749998577524088, 0.749998577524088, 1.18510576695807,
0.858171614262325, 0.858171614262325, 0.858171614262325, 0.749998577524088,
0.749998577524088, 1.20792880820659, 0.920717960779785, 0.920717960779785,
0.920717960779785, 0.852052206426261, 0.852052206426261, 1.24450960892428,
0.988946105418831, 0.988946105418831, 0.988946105418831, 0.919376727690091,
0.919376727690091, 1.26972967056357, 1.03454490825271, 1.03454490825271,
1.03454490825271, 0.954911145157684, 0.954911145157684, 1.13733607943079,
1.04483735235382, 1.04483735235382, 1.04483735235382, 0.988385269358649,
0.988385269358649, 1.09128664878198, 1.01326552242001, 1.01326552242001,
1.01326552242001, 0.948730052972523, 0.948730052972523, 1.09407460180658,
1.05197207574662, 1.05197207574662, 1.05197207574662, 0.978397483932468,
0.978397483932468),
cure = c(0.114600739520109, 0.0537234127141967,
0.0537234127141967, 0.0537234127141967, 0.0308512461181499, 0.0308512461181499,
0.114600739520109, 0.0537234127141967, 0.0537234127141967, 0.0537234127141967,
0.0308512461181499, 0.0308512461181499, 0.0862550079726627, 0.0399447271614142,
0.0399447271614142, 0.0399447271614142, 0.0259518317972345, 0.0259518317972345,
0.0823122881511086, 0.0355479853523502, 0.0355479853523502, 0.0355479853523502,
0.0190476453614591, 0.0190476453614591, 0.0816458941897609, 0.0412696810847696,
0.0412696810847696, 0.0412696810847696, 0.0319927106467962, 0.0319927106467962,
0.153015363026187, 0.0974421665330981, 0.0974421665330981, 0.0974421665330981,
0.0563274705516883, 0.0563274705516883, 0.2051853623161, 0.169506663277659,
0.169506663277659, 0.169506663277659, 0.0887344403643672, 0.0887344403643672,
0.295865002768107, 0.24644310712191, 0.24644310712191, 0.24644310712191,
0.151125704749819, 0.151125704749819),
scale = c(6.58183283712968,
3.56486060274202, 0.96113775302889, 0.784828252112451, 0.937020402040802,
1.59551751310995, 6.58183283712968, 3.56486060274202, 0.96113775302889,
0.784828252112451, 0.937020402040802, 1.59551751310995, 5.47461832642227,
4.85185826844681, 1.19197420365362, 0.960420044179229, 1.12429005910723,
1.6320163836042, 4.5256676489079, 4.92695557504941, 1.37335505106699,
1.08828035047345, 1.16345156553308, 1.88292622499511, 3.88781952822644,
4.07180517697073, 1.814341631623, 1.10397958554216, 1.27491455185655,
2.0730535296673, 1.4543783566239, 1.01084363142901, 1.01084363142901,
1.01084363142901, 0.644940081254986, 0.644940081254986, 1.4079170674224,
1.03720449704243, 1.03720449704243, 1.03720449704243, 0.643101192871478,
0.643101192871478, 1.29524623033074, 1.02176143186057, 1.02176143186057,
1.02176143186057, 0.673175882772333, 0.673175882772333))
#' @title ageStandards
#' @description Age standardisation weights
#' @format A data frame with 18 rows and 3 variables:
#' \describe{
#' \item{\code{Age}}{integer Left-hand age for five-year age groups to 85+}
#' \item{\code{Sweden2000}}{double Sweden 2000 standard weights}
#' \item{\code{World}}{double World standard weights}
#'}
#' @details Age groups
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname ageStandards
#' @export
"ageStandards"
ageStandards <- data.frame(Age = cut(seq(0, 85, 5),
breaks = c(seq(0,85,5), Inf),
right = FALSE),
Sweden2000 = c(5.3, 6.9, 6.4, 5.7, 5.9, 6.7, 7.2, 6.9,
6.6, 6.6, 7.6, 6.3, 4.9, 4.3, 4.1, 3.9,
2.6, 2.3),
World = c(12.0, 10.0, 9.0, 9.0, 8.0, 8.0, 6.0, 6.0, 6.0,
6.0, 5.0, 4.0, 4.0, 3.0, 2.0, 1.0, 0.5, 0.5))
#' @title dre
#' @description Sensitivity and specificity for digital rectal examination (DRE)
#' @format A data frame with 6 rows and 3 variables:
#' \describe{
#' \item{\code{psa_low}}{double Left-hand interval for PSA}
#' \item{\code{se}}{double sensitivity}
#' \item{\code{sp}}{double specificity}
#'}
#' @details See https://doi.org/10.1093/jnci/90.23.1817, table 2.
#' @rdname dre
#' @export
dre = data.frame(psa_low=c(0:4,10),
psa_high=c(1:4,10,1e6),
dre=c(1702,3305,1314,734,1095,217),
appa=c(19,119,97,89,254,128), # predicted!!
bx=c(109,296,128,106,249,82),
pc=c(4,29,14,35,115,67))
dre = transform(dre, sensitivity=pc/appa, specificity=(dre-appa-(bx-pc))/(dre-appa))
.xtabsDFCheck <- function(df, dim_cols) {
stopifnot(is.data.frame(df))
stopifnot(all(dim_cols %in% names(df)))
stopifnot(prod(sapply(dim_cols, function(name) length(unique(df[[name]])))) == nrow(df))
}
## fix pradt: repeated values in 2004
fhcrcData$pradt <- unique(fhcrcData$pradt)
#' @title Define and run simulation
#' @description Function to specify and run the microsimulation. A large number
#' of simulated men, \code{n}, will cause the simulation to take longer time
#' but will reduce the stochastic variation particularly for rare events.
#' @param n Integer number of men to simulate. Default: 10
#' @param screen String with one of the following simulated screening scenarios:
#' \describe{
#' \item{\code{noScreening}}{no screening test, only diagnosis from symptoms}
#' \item{\code{randomScreen50to70}}{TBA}
#' \item{\code{twoYearlyScreen50to70}}{two-yearly screening from age 50 to 70}
#' \item{\code{fourYearlyScreen50to70}}{four-yearly screening from age 50 to 70}
#' \item{\code{screen50}}{one screen at age 50}
#' \item{\code{screen60}}{one screen at age 60}
#' \item{\code{screen70}}{one screen at age 70}
#' \item{\code{screenUptake}}{current testing pattern in Sweden}
#' \item{\code{stockholm3_goteborg}}{TBA}
#' \item{\code{stockholm3_risk_stratified}}{TBA}
#' \item{\code{goteborg}}{risk stratified re-screening 2+4 from age 50 to 70}
#' \item{\code{risk_stratified}}{risk stratified re-screening 4+8 from age 50 to 70}
#' \item{\code{mixed_screening}}{risk stratified re-screening 2+4 from age 50 to 70 & opportunistic testing for other ages}
#' \item{\code{regular_screen}}{TBA}
#' \item{\code{single_screen}}{TBA}
#' \item{\code{introduced_screening_only}}{TBA}
#' \item{\code{introduced_screening_preference}}{TBA}
#' \item{\code{introduced_screening}}{TBA}
#' \item{\code{stopped_screening}}{TBA}
#' } . Default: 'noScreening'
#'
#' @param nLifeHistories Integer with number of men for all events should be
#' recorded, Default: 10
#' @param seed Integer random number seed, Default: 12345
#' @param panel Boolean to use the Stockholm3 biomarker panel test
#' characteristics, otherwise PSA is used, Default: FALSE
#' @param flatPop Boolean to set all birth cohorts of equal size, Default: FALSE
#' @param pop Data.frame with two integer columns \code{cohort} with year of the
#' birth cohorts and \code{pop} with the size of the corresponding birth
#' cohorts, Default: pop1
#' @param tables List of data.frames to update the tables in fhcrcData, Default:
#' IHE (NB: fhcrcData$ptrx gets changed below to stockholmTreatment if parameters$stockholmTreatment = TRUE -- which is the default! This implies that IHE is the default when parameters$stockholmTreatment = FALSE)
#' @param debug Boolean to print debugging, Default: FALSE
#' @param parms List to update FhcrcParameters, Default: NULL
#' @param mc.cores Integer with the number of cores to use for the computation,
#' Default: 1
#' @param cl a cluster object, created by the \code{parallel} package or by the \code{snow}
#' package. If NULL, use \code{mclapply} (i.e. does not use the registered
#' default cluster).
#' @param print.timing Boolean should the required time be printed after the
#' simulation run, Default: TRUE
#' @param ... TBA
#' @return A fhcrc object
#' @details TBA
#' @examples
#' \dontrun{
#' if(interactive()){
#' library(prostata)
#' sim1 <- callFhcrc(n=1e6, screen="screenUptake", mc.cores=3)
#' }
#' }
#' @seealso
#' \code{\link[parallel]{mclapply}}
#' @rdname callFhcrc
#' @export
#' @importFrom parallel mclapply
callFhcrc <- function(n=10, screen= "noScreening", nLifeHistories=10,
seed=12345, panel=FALSE, flatPop = FALSE, pop = pop1,
tables = IHE, debug=FALSE, parms = NULL, mc.cores = 1,
cl = NULL,
print.timing = TRUE,...) {
## save the random number state for resetting later
state <- RNGstate(); on.exit(state$reset())
## yes, we use the user-defined RNG
RNGkind("user")
set.user.Random.seed(seed)
call <- match.call()
## birth cohorts that should give approximately the number of men alive in Stockholm in 2012
## check the input arguments
screenT <- c("noScreening", "randomScreen50to70",
"twoYearlyScreen50to70", "fourYearlyScreen50to70",
"screen50", "screen60", "screen70", "screenUptake",
"stockholm3_goteborg", "stockholm3_risk_stratified",
"goteborg", "risk_stratified", "mixed_screening",
"regular_screen", "single_screen",
"introduced_screening_only", "introduced_screening_preference",
"introduced_screening", "stopped_screening",
"cap_control", "cap_study", "sthlm3_mri_arm", "grs_stratified", "grs_stratified_age",
"germany_2018")
screen <- match.arg(screen, screenT)
stopifnot(is.na(n) || is.integer(as.integer(n)))
stopifnot(is.integer(as.integer(nLifeHistories)))
## these enum strings should be moved to C++
stateT <- c("Healthy","Localised","Metastatic")
ext_stateT <- c("Healthy","T1_T2","T3plus","Metastatic")
gradeT <- c("Gleason_le_6","Gleason_7","Gleason_ge_8","Healthy")
eventT <- c("toLocalised","toMetastatic","toClinicalDiagnosis", "toCancerDeath",
"toOtherDeath","toScreen","toBiopsyFollowUpScreen",
"toScreenInitiatedBiopsy","toClinicalDiagnosticBiopsy",
"toScreenDiagnosis", "toOverDiagnosis", "toOrganised","toTreatment",
"toCM","toRP", "toRT","toADT","toUtilityChange","toUtilityRemove",
"toSTHLM3", "toOpportunistic","toT3plus", "toCancelScreens",
"toYearlyActiveSurveillance", "toYearlyPostTxFollowUp", "toMRI",
"toPalliative", "toTerminal", "toDRE")
diagnosisT <- c("NotDiagnosed","ClinicalDiagnosis","ScreenDiagnosis")
treatmentT <- c("no_treatment","CM","RP","RT")
psaT <- c("PSA<3","PSA>=3") # not sure where to put this...
screenIndex <- which(screen == screenT) - 1
timingfunction <- if (print.timing) function(x) print(system.time(x)) else function(x) x
## NB: sample() calls the random number generator (!)
if (is.vector(pop)) {
flatPop <- TRUE
pop <- data.frame(cohort=pop,pop=1)
}
if (is.na(n)) {
cohort <- pop$cohort[rep.int(1:nrow(pop),times=pop$pop)]
n <- length(cohort)
} else {
if (flatPop) {
cohort <- rep(pop$cohort,each=ceiling(n/nrow(pop))) #Need ceiling so int n=!0
n <- ceiling(n/nrow(pop)) * nrow(pop) #To get the chunks right
} else
cohort <- sample(pop$cohort,n,prob=pop$pop/sum(pop$pop),replace=TRUE)
}
cohort <- sort(cohort)
## now separate the data into chunks and set the initial random numbers
if (!is.null(cl)) mc.cores <- length(cl) # HACK
currentSeed <- user.Random.seed()
if (mc.cores==1) {
chunks <- list(cohort)
ns <- 0
initialSeeds <- list(currentSeed)
} else {
ns <- floor((0:mc.cores)/mc.cores*n)
chunks <- lapply(1:mc.cores, function(i) cohort[(ns[i]+1):ns[i+1]])
initialSeeds <- c(list(currentSeed), lapply(ns[-c(1,mc.cores+1)], function(i) advance.substream(currentSeed, i)))
ns <- ns[-length(ns)]
}
## Minor changes to fhcrcData
fhcrcData$biopsyOpportunisticComplianceTable <- swedenOpportunisticBiopsyCompliance
fhcrcData$biopsyFormalComplianceTable <- swedenFormalBiopsyCompliance
fhcrcData$secularTrendTreatment2008OR <- secularTrendTreatment2008OR
fhcrcData$rescreening <- rescreening
fhcrcData$newdre = with(dre, data.frame(psa=psa_low, sensitivity, specificity))
if (!is.null(tables))
for (name in names(tables))
fhcrcData[[name]] <- tables[[name]]
updateParameters <- parms
updateParameters$nLifeHistories <- as.integer(nLifeHistories)
updateParameters$screen <- as.integer(screenIndex)
parameter <- FhcrcParameters
for (name in names(updateParameters)){
parameter[[name]] <- updateParameters[[name]]
}
parameter$g0 <- parameter$g0 / parameter$susceptible
parameter$cost_parameters <- parameter$currency_rate * parameter$cost_parameters
parameter$production <- data.frame(ages = parameter$production$ages,
values = parameter$currency_rate * parameter$production$values)
if (parameter$stockholmTreatment)
fhcrcData$prtx <- stockholmTreatment
if (any(.intersection <- (names(parameter) %in% names(fhcrcData))))
warning("The following tables are being over-written by values from parms: ",
paste0(names(parameter)[.intersection], collapse=", "))
temp <- parameter
parameter <- fhcrcData
for (name in names(temp))
parameter[[name]] <- temp[[name]]
## table checks
parameter$rescreening$total <- parameter$rescreening$total_cat
parameter$prtx$Age <- as.double(parameter$prtx$Age)
parameter$prtx$DxY <- as.double(parameter$prtx$DxY)
parameter$prtx$G <- parameter$prtx$G - 1L
parameter$pradt$Grade <- parameter$pradt$Grade - 1L
parameter$biopsy_sensitivity$Year <- as.double(parameter$biopsy_sensitivity$Year)
parameter$neg_biopsy_to_psa$age <- as.double(parameter$neg_biopsy_to_psa$age)
parameter$neg_biopsy_to_biopsy$age <- as.double(parameter$neg_biopsy_to_biopsy$age)
parameter$cure_m_CM_to_RP$age <- as.double(parameter$cure_m_CM_to_RP$age)
parameter$cure_m_CM_to_RT$age <- as.double(parameter$cure_m_CM_to_RT$age)
parameter$pradt$Age <- as.double(parameter$pradt$Age)
parameter$pradt$DxY <- as.double(parameter$pradt$DxY)
## parameter$biopsyComplianceTable <-
## data.frame(expand.grid(psa=c(4,7,10),age=seq(55,75,by=5)),
## compliance=unlist(parameter$biopsy_frequency[,-(1:2),]))
parameter$survival_local <-
with(parameter$survival_local,
data.frame(Age=as.double(AgeLow),Grade=Grade,Time=as.double(Time),
Survival=Survival))
parameter$survival_dist <-
with(parameter$survival_dist,
data.frame(Grade=Grade,Time=as.double(Time),
Survival=Survival))
.xtabsDFCheck(parameter$prtx,c("Age","DxY","G"))
.xtabsDFCheck(parameter$pradt,c("Tx","Age","DxY","Grade"))
.xtabsDFCheck(parameter$hr_locoregional,c("age","ext_grade","psa10")) # should this be in fhcrcData?
.xtabsDFCheck(parameter$biopsyOpportunisticComplianceTable, c("psa","age"))
.xtabsDFCheck(parameter$rescreening,c("age5","total_cat"))
.xtabsDFCheck(parameter$survival_dist,c("Grade","Time"))
.xtabsDFCheck(parameter$survival_local,c("Age","Grade","Time"))
## .xtabsDFCheck(parameter$dre,c("psa_low","sensitivity","specificity"))
## parameter <- modifyList(fhcrcData, parameter)
pind <- sapply(parameter,class)=="numeric" & sapply(parameter,length)==1
bInd <- sapply(parameter,class)=="logical" & sapply(parameter,length)==1
## check some parameters for sanity
if (panel && parameter$rTPF>1) stop("Panel: rTPF>1 (not currently implemented)")
if (panel && parameter$rFPF>1) stop("Panel: rFPF>1 (not currently implemented)")
if (parameter$Andreas && parameter$MRI_screen)
stop("For 'MRI_screen=TRUE', use 'parms=c(prostata:::ShuangParameters, MRI_screen=TRUE)'")
## if (panel && parameter$MRI_screen)
## stop("Scenarios for 'panel=TRUE' and 'MRI_screen=TRUE' have not been defined")
if (parameter$MRI_clinical && !parameter$MRI_screen)
stop("Scenarios for 'MRI_clinical=TRUE' and 'MRI_screen=FALSE' have not been defined")
## now run the chunks separately
step <- function(i) {
chunk <- chunks[[i]]
set.user.Random.seed(initialSeeds[[i]])
.Call("callFhcrc",
parms=list(n=as.integer(length(chunk)),
firstId=ns[i],
panel=panel, # bool
debug=debug, # bool
cohort=as.double(chunk),
parameter=unlist(parameter[pind]),
bparameter=unlist(parameter[bInd]),
otherParameters=parameter[!pind & !bInd]),
PACKAGE="prostata")
}
if (is.null(cl)) {
timingfunction(out <- parallel::mclapply(1:mc.cores, step, mc.cores=mc.cores))
} else {
clusterEvalQ(cl, {library(prostata); RNGkind("user")})
clusterExport(cl, c("chunks", "initialSeeds", "ns", "panel", "debug", "pind",
"bInd", "parameter"), envir=environment())
timingfunction(out <- parallel::parLapply(cl, 1:length(cl), step))
}
## Apologies: we now need to massage the chunks from C++
## reader <- function(obj) {
## out <- cbind(data.frame(state=enum(obj$state[[1]],stateT),
## dx=enum(obj$state[[2]],diagnosisT),
## psa=enum(obj$state[[3]],psaT),
## cohort=obj$state[[4]]),
## data.frame(obj[-1]))
## out$year <- out$cohort + out$age
## out
## }
ext_state2state <- function(obj) # collapses T-stages to localised
`levels<-`(factor(obj),list(Healthy="Healthy",Localised=list("T1_T2","T3plus"),Metastatic="Metastatic"))
cbindList <- function(obj) # recursive (not used)
if (is.list(obj)) do.call("cbind",lapply(obj,cbindList)) else data.frame(obj)
rbindList <- function(obj) # recursive (not used)
if (is.list(obj)) do.call("rbind",lapply(obj,rbindList)) else data.frame(obj)
rbindExtract <- function(obj,name)
do.call("rbind",lapply(obj, function(obji) data.frame(obji[[name]])))
reader <- function(obj) {
cbind(data.frame(state=ext_state2state(enum(obj[[1]],ext_stateT)),
ext_state=enum(obj[[1]],ext_stateT),
grade=enum(obj[[2]],gradeT),
dx=enum(obj[[3]],diagnosisT),
psa=enum(obj[[4]],psaT),
cohort=obj[[5]]),
data.frame(obj[,-(1:5)]))
}
## grab all of the pt, prev, ut, events from summary
## pt <- lapply(out, function(obj) obj$summary$pt)
if (length(out[[1]]$summary) > 0) {
summary <- lapply(names(out[[1]]$summary),
function(name) do.call(rbind,
lapply(out, function(obj) reader(obj$summary[[name]]))))
names(summary) <- names(out[[1]]$summary)
states <- c("state","ext_state","grade","dx","psa","cohort")
names(summary$prev) <- c(states,"age","count")
names(summary$pt) <- c(states,"age","pt")
names(summary$ut) <- c(states,"age","ut")
names(summary$events) <- c(states,"event","age","n")
if(FALSE) age <- NULL # To pass false-positive check note
summary <- lapply(summary,function(obj) within(obj,year <- cohort+age))
enum(summary$events$event) <- eventT
}
else if (length(out[[1]]$shortSummary) > 0) {
## use shortSummary
summary <- lapply(names(out[[1]]$shortSummary),
function(name) do.call(rbind,
lapply(out, function(obj) obj$shortSummary[[name]])))
names(summary) <- names(out[[1]]$shortSummary)
states <- "state"
names(summary$prev) <- c(states,"age","count")
names(summary$pt) <- c(states,"age","pt")
names(summary$ut) <- c(states,"age","ut")
names(summary$events) <- c(states,"event","age","n")
if(FALSE) age <- NULL # To pass false-positive check note
if ("cohort" %in% names(summary))
summary <- lapply(summary,function(obj) within(obj,year <- cohort+age))
enum(summary$events$event) <- eventT
}
else {
summary <- list()
}
## lifeHistories <- do.call("rbind",lapply(out,function(obj) data.frame(obj$lifeHistories)))
## psarecord <- do.call("rbind",lapply(out,function(obj) data.frame(obj$psarecord)))
## diagnoses <- do.call("rbind",lapply(out,function(obj) data.frame(obj$diagnoses)))
## falsePositives <- do.call("rbind",lapply(out,function(obj) data.frame(obj$falsePositives)))
## parameters <- do.call("rbind",lapply(out,function(obj) data.frame(obj$parameters)))
lifeHistories <- rbindExtract(out,"lifeHistories")
psarecord <- rbindExtract(out,"psarecord")
bxrecord <- rbindExtract(out,"bxrecord")
diagnoses <- rbindExtract(out,"diagnoses")
falsePositives <- rbindExtract(out,"falsePositives")
parameters <- rbindExtract(out,"parameters")
indiv_costs <- do.call(c, lapply(out, "[[", "indiv_costs"))
indiv_utilities <- do.call(c, lapply(out, "[[", "indiv_utilities"))
combine.Means <- function(df) {
n <- sum(df$n)
sum <- sum(df$sum)
sumsq <- sum(df$sumsq)
mean <- sum/n
var <- n/(n-1)*(sumsq/n-mean*mean)
sd <- sqrt(var)
se <- sd/sqrt(n)
data.frame(n,mean,var,sd,se,sum,sumsq)
}
mean_utilities <- combine.Means(df=rbindExtract(out, "mean_utilities"))
mean_costs <- combine.Means(df=rbindExtract(out, "mean_costs"))
appendMeans <- function(x) c(x,
mean.sum = x[["sum"]] / x[["n"]],
mean.sumsq = x[["sumsq"]] / x[["n"]])
natural.history.summary <- data.frame(tmc_minus_t0 = appendMeans(sapply(rbindExtract(out,"tmc_minus_t0"), sum)))
## Identifying elements without name which also need to be rbind:ed
societal.costs <- do.call("rbind",lapply(out,function(obj) data.frame(obj$costs))) #split in sociatal and healthcare perspective
## names(costs) <- c("type","item","cohort","age","costs")
names(societal.costs) <- c("type","item","age","costs")
societal.costs$type <- factor(ifelse(societal.costs$type,
"Productivity loss",
"Health sector cost")) # societal perspective
healthsector.costs <- societal.costs[societal.costs["type"] == "Health sector cost", c("item", "age", "costs")] # healthcare perspective
names(lifeHistories) <- c("id", "ext_state", "ext_grade", "dx", "event", "begin", "end", "year", "psa", "utility")
enum(lifeHistories$ext_state) <- ext_stateT
lifeHistories$state <- ext_state2state(lifeHistories$ext_state)
lifeHistories <- lifeHistories[c(names(lifeHistories)[1], "state", names(lifeHistories)[-1])] # shift col order
enum(lifeHistories$dx) <- diagnosisT
enum(lifeHistories$event) <- eventT
enum(diagnoses$ext_state) <- ext_stateT
diagnoses$state <- ext_state2state(diagnoses$ext_state)
enum(diagnoses$ext_grade) <- gradeT
enum(diagnoses$dx) <- diagnosisT
enum(diagnoses$tx) <- treatmentT
enum <- list(stateT = stateT, ext_stateT = ext_stateT, eventT = eventT, screenT = screenT,
diagnosisT = diagnosisT, psaT = psaT)
out <- list(n=n,screen=screen,enum=enum,lifeHistories=lifeHistories,
parameters=parameters, summary=summary,
healthsector.costs=healthsector.costs, societal.costs=societal.costs,
psarecord=psarecord, diagnoses=diagnoses, bxrecord=bxrecord,
cohort=data.frame(table(cohort)),simulation.parameters=parameter,
falsePositives=falsePositives, panel=panel, call = call,
natural.history.summary=natural.history.summary,
indiv_costs=indiv_costs,
indiv_utilities=indiv_utilities,
mean_utilities=mean_utilities,
mean_costs=mean_costs)
class(out) <- "fhcrc"
out
}
## R --slave -e "options(width=200); require(microsimulation); callFhcrc(100,nLifeHistories=1e5,screen=\"screen50\")[[\"parameters\"]]"
#' @title Summarise simulation results
#' @description FUNCTION_DESCRIPTION
#' @param object PARAM_DESCRIPTION
#' @param from Age from which to report
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname summary.fhcrc
#' @export
summary.fhcrc <- function(object, from=0, ...) {
if (from<0) {
warning("from argument should be non-negative - changed to 0")
from <- 0
}
if (from != round(from)) {
warning("from argument should be integer - changed to round(from)")
from <- round(from)
}
if (from>0) {
adjust.costs <- (1+object$simulation.parameters$discountRate.costs)^from
adjust.effectiveness <- (1+object$simulation.parameters$discountRate.effectiveness)^from
object$n <- sum(subset(object$summary$prev,age==from)$count)
object$summary$ut <- subset(object$summary$ut, age>=from)
object$summary$pt <- subset(object$summary$pt, age>=from)
object$summary$prev <- subset(object$summary$prev, age>=from)
object$summary$events <- subset(object$summary$events, age>=from)
object$societal.costs <- subset(object$societal.costs, age>=from)
object$healthsector.costs <- subset(object$healthsector.costs, age>=from)
} else {
adjust.costs <- adjust.effectiveness <- 1
}
newobj <- object[c("n","screen")]
with(object,
structure(.Data=c(newobj,
with(object, list(
from=from,
discountRate.costs = simulation.parameters$discountRate.costs,
discountRate.effectiveness = simulation.parameters$discountRate.effectiveness,
screening.tests = with(summary$events,
sum(n[event == "toScreen"])) / n,
biopsies = with(summary$events,
sum(n[event %in% c("toScreenInitiatedBiopsy",
"toClinicalDiagnosticBiopsy")])) / n,
biopsies.screen.initiated = with(summary$events,
sum(n[event %in% c("toScreenInitiatedBiopsy")])) / n,
biopsies.clinical.diagnostic = with(summary$events,
sum(n[event %in% c("toClinicalDiagnosticBiopsy")])) / n,
negative.biopsies = with(summary$events,
sum(n[event %in% c("toScreenInitiatedBiopsy",
"toClinicalDiagnosticBiopsy")])
- sum(n[event %in% c("toScreenDiagnosis",
"toClinicalDiagnosis")])) / n,
screen.diagnosis = with(summary$events,
sum(n[event == "toScreenDiagnosis"])) / n,
clinical.diagnosis = with(summary$events,
sum(n[event == "toClinicalDiagnosis"])) / n,
over.diagnosis = with(summary$events,
sum(n[event == "toOverDiagnosis"])) / n,
diagnosis = with(summary$events,
sum(n[event %in% c("toClinicalDiagnosis",
"toScreenDiagnosis")])) / n,
cancer.deaths = with(summary$events,
sum(n[event == "toCancerDeath"])) / n,
LE = sum(summary$pt$pt) / n,
QALE = sum(summary$ut$ut) / n * adjust.effectiveness,
healthsector.costs = sum(healthsector.costs$costs) / n * adjust.costs,
productivity.loss = tapply(societal.costs$costs,
societal.costs$type=="Productivity loss",
sum)[["TRUE"]] / n * adjust.costs,
societal.costs = sum(societal.costs$costs) / n * adjust.costs,
healthsector.cost.per.qaly = sum(healthsector.costs$costs) /
sum(summary$ut$ut),
productivity.loss.per.qaly = tapply(societal.costs$costs,
societal.costs$type=="Productivity loss",
sum)[["TRUE"]] /
sum(summary$ut$ut),
societal.cost.per.qaly = sum(societal.costs$costs) /
sum(summary$ut$ut)))),
class="summary.fhcrc"))
}
#' @title Minus operator method for \code{summary.fhcrc} objects
#' @description Calculates the difference between two \code{summary.fhcrc}
#' objects.
#' @param object1 The positive reference \code{summary.fhcrc} object.
#' @param object2 The negative \code{summary.fhcrc} object which is subtracted
#' from \code{object1}.
#' @return A \code{summary.fhcrc} object representing the difference between the
#' two compared \code{summary.fhcrc}.
#' @details Calculates the difference between two \code{summary.fhcrc} objects,
#' e.g. the difference in life expectancy, quality of life etc. The
#' exceptions are \code{discountRate.costs} and
#' \code{discountRate.effectiveness} which are required two be the same for
#' the two scenarios and the original value is returned. If \code{n} has the
#' same value the original value is returned if they differ a list with both
#' is returned.
#' @examples
#' \dontrun{
#' if(interactive()){
#' scenario1 <- summary(callFhcrc(screen = "screenUptake"))
#' scenario2 <-summary(callFhcrc(screen = "noScreening"))
#' scenario1 - scenario2
#' }
#' }
#' @rdname summary.fhcrc
#' @export
'-.summary.fhcrc' <- function(object1, object2) {
if(object1$from != object2$from)
stop("The two scenarios have different 'from' arguments.")
if(object1$discountRate.costs != object2$discountRate.costs)
stop("The two scenarios have different discount rates for the costs.")
if(object1$discountRate.effectiveness != object2$discountRate.effectiveness)
stop("The two scenarios have different discount rates for the effectiveness.")
if(object1$n != object2$n)
warning("The two scenarios have different number of men.")
structure(.Data= c(list(n = if(object1$n == object2$n) {object1$n} else {list(object1$n, object2$n)},
screen = sprintf("%s-%s", object1$screen, object2$screen),
discountRate.costs = object1$discountRate.costs,
discountRate.effectiveness = object1$discountRate.effectiveness),
## exempt special elements calc difference on all others
sapply(names(object1)[!names(object1) %in% c("n", "screen",
"discountRate.costs",
"discountRate.effectiveness")],
function(x) object1[[x]] - object2[[x]],
simplify = FALSE, USE.NAMES = TRUE)),
class="summary.fhcrc")
}
#' @title Multiplaction with scalar operator method for \code{summary.fhcrc} objects
#' @description Calculates the results per a number of persons of a
#' \code{summary.fhcrc} object.
#' @param obj The \code{summary.fhcrc} object.
#' @param perpersons a scalar to multiply the summary results with.
#' @return A \code{summary.fhcrc} object scaled for a number of persons.
#' @details N.b. this is unfortunatly not a commutative
#' operator. Therefor the first parameter must be the
#' \code{summary.fhcrc} object and the second parameter a nummeric
#' scalar.
#' @examples
#' \dontrun{
#' if(interactive()){
#' ## Example 1
#' scenario1 <- summary(callFhcrc(screen = "screenUptake"))
#' scenario1 * 1000 # the summery per thousand persons
#'
#' ## Example 2
#' scenario2 <-summary(callFhcrc(screen = "noScreening"))
#' (scenario1 - scenario2) * 1000 # the difference per thousand persons
#' }
#' }
#' @rdname summary.fhcrc
#' @export
'*.summary.fhcrc' <- function(obj, perpersons) {
if(!is.numeric(perpersons))
stop("The 'perpersons' variable needs to be a numeric.")
structure(.Data= c(list(n = ifelse(is.list(obj$n),
list(lapply(obj$n, function(x) x * perpersons)),
obj$n * perpersons),
screen = sprintf("(%s)*%d", obj$screen, perpersons),
discountRate.costs = obj$discountRate.costs,
discountRate.effectiveness = obj$discountRate.effectiveness),
## exempt special elements scale all others
lapply(obj[!names(obj) %in% c("n", "screen",
"discountRate.costs",
"discountRate.effectiveness")],
function(x) x * perpersons)),
class="summary.fhcrc")
}
#' @title Divide with scalar operator method for \code{summary.fhcrc} objects
#' @description Calculates the results divided by a denominator \code{summary.fhcrc} object.
#' @param obj The \code{summary.fhcrc} object.
#' @param denominator a scalar to divide the summary results with.
#' @return A \code{summary.fhcrc} object scaled by a denominator.
#' @details N.b. this is unfortunatly not a commutative
#' operator. Therefor the first parameter must be the
#' \code{summary.fhcrc} object and the second parameter a nummeric
#' scalar.
#' @examples
#' \dontrun{
#' if(interactive()){
#' }
#' }
#' @rdname summary.fhcrc
#' @export
'/.summary.fhcrc' <- function(obj, denominator) {
if(!is.numeric(denominator))
stop("The 'denominator' variable needs to be a numeric.")
structure(.Data= c(list(n = ifelse(is.list(obj$n),
list(lapply(obj$n, function(x) x / denominator)),
obj$n * denominator),
screen = sprintf("(%s)/%d", obj$screen, denominator),
discountRate.costs = obj$discountRate.costs,
discountRate.effectiveness = obj$discountRate.effectiveness),
## exempt special elements scale all others
lapply(obj[!names(obj) %in% c("n", "screen",
"discountRate.costs",
"discountRate.effectiveness")],
function(x) x / denominator)),
class="summary.fhcrc")
}
#' @title Print a selection of the summarised simulation
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname print.summary.fhcrc
#' @export
print.summary.fhcrc <- function(x, ...) {
obj <- x
cat(sprintf(
"Screening scenario: %s
Life expectancy: %f
Discounted QALE: %f
Discounted health sector costs: %f
Discounted societal costs: %f
Discounted rate (effect.): %f
Discounted rate (costs): %f
", obj$screen, obj$LE, obj$QALE,
obj$healthsector.costs, obj$societal.costs,
obj$discountRate.effectiveness,
obj$discountRate.costs))
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param object1 PARAM_DESCRIPTION
#' @param object2 PARAM_DESCRIPTION
#' @param perspective PARAM_DESCRIPTION, Default: c("societal.costs", "healthsector.costs")
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname ICER.fhcrc
#' @export
ICER.fhcrc <- function(object1, object2,
perspective = c("societal.costs",
"healthsector.costs"), from=0, ...) {
perspective <- match.arg(perspective)
p1 <- object1$simulation.parameters
p2 <- object2$simulation.parameters
stopifnot(p1$discountRate.costs == p2$discountRate.costs)
stopifnot(p1$discountRate.effectiveness == p2$discountRate.effectiveness)
summary1 <- summary(object1, from=from, ...)
summary2 <- summary(object2, from=from, ...)
out <- list(ICER.QALE=(summary1[[perspective]] - summary2[[perspective]]) /
(summary1$QALE - summary2$QALE),
delta.QALE=summary1$QALE - summary2$QALE,
delta.costs=summary1[[perspective]] - summary2[[perspective]],
from=from, perspective=perspective)
if (p1$discountRate.costs == 0 && p2$discountRate.costs == 0 &&
p1$discountRate.effectiveness == 0 && p2$discountRate.effectiveness == 0) {
out <- c(out,
list(ICER.LE = (summary1[[perspective]] - summary2[[perspective]]) /
(summary1$LE - summary2$LE),
delta.LE = summary1$LE - summary2$LE))
}
out
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname print.fhcrc
#' @export
print.fhcrc <- function(x, ...)
cat(sprintf("FHCRC prostate cancer model with %i individual(s) under scenario '%s'.\n",
x$n, x$screen),
...)
## TODO: better solve issue below with testing.rate for noScreening scenario
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param object PARAM_DESCRIPTION
#' @param scenarios PARAM_DESCRIPTION, Default: NULL
#' @param type PARAM_DESCRIPTION, Default: 'incidence.rate'
#' @param group PARAM_DESCRIPTION, Default: 'age'
#' @param age.breaks PARAM_DESCRIPTION, Default: NULL
#' @param year.breaks PARAM_DESCRIPTION, Default: NULL
#' @param cohort.breaks PARAM_DESCRIPTION, Default: NULL
#' @param age.weights PARAM_DESCRIPTION, Default: NULL
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @importFrom stats predict
#' @rdname predict.fhcrc
#' @export
predict.fhcrc <- function(object, scenarios=NULL, type = "incidence.rate",
group = "age", group.pt = group, age.breaks = NULL, year.breaks = NULL,
cohort.breaks = NULL, age.weights = NULL, ...) {
if(!inherits(object,"fhcrc")) stop("Expecting object to be an fhcrc object")
if(!(is.null(scenarios) || all(sapply(scenarios,inherits,"fhcrc")) || inherits(object,"fhcrc")))
stop("Expecting scenarios is NULL, a fhcrc object or a list of fhcrc objects")
if(is.null(scenarios) && (grepl(".?rate.?ratio$|.?rr$", type, ignore.case = TRUE)))
stop("More than one simulation object is needed to calculate a rate-ratio")
## Stripping of potential rate ratio option before matching
abbr_type <- match.arg(sub(".?rr$|.?rate.?ratio$",
"", type, ignore.case = TRUE),
c("incidence.rate", "symptomatic.incidence.rate",
"screen.incidence.rate", "overdiagnosis.rate",
"testing.rate", "biopsy.rate", "metastasis.rate",
"pc.mortality.rate", "allcause.mortality.rate",
"prevalence", "ly", "life.years"))
## Allowing for several groups
group <- match.arg(group,
c("state", "ext_state", "grade", "dx", "psa", "age", "year", "cohort"),
several.ok = TRUE)
event_types <- switch(abbr_type,
incidence.rate = c("toClinicalDiagnosis", "toScreenDiagnosis"),
symptomatic.incidence.rate = "toClinicalDiagnosis",
screen.incidence.rate = "toScreenDiagnosis",
overdiagnosis.rate = "toOverDiagnosis",
testing.rate = "toScreen",
biopsy.rate = c("toClinicalDiagnosticBiopsy", "toScreenInitiatedBiopsy"),
metastasis.rate = "toMetastatic",
pc.mortality.rate = "toCancerDeath",
allcause.mortality.rate = c("toCancerDeath", "toOtherDeath"))
interval2break <- function(interval_vector) {
as.numeric(unique(unlist(regmatches(interval_vector,
gregexpr("[0-9]+|Inf", interval_vector)))))
}
if(FALSE) {age <- year <- cohort <- Freq <- event <- Freq.y <- Freq.x <-
scenario.x <- rate.x <- rate.y <- NULL} # To pass false-positive check note
## Manipulate input, make sure age.breaks exists from age standardisation
if(!is.null(age.weights)) age.breaks <- interval2break(age.weights[,1])
## Manipulate input, make sure time group exists for specified time interval
for (scale in c("age", "year", "cohort")) {
if(!is.null(eval(parse(text = paste(scale, "breaks",sep="."))))) {
if(!(scale %in% group)) group <- c(group, scale)
if(!(scale %in% group.pt)) group.pt <- c(group.pt, scale)
}
}
## fast operations by group using base-R
## http://stackoverflow.com/questions/3685492/r-speeding-up-group-by-operations
grp_apply = function(XS, INDEX, FUN, ..., simplify=T) {
FUN = match.fun(FUN)
if (!is.list(XS))
XS = list(XS)
as.data.frame(as.table(tapply(1:length(XS[[1L]]), INDEX, function(s, ...)
do.call(FUN, c(lapply(XS, `[`, s), list(...))), ..., simplify=simplify)))
}
## Fixes colnames after group operation
name_grp <- function(x, grp = group) {names(x)[grep("^Var[0-9]+$", names(x))] <- grp; x}
## After a group operation the group will become a factor, if the time scale
## is not reported as an interval we convert it back to numeric. N.b. This
## needs to be called after name_grp().
numeric_time_scale <- function(df, group) {
within(df, {
if("age" %in% group && is.null(age.breaks)) age <- as.numeric(levels(age))[age] #important factor conversion
if("year" %in% group && is.null(year.breaks)) year <- as.numeric(levels(year))[year] #important factor conversion
if("cohort" %in% group && is.null(cohort.breaks)) cohort <- as.numeric(levels(cohort))[cohort] #important factor conversion
})
}
categorise_time <- function(df, age.breaks, year.breaks, cohort.breaks) {
cut_t <- function(time, time.breaks) {
## using first column or single vector as breaks
as.factor(cut(time, breaks = data.frame(time.breaks)[,1],
right = FALSE, dig.lab=4, ordered_result = TRUE))
}
if(!is.null(age.breaks)) df <- transform(df, age = cut_t(age, age.breaks))
if(!is.null(year.breaks)) df <- transform(df, year = cut_t(year, year.breaks))
if(!is.null(cohort.breaks)) df <- transform(df, cohort = cut_t(cohort, cohort.breaks))
df
}
## For standardising over a time scale e.g. ages
standardise_time <- function(df, timestr = "age", parstr = "rate",
time.weights = age.weights) {
if(is.null(time.weights)) {
return(df)
} else {
collapsed_grps <- c(group[!group == timestr], "scenario")
## Scale with weight
weight <- function(time, par, time.weights) {
par * sapply(t(time), {function(x) time.weights[x == time.weights[,1],2]})
}
## Sum over the groups (at least scenario)
collapse <- function(df, collapsed_grps, parstr) {
within(with(df, grp_apply(eval(parse(text = parstr)),
lapply(as.list(collapsed_grps),
{function(x) eval(parse(text = x))}),
sum, na.rm = TRUE)) ,{
assign(eval(substitute(parstr)),
Freq / sum(time.weights[,2]))
rm(Freq)
})
}
## Standardise and rename
return(numeric_time_scale(name_grp(collapse(
within(df,
{assign(eval(substitute(parstr)), weight(eval(parse(text =timestr)),
eval(parse(text = parstr)),
time.weights))}),
collapsed_grps, parstr), collapsed_grps), collapsed_grps))
}
}
## Calculates rates of specific events by specified groups
calc_rate <- function(object, event_types, group){
pt <- with(categorise_time(object$summary$pt, age.breaks, year.breaks,
cohort.breaks),
numeric_time_scale(name_grp(grp_apply(pt,
lapply(as.list(group.pt), function(x) eval(parse(text = x))),
sum), grp = group.pt), group.pt))
## TODO: if subset has no dim replace with zeros, temp fix below:
if(!any(object$summary$event$event %in% event_types)) {
stop(paste("The event(s)", paste(event_types, collapse = ", "),
"was not found in the", object$screen, "scenario"))
}
events <- with(subset(
categorise_time(object$summary$events, age.breaks, year.breaks,
cohort.breaks), event %in% event_types),
numeric_time_scale(name_grp(grp_apply(n,
lapply(as.list(group),
{function(x) eval(parse(text = x))}),
sum)), group))
within(merge(pt, events, by = group.pt, all.y = TRUE),{
rate <- ifelse(is.na(Freq.y) & !is.na(Freq.x), 0, Freq.y/Freq.x) #no events but some pt -> 0
n <- Freq.y
pt <- Freq.x
rm(Freq.x,Freq.y)})
}
## Calculate prevalences by specified groups
calc_prev <- function(object, group){
within(with(categorise_time(object$summary$prev, age.breaks,
year.breaks, cohort.breaks),
numeric_time_scale(name_grp(grp_apply(count,
lapply(as.list(group),
function(x) eval(parse(text = x))), sum)), group)), {
prevalence <- Freq/object$n
rm("Freq")})
}
## Calculate life-years by specified groups
calc_ly <- function(object, group){
within(with(categorise_time(object$summary$pt, age.breaks,
year.breaks, cohort.breaks),
numeric_time_scale(name_grp(grp_apply(pt,
lapply(as.list(group),
function(x) eval(parse(text = x))), sum)), group)), {
ly <- Freq/object$n
rm("Freq")})
}
## Calculates the outcome in the passed function for all
## simulation objects in the 'scenarios' list. Then the object
## outcome (e.g. rates or prev) are for the scenarios are added as
## rows and the scenario name as a column.
predict_scenarios <- function(scenarios, calc_outcome, ...) {
do.call(rbind, lapply(scenarios,
{function(object, ...)
cbind(calc_outcome(object, ...), scenario = object$screen)}, ...))
}
## Input checks allow for scenarios to be a single fhcrc object or
## list of fhcrc objects. Now make sure scenarios is a list.
if(inherits(scenarios, "fhcrc")) {scenarios <- list(scenarios)}
## Rate-ratio if type ends with rate.ratio or RR
if(grepl(".?rate.?ratio$|.?rr$", type, ignore.case = TRUE)){
scenario_rates <- standardise_time(predict_scenarios(unique(scenarios),
calc_rate,
event_types, group),
timestr = "age",
parstr = "rate",
time.weights = age.weights)
reference_rate <- standardise_time(predict_scenarios(list(object),
calc_rate,
event_types, group),
timestr = "age",
parstr = "rate",
time.weights = age.weights)
## Alt. check if it is part of the rates
if(!is.null(age.weights)) group <- group[!group == "age"]
within(merge(scenario_rates, reference_rate, by = group),{
scenario <- scenario.x
rate.ratio <- rate.x/rate.y
rate.ratio[!is.finite(rate.ratio)] <- NaN
rm(list=ls(pattern=".x$|.y$"))})
## Prevalence if type ends with rate.ratio or RR
} else if(grepl(".?prev$|.?prevalence$", type, ignore.case = TRUE)){
standardise_time(predict_scenarios(unique(c(list(object),scenarios)), calc_prev, group),
timestr = "age", parstr = "prevalence", time.weights = age.weights)
## Life-years if type is ly or life*years
} else if(grepl("^ly$|^life.?years$", type, ignore.case = TRUE)){
standardise_time(predict_scenarios(unique(c(list(object),scenarios)), calc_ly, group),
timestr = "age", parstr = "LY", time.weights = age.weights)
## Defaults to plain rates. If reference object exist add it to
## scenario list and remove duplicates.
} else {
standardise_time(predict_scenarios(unique(c(list(object),scenarios)), calc_rate, event_types, group),
timestr = "age", parstr = "rate", time.weights = age.weights)
}
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param type PARAM_DESCRIPTION, Default: c("incidence.rate", "testing.rate", "biopsy.rate", "metastasis.rate",
#' "pc.mortality.rate", "allcause.mortality.rate")
#' @param plot.type PARAM_DESCRIPTION, Default: 'l'
#' @param add PARAM_DESCRIPTION, Default: FALSE
#' @param xlab PARAM_DESCRIPTION, Default: 'Age (years)'
#' @param ylab PARAM_DESCRIPTION, Default: NULL
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @importFrom graphics plot
#' @rdname plot.fhcrc
#' @export
plot.fhcrc <- function(x, type=c("incidence.rate", "testing.rate",
"biopsy.rate", "metastasis.rate",
"pc.mortality.rate",
"allcause.mortality.rate"),
plot.type="l", add=FALSE, xlab="Age (years)",
ylab=NULL, ...) {
type <- match.arg(type)
if (is.null(ylab)) {ylab <- switch(type,
incidence.rate="Prostate cancer incidence rates per 100,000",
testing.rate="PSA rates per 1000",
biopsy.rate="Biopsies per 1000",
metastasis.rate="Metastatic onset per 100,000",
pc.mortality.rate="Cancer mortality rates per 100,000",
allcause.mortality.rate="All cause mortality rates per 100,000")}
rates <- predict(object = x, type = type)
rates$rate = rates$rate*switch(type, testing.rate=1000, biopsy.rate=1000, incidence.rate=1e5, metastasis.rate=1e5,pc.mortality.rate=1e5,allcause.mortality.rate=1e5)
if (!add) plot(rate~age, data=rates, type=plot.type, xlab=xlab, ylab=ylab, ...) else lines(rate~age, data=rates, ...)
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param x PARAM_DESCRIPTION
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname lines.fhcrc
#' @importFrom graphics lines
#' @export
lines.fhcrc <- function(x,...) {
plot(x, ..., add=TRUE)
}
#' @title FUNCTION_TITLE
#' @description FUNCTION_DESCRIPTION
#' @param object1 PARAM_DESCRIPTION
#' @param object2 PARAM_DESCRIPTION
#' @param startAge PARAM_DESCRIPTION, Default: 50
#' @param stopAge PARAM_DESCRIPTION, Default: Inf
#' @param ... PARAM_DESCRIPTION
#' @return OUTPUT_DESCRIPTION
#' @details DETAILS
#' @examples
#' \dontrun{
#' if(interactive()){
#' #EXAMPLE1
#' }
#' }
#' @rdname nn.fhcrc
#' @export
nn <- function(object1, object2, startAge = 50, stopAge = Inf, ...){
UseMethod("nn")
}
#' @export
nn.fhcrc <- function(object1, object2, startAge = 50, stopAge = Inf, ...) {
pNNS <- function(thisScenario) {
with(thisScenario$summary$events,
sum(n[event=="toCancerDeath" & age>=startAge
& age<stopAge])) / # divided by
with(thisScenario$summary$prev,
sum(count[age==round(startAge)]))
}
pNND <- function(thisScenario) {
with(thisScenario$summary$events,
sum(n[event=="toCancerDeath" & age>=startAge & age<stopAge])) / # divided by
with(thisScenario$summary$events,
sum(n[event %in% c("toScreenDiagnosis","toClinicalDiagnosis")
& age>=startAge & age<stopAge]))
}
NNS <- 1 / (pNNS(object2) - pNNS(object1)) #number needed to screen to prevent 1 PCa death
NND <- 1 / (pNND(object2) - pNND(object1)) #number needed to detect to prevent 1 PCa death
## Include additional number needed to treat (NNT) [Gulati 2011] to show overdiagnosis?
structure(.Data = list(NNS=NNS,NND=NND), class = "nn.fhcrc")
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.