Operating System: Windows 7, brms Version: 2.9.0
We try to analyse an effect developing in traumatized rats over time (saturating around 120 min). The problem is that effect could grow to be negative (in some rats) and positive (in other rats). While we analyse this direction as a separate (bernoulli) model, we want to assert the effect’s presence by modelling the absolute value of the effect (mpa variable). And this is where our problems started:
-
Gaussian model mpa ~ Trt*Time converges nicely, but some of 95% HPD intervals span negative numbers. Sure one could cover this ugliness with “bar + only upper part of SEM” plots, but I really hope someone here could advise a better way.
-
Truncated gaussian model mpa | trunc(lb=0) ~ Trt*Time simply NEVER converges. Tightening max_treedepth=13, adapt_delta=0.99 and iter=1e4 do not help (Rhat >= 8830, Eff.Sample = 2). Below is the reproducible example and EVERY truncated gaussian model I faked by playing with number of rats, effect sizes, etc. suffer from these convergence problems.
-
I tried to use “inherently positive” distribution families, like LogNormal. Unfortunately our data is measured by rounding to the nearest 0.5mm and then averaged over 6 replicates, creating many “strict zeros” and breaking the positive distribution.
This is my fake example:
require("tidyverse")
require("tidybayes")
require("brms")
require("emmeans")
set.seed(1)
nRats <- 20
rats <- data.frame(ratNo=1:nRats, mu=rnorm(nRats)/3, sigma=rlnorm(nRats),
TrtEffect=rnorm(nRats,mean=1)) # each rat is special
trueParams <- expand.grid(replicateNo=1:6, ratNo=1:nRats, Time=c(0,30,60)) %>%
mutate(Trt = ifelse(ratNo <= nRats/2, "Sham", "TBI"),
mu=rats$mu[ratNo]+ifelse(Trt=="Sham", 0, rats$TrtEffect[ratNo]*Time/10),
sigma = rats$sigma[ratNo] + 1 + mu/3)
data.measured <- trueParams %>%
mutate( pa = rnorm(n(), mean=mu, sd=sigma), # the "real" effect
mpa = abs(round(pa*2)/2) ) # and its rounded measurement, producing many zeroes
data.recorded <- data.measured %>% ungroup %>%
mutate(Trt = as.factor(Trt), Time = as.factor(Time), ratNo = as.factor(ratNo) ) %>%
group_by(Trt, Time, ratNo) %>% summarize(mpa = mean(mpa))
m.mpa <- brm( data = data.recorded, family = gaussian,
# mpa ~ Trt*Time, # converges, but fails to produce strictly positive 95% HPDI
mpa | trunc(lb=0) ~ Trt*Time, # convergence fails horribly
control=list(max_treedepth=13, adapt_delta=0.99),
seed = 1, chains=3, cores=3, iter=1e4)
Of course, I have many other questions, like (4) how do I set reasonable priors if I know that untreated rats should have effect <1mm ( -1 <= pa <= 1), or (5) if it is worth to make Time a continuous factor, which will force us into non-linear models to account for saturation.