Using {projpred} latent projection with {brms} Weibull family models

Interfaces brms
, , ,
15

Thanks for this, @fweber144.

I had a go at writing some pseudocode for latent_ll_oscale and latent_ppd_oscale based on how brms:::log_lik_censor and brms:::rcontinuous do this for censored models, respectively (assuming I’ve understood how both of these work correctly). I’ve done this for the log-normal family just to make it easier to track how dis fits into it, rather than the aforementioned ad hoc use of the Weibull model’s shape draws. I assume latent_ilink is unaffected by censoring given that the inverse link function is the same for censored or uncensored observations? I’m by no means an expert R programmer, so I imagine there are faster ways for indexing/structuring the code. Note that idxs_cens is assumed to be a vector of column indices.

latent_ll_oscale_cens <- function(
    ilpreds 
    ,dis
    ,y_oscale
    ,wobs = rep(1, length(y_oscale))
    ,cl_ref
    ,idxs_cens
    ,wdraws_ref = rep(1, length(cl_ref))
) {
  
  y_oscale_mat <- matrix(
    y_oscale
    ,nrow = nrow(ilpreds)
    ,ncol = ncol(ilpreds)
    ,byrow = TRUE
  )
  
  wobs_mat <- matrix(
    wobs
    ,nrow = nrow(ilpreds)
    ,ncol = ncol(ilpreds)
    ,byrow = TRUE
  )
  
  ll_unw <- matrix(NA, nrow = nrow(ilpreds), ncol = ncol(ilpreds), byrow = TRUE)
  
  # y_oscale_mat and ilpreds are S x N matrices, so subset by column for now
  # Censored observations (based on brms:::log_lik_censor)
  ll_unw[, idxs_cens] <- plnorm(
    y_oscale_mat[, idxs_cens]
    ,meanlog = ilpreds[, idxs_cens]
    ,sdlog = dis
    ,lower.tail = FALSE
    ,log.p = TRUE
  )
  
  # Uncensored observations
  ll_unw[, -idxs_cens] <- dlnorm(
    y_oscale_mat[, -idxs_cens]
    ,meanlog = ilpreds[, -idxs_cens]
    ,sdlog = dis
    ,log = TRUE
  )
  
  return(wobs_mat * ll_unw)
  
}

latent_ppd_oscale_cens <- function(
    ilpreds_resamp
    ,dis_resamp
    ,wobs
    ,cl_ref
    ,wdraws_ref = rep(1, length(cl_ref))
    ,idxs_prjdraws
    ,idxs_cens
    ,newdata
) {
  
  if (is.null(newdata)) {
    newdata <- # <original dataset>
  }
  observed_times <- newdata[, <column for observed timepoints>]
  
  # Initialise S x N matrix
  ppd <- matrix(NA, nrow = nrow(ilpreds_resamp), ncol = ncol(ilpreds_resamp))
    
  for (i in <all observations in newdata>) {
    
    if (i %in% idxs_cens) {
     
      # Right-censored observations (based on brms:::rcontinuous)
      plb <- plnorm(-Inf, meanlog = ilpreds_resamp[, i], sdlog = dis_resamp)
      pub <- plnorm(observed_times[i], meanlog = ilpreds_resamp[, i], sdlog = dis_resamp)
      ppd[, i] <- runif(nrow(ilpreds_resamp), min = plb, max = pub) |> 
        qlnorm(meanlog = ilpreds_resamp[, i], sdlog = dis_resamp)
      
    } else {
      
      # Uncensored observations
      ppd[, i] <- rlnorm(
        nrow(ilpreds_resamp)
        ,meanlog = ilpreds_resamp[, i]
        ,sdlog = dis_resamp
      )
      
    }
    
  }
  
  return(ppd)
  
}

@avehtari, it would be great to know if I’m on the right track here - thanks in advance!

Edit: grammar