Gianluca Baio
Department of Statistical Science | University College London
https://gianluca.statistica.it
https://egon.stats.ucl.ac.uk/research/statistics-health-economics
https://github.com/giabaio https://github.com/StatisticsHealthEconomics
@gianlubaio@mas.to @gianlubaio
EsCHER’s R-evolution: Mastering R Package Development for Health Economics
Erasmus University Rotterdam
23 April 2024
Check out our departmental podcast “Random Talks” on Soundcloud!
Follow our departmental social media accounts + magazine “Sample Space”
R packages
R packages in statistics and beyond:: and ::: operators in R: Usage, importance and differencespackrat and/or renv: functionality, strengths, and best usage scenariosR packages from sources like CRAN, GitHub, and local/zip files, ensuring a wide accessibility(Kind of stolen from here)
S
Fortran libraries to do internal “statistical analysis”C (v. 3)S to StatSci (now Insightful Corp)
S-Plus, building in a bunch of fancy featuresIn the meantime…
R using syntax that is very similar to S
R is a spinoff from S — but there are some fundamental distinctions among them…R free softwareR-help and R-devel) + R-Core Group formed
RRSAS)…R installation comes with a basic core of functions and functionalities (base)
R package is a collection of functions that are bundled together in a way that lets them be easily sharedhttps://r-pkgs.org/introduction.html
Californian band that was big for like 20 minutes in the very early 2000s
2014 Canadian crime thriller film adapted from the 2008 novel of the same name by Michael Redhill
## Runs the C/E analysis, given matrices of costs and effects from BUGS model
## Differential effects and costs
delta.e <- e[,2] - e[,1]
delta.c <- c[,2] - c[,1]
## Plots the C/E plane
m.e <- range(delta.e)[1]
M.e <- range(delta.e)[2]
m.c <- range(delta.c)[1]
M.c <- range(delta.c)[2]
step <- (M.e-m.e)/10
m.e <- ifelse(m.e<0,m.e,-m.e)
m.c <- ifelse(m.c<0,m.c,-m.c)
x.pt <- .95*m.e
y.pt <- ifelse(x.pt*wtp<m.c,m.c,x.pt*wtp)
xx <- seq(100*m.c/wtp,100*M.c/wtp,step)
yy <- xx*wtp
xx[1] <- ifelse(min(xx)<m.e,xx[1],2*m.e)
yy[1] <- ifelse(min(yy)<m.c,yy[1],2*m.c)
xx[length(xx)] <- ifelse(xx[length(xx)]<M.e,1.5*M.e,xx[length(xx)])
if(!is.null(xlim)) {m.e <- xlim[1]; M.e <- xlim[2]}
if(!is.null(ylim)) {m.c <- ylim[1]; M.c <- ylim[2]}
plot(xx,yy,col="white",xlim=c(m.e,M.e),ylim=c(m.c,M.c),xlab=xlab,ylab=ylab,main=title,axes=F)
...
That’s when you know you need to turn it all into a package!
https://r4ds.had.co.nz/functions.html#when-should-you-write-a-function
devtools & roxygen✔ Creating 'testPackage/'
✔ Setting active project to '/home/gianluca/erasmus/testPackage'
✔ Creating 'R/'
✔ Writing 'DESCRIPTION'
Package: testPackage
Title: What the Package Does (One Line, Title Case)
Version: 0.0.0.9000
Authors@R (parsed):
* First Last <first.last@example.com> [aut, cre] (YOUR-ORCID-ID)
Description: What the package does (one paragraph).
License: `use_mit_license()`, `use_gpl3_license()` or friends to
pick a license
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.3.1
✔ Writing 'NAMESPACE'
✔ Writing 'testPackage.Rproj'
✔ Adding '^testPackage\\.Rproj$' to '.Rbuildignore'
✔ Adding '.Rproj.user' to '.gitignore'
✔ Adding '^\\.Rproj\\.user$' to '.Rbuildignore'
✔ Opening '/home/gianluca/erasmus/testPackage/' in new RStudio session
✔ Setting active project to '<no active project>'testPackage
DESCRIPTIONNAMESPACER folder:: vs ::: vs library()…NAMESPACE etc
Functions in a package can be accessed in different ways (and giving/using different levels of access)
NAMESPACEdevtools can do a lot of the heavy-lifting for you…Example: BCEA NAMESPACE file
# Generated by roxygen2: do not edit by hand
S3method("CEriskav<-",bcea)
S3method("CEriskav<-",default)
S3method("mixedAn<-",bcea)
S3method("mixedAn<-",default)
S3method("setComparisons<-",bcea)
S3method("setComparisons<-",default)
S3method("setKmax<-",bcea)
S3method("setKmax<-",default)
S3method("setReferenceGroup<-",bcea)
S3method("setReferenceGroup<-",default)
S3method(bcea,bugs)
S3method(bcea,default)
S3method(bcea,rjags)
S3method(bcea,rstan)
S3method(ceac.plot,bcea)
S3method(ceaf.plot,pairwise)
S3method(ceef.plot,bcea)
S3method(ceplane.plot,bcea)
S3method(contour,bcea)
S3method(contour2,bcea)
S3method(createInputs,bugs)
S3method(createInputs,data.frame)
S3method(createInputs,default)
S3method(createInputs,numeric)
S3method(createInputs,rjags)
S3method(createInputs,stanfit)
S3method(eib.plot,bcea)
S3method(evi.plot,bcea)
S3method(evi.plot,mixedAn)
S3method(evppi,bcea)
S3method(evppi,default)
S3method(ib.plot,bcea)
S3method(info.rank,bcea)
S3method(line_labels,default)
S3method(line_labels,pairwise)
S3method(multi.ce,bcea)
S3method(plot,CEriskav)
S3method(plot,bcea)
S3method(plot,evppi)
S3method(plot,mixedAn)
S3method(print,bcea)
S3method(sim_table,bcea)
S3method(sim_table,default)
S3method(summary,bcea)
S3method(summary,mixedAn)
S3method(summary,pairwise)
export("CEriskav<-")
export("mixedAn<-")
export("setComparisons<-")
export("setKmax<-")
export("setReferenceGroup<-")
export(bcea)
export(best_interv_given_k)
export(ce_table)
export(ceac.plot)
export(ceaf.plot)
export(ceef.plot)
export(ceplane.plot)
export(compute_CEAC)
export(compute_EIB)
export(compute_EVI)
export(compute_IB)
export(compute_ICER)
export(compute_U)
export(compute_Ustar)
export(compute_kstar)
export(compute_ol)
export(compute_vi)
export(contour)
export(contour2)
export(createInputs)
export(diag.evppi)
export(eib.plot)
export(evi.plot)
export(evppi)
export(ib.plot)
export(info.rank)
export(info_rank_plotly)
export(is.bcea)
export(line_labels)
export(make.report)
export(make_legend_plotly)
export(mce.plot)
export(multi.ce)
export(new_bcea)
export(prep_ceplane_params)
export(setComparisons)
export(sim_table)
export(struct.psa)
export(tabulate_means)
export(validate_bcea)
import(dplyr)
import(ggplot2)
import(grid)
import(reshape2)
importFrom(MASS,kde2d)
importFrom(MCMCvis,MCMCchains)
importFrom(Matrix,Matrix)
importFrom(Matrix,expm)
importFrom(Rdpack,reprompt)
importFrom(cli,cli_alert_warning)
importFrom(dplyr,arrange)
importFrom(dplyr,desc)
importFrom(dplyr,filter)
importFrom(dplyr,mutate)
importFrom(dplyr,slice)
importFrom(grDevices,colors)
importFrom(grDevices,colours)
importFrom(grDevices,dev.new)
importFrom(grDevices,devAskNewPage)
importFrom(grDevices,grey.colors)
importFrom(grDevices,pdf.options)
importFrom(grDevices,ps.options)
importFrom(graphics,abline)
importFrom(graphics,axis)
importFrom(graphics,barplot)
importFrom(graphics,box)
importFrom(graphics,contour)
importFrom(graphics,legend)
importFrom(graphics,lines)
importFrom(graphics,matlines)
importFrom(graphics,matplot)
importFrom(graphics,par)
importFrom(graphics,points)
importFrom(graphics,polygon)
importFrom(graphics,rect)
importFrom(graphics,text)
importFrom(grid,unit)
importFrom(gridExtra,grid.arrange)
importFrom(purrr,keep)
importFrom(purrr,list_cbind)
importFrom(purrr,map)
importFrom(purrr,map_int)
importFrom(purrr,map_lgl)
importFrom(purrr,pluck)
importFrom(reshape2,melt)
importFrom(rlang,.data)
importFrom(rlang,int)
importFrom(rstan,extract)
importFrom(scales,label_dollar)
importFrom(stats,as.formula)
importFrom(stats,cov)
importFrom(stats,density)
importFrom(stats,lm)
importFrom(stats,qnorm)
importFrom(stats,qqline)
importFrom(stats,qqnorm)
importFrom(stats,quantile)
importFrom(stats,reorder)
importFrom(stats,rnorm)
importFrom(stats,runif)
importFrom(stats,sd)
importFrom(stats,setNames)
importFrom(utils,methods)
importFrom(utils,modifyList)
importFrom(utils,str)
importFrom(utils,tail)
Functions that are exported are directly available when you load a package onto your workspace, using the library(package) command
So you can do something like
and bcea is available for you to use…
OR, even without loading the package (assuming it’s installed in your R installation), you can still access the exposed functions using the “::” notation
DESCRIPTION file (to add your details and more info on the project/package)
DESCRIPTION file (to add your details and more info on the project/package)R folder, start writing your functions is suitable .R files#' Helper function to make the Kaplan-Meier analysis of the underlying data
#' for a given formula and dataset
#'
#' @param formula a formula specifying the model to be used, in the form
#' \code{Surv(time,event)~treatment[+covariates]} in flexsurv terms, or
#' \code{inla.surv(time,event)~treatment[+covariates]} in INLA terms.
#' @param method A string specifying the inferential method (\code{'mle'},
#' \code{'inla'} or \code{'hmc'}). If \code{method} is set to \code{'hmc'},
#' then \code{survHE} will write suitable model code in the Stan language
#' (according to the specified distribution), prepare data and initial values
#' and then run the model.
#' @param data A data frame containing the data to be used for the analysis.
#' This must contain data for the 'event' variable. In case there is no
#' censoring, then \code{event} is a column of 1s.
#' @return \item{ObjSurvfit}{A 'rms::npsurv' estimate of the KM curves}.
#' @author Gianluca Baio
#' @seealso fit.models
#' @template refs
#' @keywords Kaplan-Meier estimate
#' @noRd
#'
make_KM <- function(formula,data) {
km.formula <- as.formula(gsub("inla.surv","Surv",deparse(formula)))
# Computes the Kaplan Meier curve using the package "rms"
ObjSurvfit <- rms::npsurv( # Uses the function "npsurv" from the package "rms"
formula = km.formula, # to fit the model specified in the "formula" object
data = data # to the dataset named "data"
)
return(ObjSurvfit)
}In the file man-roxygen/refs.R…
NB: You can also add a .bib file with all the relevant references you want to use (in bibtex format) in the folder /inst/REFERENCES.bib and then in your documentation
for instance…
devtools and roxygendevtools and roxygen
devtools and roxygen
devtools and roxygen
R and packages is to deal with dependencies
boot is installed on your computer (and in general the user’s computer)!Warning
library(xxx) within a package — you would fail syntax controls when submitting to CRAN (more on this later!):: instead
devtools::install_dependencies to try and automate the installation of all the required packages to run yours…R Archive Network
devtools et al have made the background preparation much easierCRAN first time around and you need to fix several quirks…There are other ways of sharing your package and making it available with others
There are other ways of sharing your package and making it available with others
CRAN…
package_version_2_2024_04_22_final_new_GB.R file diaspora…R package folder a .git repositorygit syntax and philosophy (some entry barrier…)
GitHub allows you to respond very quickly to requests for changes/bug alerts etc
main is aligned with CRAN but devel is one (or several…) step(s) aheaddevel to main and re-submit to CRAN…CRAN on GitHubCRAN have their own official GitHub repo: https://github.com/cran/
R2OpenBUGS vs R2jags
.Rmd/.qmd file with the code/text of the vignette and store it in the ./vignettes folder in your packagepkgdown this is parsed and suitable html vignettes are created for your users!---
title: "Risk Aversion Analysis"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Risk Aversion Analysis}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
editor_options:
chunk_output_type: console
---
```{r}
#| include: false
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 6
)
```
```{r}
#| label: setup-example
#| results: 'hide'
#| message: false
#| warning: false
#| echo: false
library(BCEA)
library(dplyr)
library(reshape2)
library(ggplot2)
library(purrr)
```
Set-up analysis using smoking cessation data set.
```{r}
data(Smoking)
treats <- c("No intervention", "Self-help", "Individual counselling", "Group counselling")
bcea_smoke <- bcea(eff, cost, ref = 4, interventions = treats, Kmax = 500)
```
Run the risk aversion analysis straight away with both the base R and ggplot2 versions of plots.
```{r}
r <- c(0, 0.005, 0.020, 0.035)
CEriskav(bcea_smoke) <- r
plot(bcea_smoke)
plot(bcea_smoke, graph = "ggplot")
```
Notice that the first value is asymptotically zero but the function handles that for us.
Previously, you had to use something like 1e-10 which was a bit awkward.
Now we modify the comparison group so that it doesn't contain 2 ("self-help") anymore.
We still keep the same first comparison though "no intervention".
...R-HTA consortium© Gianluca Baio (UCL) | | R Packages in research | EsCHER workshop | 23 Apr 2024 | Slides available at https://giabaio.github.io/erasmus