Skip to content

Commit

Permalink
Try to tidy basics and test-basics
Browse files Browse the repository at this point in the history
  • Loading branch information
Dominic Muston committed Sep 9, 2023
1 parent 650e126 commit 44c7e06
Show file tree
Hide file tree
Showing 9 changed files with 337 additions and 289 deletions.
8 changes: 8 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,19 @@

export(calc_allrmds)
export(calc_allrmds_boot)
export(calc_haz)
export(calc_haz_par)
export(calc_haz_psm)
export(calc_ibs)
export(calc_likes)
export(calc_likes_psm_complex)
export(calc_likes_psm_simple)
export(calc_likes_stm_cf)
export(calc_likes_stm_cr)
export(calc_pdist)
export(calc_pdist_par)
export(calc_rmd)
export(calc_rmd_par)
export(calc_surv)
export(calc_surv_psmpps)
export(check_pfsos_consistent)
Expand All @@ -18,6 +24,8 @@ export(find_bestfit_par)
export(find_bestfit_spl)
export(fit_ends_mods_par)
export(fit_ends_mods_spl)
export(give_noparams)
export(give_noparams_par)
export(graph_survs)
export(prob_os_psm)
export(prob_os_stm_cf)
Expand Down
117 changes: 59 additions & 58 deletions R/basics.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,32 +36,32 @@
#' - Generalized Gamma requires the mu, sigma and Q parameters if using the standard parameterization (`gengamma`) or shape, scale and k parameters if using the original parameterization (`gengamma.orig`).
#' @seealso [flexsurv::flexsurvreg]
#' @return The value of the distribution function, a numeric value.
#' @export
#' @examples
#' calc_pdist(10, "exp", 0.01)
#' calc_pdist(5, "lnorm", c(3, 1))
calc_pdist_par <- function(time, dist, pars) {
if (dist=="exp") {
cdf <- stats::pexp(time, rate=pars[1])
stats::pexp(time, rate=pars[1])
} else if (dist=="weibullPH") {
cdf <- flexsurv::pweibullPH(time, shape=pars[1], scale=pars[2])
flexsurv::pweibullPH(time, shape=pars[1], scale=pars[2])
} else if (dist=="weibull" | dist=="weibull.quiet") {
cdf <- stats::pweibull(time, shape=pars[1], scale=pars[2])
stats::pweibull(time, shape=pars[1], scale=pars[2])
} else if (dist=="llogis") {
cdf <- flexsurv::pllogis(time, shape=pars[1], scale=pars[2])
flexsurv::pllogis(time, shape=pars[1], scale=pars[2])
} else if (dist=="lnorm") {
cdf <- stats::plnorm(time, meanlog=pars[1], sdlog=pars[2])
stats::plnorm(time, meanlog=pars[1], sdlog=pars[2])
} else if (dist=="gamma") {
cdf <- stats::pgamma(time, shape=pars[1], rate=pars[2])
stats::pgamma(time, shape=pars[1], rate=pars[2])
} else if (dist=="gompertz") {
cdf <- flexsurv::pgompertz(time, shape=pars[1], rate=pars[2])
flexsurv::pgompertz(time, shape=pars[1], rate=pars[2])
} else if (dist=="gengamma") {
cdf <- flexsurv::pgengamma(time, mu=pars[1], sigma=pars[2], Q=pars[3])
flexsurv::pgengamma(time, mu=pars[1], sigma=pars[2], Q=pars[3])
} else if (dist=="gengamma.orig") {
cdf <- flexsurv::pgengamma.orig(time, shape=pars[1], scale=pars[2], k=pars[3])
flexsurv::pgengamma.orig(time, shape=pars[1], scale=pars[2], k=pars[3])
} else {
cdf <- NA
NA
}
cdf
}

