You’re welcome. New target date is 2023-09-25.
Version 1.0-0 binary packages are now available for Mac at howto . Priors for proportional odds parts of the model are now specified with pcontrast=list(...)
as documented here. Windows and Linux versions are available but I’ll be updating the Windows binary a little by tomorrow.
The new version is now on CRAN.
Thanks Frank. Bit of a dumb question, but how would I use pcontrast
in the example I gave above to specify a prior SD of, say, 10, for the proportional odds parts of the model with one covariate? (I am assuming that the default is SD = 100 still for the PO components)
For the PO components, predictors not named in pcontrast
have uninformative priors. For departures from PO you still need keepsep
and priorsdppo
or something similar to that in name. In the future I’ll move the non-PO part also under pcontrast
when I figure out how to do it, e.g. .(tx=‘B’, ycut=3), .(tx=‘B’, ycut=4)
.
Right, so if one is just using a fully PO model, I’m struggling to understand how to change the prior SD for the PO components using pcontrast
in blrm
, e.g. something like the example I gave above but with no PPO components:
mod <- blrm(y~x, data=data, seed = 1234, iter = 2000, chains = 4)
One of the examples I linked to should give you enough to work this out. It comes down to forming a comparison such as .(tx=‘B’), .(tx=‘A’)
where you define . <- function(…) list(…)
and the B-A example contrast is replaced with your aim comparison of interest, and you specify an SD for the prior for that.
Can you please point me to a vignette that shows how to fit a full joint Bayesian model of covariates (with missing data) and outcome? The outcome is ordinal (with absorbing state), and the time-varying covariates are also ordinal, but with considerable missing data due to missed visits, etc. I saw some of the examples that you provided for ordinal Markov longitudinal modeling, but I am not sure how to set it up so that the longitudinal predictors are also modeled. Thank you.
I’ve been looking for that. If you or someone else finds that please post it here.
Does that mean it’s not possible with the rmsb package, or you just haven’t made a vignette for it?
It has not been implemented in rmsb
. This probably requires a custom Stan
model (which really contains a set of models) under your control.
I’m having trouble with ‘predict’ and type = “fitted” - maybe someone has a solution to this
Here is a minimal reproducible example:
library(rmsb)
getHdata(support)
dd <- datadist(support); options(datadist='dd')
bsup <- blrm(hospdead ~ dzgroup + rcs(crea, 5) + rcs(meanbp, 5), data=support)
newdata = eval(bsup$call$data)[all.vars(bsup$sformula)] %>% data.frame
newdata <- newdata[complete.cases(newdata), ]
newdatapred_binary <- predict(bsup, newdata, type="fitted", fun=plogis, funint=FALSE, posterior.summary="all")
Error in ints + betas %*% t(X) : non-conformable arrays
Rversion: “R version 4.3.2 (2023-10-31)”
SessionInfo:
R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2.1
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] ggpubr_0.6.0 viridis_0.6.5 viridisLite_0.4.2 conflicted_1.2.0 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
[8] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 tidyverse_2.0.0 qreport_1.0-0
[15] ggplot2_3.4.4 data.table_1.15.0 magrittr_2.0.3 modelsummary_1.4.5 flextable_0.9.4 huxtable_5.5.6 texreg_1.39.3
[22] rvest_1.0.4 kableExtra_1.4.0 knitr_1.45 rmsb_1.1-1 rms_6.8-0 Hmisc_5.1-1 targets_1.5.1
Your best bet is probably MICE + multiple runs of the model and then you just combine the posteriors. (https://cran.r-project.org/web/packages/brms/vignettes/brms_missings.html).
It’s maybe not quite as elegant but much more flexible. Note also that since rmsb uses Stan even if the was a set-up for imputation on the fly I do not believe it could impute any discrete variables (eg, previous event counts, baseline variables like sex or race etc) since this can’t be done with Stan.
And as an aside just do predict(fit, complete.cases(support), …)
The following solution worked for me.
in predict.blrm()
instead of:
if(ns == 1) return(cumprob(ints + betas * t(X))) # binary logistic model
add intercept:
if(ns == 1) {
lp ← betas * t(X)
lp ← sweep(lp, 1, ints[,1], ‘+’)
return(cumprob(lp))
}
in removeFormulaTerms()
replace
form ← as.formula(gsub(‘offset(’, ‘.off.(’, form, fixed=TRUE))
with
form ← as.formula(paste(gsub(‘offset(’, ‘.off.(’, form, fixed=TRUE), collapse=“”))
Thanks for the cumprob() fix which will be in the next release. For the as.formula(), form can only be a single character string so that fix should not be needed.