Categorical interaction modelling deviation from expected average

I have a model design in mind which I don’t know how to implement using brms. Suppose I have a dataset which represents height averages across gender and country. This would yield a dataset with a continuous response (height) and two categorical predictors: gender (male and female) and country (A, B and C). Now I’d like to model the interactions of country and gender as deviations from the expected average. E.g. men are taller than women globally and each country has its own average height. Suppose country A exhibits a fairly moderate average height but there, the women are taller than men hence they deviate from both their expected height associated with their gender and given that the expected height in country A is moderate, this makes women’s height stand out. However, in country B, women are also taller than the global average but men are even taller hence despite deviating from the group average associated with gender, this is not remarkable as women’s height doesn’t deviate much from the expected average associated with country. In essence, I’d like the model to reflect these deviations.

However, I don’t know how to do it and I’ve tried various categorical interaction models with or without varying intercepts but so far, none seems to do what I want it to model:

brm(height ~ 0 + country + gender, data = mydat)
brm(height ~ 0 + country:gender, data = mydat)
brm(height ~ 0 + country + gender + country:gender, data = mydat)
brm(height ~ 0 + country:gender + (1|country) + (1|gender), data = mydat)

Also I don’t know whether including intercepts (e.g. 1 + variable) plays a role here but so far it doesn’t change anything except the coefficients’ interpretation. The closest I’ve come to a solution is the idea to fix the intercept for country and gender to the respective averages and then running an interaction model but I’m not sure how this could work in brms.

Do you have suggestions regarding how to implement such a model?

What I’m trying to achieve is basically a large and more universal Bayesian version of the same notion that is, albeit loosely, connected to the idea behind a frequentist chi-square contingency table test.

Perhaps what is you is to use use effect/sum-coding. See ?contr.sum for documentation. There are lots of tutorials out there that explain coding of variables in regression models (not just brms). Perhaps, you can then just use

height ~ country * gender

Thanks, I’ll look into effect/sum-coding! I’ve since found a way of approximating what I’m looking for by centring the data group-wise before running a regression. It goes like this:
library (plyr) = ddply(mydata, c("gender"), transform, soundratio = scale(height, scale = F)) = ddply(, c("country"), transform, soundratio = scale(height, scale = F))
And then run the following model:
brm(height ~ 0 + country:gender, data =
However, doing this is just hierarchically group-centering the height variable. What would be more elegant is for a model to do this transformation based on estimated (posterior) group means.