Ordinal response with monotonic predictor

Strange. Can you paste us the output from this:

get_prior(bodn ~ mo(chf) + (1 | id), family = cumulative, data = bodn)
getprior
                 prior     class               coef 
1                              b                                               
2                              b    mochf                            
3  student_t(3, 0, 10) Intercept                                               
4                      Intercept                  1                            
5                      Intercept                 10                            
6                      Intercept                 11                            
7                      Intercept                 12                            
8                      Intercept                 13                            
9                      Intercept                 14                            
10                     Intercept                  2                            
11                     Intercept                  3                            
12                     Intercept                  4                            
13                     Intercept                  5                            
14                     Intercept                  6                            
15                     Intercept                  7                            
16                     Intercept                  8                            
17                     Intercept                  9                            
18 student_t(3, 0, 10) sd                                               
19                               sd                          idno                      
20                               sd          Intercept  idno                      
21                               simo      mochf1 
p <- get_prior(bodn ~ mo(chf) + (1 | id), family = cumulative, data = bodn)
p$prior[1] <- "normal(0,1)" # sets a N(0,1) on \beta
p$prior[3] <- "normal(0,2.5)" # for each cutpoint
p$prior[18] <- "cauchy(0,2.5)" # for sd
p$prior[21] <- "dirichlet(2,2)" # sets weakly regularizing dirichlet for 1-3 Likert scale response

Then use prior = p when you call brms. Of course, at least do pp_check() when you sample with sample_prior="only" to see that the priors are sane.

I try the model again and get back here to update.
-Thanks

The problem is I do not get the similar pattern as you wrote:
The response is as blow:

 p$prior[1]
[1] ""
> p$prior[3] 
[1] "student_t(3, 0, 10)"
> p$prior[1] 
[1] ""
> p$prior[3] 
[1] "student_t(3, 0, 10)"
> p$prior[18]
[1] "student_t(3, 0, 10)"
> p$prior[21] 
[1] ""

:(
I can send some of the data, if it helps to understand why!

I set the first row in p to “normal(0,1)”. What happens when you cut and paste exactly what I have written above (I changed things to better fit your scenario) and then simply write p in the terminal?

I ran the model using the above priors. The good thing is that I got simplex results. However, the Rhats are not good.
here is the results

Data: bodi (Number of observations: 2622) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~idno (Number of levels: 898) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     4.13      0.14     3.86     4.41 1.01      457     1116

Population-Level Effects: 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1]        -11.12      0.50   -12.11   -10.17 1.01      949     1834
Intercept[2]         -8.07      0.30    -8.67    -7.51 1.01      703     1816
Intercept[3]         -5.67      0.21    -6.10    -5.28 1.02      422     1057
Intercept[4]         -3.56      0.18    -3.92    -3.20 1.03      262      803
Intercept[5]         -1.82      0.17    -2.15    -1.50 1.04      219      514
Intercept[6]         -0.06      0.16    -0.38     0.24 1.04      125      498
Intercept[7]          1.40      0.16     1.08     1.71 1.04      130      560
Intercept[8]          2.97      0.17     2.64     3.30 1.03      172      687
Intercept[9]          4.49      0.19     4.13     4.87 1.03      259      983
Intercept[10]         5.97      0.21     5.57     6.40 1.02      373     1267
Intercept[11]         7.60      0.25     7.14     8.08 1.02      554     1719
Intercept[12]         9.49      0.31     8.92    10.10 1.01      650     2034
Intercept[13]        11.21      0.41    10.43    12.04 1.01     1412     2826
Intercept[14]        13.22      0.66    12.04    14.59 1.00     2611     2885
       mochf          0.50      0.08     0.34     0.65 1.00      915     1622

Simplex Parameters: 
                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
      mochf1[1]        0.16      0.09     0.03     0.36 1.00     2083     2314
      mochf1[2]        0.84      0.09     0.64     0.97 1.00     2083     2314

Samples were drawn using sampling(NUTS). 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).

The good news is I got the simplex result but not good Rhat
Should I change the priors’ values?
Thanks

I modified the prior for cut-points to (0,5) and reran the syntax and it worked well with reasonable Rhat and convergence.

Group-Level Effects: 
~idno (Number of levels: 898) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     4.32      0.15     4.04     4.62 1.00      711     1152

