brms
brms copied to clipboard
Published
20 hours ago •
paul-buerkner
paul-buerkner
feature request: alternative approach for censored data
Context:
I often work with data that are censored. The approach for censored data implemented in brms is really useful, but I have found it to be slow - especially when a large part of the data are censored and with distributions where it is costly to integrate out censored data.
The approach taken by brms is explained in the stan user guide: https://mc-stan.org/docs/stan-users-guide/truncation-censoring.html#integrating-out-censored-values
This involves integrating out the PDF between the censoring bounds and therefore requires the CDF and/or CCDF functions depending on left, right or interval censoring.
Alternative:
The stan user guide also mentions an alternative approach: https://mc-stan.org/docs/stan-users-guide/truncation-censoring.html#estimating-censored-values
I have found this approach to be much faster (about 30 times faster in below reprex) and it recovers the same parameter estimates. In this approach, the censored data are declared as parameters to be estimated. I think this approach is called a data augmentation (imputation) approach and is discussed in https://academic.oup.com/jrsssb/article-abstract/55/1/3/7035892 (section 6.3). It only requires the PDF function.
Notes:
The below reprex is for a beta distribution. I have extracted the brms generated stan code and changed this according to the User's guide, but without changing how brms requires the data to be formatted (so standata() is the same in both implementations). I am not very familiar with programming in stan, so some further efficiency gains are probably possible.
I have also tested this with a Gamma distribution and this worked too, but the speed-up is not that spectacular. I can imagine that for a Gaussian distribution there would be no speed-up.
Feature request:
Would you be willing to consider implementing this approach in brms? Maybe not as a replacement to the current implementation, but as an alternative. Perhaps an extra argument in the cens() function to switch between both alternatives?
I think this approach might not work for discrete distributions because stan cannot sample discrete parameters? Although I have a feeling this might be related: https://users.aalto.fi/~johnsoa2/notebooks/CopulaIntro.html#bayesian-estimation-data-augmentation. (haven't thought this through)
Reproducible example:
library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.20.4). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#>
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#>
#> ar
library(tidyverse)
## generate data from beta distribution
set.seed(2154)
mu <- 0.5
phi <- 2
true_data <- tibble(
y = rbeta(100, shape1 = mu * phi, shape2 = (1 - mu) * phi)
)
## Censoring
censored_data <- true_data %>%
mutate(
y_left = floor(y * 10) / 10,
y_right = ceiling(y * 10) / 10,
y_left = ifelse(y_left == 0, 1e-7, y_left),
y_right = ifelse(y_right == 1, 1 - 1e-7, y_right),
censoring = case_when(
y_right == 0.1 ~ "left",
y_left == 0.9 ~ "right",
TRUE ~ "interval"
),
y_cens = case_when(
censoring == "left" ~ y_right,
censoring == "right" ~ y_left,
TRUE ~ y_left
)
)
## plot censored data
censored_data %>%
arrange(y) %>%
mutate(
row_number = row_number()
) %>%
ggplot(
aes(x = row_number, colour = censoring)) +
geom_point(
aes(y = y),
position = position_jitter(width = 0.5, height = 0)) +
geom_linerange(
aes(y = y, ymin = y_left, ymax = y_right),
alpha = 0.3,
position = position_jitter(width = 0.5, height = 0)
)
# fit brms model without sampling
beta_censored <- brm(
y_cens | cens(censoring, y_right) ~ 1,
family = brms::Beta(),
data = censored_data,
iter = 0
)
#> Compiling Stan program...
#> Start sampling
#> Error in config_argss(chains = chains, iter = iter, warmup = warmup, thin = thin, :
#> parameter 'iter' should be a positive integer
#> error in specifying arguments; sampling not done
# extract stan code and data
beta_stan_code <- brms::stancode(beta_censored)
beta_stan_data <- brms::standata(beta_censored)
# fit using cmdstanr (for easier comparison of timing for sampling)
writeLines(beta_stan_code, "beta_censored_brms.stan")
beta_censored_mod <- cmdstanr::cmdstan_model("beta_censored_brms.stan")
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:356,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:14,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a81bb744e7.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> In file included from C:/rtools43/ucrt64/include/c++/13.2.0/algorithm:61,
#> from stan/lib/stan_math/lib/eigen_3.4.0/Eigen/Core:95,
#> from stan/lib/stan_math/lib/eigen_3.4.0/Eigen/Dense:1,
#> from stan/lib/stan_math/stan/math/prim/fun/Eigen.hpp:22,
#> from stan/lib/stan_math/stan/math/rev.hpp:4:
#> In function 'void std::__final_insertion_sort(_RandomAccessIterator, _RandomAccessIterator, _Compare) [with _RandomAccessIterator = unsigned int*; _Compare = __gnu_cxx::__ops::_Iter_comp_iter<less<long double> >]',
#> inlined from 'void std::__sort(_RandomAccessIterator, _RandomAccessIterator, _Compare) [with _RandomAccessIterator = unsigned int*; _Compare = __gnu_cxx::__ops::_Iter_comp_iter<less<long double> >]' at C:/rtools43/ucrt64/include/c++/13.2.0/bits/stl_algo.h:1950:31,
#> inlined from 'void std::sort(_RAIter, _RAIter, _Compare) [with _RAIter = unsigned int*; _Compare = less<long double>]' at C:/rtools43/ucrt64/include/c++/13.2.0/bits/stl_algo.h:4894:18,
#> inlined from 'unsigned int boost::math::detail::set_crossover_locations(const Seq&, const Seq&, const Real&, unsigned int*) [with Seq = std::vector<double, std::allocator<double> >; Real = long double]' at stan/lib/stan_math/lib/boost_1.78.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp:113:21,
#> inlined from 'std::pair<T, T> boost::math::detail::hypergeometric_pFq_checked_series_impl(const Seq&, const Seq&, const Real&, const Policy&, const Terminal&, long long int&) [with Seq = std::vector<double, std::allocator<double> >; Real = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy
#> , boost::math::policies::default_policy, boost::math::policies::default_policy>; Terminal = iteration_terminator]' at stan/lib/stan_math/lib/boost_1.78.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp:260:56:
#> C:/rtools43/ucrt64/include/c++/13.2.0/bits/stl_algo.h:1859:32: warning: array subscript 16 is outside array bounds of 'unsigned int [5]' [-Warray-bounds=]
#> 1859 | std::__insertion_sort(__first, __first + int(_S_threshold), __comp);
#> | ~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#> In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/math/special_functions/hypergeometric_pFq.hpp:11,
#> from stan/lib/stan_math/stan/math/prim/fun/hypergeometric_pFq.hpp:7,
#> from stan/lib/stan_math/stan/math/prim/fun/hypergeometric_2F1.hpp:18,
#> from stan/lib/stan_math/stan/math/prim/fun/grad_2F1.hpp:14,
#> from stan/lib/stan_math/stan/math/prim/fun.hpp:129,
#> from stan/lib/stan_math/stan/math/rev/fun/multiply.hpp:7,
#> from stan/lib/stan_math/stan/math/rev/fun/elt_multiply.hpp:9,
#> from stan/lib/stan_math/stan/math/rev/fun.hpp:56,
#> from stan/lib/stan_math/stan/math/rev.hpp:10:
#> stan/lib/stan_math/lib/boost_1.78.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp: In function 'std::pair<T, T> boost::math::detail::hypergeometric_pFq_checked_series_impl(const Seq&, const Seq&, const Real&, const Policy&, const Terminal&, long long int&) [with Seq = std::vector<double, std::allocator<double> >; Real = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>; Terminal = iteration_terminator]':
#> stan/lib/stan_math/lib/boost_1.78.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp:258:18: note: at offset 64 into object 'crossover_locations' of size 20
#> 258 | unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS];
#> | ^~~~~~~~~~~~~~~~~~~
beta_integrated_out <- beta_censored_mod$sample(
data = beta_stan_data,
seed = 123,
chains = 4,
parallel_chains = 1
)
#> Running MCMC with 4 sequential chains...
#>
#> Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: beta_lccdf: Second shape parameter is 0, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a81bb744e7.stan', line 34, column 8 to column 68)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1 finished in 64.4 seconds.
#> Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a81bb744e7.stan', line 20, column 2 to column 41)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2 finished in 73.5 seconds.
#> Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 3 Exception: gamma_lpdf: Random variable is 0, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a81bb744e7.stan', line 20, column 2 to column 41)
#> Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 3
#> Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 3 Exception: beta_lccdf: Second shape parameter is 0, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a81bb744e7.stan', line 34, column 8 to column 68)
#> Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 3
#> Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3 finished in 56.5 seconds.
#> Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 4 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a81bb744e7.stan', line 20, column 2 to column 41)
#> Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 4
#> Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4 finished in 67.5 seconds.
#>
#> All 4 chains finished successfully.
#> Mean chain execution time: 65.5 seconds.
#> Total execution time: 262.7 seconds.
# changed brms generated stan code to use data augmentation approach
# rather than integrating out (i.e. avoid use of CDF or CCDF)
stanlines <- "
// generated with brms 2.20.4
functions {
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
array[N] int<lower=-1,upper=2> cens; // indicates censoring
vector[N] rcens; // right censor points for interval censoring
int prior_only; // should the likelihood be ignored?
}
transformed data {
// count uncensored data points
int Nobs = 0;
for (i in 1:N) {
if (cens[i] == 0) {
Nobs += 1;
}
}
// extract uncensored response values
array[Nobs] real Yobs;
int index = 1; // index to keep track of the position
// Populate with the values from Y
for (i in 1:N) {
if (cens[i] == 0) {
Yobs[index] = Y[i];
index += 1;
}
}
// count right censored data points
int Nright = 0;
for (i in 1:N) {
if (cens[i] == 1) {
Nright += 1;
}
}
// extract right censored lower bounds
array[Nright] real LBright;
int index2 = 1; // index to keep track of the position
// Populate with the values from Y
for (i in 1:N) {
if (cens[i] == 1) {
LBright[index2] = Y[i];
index2 += 1;
}
}
// count left censored data points
int Nleft = 0;
for (i in 1:N) {
if (cens[i] == -1) {
Nleft += 1;
}
}
// extract left censored upper bounds
array[Nleft] real UBleft;
int index3 = 1; // index to keep track of the position
// Populate with the values from Y
for (i in 1:N) {
if (cens[i] == -1) {
UBleft[index3] = Y[i];
index3 += 1;
}
}
// count interval censored data points
int Nint = 0;
for (i in 1:N) {
if (cens[i] == 2) {
Nint += 1;
}
}
// extract interval censored lower and upper bounds
array[Nint] real LBint;
array[Nint] real UBint;
int index4 = 1; // index to keep track of the position
// Populate with the values from Y
for (i in 1:N) {
if (cens[i] == 2) {
LBint[index4] = Y[i];
UBint[index4] = rcens[i];
index4 += 1;
}
}
}
parameters {
real Intercept; // temporary intercept for centered predictors
real<lower=0> phi; // precision parameter
array[Nright] real<lower=LBright, upper=1> Yright;
array[Nleft] real<lower=0, upper=UBleft> Yleft;
array[Nint] real<lower=LBint, upper=UBint> Yint;
}
transformed parameters {
real lprior = 0; // prior contributions to the log posterior
lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
lprior += gamma_lpdf(phi | 0.01, 0.01);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[Nobs] muobs = rep_vector(0.0, Nobs);
muobs += Intercept;
muobs = inv_logit(muobs);
for (n in 1:Nobs) {
target += beta_lpdf(Yobs[n] | muobs[n] * phi, (1 - muobs[n]) * phi);
}
vector[Nleft] muleft = rep_vector(0.0, Nleft);
muleft += Intercept;
muleft = inv_logit(muleft);
for (n in 1:Nleft) {
target += beta_lpdf(Yleft[n] | muleft[n] * phi, (1 - muleft[n]) * phi);
}
vector[Nright] muright = rep_vector(0.0, Nright);
muright += Intercept;
muright = inv_logit(muright);
for (n in 1:Nright) {
target += beta_lpdf(Yright[n] | muright[n] * phi, (1 - muright[n]) * phi);
}
vector[Nint] muint = rep_vector(0.0, Nint);
muint += Intercept;
muint = inv_logit(muint);
for (n in 1:Nint) {
target += beta_lpdf(Yint[n] | muint[n] * phi, (1 - muint[n]) * phi);
}
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept;
}
"
writeLines(stanlines, "bda.stan")
# fit data augmentation model with cmdstanr
bda <- cmdstanr::cmdstan_model("bda.stan")
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:356,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:14,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a854246415.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
beta_data_augmentation <- bda$sample(
data = beta_stan_data,
seed = 123,
chains = 4,
parallel_chains = 1
)
#> Running MCMC with 4 sequential chains...
#>
#> Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: beta_lpdf: Second shape parameter is 0, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a854246415.stan', line 113, column 6 to column 78)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1 finished in 2.0 seconds.
#> Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2 finished in 2.0 seconds.
#> Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 3 Exception: beta_lpdf: Second shape parameter is 0, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a854246415.stan', line 113, column 6 to column 78)
#> Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 3
#> Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3 finished in 2.0 seconds.
#> Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 4 Exception: beta_lpdf: Second shape parameter is 0, but must be positive finite! (in 'C:/Users/HANS_V~1/AppData/Local/Temp/Rtmp4yUrzl/model-34a854246415.stan', line 113, column 6 to column 78)
#> Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 4
#> Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
#> Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
#> Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
#> Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
#> Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
#> Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
#> Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
#> Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
#> Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
#> Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
#> Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4 finished in 2.0 seconds.
#>
#> All 4 chains finished successfully.
#> Mean chain execution time: 2.0 seconds.
#> Total execution time: 8.6 seconds.
# compare timings
beta_integrated_out$time()
#> $total
#> [1] 262.7254
#>
#> $chains
#> chain_id warmup sampling total
#> 1 1 28.933 35.458 64.391
#> 2 2 30.083 43.453 73.536
#> 3 3 28.928 27.539 56.467
#> 4 4 30.996 36.493 67.489
beta_data_augmentation$time()
#> $total
#> [1] 8.63974
#>
#> $chains
#> chain_id warmup sampling total
#> 1 1 1.091 0.907 1.998
#> 2 2 1.126 0.914 2.040
#> 3 3 1.102 0.909 2.011
#> 4 4 1.098 0.919 2.017
beta_data_augmentation$time()$total / beta_integrated_out$time()$total
#> [1] 0.03288506
# compare estimated parameters mu and phi
beta_integrated_out$summary(variables = c("b_Intercept", "phi"))
#> # A tibble: 2 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 b_Intercept 0.190 0.186 0.112 0.110 0.00636 0.378 1.00 2829. 2563.
#> 2 phi 2.42 2.40 0.335 0.342 1.91 2.98 1.00 3389. 2448.
beta_data_augmentation$summary(variables = c("b_Intercept", "phi"))
#> # A tibble: 2 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 b_Intercept 0.187 0.186 0.107 0.105 0.00940 0.363 1.00 6905. 2961.
#> 2 phi 2.42 2.40 0.328 0.334 1.92 2.99 1.00 6211. 3134.
# graphical comparison between both approaches
drw_bio <- posterior::as_draws_rvars(beta_integrated_out)
drw_bda <- posterior::as_draws_rvars(beta_data_augmentation)
bda_mu_phi <- drw_bda %>%
as_draws_df() %>%
select(Intercept, phi) %>%
mutate(
mu = plogis(Intercept),
model = "data augmentation"
)
#> Warning: Dropping 'draws_df' class as required metadata was removed.
bio_mu_phi <- drw_bio %>%
as_draws_df() %>%
select(Intercept, phi) %>%
mutate(
mu = plogis(Intercept),
model = "integrate out (brms implementation)"
)
#> Warning: Dropping 'draws_df' class as required metadata was removed.
bind_rows(bio_mu_phi, bda_mu_phi) %>%
pivot_longer(cols = c("phi", "mu")) %>%
ggplot() +
geom_density(aes(x = value, colour = model)) +
facet_wrap(~name, scales = "free")
# visualise estimated (imputed / data augmented) parameters
imputed_left <- drw_bda$Yleft %>%
as_draws_df() %>%
pivot_longer(cols = starts_with("x")) %>%
group_by(name) %>%
mutate(
censoring = "left",
interval_mean = round(mean(value), 2)
)
#> Warning: Dropping 'draws_df' class as required metadata was removed.
imputed_right <- drw_bda$Yright %>%
as_draws_df() %>%
pivot_longer(cols = starts_with("x")) %>%
group_by(name) %>%
mutate(
censoring = "right",
interval_mean = round(mean(value), 2)
)
#> Warning: Dropping 'draws_df' class as required metadata was removed.
imputed_interval <- drw_bda$Yint %>%
as_draws_df() %>%
pivot_longer(cols = starts_with("x")) %>%
group_by(name) %>%
mutate(
censoring = "interval",
interval_mean = round(mean(value), 2)
)
#> Warning: Dropping 'draws_df' class as required metadata was removed.
bind_rows(
imputed_right, imputed_left, imputed_interval
) %>%
ggplot() +
geom_density(aes(x = value, colour = name)) +
facet_wrap(interval_mean ~ censoring, scales = "free_x") +
theme(legend.position = "none") +
labs(title = "data augmentation / imputed values")
Created on 2024-05-19 with reprex v2.1.0
Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.3.2 (2023-10-31 ucrt)
#> os Windows 10 x64 (build 19045)
#> system x86_64, mingw32
#> ui RTerm
#> language (EN)
#> collate Dutch_Belgium.utf8
#> ctype Dutch_Belgium.utf8
#> tz Europe/Brussels
#> date 2024-05-19
#> pandoc 3.1.12.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> ! package * version date (UTC) lib source
#> abind 1.4-5 2016-07-21 [1] CRAN (R 4.3.0)
#> backports 1.4.1 2021-12-13 [1] CRAN (R 4.3.0)
#> base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.3.0)
#> bayesplot 1.11.1 2024-02-15 [1] CRAN (R 4.3.2)
#> bridgesampling 1.1-2 2021-04-16 [1] CRAN (R 4.3.0)
#> brms * 2.20.4 2023-09-25 [1] CRAN (R 4.3.1)
#> Brobdingnag 1.2-9 2022-10-19 [1] CRAN (R 4.3.0)
#> checkmate 2.3.1 2023-12-04 [1] CRAN (R 4.3.2)
#> cli 3.6.2 2023-12-11 [1] CRAN (R 4.3.2)
#> cmdstanr 0.6.1 2023-09-12 [1] local
#> coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.3.2)
#> codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.2)
#> colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
#> colourpicker 1.3.0 2023-08-21 [1] CRAN (R 4.3.1)
#> crosstalk 1.2.1 2023-11-23 [1] CRAN (R 4.3.2)
#> curl 5.2.0 2023-12-08 [1] CRAN (R 4.3.2)
#> data.table 1.15.0 2024-01-30 [1] CRAN (R 4.3.2)
#> digest 0.6.35 2024-03-11 [1] CRAN (R 4.3.3)
#> distributional 0.4.0 2024-02-07 [1] CRAN (R 4.3.2)
#> dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.2)
#> DT 0.32 2024-02-19 [1] CRAN (R 4.3.2)
#> dygraphs 1.1.1.6 2018-07-11 [1] CRAN (R 4.3.0)
#> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
#> emmeans 1.10.0 2024-01-23 [1] CRAN (R 4.3.2)
#> estimability 1.5 2024-02-20 [1] CRAN (R 4.3.2)
#> evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
#> fansi 1.0.6 2023-12-08 [1] CRAN (R 4.3.2)
#> farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
#> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
#> forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.0)
#> fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.1)
#> generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
#> ggplot2 * 3.5.0 2024-02-23 [1] CRAN (R 4.3.2)
#> glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.2)
#> gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
#> gtable 0.3.4 2023-08-21 [1] CRAN (R 4.3.1)
#> gtools 3.9.5 2023-11-20 [1] CRAN (R 4.3.2)
#> highr 0.10 2022-12-22 [1] CRAN (R 4.3.0)
#> hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
#> htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.2)
#> htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.2)
#> httpuv 1.6.14 2024-01-26 [1] CRAN (R 4.3.2)
#> igraph 2.0.2 2024-02-17 [1] CRAN (R 4.3.3)
#> inline 0.3.19 2021-05-31 [1] CRAN (R 4.3.0)
#> jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.2)
#> knitr 1.45 2023-10-30 [1] CRAN (R 4.3.2)
#> labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.1)
#> later 1.3.2 2023-12-06 [1] CRAN (R 4.3.2)
#> lattice 0.22-5 2023-10-24 [1] CRAN (R 4.3.2)
#> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.2)
#> loo 2.7.0 2024-02-24 [1] CRAN (R 4.3.2)
#> lubridate * 1.9.3.9000 2024-01-29 [1] https://inbo.r-universe.dev (R 4.3.2)
#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
#> markdown 1.12 2023-12-06 [1] CRAN (R 4.3.2)
#> MASS 7.3-60.0.1 2024-01-13 [2] CRAN (R 4.3.2)
#> Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.2)
#> matrixStats 1.2.0 2023-12-11 [1] CRAN (R 4.3.2)
#> mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
#> miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.0)
#> multcomp 1.4-25 2023-06-20 [1] CRAN (R 4.3.1)
#> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
#> mvtnorm 1.2-4 2023-11-27 [1] CRAN (R 4.3.2)
#> nlme 3.1-164 2023-11-27 [2] CRAN (R 4.3.2)
#> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
#> pkgbuild 1.4.4 2024-03-17 [1] CRAN (R 4.3.3)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
#> plyr 1.8.9 2023-10-02 [1] CRAN (R 4.3.1)
#> posterior 1.5.0 2023-10-31 [1] CRAN (R 4.3.2)
#> processx 3.8.4 2024-03-16 [1] CRAN (R 4.3.3)
#> promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.1)
#> ps 1.7.6 2024-01-18 [1] CRAN (R 4.3.2)
#> purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.1)
#> QuickJSR 1.1.3 2024-01-31 [1] CRAN (R 4.3.2)
#> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.3.0)
#> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0)
#> R.oo 1.26.0 2024-01-24 [1] CRAN (R 4.3.2)
#> R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.3.2)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
#> Rcpp * 1.0.12 2024-01-09 [1] CRAN (R 4.3.2)
#> D RcppParallel 5.1.7 2023-02-27 [1] CRAN (R 4.3.0)
#> readr * 2.1.5 2024-01-10 [1] CRAN (R 4.3.2)
#> reprex 2.1.0 2024-01-11 [1] CRAN (R 4.3.2)
#> reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.3.0)
#> rlang 1.1.3 2024-01-10 [1] CRAN (R 4.3.2)
#> rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1)
#> rstan 2.32.6 2024-03-05 [1] CRAN (R 4.3.3)
#> rstantools 2.4.0 2024-01-31 [1] CRAN (R 4.3.2)
#> rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.1)
#> sandwich 3.1-0 2023-12-11 [1] CRAN (R 4.3.2)
#> scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.2)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
#> shiny 1.8.0 2023-11-17 [1] CRAN (R 4.3.2)
#> shinyjs 2.1.0 2021-12-23 [1] CRAN (R 4.3.0)
#> shinystan 2.6.0 2022-03-03 [1] CRAN (R 4.3.0)
#> shinythemes 1.2.0 2021-01-25 [1] CRAN (R 4.3.0)
#> StanHeaders 2.32.6 2024-03-01 [1] CRAN (R 4.3.3)
#> stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.2)
#> stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.3.2)
#> styler 1.10.2 2023-08-29 [1] CRAN (R 4.3.1)
#> survival 3.5-8 2024-02-14 [2] CRAN (R 4.3.2)
#> tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.3.2)
#> TH.data 1.1-2 2023-04-17 [1] CRAN (R 4.3.1)
#> threejs 0.3.3 2020-01-21 [1] CRAN (R 4.3.0)
#> tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
#> tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.3.2)
#> tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.3)
#> tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.0)
#> timechange 0.3.0 2024-01-18 [1] CRAN (R 4.3.2)
#> tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.1)
#> utf8 1.2.4 2023-10-22 [1] CRAN (R 4.3.2)
#> V8 4.4.2 2024-02-15 [1] CRAN (R 4.3.2)
#> vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.2)
#> withr 3.0.0 2024-01-16 [1] CRAN (R 4.3.2)
#> xfun 0.41 2023-11-01 [1] CRAN (R 4.3.2)
#> xml2 1.3.6 2023-12-04 [1] CRAN (R 4.3.2)
#> xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0)
#> xts 0.13.2 2024-01-21 [1] CRAN (R 4.3.2)
#> yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.2)
#> zoo 1.8-12 2023-04-13 [1] CRAN (R 4.3.0)
#>
#> [1] C:/R/library
#> [2] C:/R/R-4.3.2/library
#>
#> D ── DLL MD5 mismatch, broken installation.
#>
#> ──────────────────────────────────────────────────────────────────────────────
