Hi all,
I’m trying to fit a left-censored model using brms, but the estimates seem to be off for both the intercept and sigma parameter. The censReg model estimates the parameters just fine.
See below for an RME.
library(broom.mixed)
seed <- 1511
set.seed(seed)
generate_data <- function(sample_size, mean, sd, perc_missing) {
df <- tibble::tibble(y = rlnorm(sample_size, meanlog = mean, sdlog = sd),
lod = quantile(y, perc_missing, names = FALSE), # limit of detection
cens = ifelse(y < lod, 1, 0),
log_y1 = ifelse(cens, log(lod), log(y)))
return(df)
}
df <- generate_data(2000, 10, 2, 0.3)
m_1 <- censReg::censReg(log_y1 ~ 1, data = df,
method = "BHHH",
left = min(df$log_y1), right = Inf)
m_1 |> broom::tidy()
# A tibble: 2 × 5
# term estimate std.error statistic p.value
# <chr> <dbl> <dbl> <dbl> <dbl>
# 1 (Intercept) 9.99 0.0478 209. 0
# 2 logSigma 0.682 0.0192 35.6 2.07e-277
m_2 <- brms::brm(log_y1 | cens(cens) ~ 1,
data = df, seed = seed, backend = 'cmdstanr')
m_2 |> broom::tidy()
# A tibble: 2 × 8
# effect component group term estimate std.error conf.low conf.high
# <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
# 1 fixed cond NA (Intercept) 11.0 0.0353 11.0 11.1
# 2 ran_pars cond Residual sd__Observa… 1.36 0.0244 1.32 1.41
> sessionInfo()
# R version 4.3.1 (2023-06-16)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Ubuntu 22.04.1 LTS
#
# Matrix products: default
# BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
# LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#
# locale:
# [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
# [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
# [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
# [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
# [9] LC_ADDRESS=C LC_TELEPHONE=C
# [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#
# time zone: UTC
# tzcode source: system (glibc)
#
# attached base packages:
# [1] stats graphics grDevices datasets utils methods base
#
# other attached packages:
# [1] broom.mixed_0.2.9.4
#
# loaded via a namespace (and not attached):
# [1] Rdpack_2.6 gridExtra_2.3 inline_0.3.19
# [4] sandwich_3.0-2 rlang_1.1.1 magrittr_2.0.3
# [7] censReg_0.5-36 furrr_0.3.1 matrixStats_1.0.0
# [10] compiler_4.3.1 glmmML_1.1.5 loo_2.6.0
# [13] callr_3.7.3 vctrs_0.6.3 reshape2_1.4.4
# [16] stringr_1.5.0 pkgconfig_2.0.3 crayon_1.5.2
# [19] fastmap_1.1.1 backports_1.4.1 ellipsis_0.3.2
# [22] utf8_1.2.3 cmdstanr_0.5.3 threejs_0.3.3
# [25] promises_1.2.0.1 markdown_1.7 ps_1.7.5
# [28] miscTools_0.6-28 purrr_1.0.1 xfun_0.39
# [31] jsonlite_1.8.7 collapse_2.0.6 later_1.3.1
# [34] broom_1.0.5 parallel_4.3.1 prettyunits_1.1.1
# [37] R6_2.5.1 dygraphs_1.1.1.6 StanHeaders_2.26.27
# [40] stringi_1.7.12 parallelly_1.36.0 lmtest_0.9-40
# [43] knitr_1.43 Rcpp_1.0.10 rstan_2.21.8
# [46] zoo_1.8-12 base64enc_0.1-3 bayesplot_1.10.0
# [49] httpuv_1.6.11 Matrix_1.5-4.1 splines_4.3.1
# [52] igraph_1.5.0 tidyselect_1.2.0 abind_1.4-5
# [55] maxLik_1.5-2 codetools_0.2-19 miniUI_0.1.1.1
# [58] processx_3.8.2 listenv_0.9.0 pkgbuild_1.4.2
# [61] lattice_0.21-8 tibble_3.2.1 plyr_1.8.8
# [64] withr_2.5.0 shiny_1.7.4 bridgesampling_1.1-2
# [67] posterior_1.4.1 coda_0.19-4 future_1.33.0
# [70] RcppParallel_5.1.7 xts_0.13.1 pillar_1.9.0
# [73] tensorA_0.36.2 checkmate_2.2.0 DT_0.28
# [76] stats4_4.3.1 shinyjs_2.1.0 distributional_0.3.2
# [79] generics_0.1.3 ggplot2_3.4.2 rstantools_2.3.1
# [82] munsell_0.5.0 scales_1.2.1 globals_0.16.2
# [85] gtools_3.9.4 xtable_1.8-4 bspm_0.5.2
# [88] glue_1.6.2 tools_4.3.1 shinystan_2.6.0
# [91] data.table_1.14.8 colourpicker_1.2.0 forcats_1.0.0
# [94] mvtnorm_1.2-2 grid_4.3.1 tidyr_1.3.0
# [97] rbibutils_2.2.16 crosstalk_1.2.0 bdsmatrix_1.3-6
# [100] colorspace_2.1-0 nlme_3.1-162 Formula_1.2-5
# [103] cli_3.6.1 plm_2.6-3 fansi_1.0.4
# [106] Brobdingnag_1.2-9 dplyr_1.1.2 gtable_0.3.3
# [109] digest_0.6.32 brms_2.19.0 htmlwidgets_1.6.2
# [112] farver_2.1.1 htmltools_0.5.5 lifecycle_1.0.3
# [115] mime_0.12 shinythemes_1.2.0 MASS_7.3-60
Any ideas? Thanks!