Strange. Can you paste us the output from this:
get_prior(bodn ~ mo(chf) + (1 | id), family = cumulative, data = bodn)
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,
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).
Blockquote I’d still be careful. Have you run the model with
sample_prior = "only"
and usedpp_check()
to see how the priors look like?
Blockquote
We expect the weak priors and therefore a weak posterior coefficient, right?
Blockquote brms
](Profile - ssalimi - The Stan Forums)
torkar
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
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
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
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
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
Hi,
please post the output of this:
pp_check(fit_bodn, type = "bars")
please post the output of this:
pp_check(fit_bodn, type = "bars")
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!
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)