## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(fssg)
library(flexsurv)

## -----------------------------------------------------------------------------
# 'Density' function
dgamgomp(x=1, b=1, sigma=1, beta=1)

# 'Distribution' function
pgamgomp(q=1, b=1, sigma=1, beta=1)


# Constructing our distribution object
fssg_gamgomp <- fssg_dist(
    name='gamgomp',                       # a simple name
    pars=c('b','sigma','beta'),           # parameters are scale, shape, shape respectively
    location='b',                         # name of the parameter we want to vary 
    transforms=c(log,log,log),            # transformation functions for each of our parameters (all must be >0)
    inv.transforms=c(exp,exp,exp),        # back-transformation functions for each of our parameters
    inits=function(t){c(1/max(t), 1, 1)}, # function that constructs the initial parameter estimates for optimization
    d = dgamgomp,                         # density function
    p = pgamgomp,                         # distribution function.
    q = quantilify(pgamgomp),             # quantile function (optional) 
    h = hazardify(dgamgomp, pgamgomp),    # hazard function (optional)
    H = cumhazardify(pgamgomp),           # cumulative hazard function (optional)
    fullname='gamma_gompertz'             # secondary name(s) which is used in the back-end of `fssg` to better label outputs
  )

## ----message=F----------------------------------------------------------------
# Sample data set provided in the package
data('pseudo', package='fssg')  

flexsurvreg(
  formula = Surv(time, death) ~ 1, 
  data = pseudo, 
  dist = fssg_gamgomp, 
  dfns = list(fssg_gamgomp$d, fssg_gamgomp$p)
) -> f1

print(f1)
plot(f1)

## -----------------------------------------------------------------------------
get_fssg_dist('gamma_gompertz')

## -----------------------------------------------------------------------------
fssg(
  formula = Surv(time, death) ~ 1, 
  data = pseudo, 
  model='gamma_gompertz'
) -> model1

model1$models$gamma_gompertz

## -----------------------------------------------------------------------------
fssg(
  survival::Surv(time, status) ~ x, 
  data = survival::aml
) -> aml_model

# Summary Table
aml_model$summary

# List of Models
head(aml_model$models)

## -----------------------------------------------------------------------------
fssg(
  survival::Surv(time, death) ~ gender, 
  data = pseudo, 
  models = c('gamma_gompertz','gompertz','gumbel','gumbel_T2')
)$summary

## ----message=F----------------------------------------------------------------
check_inits(times = pseudo$time,  get_fssg_dist('gumbel'))
check_inits(times = pseudo$time,  get_fssg_dist('gumbel_T2'))
check_inits(times = pseudo$time,  get_fssg_dist('lomax'))

## -----------------------------------------------------------------------------
fssg(survival::Surv(time, status) ~ age + sex, 
  data = survival::cancer, 
  models = c('gamma_gompertz','gompertz','gumbel','gumbel_T2'), progress = F) -> output

output$models$gumbel

## -----------------------------------------------------------------------------
plot(
  output$models$gumbel, 
  newdata = data.frame(sex=c(1,2), age=70), col=c('blue','red'),
  type='survival',
  lwd=2,
)
legend(
  'topright',
  legend = c('Male','Female'),
  col =c('blue','red'),
  lwd=2
)

## -----------------------------------------------------------------------------
plot(
  output$models$gumbel, 
  newdata = data.frame(sex=c(1), age=c(65,80)), col=c('blue','red'),
  type='survival',
  lwd=2
)
legend(
  'topright',
  legend = c('Age 65','Age 80'),
  col =c('blue','red'),
  lwd=2
)

## -----------------------------------------------------------------------------
fssg(
  survival::Surv(time, status) ~ x, 
  data = survival::aml, 
  models = c('weibull'), 
  spline = c('rp','wy'), # include both methods
  max_knots = 3          # attempt up to 3 knots
) -> spline_models

spline_models$summary

# Best model has 1 knot
# plot(spline_models$models$spline_wy_normal_1)

# Model with 3 knots
# plot(spline_models$models$spline_wy_normal_3)

## -----------------------------------------------------------------------------
# Example using our models on the aml dataset
head(aml_model$summary)

## -----------------------------------------------------------------------------
get_fit_stats(aml_model$models$inv_lind, ibs=T) -> fitstats

fitstats