"acat" and "cumulative" outputting different predictions?

[Peter Phalen]

I seem to be getting very different predictions when I fit a multilevel model using the "acat" family rather than the "cumulative" family.

I'm trying to fit a multilevel model where the outcome is a 3-level ordinal category. I'm estimating whether higher scores on a measure of a theoretical construct predict higher levels of the outcome variable, and I'm particularly interested in variability by subscale.

The code I'm using:

model <- brm(ordered_outcome ~ score + (1 + score | measure_subscale),

                            data = dt.long, iter=4000, cores=4, chains=4, control = list(adapt_delta = .975), family="cumulative")

predictions <- predict(model) # get predictions for the input data

predictions <- apply(model,1,which.max) # get category with the highest probability for each case

The above code gives sensible predictions that roughly match what I get by fitting a series of unpooled proportional odds logistic regressions for each of the different subscales. However, when I run the same code with family set to "acat" (rather than "cumulative") I get very different predictions that don't really match the data. For example, the majority of predictions end up in favor of the category that has the lowest observed count...

Perhaps I'm doing something wrong in processing the predictions?

Thank you!

[Paul Buerkner]

Hmm this sounds strange. Maybe there is a problem with the prediction of 'acat'. Would you mind writing a short reproducible example of your problem?

--
You received this message because you are subscribed to the Google Groups "brms-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.
To post to this group, send email to brms-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/c1a5d283-af46-4cde-ab5e-93f6802b8def%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Peter Phalen]

Sure. This is kind of a janky way of simulating the data but it illustrates the problem:

##### simulate

set.seed(1)

# the majority of outcomes are category "1"

ordered_outcome <- c(rep(1,300),rep(2,150),rep(3,75))

measure_subscale <- c(rep(1,100),rep(2,100),rep(3,100),

                      rep(1,50),rep(2,50),rep(3,50),

                      rep(1,25),rep(2,25),rep(3,25))

# crude way of simulating a relationship between subscales and outcome

score <- rnorm(length(ordered_outcome), ordered_outcome * measure_subscale, .5)

ordered_outcome <- as.ordered(ordered_outcome)

measure_subscale <- factor(measure_subscale, levels=1:3, labels=c("A", "B", "C"))

dt.lng <- data.frame(ordered_outcome, measure_subscale, score)

## Draw a histogram of the density of scores on each measure subscale by outcome

densplot <- dt.lng %>% 

  ggplot() + 

  aes(x = score, fill=ordered_outcome) +

  geom_density(show.legend = T) + 

  xlab("score on subscales A, B, and C") + ylab("frequency")

densplot + facet_wrap(~measure_subscale)

cumulative_model <- brm(ordered_outcome ~ score + (1 + score | measure_subscale),

             data = dt.lng, iter=500, cores=4, chains=4, family="cumulative")

cumulative_predictions <- predict(cumulative_model) # get predictions for the input data

cumulative_predictions <- apply(cumulative_predictions,1,which.max) # get category with the highest probability for each case

dt.lng <- cbind(dt.lng, cumulative_predictions)

p <- ggplot(dt.lng) +

  aes(x = cumulative_predictions, y = ordered_outcome, col=score) +

  geom_jitter(height=.1,width=.1,alpha=.6) +

  scale_colour_gradient(low="blue", high="red") +

  ylab("Outcome; as observed") +

  labs(col="Observed\nscore", title="family=cumulative (this model looks appropriate)") +

  xlab("Outcome; as predicted")

p + facet_wrap(~measure_subscale) 

acat_model <- brm(ordered_outcome ~ score + (1 + score | measure_subscale),

                        data = dt.lng, iter=500, cores=4, chains=4, family="acat")

acat_predictions <- predict(acat_model) # get predictions for the input data

acat_predictions <- apply(acat_predictions,1,which.max) # get category with the highest probability for each case

dt.lng <- cbind(dt.lng, acat_predictions)

p <- ggplot(dt.lng) +

  aes(x = acat_predictions, y = ordered_outcome, col=score) +

  geom_jitter(height=.1,width=.1,alpha=.6) +

  scale_colour_gradient(low="blue", high="red") +

  ylab("Outcome; as observed") +

  labs(col="Observed\nscore", labs="family=acat (this doesn't look right)") +

  xlab("Outcome; as predicted")

p + facet_wrap(~measure_subscale) 

# everyone gets predicted level "3"

# in my 'real' dataset there is some variation in predictions (not everyone 

# is predicted to fall into the same category) but the predictions are about 

# as unreasonable

table(dt.lng$acat_predictions)

On Saturday, December 16, 2017 at 6:51:22 PM UTC-5, Paul Buerkner wrote:

Hmm this sounds strange. Maybe there is a problem with the prediction of 'acat'. Would you mind writing a short reproducible example of your problem?

2017-12-16 22:42 GMT+01:00 Peter Phalen <peter....@gmail.com>:

I seem to be getting very different predictions when I fit a multilevel model using the "acat" family rather than the "cumulative" family.

I'm trying to fit a multilevel model where the outcome is a 3-level ordinal category. I'm estimating whether higher scores on a measure of a theoretical construct predict higher levels of the outcome variable, and I'm particularly interested in variability by subscale.

The code I'm using:

model <- brm(ordered_outcome ~ score + (1 + score | measure_subscale),

                            data = dt.long, iter=4000, cores=4, chains=4, control = list(adapt_delta = .975), family="cumulative")

predictions <- predict(model) # get predictions for the input data

predictions <- apply(model,1,which.max) # get category with the highest probability for each case

The above code gives sensible predictions that roughly match what I get by fitting a series of unpooled proportional odds logistic regressions for each of the different subscales. However, when I run the same code with family set to "acat" (rather than "cumulative") I get very different predictions that don't really match the data. For example, the majority of predictions end up in favor of the category that has the lowest observed count...

Perhaps I'm doing something wrong in processing the predictions?

Thank you!

--
You received this message because you are subscribed to the Google Groups "brms-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+...@googlegroups.com.

[Paul Buerkner]

Many thanks for taking the time to give me a reproducible example!

I can confirm it is a bug (see https://github.com/paul-buerkner/brms/issues/302).

I have already fixed it in the dev version of brms to be installed via

if (!require("devtools")) {
  install.packages("devtools")
}
devtools::install_github("paul-buerkner/brms", dependencies = TRUE)

For the predicitons to work correctly, you have to refit your acat model.

To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.

To post to this group, send email to brms-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/d1a1169a-0fe2-470c-92e6-44b31836040f%40googlegroups.com.

[Peter Phalen]

So fast! Thank you very much. 

Incidentally, what's the substantive difference between these three different ordered logistic regression models? Is one generally preferable to another?

[Paul Buerkner]

With regard to predictions and model fit, they are usually very similar. The main difference is their theoretical derivation.

The cumulative model assumes a continuous latent process that produces the observered categories through categorization.

The sequential model(s) assummes a sequential process the in the end provided the obserserved categories.

The adjacent category model has (in my opinion) no such derivation and is usually used for mathematical convenience,

Here is what I write about them in some of my workshop slides:

## Which Ordinal Model Should I Use?

Is the assumption of a latent continuous variable plausible?

- If yes: use the cumulative model

Is the assumption of a sequential process plausible?

- If yes: use the sequential model

Does sampling speed matter?

- If yes: use the cumulative model

Do I want to model category specific effects?

- If yes: do not use the cumulative model