parameters for random need to reset the index

[DominiqueMakowski]

If I run this model for more than one participant with participants as random intercept, I get the following:

I think there's a missing row.names() <- NULL somwhere

reprex

``` r library(brms) #> Loading required package: Rcpp #> Loading 'brms' package (version 2.15.0). 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 # Stan family ------------------------------------------------------------- rtmix <- brms::custom_family( "rtmix", dpars = c("mu", "sigma", "mix", "shiftprop"), # Those will be estimated links = c("identity", "log", "logit", "logit"), type = "real", lb = c(NA, 0, 0, 0), # bounds for the parameters ub = c(NA, NA, 1, 1), vars = c("vreal1[n]", "vreal2[n]") # Data for max_shift and upper (known) ) # Stan model ------------------------------------------------------------------- stan_functions <- brms::stanvar(block = "functions", scode = " real rtmix_lpdf(real y, real mu, real sigma, real mix, real shiftprop, real max_shift, real upper) { real shift = shiftprop * max_shift; if(y <= shift) { // Could only be created by the contamination return log(mix) + uniform_lpdf(y | 0, upper); } else if(y >= upper) { // Could only come from the lognormal return log1m(mix) + lognormal_lpdf(y - shift | mu, sigma); } else { // Actually mixing real lognormal_llh = lognormal_lpdf(y - shift | mu, sigma); real uniform_llh = uniform_lpdf(y | 0, upper); return log_mix(mix, uniform_llh, lognormal_llh); } } ") + brms::stanvar(block = "functions", scode = " real rtmix_lcdf(real y, real mu, real sigma, real mix, real shiftprop, real max_shift, real upper) { real shift = shiftprop * max_shift; if(y <= shift) { return log(mix) + uniform_lcdf(y | 0, upper); } else if(y >= upper) { // The whole uniform part is below, so the mixture part is log(1) = 0 return log_mix(mix, 0, lognormal_lcdf(y - shift | mu, sigma)); } else { real lognormal_llh = lognormal_lcdf(y - shift | mu, sigma); real uniform_llh = uniform_lcdf(y | 0, upper); return log_mix(mix, uniform_llh, lognormal_llh); } } real rtmix_lccdf(real y, real mu, real sigma, real mix, real shiftprop, real max_shift, real upper) { real shift = shiftprop * max_shift; if(y <= shift) { // The whole lognormal part is above, so the mixture part is log(1) = 0 return log_mix(mix, uniform_lccdf(y | 0, upper), 0); } else if(y >= upper) { return log1m(mix) + lognormal_lccdf(y - shift | mu, sigma); } else { real lognormal_llh = lognormal_lccdf(y - shift | mu, sigma); real uniform_llh = uniform_lccdf(y | 0, upper); return log_mix(mix, uniform_llh, lognormal_llh); } } ") # Model components -------------------------------------------------------- posterior_predict_rtmix <- function(i, prep, ...) { if ((!is.null(prep$data$lb) && prep$data$lb[i] > 0) || (!is.null(prep$data$ub) && prep$data$ub[i] < Inf)) { stop("Predictions for truncated distributions not supported") } mu <- brms:::get_dpar(prep, "mu", i = i) sigma <- brms:::get_dpar(prep, "sigma", i = i) mix <- brms:::get_dpar(prep, "mix", i = i) shiftprop <- brms:::get_dpar(prep, "shiftprop", i = i) max_shift <- prep$data$vreal1[i] upper <- prep$data$vreal2[i] shift <- shiftprop * max_shift rtmix(prep$nsamples, meanlog = mu, sdlog = sigma, mix = mix, shift = shift, upper = upper ) } # Needed for numerical stability # from http://tr.im/hH5A logsumexp <- function(x) { y <- max(x) y + log(sum(exp(x - y))) } rtmix_lpdf <- function(y, meanlog, sdlog, mix, shift, upper) { unif_llh <- dunif(y, min = 0, max = upper, log = TRUE) lognormal_llh <- dlnorm(y - shift, meanlog = meanlog, sdlog = sdlog, log = TRUE) - plnorm(upper - shift, meanlog = meanlog, sdlog = sdlog, log.p = TRUE) # Computing logsumexp(log(mix) + unif_llh, log1p(-mix) + lognormal_llh) # but vectorized llh_matrix <- array(NA_real_, dim = c(2, max(length(unif_llh), length(lognormal_llh)))) llh_matrix[1, ] <- log(mix) + unif_llh llh_matrix[2, ] <- log1p(-mix) + lognormal_llh apply(llh_matrix, MARGIN = 2, FUN = logsumexp) } log_lik_rtmix <- function(i, draws) { mu <- brms:::get_dpar(draws, "mu", i = i) sigma <- brms:::get_dpar(draws, "sigma", i = i) mix <- brms:::get_dpar(draws, "mix", i = i) shiftprop <- brms:::get_dpar(draws, "shiftprop", i = i) max_shift <- draws$data$vreal1[i] upper <- draws$data$vreal2[i] shift <- shiftprop * max_shift y <- draws$data$Y[i] rtmix_lpdf(y, meanlog = mu, sdlog = sigma, mix = mix, shift = shift, upper = upper ) } # Test -------------------------------------------------------------------- library(dplyr) #> #> Attaching package: 'dplyr' #> The following objects are masked from 'package:stats': #> #> filter, lag #> The following objects are masked from 'package:base': #> #> intersect, setdiff, setequal, union library(ggplot2) # Generate Data generate_RT <- function(n, meanlog, sdlog, mix, shift, upper) { ifelse(runif(n) < mix, runif(n, 0, upper), shift + rlnorm(n, meanlog = meanlog, sdlog = sdlog) ) } # Parameters n <- 5 n_obs <- 500 max_shift <- runif(n, 0.25, 0.5) shift <- runif(n) * max_shift upper <- runif(n, 8, 10) mix <- runif(n, 0, 0.2) intercept <- runif(n, 0.2, 1) beta <- abs(rnorm(n, 0.5, 0.5)) sigma <- abs(rnorm(n, 0.5, 0.2)) df <- data.frame() for(i in 1:n){ X <- rnorm(n_obs) mu <- rep(intercept[i], n_obs) + beta[i] * X df <- rbind(df, data.frame( RT = generate_RT(n = n_obs, meanlog = mu, sdlog = sigma[i], mix = mix[i], shift = shift[i], upper = upper[i]), x = X, max_shift = max_shift[i], upper = upper[i], participant = as.factor(i))) df$RT[df$RT > 10] <- NA } ggplot(df, aes(x = RT, color = participant)) + geom_density() #> Warning: Removed 33 rows containing non-finite values (stat_density). ```  ``` r f <- brmsformula( # RT | vreal(max_shift, upper) ~ x RT | vreal(max_shift, upper) ~ x + (1|participant) # beta ~ 1 + (1|G|Participant), # beta is the "tau" # sigma ~ 1 + (1|G|Participant) ) model <- brms::brm(f, data = df, family = rtmix, stanvars = stan_functions, refresh = 0, iter = 500, prior = c(brms::prior(beta(1, 5), class = "mix"))) #> Warning: Rows containing NAs were excluded from the model. #> Compiling Stan program... #> Start sampling #> Warning: There were 3 divergent transitions after warmup. See #> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup #> to find out why this is a problem and how to eliminate them. #> Warning: Examine the pairs() plot to diagnose sampling problems #> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#bulk-ess #> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#tail-ess # model # insight::get_priors(model) # bayestestR::describe_prior(model) # parameters::parameters(model) rand <- as.data.frame(parameters::parameters(model, effects = "random")) #> Warning: Note that the default rope range for linear models might change in #> future versions (see https://github.com/easystats/bayestestR/issues/364).Please #> set it explicitly to preserve current results. #> Warning: Could not estimate a good default ROPE range. Using 'c(-0.1, 0.1)'. rand #> Parameter Median CI CI_low CI_high pd #> 3 r_participant.1.Intercept. -0.14191265 0.89 -0.41190957 0.12318264 0.853 #> 4 r_participant.2.Intercept. -0.32425069 0.89 -0.61386007 -0.07668998 0.973 #> 5 r_participant.3.Intercept. 0.04945409 0.89 -0.24066854 0.29530996 0.635 #> 6 r_participant.4.Intercept. 0.05624344 0.89 -0.18757755 0.34616650 0.647 #> 7 r_participant.5.Intercept. 0.36355684 0.89 0.08977932 0.63979605 0.970 #> 8 sd_participant__Intercept 0.31869009 0.89 0.15314769 0.58594382 1.000 #> 1 b_Intercept NA 0.89 NA NA NA #> 2 b_x NA 0.89 NA NA NA #> ROPE_Percentage Rhat ESS Prior_Distribution Prior_df Prior_Location #> 3 0.310 1.032571 266.5525 NA NA #> 4 0.052 1.032842 263.6838 NA NA #> 5 0.468 1.031753 262.2402 NA NA #> 6 0.474 1.031781 269.0068 NA NA #> 7 0.057 1.031467 266.2673 NA NA #> 8 0.000 1.024103 110.6630 NA NA #> 1 NA NA NA student_t 3 2.1 #> 2 NA NA NA uniform NA NA #> Prior_Scale #> 3 NA #> 4 NA #> 5 NA #> 6 NA #> 7 NA #> 8 NA #> 1 2.5 #> 2 NA ``` ^{Created on 2021-04-03 by the [reprex package](https://reprex.tidyverse.org) (v1.0.0)}

[strengejacke]

Can you check if this is still an issue? At least rownames should be in correct order.

[DominiqueMakowski]

It doesn't return the two last rows, and the row names are correct :)