Population-Level Effects: 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1]        -11.74      0.54   -12.87   -10.72 1.00     1560     2403
Intercept[2]         -8.35      0.31    -8.97    -7.74 1.00     1307     2258
Intercept[3]         -5.80      0.23    -6.25    -5.36 1.00      859     1851
Intercept[4]         -3.61      0.19    -3.98    -3.23 1.01      639     1427
Intercept[5]         -1.81      0.18    -2.16    -1.46 1.01      560     1154
Intercept[6]         -0.01      0.17    -0.34     0.33 1.01      528     1297
Intercept[7]          1.50      0.17     1.17     1.85 1.01      529     1133
Intercept[8]          3.12      0.18     2.78     3.49 1.01      578     1077
Intercept[9]          4.70      0.20     4.32     5.08 1.01      677     1627
Intercept[10]         6.25      0.22     5.82     6.69 1.01      830     1734
Intercept[11]         7.98      0.26     7.48     8.50 1.01     1071     2221
Intercept[12]        10.02      0.33     9.38    10.67 1.00     1736     2649
Intercept[13]        12.00      0.45    11.16    12.90 1.00     2733     2877
Intercept[14]        14.64      0.87    13.10    16.56 1.00     4483     2673
        mochf         0.52      0.08     0.35     0.68 1.00     2247     2494

Simplex Parameters: 
                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
          mochf1[1]     0.17      0.09     0.03     0.36 1.00     3518     3121
          mochf1[2]     0.83      0.09     0.64     0.97 1.00     3518     3121

Thanks for your help and guidance,

3 Likes

Hi,

I’d still be careful. Have you run the model with sample_prior = "only" and used pp_check() to see how the priors look like?

Make sure to check that \widehat{R} < 1.01, that neff is more than 0.1-0.2, and that the trace plots look alright. It seems you have a Bulk_ESS which is close to 0.1 (if you ran with four chains each with 2000 iterations).

1 Like

Blockquote I’d still be careful. Have you run the model with sample_prior = "only" and used pp_check() to see how the priors look like?

Blockquote
We expect the weak priors and therefore a weak posterior coefficient, right?

Blockquote brms

Ordinal response with monotonic predictor

Interfacesbrms

](Profile - ssalimi - The Stan Forums)

ssalimi

torkar

18h

The problem is I do not get the similar pattern as you wrote:
The response is as blow:

 p$prior[1]
[1] ""
> p$prior[3] 
[1] "student_t(3, 0, 10)"
> p$prior[1] 
[1] ""
> p$prior[3] 
[1] "student_t(3, 0, 10)"
> p$prior[18]
[1] "student_t(3, 0, 10)"
> p$prior[21] 
[1] ""

:(
I can send some of the data, if it helps to understand why!

SolutionReply

torkar

18h

I set the first row in p to “normal(0,1)”. What happens when you cut and paste exactly what I have written above (I changed things to better fit your scenario) and then simply write p in the terminal?

SolutionReply

ssalimi

17h

I ran the model using the above priors. The good thing is that I got simplex results. However, the Rhats are not good.
here is the results

Data: bodi (Number of observations: 2622) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~idno (Number of levels: 898) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     4.13      0.14     3.86     4.41 1.01      457     1116

Population-Level Effects: 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1]        -11.12      0.50   -12.11   -10.17 1.01      949     1834
Intercept[2]         -8.07      0.30    -8.67    -7.51 1.01      703     1816
Intercept[3]         -5.67      0.21    -6.10    -5.28 1.02      422     1057
Intercept[4]         -3.56      0.18    -3.92    -3.20 1.03      262      803
Intercept[5]         -1.82      0.17    -2.15    -1.50 1.04      219      514
Intercept[6]         -0.06      0.16    -0.38     0.24 1.04      125      498
Intercept[7]          1.40      0.16     1.08     1.71 1.04      130      560
Intercept[8]          2.97      0.17     2.64     3.30 1.03      172      687
Intercept[9]          4.49      0.19     4.13     4.87 1.03      259      983
Intercept[10]         5.97      0.21     5.57     6.40 1.02      373     1267
Intercept[11]         7.60      0.25     7.14     8.08 1.02      554     1719
Intercept[12]         9.49      0.31     8.92    10.10 1.01      650     2034
Intercept[13]        11.21      0.41    10.43    12.04 1.01     1412     2826
Intercept[14]        13.22      0.66    12.04    14.59 1.00     2611     2885
       mochf          0.50      0.08     0.34     0.65 1.00      915     1622

Simplex Parameters: 
                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
      mochf1[1]        0.16      0.09     0.03     0.36 1.00     2083     2314
      mochf1[2]        0.84      0.09     0.64     0.97 1.00     2083     2314

