#+TITLE: Prostata package for R

#+OPTIONS: toc:nil #+OPTIONS: num:nil #+OPTIONS: html-postamble:nil

Babel settings

+PROPERTY: cache yes

+PROPERTY: results output graphics

+PROPERTY: exports both

+PROPERTY: tangle yes

+PROPERTY: exports both

[[http://www.gnu.org/licenses/gpl-3.0.html][http://img.shields.io/:license-gpl3-blue.svg]]

• 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.

• Installing prostata

• 1 Dependencies: :: A convenient, but not required, way of installing github-packages in R is to use [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]]. Since [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]] is available on [[https://cran.r-project.org/][CRAN]] just run the following in R. #+BEGIN_SRC R :session r-prostata-readme :exports code :eval never install.packages("devtools") #+END_SRC

• 2a Installation with devtools: :: To install the microsimulation using [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]] just run the following in R: #+BEGIN_SRC R :session r-prostata-readme :exports code :eval never require(devtools) install_github("mclements/prostata") #+END_SRC

• 2b Alternative installation from shell: ::

Some thing OS-specific?

If you prefer the shell over [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]], just run the following to download the microsimulation R-package: #+BEGIN_SRC shell :exports code :eval never git clone https://github.com/mclements/prostata.git #+END_SRC

To install the prostata R-package run this in your shell: #+BEGIN_SRC shell :exports code :eval never R CMD INSTALL path_to_prostata #+END_SRC

• Running the simulation

#+HEADERS: :var reRunSimulation = 0 #+BEGIN_SRC R :session r-prostata-readme :exports none require(prostata) myFile <- file.path("~/Dropbox/microsimulation_runs","README_sim.RData")

if (reRunSimulation || !file.exists(myFile)){ sim1 <- callFhcrc(n=1e6, mc.cores=3, screen="screenUptake") sim2 <- callFhcrc(n=1e6, mc.cores=3, screen="noScreening") save(sim1, sim2, file=myFile) } else { load(file=myFile) } #+END_SRC

** Available screening scenarios There are a number of available testing scenarios. They determine testing frequencies and re-testing intervals over calendar period and ages.

• =noScreening= - no screening test, only diagnosis from symptoms
• =twoYearlyScreen50to70= - two-yearly screening from age 50 to 70
• =fourYearlyScreen50to70= - four-yearly screening from age 50 to 70
• =screen50= - one screen at age 50
• =screen60= - one screen at age 60
• =screen70= - one screen at age 70
• =screenUptake= - current testing pattern in Sweden
• =goteborg= - risk stratified re-screening 2+4 from age 50 to 70
• =risk_stratified= - risk stratified re-screening 4+8 from age 50 to 70
• =mixed_screening= - risk stratified re-screening 2+4 from age 50 to 70 & opportunistic testing for other ages

+ =randomScreen50to70=

+ =stockholm3_goteborg=

+ =stockholm3_risk_stratified=

+ =regular_screen=

+ =single_screen=

#+name: commentify #+begin_src emacs-lisp :var result="" :exports none (concat "#> "(mapconcat 'identity (split-string result "\n") "\n#> ")) #+end_src

#+BEGIN_SRC R :session r-prostata-readme :post commentify(this) :results output :exports both :eval never-export require(prostata) sim1 <- callFhcrc(n=1e6, mc.cores=3, screen="screenUptake") #+END_SRC

#+RESULTS: : #> user system elapsed : #> 144.180 0.180 91.173

• Results ** Type of outcome: prevalence, event rates or rate ratios Some of the more commonly used outcomes are provided through a =plot= and a =predict= function. The available outcomes are:

• =prevalence= - proportion of a population in the groups described below
• =incidence.rate= - rate of /clinical diagnosis/ & /screen initiated diagnosis/
• =symptomatic.incidence.rate= - rate of /clinical diagnosis/
• =screen.incidence.rate= - rate of /screen initiated diagnosis/
• =overdiagnosis.rate= - rate of diagnosed men who would have died from other causes before having symptoms
• =testing.rate= - rate of /screening tests/ (e.g. psa tests)
• =biopsy.rate= - rate of /clinical diagnostic biopsies/ & /screen initiated biopsies/
• =metastasis.rate= - rate of natural history transitions to /metastatic/ cancer
• =pc.mortality.rate= - rate of /cancer deaths/
• =allcause.mortality.rate= - rate of /cancer deaths/ & /other deaths/
• =incidence.rr= - rate ratio of /clinical diagnosis/ & /screen initiated diagnosis/
• =testing.rr= - rate ratio of /screening tests/ (e.g. psa tests)
• =biopsy.rr= - rate ratio of /clinical diagnostic biopsies/ & /screen initiated biopsies/
• =metastasis.rr= - rate ratio of natural history transitions to /metastatic/ cancer
• =pc.mortality.rr= - rate ratio of /cancer deaths/
• =allcause.mortality.rr= - rate ratio of /cancer deaths/ & /other deaths/ To construct an outcome not listed above you can use the objects in ~sim1$summary~ to construct them.