#' Calculate the value of the distribution function
Expand All @@ -84,6 +84,7 @@ calc_pdist_par <- function(time, dist, pars) {
#' ** scale - "hazard", "odds", or "normal", as described in [flexsurv::flexsurvspline]. With the default of no knots in addition to the boundaries, this model reduces to the Weibull, log-logistic and log-normal respectively. The scale must be common to all times.
#' @seealso [flexsurv::flexsurvspline] and [flexsurv::flexsurvreg]
#' @return The value of the distribution function, a numeric value.
#' @export
#' @examples
#' calc_pdist(time=1:5,
#' type="spl",
Expand All @@ -96,7 +97,7 @@ calc_pdist_par <- function(time, dist, pars) {
calc_pdist <- function(time, type, spec){
type <- tolower(substr(type, 1, 3))
if (type=="spl") {
cdf <- flexsurv::psurvspline(q = time,
flexsurv::psurvspline(q = time,
gamma = spec$gamma,
beta = 0,
X = 0,
Expand All @@ -107,11 +108,10 @@ calc_pdist <- function(time, type, spec){
lower.tail = TRUE,
log.p = FALSE)
} else if (type=="par") {
cdf <- calc_pdist_par(time, dist=spec$dist, pars=spec$pars)
calc_pdist_par(time, dist=spec$dist, pars=spec$pars)
} else {
cdf <- NA
NA
}
cdf
}

#' Calculate the value of the survival function
Expand All @@ -138,38 +138,39 @@ calc_surv <- function(time, type, spec) {
#' @inheritParams calc_pdist_par
#' @inherit calc_pdist_par seealso
#' @return The value of a hazard function, a numeric value.
#' @export
#' @examples
#' calc_haz_par(10, "exp", 0.01)
#' calc_haz_par(5, "lnorm", c(3, 1))
calc_haz_par <- function(time, dist, pars) {
if (dist=="exp") {
haz <- flexsurv::hexp(time, rate=pars[1])
flexsurv::hexp(time, rate=pars[1])
} else if (dist=="weibullPH") {
haz <- flexsurv::hweibullPH(time, shape=pars[1], scale=pars[2])
flexsurv::hweibullPH(time, shape=pars[1], scale=pars[2])
} else if (dist=="weibull" | dist=="weibull.quiet") {
haz <- flexsurv::hweibull(time, shape=pars[1], scale=pars[2])
flexsurv::hweibull(time, shape=pars[1], scale=pars[2])
} else if (dist=="llogis") {
haz <- flexsurv::hllogis(time, shape=pars[1], scale=pars[2])
flexsurv::hllogis(time, shape=pars[1], scale=pars[2])
} else if (dist=="lnorm") {
haz <- flexsurv::hlnorm(time, meanlog=pars[1], sdlog=pars[2])
flexsurv::hlnorm(time, meanlog=pars[1], sdlog=pars[2])
} else if(dist=="gamma") {
haz <- flexsurv::hgamma(time, shape=pars[1], rate=pars[2])
flexsurv::hgamma(time, shape=pars[1], rate=pars[2])
} else if(dist=="gompertz") {
haz <- flexsurv::hgompertz(time, shape=pars[1], rate=pars[2])
flexsurv::hgompertz(time, shape=pars[1], rate=pars[2])
} else if(dist=="gengamma") {
haz <- flexsurv::hgengamma(time, mu=pars[1], sigma=pars[2], Q=pars[3])
flexsurv::hgengamma(time, mu=pars[1], sigma=pars[2], Q=pars[3])
} else if(dist=="gengamma.orig") {
haz <- flexsurv::hgengamma.orig(time, shape=pars[1], scale=pars[2], k=pars[3])
flexsurv::hgengamma.orig(time, shape=pars[1], scale=pars[2], k=pars[3])
} else {
haz <- NA
NA
}
haz
}

#' Calculate the value of the hazard function
#' @description Calculate the value of the hazard function, given specification as either parametric or Royston-Parmar splines model
#' @inheritParams calc_surv
#' @inherit calc_haz_par return
#' @export
#' @examples
#' calc_haz(time=1:5,
#' type="spl",
Expand All @@ -182,7 +183,7 @@ calc_haz_par <- function(time, dist, pars) {
calc_haz <- function(time, type, spec){
type <- tolower(substr(type, 1, 3))
if (type=="spl") {
haz <- flexsurv::hsurvspline(x = time,
flexsurv::hsurvspline(x = time,
gamma = spec$gamma,
beta = 0,
X = 0,
Expand All @@ -191,11 +192,10 @@ calc_haz <- function(time, type, spec){
timescale = "log",
offset = 0)
} else if (type=="par") {
haz <- calc_haz_par(time, dist=spec$dist, pars=spec$pars)
calc_haz_par(time, dist=spec$dist, pars=spec$pars)
} else {
haz <- NA
NA
}
haz
}

#' Calculate restricted mean durations (parametric form)
Expand All @@ -205,32 +205,32 @@ calc_haz <- function(time, type, spec){
#' @param Tw is the time period over which the restricted mean is calculated
#' @inherit calc_pdist seealso
#' @return the restricted mean duration, a numeric value.
#' @export
#' @examples
#' calc_rmd_par("exp", 0.2, 20)
#' calc_rmd_par("lnorm", c(3, 1), 10)
calc_rmd_par <- function(Tw, dist, pars) {
if(dist=="exp") {
rmean <- flexsurv::rmst_exp(Tw, rate=pars[1], start=0)
flexsurv::rmst_exp(Tw, rate=pars[1], start=0)
} else if (dist=="weibullPH") {
rmean <- flexsurv::rmst_weibullPH(Tw, shape=pars[1], scale=pars[2], start=0)
flexsurv::rmst_weibullPH(Tw, shape=pars[1], scale=pars[2], start=0)
} else if (dist=="weibull" | dist=="weibull.quiet") {
rmean <- flexsurv::rmst_weibull(Tw, shape=pars[1], scale=pars[2], start=0)
flexsurv::rmst_weibull(Tw, shape=pars[1], scale=pars[2], start=0)
} else if (dist=="llogis") {
rmean <- flexsurv::rmst_llogis(Tw, shape=pars[1], scale=pars[2], start=0)
flexsurv::rmst_llogis(Tw, shape=pars[1], scale=pars[2], start=0)
} else if (dist=="lnorm") {
rmean <- flexsurv::rmst_lnorm(Tw, meanlog=pars[1], sdlog=pars[2], start=0)
flexsurv::rmst_lnorm(Tw, meanlog=pars[1], sdlog=pars[2], start=0)
} else if (dist=="gamma") {
rmean <- flexsurv::rmst_gamma(Tw, shape=pars[1], rate=pars[2], start=0)
flexsurv::rmst_gamma(Tw, shape=pars[1], rate=pars[2], start=0)
} else if (dist=="gompertz") {
rmean <- flexsurv::rmst_gompertz(Tw, shape=pars[1], rate=pars[2], start=0)
flexsurv::rmst_gompertz(Tw, shape=pars[1], rate=pars[2], start=0)
} else if (dist=="gengamma") {
rmean <- flexsurv::rmst_gengamma(Tw, mu=pars[1], sigma=pars[2], Q=pars[3], start=0)
flexsurv::rmst_gengamma(Tw, mu=pars[1], sigma=pars[2], Q=pars[3], start=0)
} else if (dist=="gengamma.orig") {
rmean <- flexsurv::rmst_gengamma.orig(Tw, shape=pars[1], scale=pars[2], k=pars[3], start=0)
flexsurv::rmst_gengamma.orig(Tw, shape=pars[1], scale=pars[2], k=pars[3], start=0)
} else {
rmean <- NA
NA
}
rmean
}

#' Calculate restricted mean durations
Expand All @@ -253,6 +253,7 @@ calc_rmd_par <- function(Tw, dist, pars) {
#' ** scale - "hazard", "odds", or "normal", as described in [flexsurv::flexsurvspline]. With the default of no knots in addition to the boundaries, this model reduces to the Weibull, log-logistic and log-normal respectively. The scale must be common to all times.
#' @inherit calc_haz_par seealso
#' @inherit calc_haz_par return
#' @export
#' @examples
#' calc_rmd(Tw=200,
#' type="spl",
Expand All @@ -265,7 +266,7 @@ calc_rmd_par <- function(Tw, dist, pars) {
calc_rmd <- function(Tw, type, spec){
type <- tolower(substr(type, 1, 3))
if (type=="spl") {
rmd <- flexsurv::rmst_survspline(t = Tw,
flexsurv::rmst_survspline(t = Tw,
gamma = spec$gamma,
beta = 0,
X = 0,
Expand All @@ -275,11 +276,10 @@ calc_rmd <- function(Tw, type, spec){
offset = 0
)
} else if (type=="par") {
rmd <- calc_rmd_par(Tw=Tw, dist=spec$dist, pars=spec$pars)
calc_rmd_par(Tw=Tw, dist=spec$dist, pars=spec$pars)
} else {
rmd <- NA
NA
}
rmd
}

#' Number of parameters used by parametric statistical distributions
Expand All @@ -288,38 +288,39 @@ calc_rmd <- function(Tw, type, spec){
#' e.g. "llogis" is the log-logistic distribution.
#' @return a numeric value
#' @inherit calc_pdist seealso
#' @export
#' @examples
#' give_noparams("llogis")
#' give_noparams(c("exp", "gengamma"))
give_noparams_par <- function(dist) {
knowndists <- c("exp", "weibull", "weibullPH", "llogis",
"lnorm", "gamma", "gompertz", "gengamma", "gengamma.orig")
if (length(dist)!=1) stop("Multiple distributions entered, may call only one.")
knownflag <- length(intersect(knowndists, dist))>0
if (knownflag==FALSE) return(NA) else {
params <- if (dist=="exp") 1 else if (dist=="gengamma" | dist=="gengamma.orig") 3 else 2
as.numeric(params)
}
dplyr::case_when(
dist=="exp" ~ 1,
dist=="weibull" | dist=="weibullPH" | dist=="llogis" | dist=="lnorm" ~ 2,
dist=="gamma" | dist=="gompertz" ~ 2,
dist=="gengamma" | dist=="gengamma.orig" ~ 3,
.default = NA
)
}

#' Number of parameters used by parametric statistical distributions
#' @description Returns the number of parameters used by the specified statistical model
#' @inheritParams calc_surv
#' @return a numeric value
#' @inherit calc_pdist seealso
#' @export
#' @examples
#' give_noparams(type="par", spec=list(dist="weibullPH"))
#' give_noparams(type="spl", spec=list(gamma=c(1.1,2.1,3.1)))
give_noparams <- function(type, spec) {
type <- tolower(substr(type, 1, 3))
type <- base::tolower(base::substr(type, 1, 3))
if (type=="spl") {
np <- length(spec$gamma)
length(spec$gamma)
}
else if (type=="par") {
np <- give_noparams_par(dist=spec$dist)
give_noparams_par(dist=spec$dist)
}
else {
np <- NA
NA
}
np
}
2 changes: 1 addition & 1 deletion R/fitting-spl.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ fit_ends_mods_spl <- function(ds,
) {
# Derive additional fields, as with regular function
ds <- create_extrafields(ds, cuttime=0)
dspps <- ds |> filter(pps.durn>0)
dspps <- ds |> dplyr::filter(pps.durn>0)
# Captures lists of fitted models to each endpoint
fits.ppd <- fit_mods_spl(durn1 = ds$tzero,
durn2 = ds$ppd.durn,
Expand Down
2 changes: 1 addition & 1 deletion R/fitting.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,7 +240,7 @@ fit_ends_mods_par <- function(ds,
expvar = NA) {
# Derive additional fields, as with regular function
ds <- create_extrafields(ds, cuttime)
dspps <- ds |> filter(pps.durn>0)
dspps <- ds |> dplyr::filter(pps.durn>0)
# Captures lists of fitted models to each endpoint
fits.ppd <- fit_mods_par(durn1 = ds$tzero,
durn2 = ds$ppd.durn,
Expand Down
Loading

0 comments on commit 44c7e06

Please sign in to comment.