How to use bam (mgcv) including crossed random effects?

I’m interested in the relationship of the arterial ‘po2’ (response) according to the total delay of analysis ‘deltot’ (main predictor) since po2 is known to increase with deltot.

More precisely, deltot is the total delay between the time of blood sampling ‘tbs’ - as informed by nurses from medical units ‘medu’ - and the end of measurement at laboratory (i.e., when the result is given by the analyser).

Purpose: I’d like to test if rounding tbs, which may false deltot, may impact the po2-to-deltot (nonlinear) relationship.

For this purpose, I’ve classified the data according to 3 conditions ‘cond’ based on tbs:

• cond = 0: patients for whom tbs have been precisely informed by nurses, which means likely reliable deltot
• cond = 1: patients for whom tbs have been rounded per half-hour by nurses, which means likely less reliable deltot
• cond = 2: patients for whom tbs have been rounded per hour by nurses, which means likely not reliable deltot

I’ve to account for potential random effects given that:

• po2 may vary between patients (e.g., blood po2 decreases with age),
• po2 may vary depending on pathophysiological contexts; I therefore assumed po2 to vary depending on medu,
• Since over years, a given patient may have one or more po2 results in one or more medu, I think crossed random effects are appropriate.

I want to use bam function (mgcv) since it is more efficient for the dataset ‘dat’ with nearly 750,000 po2 results, measured on a total of 109,000 patients hospitalized in nearly 640 medu.

Since po2 distribution is markedly right-skewed, I think ‘family = scat’ is appropriate.

Lastly, po2 and deltot are positive continuous variables, and cond and medu are non-ordered factors.

Question: what would the best model for the purpose and why?:

  1) m1 <- bam(po2 ~
                 cond + s(deltot, by = cond, bs = "cr") +
                 s(medu, bs = "re") +
                 s(medu, deltot, bs = "re"),
               data = dat, method = 'fREML', family = "scat", discrete = TRUE)
  
  2) m2 <- bam(po2 ~
                 cond + s(deltot, by = cond, bs = "cr") +
                 s(medu, bs = "re", by = cond) +
                 s(medu, deltot, bs = "re", by = cond),
               data = dat, method = 'fREML', family = "scat", discrete = TRUE)
  
  3) m3 <- bam(po2a ~
                 s(deltot, bs = "cr") +
                 s(deltot, medu, bs = "fs", m = 1),
               data = dat, method = 'fREML', family = "scat", discrete = TRUE)
  
  4) An other model ?

Thanks for help