# Confusing result of brms predict() and emmeans() for the ordinal continuation model **URL:** \ **Category:** brms\ **Tags:** emmeans\ **Created:** [January 29, 2021, 12:50pm UTC](https://discourse.mc-stan.org/t/confusing-result-of-brms-predict-and-emmeans-for-the-ordinal-continuation-model/20505 "2021-01-29T12:50:43Z")\ **Posts on this page:** 2\ **Page:** 1

**Author:** ![petrovska](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/petrovska/32/8375_2.png) [@petrovska](https://discourse.mc-stan.org/u/petrovska)\ **Post date:** [January 29, 2021, 12:50pm UTC](https://discourse.mc-stan.org/t/confusing-result-of-brms-predict-and-emmeans-for-the-ordinal-continuation-model/20505/1 "2021-01-29T12:50:43Z")

Dear all, Could someone advise if the ordinal continuation ratio model accounts for the frequency of response when estimating its conditional probability? The count (response) and its frequency are low on all plots for some species (explanatory variable), nevertheless the model assigns almost 100% probability for Acer in response 1. At the same time, it assesses lower probability for Fagus that is frequent and abundant. Moreover, the number of species in 20 cm dbh class is low and thus, the line should decline but it is going up. In contrast, Fagus has a higher probability for the response 2 and 3, which I assume is the correct interpretation of the observed data. ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/5/56cf9a904a836699d83783c59a409d0f441502f5.png) . Besides, I am puzzled with marginal mean in `emmeans`. The difference in response is huge. I cannot figure out what threshold was used in this case. I would appreciate any hints or advice. ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/e14f0a6a294590116fc5f4ae5f4a2883657b536f.png) ``` inits <- list(`Intercept[1]` = 0.19, `Intercept[2]` = 1.33, `Intercept[3]` = 2.15, `Intercept[4]` = 2.92, `Intercept[5]` = 3.69, `Intercept[6]` = 4.34, `Intercept[7]`= 5.26, `Intercept[8]`=6.7) inits_list <- list(inits, inits) p= c(prior_string("student_t(5,0,0.5)", class="Intercept"), prior(normal(0,0.5), class=b), prior(student_t(5,0.1,0.5), class=sd)) f=brmsformula(count1 ~ ba_sqrt + dbh_log + species.p + species.m + position.m + status.p + y + offset(log(count_dead)) + (y|subplot_id), center=TRUE) + lf(disc ~ 0 + status.p+y+species.p, cmc = FALSE) fit4=brm(f, data=d1, family = cratio(link="logit"), save_pars=save_pars(all = TRUE), cores = 2, iter = 800 + 800, warmup = 800, chains = 2, seed=123, inits = inits_list, sample_prior="yes", silent=TRUE, open_progress=FALSE, control = list(adapt_delta = 0.85, max_treedepth = 10)) nd <-data.frame(species.p=c("Acer platanoides", "Fagus sylvatica"), species.m=d1$species.m, position.m=d1$position.m, ba_sqrt=seq(2, 4, by=0.2), dbh_log = seq(2.3, 3, by=0.07), count_dead=1, y="2015", status.p="alive", subplot_id=d1$subplot_id) View(nd) max_iter <- 100 p=as.data.frame(predict(fit4, newdata = nd, subset = 1:max_iter, summary = T, re_formula = NULL, scale = "response")) together=cbind(nd$species.p, nd$dbh_log, p) together$gr=rep(1:368, times=1, each=11) prep=melt(together, id=c("nd$species.p", "nd$dbh_log", "gr")) prep$dbh_exp=exp(nd$dbh_log) names(prep)[names(prep) == "nd$species.p"] <- "species.p" prep %>% group_by(species.p, variable, dbh_exp, gr) %>% distinct(value) ggplot(subset(prep, variable %in% c("P(Y = 1)", "P(Y = 2)", "P(Y = 3)")), aes(dbh_exp, value, group=gr, color=species.p), alpha=0.2)+ stat_summary(fun = mean, geom = "line", aes(group=species.p))+ facet_wrap(.~variable, scale="free")+theme(legend.position = "top") ``` - Operating System: Windows 10 - brms Version: 2.14.6 ---

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\ **Post date:** [February 5, 2021, 4:49pm UTC](https://discourse.mc-stan.org/t/confusing-result-of-brms-predict-and-emmeans-for-the-ordinal-continuation-model/20505/2 "2021-02-05T16:49:30Z")

Sorry for not getting to your question earlier. It is relevant and on-topic. > [@petrovska](#): > > The count (response) and its frequency are low on all plots for some species (explanatory variable), nevertheless the model assigns almost 100% probability for Acer in response 1. At the same time, it assesses lower probability for Fagus that is frequent and abundant. That sounds correct to me - if the counts are low, then the model should assign high probability to `Y = 1`, i.e. that the count is the lowest possible. If the counts are high, than assigning low probability of `Y = 1` sounds correct. > [@petrovska](#): > > Moreover, the number of species in 20 cm dbh class is low and thus, the line should decline but it is going up. That - once again - looks correct to me: the probabilities for individual count response cannot be interpreted very directly - increase in probability that `Y = 1` combined with decrease in probability that `Y = 2` is likely and indication of decrease in average count. You might want to use `posterior_epred` to get you predictions of the posterior mean directly. Additionally, the model is quite complex and I honestly don’t fully understand what is going on. I think there are two simple explanations that also cannot be ruled out immediately: 1. The model is a poor fit to the data. This could be easily investigated with some posterior predictive checks. I would guess that things like `pp_check(fit4, type = "bars_grouped", group = "species.p")` `pp_check(fit4, type = "stat_grouped", stat = "sd", group = "species.p")` could show some potential problems. 2. There is a bug in how you handle the predictions. I am for example somewhat suspicious of the use of `distinct` n the plot preparation, but I don’t really understand what is going on there, so I am likely to be wrong. > [@petrovska](#): > > Besides, I am puzzled with marginal mean in `emmeans` . The difference in response is huge. I cannot figure out what threshold was used in this case. I would appreciate any hints or advice. Unfortunately `emmeans` isn’t part of my standard toolkit, so can’t really help there. Is this comparable to results with `conditional_effects` from `brms` (which I have more experience with and could thus help)? Best of luck with your model!