How can I deal with this heteroskedastic model?

Hi all.

I am interested in the latency it takes to flies to get attacked by predators, and the relationship between such latency and several covariates. I tested each spider more than once, so spider ID is a random factor. In traditional R language, my model looks like:

mod<- lmer(Latency~Sex+Age+Distance+Lift+Walk_back+Size+condition+Sex:Age+Lift:Walk_back+(1|ID_Spider)) 


The problem with this model is that the residual versus fitted plot looks heteroskedastic. However, the plots of residuals versus individual explanatory variables all look fine. This shows me that whatever variable is generating heteroskedasticity, is not in my data set. If I was going with a full frequentist approach, I would use the HC (also known as “Sandwich”) estimators for the standard errors, to correct for this problem (see Cleasby & Nagakawa 2011). However, I am not aware of any similar approach in the Bayesian framework. Of course, I technically could model something like:

X<-Sex+Age+Distance+Lift+Walk_back+Size+condition+Sex:Age+Lift:Walk_back+(1|ID_Spider))

brm(Latency~X, sigma~X)
  

But that looks like a very heavy model. I am almost sure convergence problems will arise.

Can anyone help me with this?

Have you tried fitting the “heavy” model in brms to see if issues arise? Would a log-linear model reduce the heteroskedasticity? Perhaps you could model the scale directly as a “distributional model”?