How to get the same output each time I run the code

[JFormoso]

Hi everybody!

I am trying to fix the output of my model, since I get a different result each time I run it and set.seed does not seem to work.
Because I want to share the script and data set with my manuscript, it would be nice to get the exact results each time or have a reasonable justification I can add to it.

This is an example of what I get:

Family: bernoulli(link = 'logit')
Iterations: 9000 from 3 chains.
Segments:
  1: accuracy_enumeration ~ numerosity
  2: accuracy_enumeration ~ 1 + (1 | id) ~ 0 + numerosity

Population-level parameters:
         name  mean lower upper Rhat n.eff
         cp_1  4.38  2.79  5.89  1.2    39
      cp_1_sd  7.83  0.38 24.98  1.1   233
        int_1  6.10  3.68  7.72  1.1    55
 numerosity_1 -0.34 -0.78  0.65  1.2    19
 numerosity_2 -1.21 -1.44 -0.96  1.0   247
#----------

Family: bernoulli(link = 'logit')
Iterations: 9000 from 3 chains.
Segments:
  1: accuracy_enumeration ~ numerosity
  2: accuracy_enumeration ~ 1 + (1 | id) ~ 0 + numerosity

Population-level parameters:
         name  mean lower upper Rhat n.eff
         cp_1  4.40  2.73  6.19  1.2    26
      cp_1_sd  9.99  0.36 31.34  1.1    28
        int_1  6.42  4.40  8.34  1.1    47
 numerosity_1 -0.41 -0.92  0.37  1.2    23
 numerosity_2 -1.23 -1.51 -0.97  1.1   101
#----------

Family: bernoulli(link = 'logit')
Iterations: 9000 from 3 chains.
Segments:
  1: accuracy_enumeration ~ numerosity
  2: accuracy_enumeration ~ 1 + (1 | id) ~ 0 + numerosity

Population-level parameters:
         name  mean lower  upper Rhat n.eff
         cp_1  4.49  2.87  5.934  1.1    36
      cp_1_sd  9.76  0.50 26.737  1.0   336
        int_1  6.41  4.53  8.064  1.0    70
 numerosity_1 -0.46 -0.85  0.085  1.2    30
 numerosity_2 -1.21 -1.46 -0.979  1.0   256
#----------

Family: bernoulli(link = 'logit')
Iterations: 9000 from 3 chains.
Segments:
  1: accuracy_enumeration ~ numerosity
  2: accuracy_enumeration ~ 1 + (1 | id) ~ 0 + numerosity

Population-level parameters:
         name  mean lower upper Rhat n.eff
         cp_1  4.36  2.82  5.80  1.1    38
      cp_1_sd  7.82  0.41 24.09  1.1   121
        int_1  6.16  4.19  7.93  1.2    65
 numerosity_1 -0.36 -0.78  0.47  1.3    17
 numerosity_2 -1.21 -1.43 -0.97  1.0   280

Any help you can give me will be welcome!
Thanks in advance!

Jesi

[lindeloev]

Thanks! This may be possible, so I'll look into it.

In the meantime, your main issue is the low effective sample rate. Inspect via plot_pars(fit) whether this is a convergence issue (increase adapt to e.g. 5000) or a low-identifiability issue (set iter to e.g. 15000). There's a good chance that this will make the results replicable at least to one digit precision.

[JFormoso]

Thank you very much! I'll try that.