To simply plot e.g. the /incidence rate/ of the simulated screening scenario the following line can be used: #+BEGIN_SRC R :session r-prostata-readme :file inst/inc.png :results output graphics :exports both plot(sim1, type = "incidence.rate", xlab="Age (years)", xlim=c(40, 90)) #+END_SRC

#+RESULTS: [[file:inst/inc.png]]

** Groups in natural and clinical history The =predict= function returns various outcomes (/rate/, /rate ratios/ or /prevalence/) as described above. It can also be used to predict outcomes by a number of subgroups. The available subgroups are two time-scales and four natural history categories:

• =age= - grouping by single /year of age/ this is the default time-scale
• =year= - grouping by single /calendar year/ as an alternative time-scale
• =cohort= - grouping by single year /birth cohort/ as an alternative time-scale
• =state= - grouping by /healthy/, /localised/ & /metastatic/
• =ext_state= - grouping by /healthy/, /T1-T2-stages/, /T3+-stages/ and /metastatic/
• =grade= - grouping by /gleason grade/ ~<=6~, ~7~ & ~>=8~
• =dx= - grouping by /not diagnosed/, /screen diagnosis/ & /clinical diagnosis/
• =psa= - grouping by psa ~<3~ & ~>=3~

To use a time increment other then single year a specific interval can be supplied in addition to the corresponding time-scale group argument:

• =age.breaks= a vector of interval breaks e.g. =c(seq(50, 85, 5), Inf)=
• =year.breaks= a vector of interval breaks e.g. =seq(1990, 2020, 10)=
• =cohort.breaks= a vector of interval breaks e.g. =seq(1900, 2000, 20)=

The predict function also support age-standardisation through the =age.weights= argument. As a starting point the prostata package supplies a data.frame, =ageStandards=, with age weights for Sweden the year 2000 and the Standard World Population. User defined weight on the same format is also supported.

Below is the PSA testing rate by calendar period predicted and displayed with =ggplot=. In addition we here also give an example on how to age-standardise with the =age.weights= argument. #+BEGIN_SRC R :session r-prostata-readme :file inst/psa.png :results output graphics :exports both require(ggplot2) ggplot(predict(sim1, group = "year", type = "testing.rate", age.weights = ageStandards[c("Age", "Sweden2000")]), aes(x=year, y=rate)) + geom_line() + xlim(1990, 2015) + ylab("PSA testing rate") + xlab("Calendar period (years)") #+END_SRC

#+RESULTS: [[file:inst/psa.png]]

The outcomes can also be predicted by several subgroups at once. Plotted below is the prevalence by age, clinical state and diagnoses. Note that since this is a natural history of disease model also the unobserved /not diagnosed/ cancers are predicted. #+BEGIN_SRC R :session r-prostata-readme :file inst/prev.png :results output graphics :exports both ggplot(predict(sim1, type = "prevalence", group=c("age", "state", "dx")), aes(x=age, y=prevalence*1e5, colour = dx)) + geom_line() + ylab("Prevalence (cases per 100,000)") + xlab("Age (years)") + facet_grid(. ~ state) #+END_SRC

#+RESULTS: [[file:inst/prev.png]]

** Comparing multiple scenarios In order to compare multiple screening scenarios the =predict= function has a second argument for simulation objects. It can be used to pass a second simulation objects if you which to compare two screening scenarios or a list of simulation objects for comparing several screening scenarios. The =type= and =group= argument works as described earlier. Below is the incidence rate with the /current uptake/ pattern compared with the hypothetical /no screening/ scenario.

#+BEGIN_SRC R :session r-prostata-readme :post commentify(this) :results output :exports both :eval never-export sim2 <- callFhcrc(n=1e6, mc.cores=3, screen="noScreening") #+END_SRC

#+RESULTS: : #> user system elapsed : #> 106.032 0.660 54.287

#+BEGIN_SRC R :session r-prostata-readme :file inst/scen.png :results output graphics :exports both ggplot(predict(sim1, sim2, group= "age", type = "incidence.rate"), aes(x=age, y=rate, colour = scenario)) + geom_line() + xlim(50, 85) + ylab("Incidence rate") + xlab("Age (years)") #+END_SRC

#+RESULTS: [[file:inst/scen.png]]

** Rate ratios and reference scenarios If