Samples were drawn using sampling(NUTS). 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).

The good news is I got the simplex result but not good Rhat
Should I change the priors’ values?
Thanks

SolutionReply

ssalimi

15h

I modified the prior for cut-points to (0,5) and reran the syntax and it worked well with reasonable Rhat and convergence.

Group-Level Effects: 
~idno (Number of levels: 898) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     4.32      0.15     4.04     4.62 1.00      711     1152

Population-Level Effects: 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1]        -11.74      0.54   -12.87   -10.72 1.00     1560     2403
Intercept[2]         -8.35      0.31    -8.97    -7.74 1.00     1307     2258
Intercept[3]         -5.80      0.23    -6.25    -5.36 1.00      859     1851
Intercept[4]         -3.61      0.19    -3.98    -3.23 1.01      639     1427
Intercept[5]         -1.81      0.18    -2.16    -1.46 1.01      560     1154
Intercept[6]         -0.01      0.17    -0.34     0.33 1.01      528     1297
Intercept[7]          1.50      0.17     1.17     1.85 1.01      529     1133
Intercept[8]          3.12      0.18     2.78     3.49 1.01      578     1077
Intercept[9]          4.70      0.20     4.32     5.08 1.01      677     1627
Intercept[10]         6.25      0.22     5.82     6.69 1.01      830     1734
Intercept[11]         7.98      0.26     7.48     8.50 1.01     1071     2221
Intercept[12]        10.02      0.33     9.38    10.67 1.00     1736     2649
Intercept[13]        12.00      0.45    11.16    12.90 1.00     2733     2877
Intercept[14]        14.64      0.87    13.10    16.56 1.00     4483     2673
        mochf         0.52      0.08     0.35     0.68 1.00     2247     2494

Simplex Parameters: 
                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
          mochf1[1]     0.17      0.09     0.03     0.36 1.00     3518     3121
          mochf1[2]     0.83      0.09     0.64     0.97 1.00     3518     3121

Thanks for your help and guidance,

Solution

2

Reply

torkar

10h

Hi,

I’d still be careful. Have you run the model with sample_prior = "only" and used pp_check() to see how the priors look like?

Make sure to check that ˆR < 1.01, that neff is more than 0.1-0.2, and that the trace plots look alright. It seems you have a Bulk_ESS which is close to 0.1 (if you ran with four chains each with 2000 iterations).

Blockquote
sample_prior density_overlay

Hi,

please post the output of this:

pp_check(fit_bodn, type = "bars")
1 Like

please post the output of this:

pp_check(fit_bodn, type = "bars")

PPcheck_bar

@torkar
just one person has bodn=15, is this a problem?

Your priors are a bit heavy on low and high numbers, while they don’t cover, e.g., category 2 very much. I’d like to see a more balance view.

I presume your pp_check() on the full model (not with priors only) is ok?

[quote=“torkar, post:18, topic:13165, full:true”]
Your priors are a bit heavy on low and high numbers, while they don’t cover, e.g., category 2 very much. I’d like to see a more balance view.

I presume your pp_check() on the full model (not with priors only) is ok?

It seems so, this is the pp_check (fit_bodn) for the full model!
What do you think? Thanks a lot for all your time. I learned many points!
bodn_chf_ppCheck

I think it looks ok :) I would now compare that model with a model with no monotonic effects,

bodn ~ chf + (1 | id)

and a general intercept model,

bodn ~ 1

using LOO: https://rdrr.io/cran/brms/man/loo_compare.brmsfit.html

That is however completely up to you :)

I will do in a few hours and post here.

Thanks

I did compare four models
fit1<-bodn ~ 1
fit2<-bodn ~ chf + (1 | id)
fit3<-bodn ~ chf + (1 | id), family = cumulative, data = bodn)
fit4<-p <bodn ~ mo(chf) + (1 | id), family = cumulative, data = bodn)
using
model_weights(fit1,fit2,fit3,fit4,weights=loo)
fit 4 had higher values than others.
But when I used “loo_compare” I got the message that I better use k-fold cross-validation with k=10 instead of LOO.
cross-validation is the next to learn!

Hi

fit1 <- add_criterion(fit, criterion = "loo")
fit2 <- add_criterion(fit, criterion = "loo")
fit3 <- add_criterion(fit, criterion = "loo")
fit4 <- add_criterion(fit, criterion = "loo")

Some of them might need reloo = TRUE. Then you do

loo_compare(fit1, fit2, fit3, fit4)