# Modeling clinical depression data in ADHD

Hi,

I am trying to find a clinically sensible way to describe group difference in clinical psychological symptoms.

For example, we had 200 individuals come to the lab for testing. These are their depression scores taken from a self-report that can range between 0 and 63 (BDI).

The mean of ADHD is ~10pts higher, which has clinical meaning. But - I feel that just talking about the mean misses an important aspect here. We have what seems to me as a “higher rate” of more severe case, and that is important. The extreme cases are the most meaningful even if they are few. I fitted a y~1+group in brms with a gaussian that seems to do a good job only at capturing the mean. Also disregarding that the scores can only be between 0 and 63. Here is a pp_check:

So, my question is how would you get parameters that can have a clear interpretation regarding the rate of high cases. What I tried so far is a beta_binomial

model<-brm(bf(y | trials(100) ~ 0 + Intercept + group,
phi ~ 0 + Intercept + group),
data   = df ,
warmup = 8000,
iter   = 10000,
cores  = 4,
chains = 4,
silent = 0,
#              prior  = myprior,
#              sample_prior = "only",
backend='cmdstan')


which gave me this pp_check and posterior estimates:

Yes - it looks better, but I can’t convey the phi differences in any clinically meaningful way. I converted this to alpha/beta medians and got approximately:
alpha_control = 0.95
beta_control = 8.5

I think I can say “adhd has 1.3/0.95” more right skewness" but I’m not even sure this is mathematically correct. I ended up sampling and calculating the precent of response for each level of symptoms for each group for from the posterior median parameters:

####beta
library(ggplot2)
library(dplyr)

# Define parameters
alpha_control <- 0.95
beta_control <- 8.5
n <- 1000000  # Number of simulations
max_score <- 63

# Simulate scores
set.seed(123)  # For reproducibility
scores_control <- rbeta(n, alpha_control, beta_control) * max_score

library(ggplot2)
library(dplyr)
library(tidyr)

# Assuming scores_adhd and scores_control have been generated as before

# Prepare the data for plotting
data <- tibble(
group = rep(c("ADHD", "Control"), each = n)
)

# Calculate the percentage of each score within each group
data_summary <- data %>%
group_by(group, score) %>%
summarise(n_bin = n(), .groups = 'drop') %>%
mutate(percentage = (n_bin / n) * 100) %>%
ungroup() %>%
# Ensure every possible score is represented for both groups
complete(group, score, fill = list(n = 0, percentage = 0)) %>%
arrange(group, score)

# Plot
ggplot(data_summary, aes(x = score, y = percentage, color = group)) +
geom_col(position = "dodge") +
ylim(0, 20) +
labs(title = "Probability Distribution of Depression Scores",
x = "Depression Score",
y = "Percentage Chance of Each Score") +
scale_x_continuous(breaks = seq(0, max_score, by = 5)) +
scale_y_continuous(
breaks = seq(0, 20, by = 2),
labels = function(x) paste0(x, "%"),
limits = c(0, 20)
) +
scale_color_manual(values = c("ADHD" = "blue", "Control" = "red")) +
theme_minimal() +
theme(
panel.grid.major.x = element_blank(), # Remove vertical major grid lines
panel.grid.minor.x = element_blank(), # Remove vertical minor grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal minor grid lines
panel.grid.major.y = element_line(color = "grey80", size = 0.5) # Customize horizontal major grid lines
)


So this now has finally more clinical strightforward interpetation I think. We can see how much more adhd response are more probable across the self-report range. But now, I lost all the nice uncertainty estimates that bayesian inferences gives.

My questions are:

1. Do you have other suggestions on how to gain parameters that will have estimates that can be clearly interpreted, specifically in regards to the heaviness of the tail.
2. What will be a good way to convey the change in probability of getting a specific score between control to adhd?
I should say that I rather not work with the clinical threshold, these are really arbitrary in the end. I also love the exgaussian, but for some reason brms/stan has a really hard time sampling for that (I guess this is for a different topic).

Thank you for reading, and I will really appreciate any thought!
Nitzan

3 Likes

Hi. Friendly neighborhood clinical psychologist, here. I really like your beta-binomial approach; it’s a much better fit than your initial Gaussian model, and I’d love to see more beta-binomial models for sum-score data out in the wild (especially in the clinical psychology literature). To my eye, you can capture the differences in the distributions by describing their means and variances. If some variable y is distributed beta-binomial, you can define the mean as

\mathbb E(y) = n\mu,

and define the variance as

\operatorname{Var}(y) = n\mu(1 - \mu) \frac{\phi + n}{\phi + 1}.

If you’re afraid us clinical psychologists won’t like variances, just take the square root to convert them to standard deviations.

2 Likes

Nice to be in a friendly neighborhood - thanks for writing :)

I tend to think that variance (or sd) won’t really say anything specific about the right-edge, and the tendency for more severe scores… But maybe I’m missing something in the math? It seems to me like the var conveys spread from on both edges.

Also, I am not sure I really understand what will be the benefit of reporting var from beta-binomial compared to Gaussian likelihood… As rough as the pp_check might be for the Gaussian - it still get the first and second moment ok I think…

I dunno - but it feels like if we are back to mean and var - we are missing in the report something important about the skewness…

To my mind, a left- and right-bounded likelihood like the beta-binomial presumes you’ll have skew the closer you get to either of the boundaries. Further, if you hold a mean constant near the lower boundary but increase the variance, you’ll necessarily increase the skew.

1 Like

You can use the marginaleffects or the brmsmargins package to calculate the mean difference between the two groups.

