Closed DominiqueMakowski closed 3 years ago
Also, should ndt_Intercept
, which is a population effect (corresponding to the shift parameter of the shifted lognormal distribution) be moved from the Fixed effects as currently defined by parameters
to another group of auxiliary parameters? @mattansb
In the above it seems that sigma_Intercept (a population effect) gets duplicated as sigma_Intercept and Intercept_sigma.
Not only that, but sigma_Intercept
is also the duplicate name of different parameters (the fixed intercept, and the sd of the random intercept).
Also, should ndt_Intercept, which is a population effect (corresponding to the shift parameter of the shifted lognormal distribution) be moved from the Fixed effects as currently defined by parameters to another group of auxiliary parameters?
I would say so, yes.
Perhaps we need a: location, scale, shape, aux families of parameters? I don't know if we can cover all of the billions of options brms
offers....
that model takes forever to run, could you attach the model as RData/RDS?
Here is what clean_parameters()
returns:
load("~/../Desktop/test.RData")
insight::clean_parameters(model)
#> Parameter Effects Component
#> 1 b_Intercept fixed conditional
#> 2 b_sigma_Intercept fixed conditional
#> 3 b_ndt_Intercept fixed conditional
#> 4 b_conditionspeed fixed conditional
#> 5 r_1.1.1. random conditional
#> 6 r_1.2.1. random conditional
#> 7 r_1.3.1. random conditional
#> 8 r_1.4.1. random conditional
#> 9 r_1.5.1. random conditional
#> 10 r_1.6.1. random conditional
#> 11 r_1.1.2. random conditional
#> 12 r_1.2.2. random conditional
#> 13 r_1.3.2. random conditional
#> 14 r_1.4.2. random conditional
#> 15 r_1.5.2. random conditional
#> 16 r_1.6.2. random conditional
#> 17 r_1.1.3. random conditional
#> 18 r_1.2.3. random conditional
#> 19 r_1.3.3. random conditional
#> 20 r_1.4.3. random conditional
#> 21 r_1.5.3. random conditional
#> 22 r_1.6.3. random conditional
#> 23 r_id.1.Intercept. random conditional
#> 24 r_id.2.Intercept. random conditional
#> 25 r_id.3.Intercept. random conditional
#> 26 r_id.4.Intercept. random conditional
#> 27 r_id.5.Intercept. random conditional
#> 28 r_id.6.Intercept. random conditional
#> 29 r_id__ndt.1.Intercept. random conditional
#> 30 r_id__ndt.2.Intercept. random conditional
#> 31 r_id__ndt.3.Intercept. random conditional
#> 32 r_id__ndt.4.Intercept. random conditional
#> 33 r_id__ndt.5.Intercept. random conditional
#> 34 r_id__ndt.6.Intercept. random conditional
#> 35 r_id__sigma.1.Intercept. random conditional
#> 36 r_id__sigma.2.Intercept. random conditional
#> 37 r_id__sigma.3.Intercept. random conditional
#> 38 r_id__sigma.4.Intercept. random conditional
#> 39 r_id__sigma.5.Intercept. random conditional
#> 40 r_id__sigma.6.Intercept. random conditional
#> 41 sd_id__Intercept random conditional
#> 42 sd_id__ndt_Intercept random conditional
#> 43 sd_id__sigma_Intercept random conditional
#> 44 cor_id__Intercept__ndt_Intercept random conditional
#> 45 cor_id__Intercept__sigma_Intercept random conditional
#> 46 cor_id__ndt_Intercept__sigma_Intercept random conditional
#> 47 b_sigma_Intercept fixed sigma
#> 48 sd_id__sigma_Intercept fixed sigma
#> 49 cor_id__Intercept__sigma_Intercept fixed sigma
#> 50 cor_id__ndt_Intercept__sigma_Intercept fixed sigma
#> 51 Intercept_sigma fixed sigma
#> 52 r_id__sigma.1.Intercept. fixed sigma
#> 53 r_id__sigma.2.Intercept. fixed sigma
#> 54 r_id__sigma.3.Intercept. fixed sigma
#> 55 r_id__sigma.4.Intercept. fixed sigma
#> 56 r_id__sigma.5.Intercept. fixed sigma
#> 57 r_id__sigma.6.Intercept. fixed sigma
#> Group Cleaned_Parameter
#> 1 (Intercept)
#> 2 sigma_Intercept
#> 3 ndt_Intercept
#> 4 conditionspeed
#> 5 1: 1 1.1
#> 6 1: 1 1.2
#> 7 1: 1 1.3
#> 8 1: 1 1.4
#> 9 1: 1 1.5
#> 10 1: 1 1.6
#> 11 2: 1 1.1
#> 12 2: 1 1.2
#> 13 2: 1 1.3
#> 14 2: 1 1.4
#> 15 2: 1 1.5
#> 16 2: 1 1.6
#> 17 3: 1 1.1
#> 18 3: 1 1.2
#> 19 3: 1 1.3
#> 20 3: 1 1.4
#> 21 3: 1 1.5
#> 22 3: 1 1.6
#> 23 Intercept: id id.1
#> 24 Intercept: id id.2
#> 25 Intercept: id id.3
#> 26 Intercept: id id.4
#> 27 Intercept: id id.5
#> 28 Intercept: id id.6
#> 29 Intercept: id__ndt id__ndt.1
#> 30 Intercept: id__ndt id__ndt.2
#> 31 Intercept: id__ndt id__ndt.3
#> 32 Intercept: id__ndt id__ndt.4
#> 33 Intercept: id__ndt id__ndt.5
#> 34 Intercept: id__ndt id__ndt.6
#> 35 Intercept: id__sigma id__sigma.1
#> 36 Intercept: id__sigma id__sigma.2
#> 37 Intercept: id__sigma id__sigma.3
#> 38 Intercept: id__sigma id__sigma.4
#> 39 Intercept: id__sigma id__sigma.5
#> 40 Intercept: id__sigma id__sigma.6
#> 41 SD/Cor: id (Intercept)
#> 42 SD/Cor: id ndt_Intercept
#> 43 SD/Cor: id sigma_Intercept
#> 44 SD/Cor: id Intercept ~ ndt_Intercept
#> 45 SD/Cor: id Intercept ~ sigma_Intercept
#> 46 SD/Cor: id ndt_Intercept ~ sigma_Intercept
#> 47 sigma_Intercept
#> 48 SD/Cor: id sigma_Intercept
#> 49 SD/Cor: id Intercept ~ sigma_Intercept
#> 50 SD/Cor: id ndt_Intercept ~ sigma_Intercept
#> 51 Intercept_sigma
#> 52 Intercept: id__sigma id__sigma.1
#> 53 Intercept: id__sigma id__sigma.2
#> 54 Intercept: id__sigma id__sigma.3
#> 55 Intercept: id__sigma id__sigma.4
#> 56 Intercept: id__sigma id__sigma.5
#> 57 Intercept: id__sigma id__sigma.6
Created on 2021-05-06 by the reprex package (v2.0.0)
The model can be created instantaneously by changing the algo:
library(brms)
library(dplyr)
data <- rtdists::speed_acc %>%
filter(rt < 2, !id %in% as.character(7:17)) %>%
mutate(rt = rt * 1000)
formula <- brms::brmsformula(
rt ~ condition + (1|G|id),
ndt ~ 1 + (1|G|id),
sigma ~ 1 + (1|G|id)
)
model <- brms::brm(formula, data=data, family = shifted_lognormal(), refresh = 0, algorithm = "fixed_param")
#> Start sampling
#> Running MCMC with 1 chain...
#>
#> Chain 1 finished in 0.1 seconds.
parameters::parameters(model)
#> Warning in stats::cor(dat): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> Warning in cor(x, y): the standard deviation is zero
#> # Fixed effects
#>
#> Parameter | Median | 89% CI | pd | % in ROPE | Rhat | ESS
#> --------------------------------------------------------------------------
#> (Intercept) | -1.56 | [-1.56, -1.56] | 100% | 100% | 0.999 | 1.00
#> ndt_Intercept | -0.20 | [-0.20, -0.20] | 100% | 100% | 0.999 | 1.00
#> conditionspeed | -0.03 | [-0.03, -0.03] | 100% | 100% | 0.999 | 1.00
#>
#> # Fixed effects sigma
#>
#> Parameter | Median | 89% CI | pd | % in ROPE | Rhat | ESS
#> ---------------------------------------------------------------------------
#> sigma_Intercept | -0.41 | [-0.41, -0.41] | 100% | 100% | 0.999 | 1.00
#> Intercept_sigma | -0.41 | [-0.41, -0.41] | 100% | 100% | |
#>
#> # Fixed effects sigma Intercept: id__sigma
#>
#> Parameter | Median | 89% CI | pd | % in ROPE | Rhat | ESS
#> ---------------------------------------------------------------------
#> id__sigma.2 | 0.24 | [0.24, 0.24] | 100% | 100% | 0.999 | 1.00
#>
#> # Random effects SD/Cor: id
#>
#> Parameter | Median | 89% CI | pd | % in ROPE | Rhat | ESS
#> -------------------------------------------------------------------------------------------
#> ndt_Intercept ~ sigma_Intercept | -0.99 | [-0.99, -0.99] | 100% | 100% | 0.999 | 1.00
#>
#> # Fixed effects sigma SD/Cor: id
#>
#> Parameter | Median | 89% CI | pd | % in ROPE | Rhat | ESS
#> -------------------------------------------------------------------------------------
#> sigma_Intercept | 0.17 | [0.17, 0.17] | 100% | 100% | 0.999 | 1.00
#> Intercept ~ sigma_Intercept | 0.95 | [0.95, 0.95] | 100% | 100% | 0.999 | 1.00
#>
#> Using highest density intervals as credible intervals.
Created on 2021-05-07 by the reprex package (v1.0.0)
brms returns all those columns, so there is no "buggy" duplication, but just all columns are shown. Maybe we should filter all columns that do not start with "b_" or so?
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
load("C:/Users/mail/Desktop/test.RData")
summary(model)
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 68 divergent transitions after warmup. Increasing
#> adapt_delta above may help. See http://mc-stan.org/misc/warnings.html#divergent-
#> transitions-after-warmup
#> Family: shifted_lognormal
#> Links: mu = identity; sigma = log; ndt = log
#> Formula: rt ~ condition + (1 | G | id)
#> ndt ~ 1 + (1 | G | id)
#> sigma ~ 1 + (1 | G | id)
#> Data: data (Number of observations: 10375)
#> Samples: 2 chains, each with iter = 500; warmup = 250; thin = 1;
#> total post-warmup samples = 500
#>
#> Group-Level Effects:
#> ~id (Number of levels: 6)
#> Estimate Est.Error l-95% CI u-95% CI Rhat
#> sd(Intercept) 0.51 0.14 0.29 0.95 1.23
#> sd(ndt_Intercept) 1.91 0.81 0.92 3.95 1.03
#> sd(sigma_Intercept) 0.38 0.17 0.20 0.85 1.09
#> cor(Intercept,ndt_Intercept) -0.46 0.25 -0.83 0.09 1.13
#> cor(Intercept,sigma_Intercept) -0.63 0.32 -0.93 0.32 1.15
#> cor(ndt_Intercept,sigma_Intercept) 0.26 0.27 -0.33 0.73 1.04
#> Bulk_ESS Tail_ESS
#> sd(Intercept) 8 13
#> sd(ndt_Intercept) 41 61
#> sd(sigma_Intercept) 17 59
#> cor(Intercept,ndt_Intercept) 30 43
#> cor(Intercept,sigma_Intercept) 13 54
#> cor(ndt_Intercept,sigma_Intercept) 33 107
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 6.10 0.20 5.58 6.42 1.70 3 22
#> sigma_Intercept -1.06 0.15 -1.28 -0.68 1.91 3 21
#> ndt_Intercept 4.47 0.66 2.75 5.49 1.27 7 68
#> conditionspeed -0.37 0.01 -0.38 -0.35 1.01 199 339
#>
#> Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
colMeans(as.data.frame(model)[1:13])
#> b_Intercept b_sigma_Intercept
#> 6.0973153 -1.0550569
#> b_ndt_Intercept b_conditionspeed
#> 4.4653645 -0.3692906
#> sd_id__Intercept sd_id__ndt_Intercept
#> 0.5149773 1.9056076
#> sd_id__sigma_Intercept cor_id__Intercept__ndt_Intercept
#> 0.3834334 -0.4649650
#> cor_id__Intercept__sigma_Intercept cor_id__ndt_Intercept__sigma_Intercept
#> -0.6265451 0.2577520
#> Intercept Intercept_sigma
#> 5.9117270 -1.0550569
#> Intercept_ndt
#> 4.4653645
Created on 2021-05-07 by the reprex package (v2.0.0)
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
load("C:/Users/mail/Desktop/test.RData")
summary(model)
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#> Warning: There were 68 divergent transitions after warmup. Increasing
#> adapt_delta above may help. See http://mc-stan.org/misc/warnings.html#divergent-
#> transitions-after-warmup
#> Family: shifted_lognormal
#> Links: mu = identity; sigma = log; ndt = log
#> Formula: rt ~ condition + (1 | G | id)
#> ndt ~ 1 + (1 | G | id)
#> sigma ~ 1 + (1 | G | id)
#> Data: data (Number of observations: 10375)
#> Samples: 2 chains, each with iter = 500; warmup = 250; thin = 1;
#> total post-warmup samples = 500
#>
#> Group-Level Effects:
#> ~id (Number of levels: 6)
#> Estimate Est.Error l-95% CI u-95% CI Rhat
#> sd(Intercept) 0.51 0.14 0.29 0.95 1.23
#> sd(ndt_Intercept) 1.91 0.81 0.92 3.95 1.03
#> sd(sigma_Intercept) 0.38 0.17 0.20 0.85 1.09
#> cor(Intercept,ndt_Intercept) -0.46 0.25 -0.83 0.09 1.13
#> cor(Intercept,sigma_Intercept) -0.63 0.32 -0.93 0.32 1.15
#> cor(ndt_Intercept,sigma_Intercept) 0.26 0.27 -0.33 0.73 1.04
#> Bulk_ESS Tail_ESS
#> sd(Intercept) 8 13
#> sd(ndt_Intercept) 41 61
#> sd(sigma_Intercept) 17 59
#> cor(Intercept,ndt_Intercept) 30 43
#> cor(Intercept,sigma_Intercept) 13 54
#> cor(ndt_Intercept,sigma_Intercept) 33 107
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 6.10 0.20 5.58 6.42 1.70 3 22
#> sigma_Intercept -1.06 0.15 -1.28 -0.68 1.91 3 21
#> ndt_Intercept 4.47 0.66 2.75 5.49 1.27 7 68
#> conditionspeed -0.37 0.01 -0.38 -0.35 1.01 199 339
#>
#> Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
parameters::model_parameters(model, effects = "all", centrality = "mean")
#> # Fixed effects
#>
#> Parameter | Mean | 89% CI | pd | % in ROPE | Rhat | ESS
#> ----------------------------------------------------------------------------
#> (Intercept) | 6.10 | [ 5.83, 6.44] | 100% | 100% | 1.572 | 3.00
#> sigma_Intercept | -1.06 | [-1.28, -0.86] | 100% | 100% | 1.799 | 2.00
#> ndt_Intercept | 4.47 | [ 3.50, 5.50] | 100% | 100% | 1.252 | 17.00
#> conditionspeed | -0.37 | [-0.38, -0.36] | 100% | 100% | 1.008 | 186.00
#>
#> # Random effects SD/Cor: id
#>
#> Parameter | Mean | 89% CI | pd | % in ROPE | Rhat | ESS
#> ---------------------------------------------------------------------------------------------
#> (Intercept) | 0.51 | [ 0.27, 0.64] | 100% | 100% | 1.180 | 26.00
#> ndt_Intercept | 1.91 | [ 0.92, 3.08] | 100% | 100% | 1.022 | 45.00
#> sigma_Intercept | 0.38 | [ 0.19, 0.59] | 100% | 100% | 1.080 | 19.00
#> Intercept ~ ndt_Intercept | -0.46 | [-0.84, -0.15] | 95.60% | 100% | 1.131 | 36.00
#> Intercept ~ sigma_Intercept | -0.63 | [-0.96, -0.13] | 93.00% | 100% | 1.164 | 21.00
#> ndt_Intercept ~ sigma_Intercept | 0.26 | [-0.10, 0.73] | 82.80% | 100% | 1.040 | 36.00
#>
#> Using highest density intervals as credible intervals.
Created on 2021-05-07 by the reprex package (v2.0.0)
In the above it seems that
sigma_Intercept
(a population effect) gets duplicated assigma_Intercept
andIntercept_sigma
.The default summaries goes like this:
Created on 2021-05-04 by the reprex package (v1.0.0)