From bb0349eb830cfacb5819c4a953f48b20a9c29646 Mon Sep 17 00:00:00 2001 From: Gianluca Baio Date: Fri, 8 May 2026 10:39:46 +0100 Subject: [PATCH 1/2] Implements the Gamma model --- DESCRIPTION | 3 +-- NEWS.md | 8 +++++++- R/runINLA.R | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++ prompt | 1 + 4 files changed, 59 insertions(+), 3 deletions(-) create mode 100644 prompt diff --git a/DESCRIPTION b/DESCRIPTION index 72b44db..6cecc90 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,6 @@ Package: survHEinla Title: Survival Analysis in Health Economic Evaluation using INLA -Version: 0.0.3 -Date: 2025-07-09 +Version: 0.0.4 Authors@R: c( person(given = "Gianluca",family = "Baio",role = c("aut", "cre"),email = "g.baio@ucl.ac.uk")) URL: https://github.com/giabaio/survHEinla, https://gianluca.statistica.it/software/survhe diff --git a/NEWS.md b/NEWS.md index 5d3a5a2..5b17db7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ -# survHEinla 0.0.1 +# survHEinla 0.0.4 (Mar 2026) +* Implements the Gamma survival model.Adds a function `lik_gam_inla()` to allow post-processing and likelihood computation. +All other changes are actually taken care of directly from the `INLA` infrastructure and facilities in `survHE`. + +# survHEinla 0.0.2 (Mar 2021) +* Implements the Gompertz survival model +# survHEinla 0.0.1 * Added a `NEWS.md` file to track changes to the package. diff --git a/R/runINLA.R b/R/runINLA.R index 0979128..632af91 100644 --- a/R/runINLA.R +++ b/R/runINLA.R @@ -337,3 +337,53 @@ lik_gom_inla <- function(model,nsim,time_max,t,d,formula,data) { npars <- 2 + (length(all.vars(formula))-1) list(logf=logf,logf.hat=logf.hat,npars=npars,f=NULL,f.bar=NULL,s=NULL,s.bar=NULL) } + +lik_gam_inla <- function(model,nsim,time_max,t,d,formula,data) { + # Gamma survival model (family="gammasurv" in INLA). + # Parameterisation: E[t_scaled] = exp(eta), shape = phi (precision hyperpar). + # Rate on scaled time: rate_scaled = phi / exp(eta) + # Back-transform intercept: +log(time_max) (same direction as lno, because + # the link is log(mean) and t_scaled = t/time_max => eta_true = eta_inla + log(time_max)) + # Rate on original scale: phi / exp(eta_true) + # Hazard: h(t) = dgamma(t, shape, rate) / pgamma(t, shape, rate, lower.tail=FALSE) + + # Sample shape (phi) from its marginal posterior + phi <- INLA::inla.rmarginal(nsim, model$marginals.hyperpar[[1]]) + # Sample fixed effects + beta <- lapply(1:nrow(model$summary.fixed), function(i) { + INLA::inla.rmarginal(nsim, model$marginals.fixed[[i]]) + }) + names(beta) <- colnames(model$model.matrix) + beta <- beta %>% bind_cols() + # Rescale intercept from [0,1] time scale to original time scale + if (grep("(Intercept)", colnames(model$model.matrix)) > 0) { + beta[, grep("(Intercept)", colnames(model$model.matrix))] <- + beta[, grep("(Intercept)", colnames(model$model.matrix))] + log(time_max) + } + phi.hat <- mean(phi) + beta.hat <- beta %>% summarise_all(mean) %>% as.numeric() + linpred <- as.matrix(beta) %*% t(model.matrix(formula, data)) + linpred.hat <- beta.hat %*% t(model.matrix(formula, data)) + + # Log-likelihood: d_i * log h(t_i) + log S(t_i) + # where h = dgamma/(1-pgamma) and log S = pgamma(..., lower.tail=FALSE, log.p=TRUE) + logf <- matrix( + unlist(lapply(1:nrow(linpred), function(i) { + rate_i <- phi[i] / exp(linpred[i, ]) + d * log(dgamma(t, shape = phi[i], rate = rate_i) / + pgamma(t, shape = phi[i], rate = rate_i, lower.tail = FALSE)) + + pgamma(t, shape = phi[i], rate = rate_i, lower.tail = FALSE, log.p = TRUE) + })), + nrow = nrow(linpred), byrow = TRUE) + logf.hat <- matrix( + { + rate.hat <- phi.hat / exp(linpred.hat) + d * log(dgamma(t, shape = phi.hat, rate = rate.hat) / + pgamma(t, shape = phi.hat, rate = rate.hat, lower.tail = FALSE)) + + pgamma(t, shape = phi.hat, rate = rate.hat, lower.tail = FALSE, log.p = TRUE) + }, + nrow = 1) + # Number of parameters: phi (shape) + fixed effects (intercept + covariates) + npars <- 1 + nrow(model$summary.fixed) + list(logf = logf, logf.hat = logf.hat, npars = npars, f = NULL, f.bar = NULL) +} diff --git a/prompt b/prompt new file mode 100644 index 0000000..609d568 --- /dev/null +++ b/prompt @@ -0,0 +1 @@ +No. I should have said, but survHEinla works in conjunction with survHE (https://github.com/giabaio/survHE), which has specific functions to handle the table of output and then all the plots. You need to modify some functions in survHE (survHE::like load_availables() and survHE:::manipulate_distributions() ) to ensure that the new Gamma distribution is available. You also need to create a new function lik_gam_inla() to go under runINLA.R to compute the log-density for the Gamma model. You can use the lik_***_inla() functions already present to see what you need to do. From 2e31eb5b52f94865ec3c03f66fa83f7373f50ccb Mon Sep 17 00:00:00 2001 From: Gianluca Baio Date: Fri, 8 May 2026 10:41:12 +0100 Subject: [PATCH 2/2] Clean up repo --- prompt | 1 - 1 file changed, 1 deletion(-) delete mode 100644 prompt diff --git a/prompt b/prompt deleted file mode 100644 index 609d568..0000000 --- a/prompt +++ /dev/null @@ -1 +0,0 @@ -No. I should have said, but survHEinla works in conjunction with survHE (https://github.com/giabaio/survHE), which has specific functions to handle the table of output and then all the plots. You need to modify some functions in survHE (survHE::like load_availables() and survHE:::manipulate_distributions() ) to ensure that the new Gamma distribution is available. You also need to create a new function lik_gam_inla() to go under runINLA.R to compute the log-density for the Gamma model. You can use the lik_***_inla() functions already present to see what you need to do.