To emphasize the longer tail, consider reporting the 80th, 90th, or 95th percentile, or you could state that X% of individuals with ADHD score above the 95th percentile of those without ADHD.

2 Likes

Would something, at least conceptually, like quantile regression help here to compare the distributions?

1 Like

Interesting thought @Richard_Erickson thanks. I actually didn’t know these were available…
I gave it a try and this is what I got. What do you think? I think the nice thing here is that we can see that for the high scores in each group (0.8 quantile) we get a difference of close to 20 pts, while the mean difference was more around 10pts. I think this is usefull since the difference in raw scores is very meangfull. There are two issues here, one that I had to use three different models, and the other is that it seems pp_check is not really usefull (since each model is capturing only a bit of the data). Is there a better way to go about?

model_0.2<-brm(bf(y  ~ 0 + Intercept + group,quantile = 0.2),
data   = df ,
warmup = 8000,
iter   = 10000,
cores  = 4,
chains = 4,
silent = 0,
family = asym_laplace(),
backend='cmdstan')

model_0.5<-brm(bf(y  ~ 0 + Intercept + group,quantile = 0.5),
data   = df ,
warmup = 8000,
iter   = 10000,
cores  = 4,
chains = 4,
silent = 0,
family = asym_laplace(),
backend='cmdstan')

model_0.8<-brm(bf(y  ~ 0 + Intercept + group,quantile = 0.8),
data   = df ,
warmup = 8000,
iter   = 10000,
cores  = 4,
chains = 4,
silent = 0,
family = asym_laplace(),
backend='cmdstan')


thank you

2 Likes

I’ve actually only used simple, non-Bayesian quantile regression a couple of times, possibly only in a stats class. I would look at the literature in your field to see how they use quantile regression. Personally, I would probably use what you have now for your paper.

It’s understandable and estimates the differences. Basically, the “typical” or “average” (near median) person with ADHD has a slightly higher depression score, but individuals with higher depression score have a much higher score compared to non-ADHD individuals.

1 Like

I’ve never had reason to dip my toes into quantile regression, which is a very rare bird in clinical psychology. To my eye, the beta-binomial model is already doing a fine job describing the data.

1 Like

You can get the group differences for quantile levels in the original scale also from the beta-binomial model by using the posterior predictions.

Let’s say PP_ADHD is a matrix of posterior predictions for the ADHD cases, where participants are in columns and iterations are in the rows. Then
Q80_ADHD = apply(pp_ADHD,1, function(x) quantile(x, prob = .8))
will produce a vector (posterior distribution) of 80% quantiles for the ADHD group.
With the same approach you can get the 80% quantile posterior distribution for the control group, which you can use properly to calculate the posterior distribution of the difference of the 80% quantile value between both groups:
Q80_delta = Q80_ADHD - Q80_CONTR.

Obviously, you can do this for any quantile of your choosing.

Having said this, I think it is easier to communicate the difference in the prevalence of extreme sum scores by saying that, I am making up numbers now: “50% (CI: 40, 60) of ADHD cases score higher than the 95% quantile of the control cases”. This can also be calculated from the posterior distribution. (because you are dealing with samples, you can add credible intervals to the made 50% value)

1 Like

Thank you everyone, this is very helpful. I now realize that the big advantage of using beta_bionomial (for me at least) is that we can ask questions regarding any bin of scores using the predictive distribution (and get uncertainty estimates). This is great because in the end I think we want to talk about “differences in symptom severity” which can be expressed in different ways (not just the mean). Here is what I ended up doing did in case its helpfull for others. Would be more then happy to get your thoughts if you have any, or if you see any thing that might be inaccurate.
Thanks!
Nitzan

Quick explanation: (A) Raw scores and posterior predictive estimates of the group mean differences. (B) posterior for the group difference between the means (mu from the beta_binomial), (C) posterior predictive distribution showing the percent of cases above the clinical cutoff for each group, (D) An estimate of the expected percent of cases for each scores, in each group (using the mean posterior mu and phi)

Text for results: “We found that ADHD were about 10 point higher on average then controls on their depression score (CI90% between 8.4 and 13.2). We then sampled from the posterior predictive distrbution and estimated the percentage of cases that scored above the clinical cutoff (21pts). We found individuals with ADHD to have 33% chance of being above the clinical cut-off (CI90% 25.4 to 41.9), compared to 3.8% in the control group (CI90% 1.5 to 7.4). Individuals in the ADHD were about 9 times more likely to score above the clinical cutoff compared with non-ADHD individuals (CI90% between 4.2 and 22.8)"


#sample
model<-brm(bf(y | trials(63) ~ 0 + Intercept + group,
phi ~ 0 + Intercept + group),
data   = df ,
warmup = 30000,
iter   = 31000,
cores  = 4,
chains = 4,
silent = 0,
prior  = myprior,
#sample_prior = "only",
backend='cmdstan')

#### percent above cutoff ----

pp_control = posterior_predict(model, newdata = (data.frame(group = rep("Control",1000))))
percent_above_control = rowSums(pp_control > 21) / ncol(pp_control) * 100

describe_posterior(percent_above_control)

as.data.frame() |>
ggplot(aes(x = percent_above, fill = group)) +
xlim(0,100)+ xlab("percent above cutoff (ADHD)") +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray")+
stat_halfeye()

3 Likes

The colors in the plot are switching: In plot A blue is Control, but in C and D blue is ADHD, which caused a moment of confusion. Otherwise the plots look great and informative!

2 Likes

Agreed. Overall, nice visualization.

1 Like

Nice plot!