Is there something similar to “weights = varIdent(form = ~1 | variable)” in brms?
Could you elaborate on what you want to do, for those less familiar with the quoted code (I assume it’s from nlme).
You can model the scale parameter of several likelihoods.
E.g.,
brmsformula(y ~ x , sigma ~ x )
Would include a (log) linear model of x on the residual standard deviation of y.
Or, with random effects:
brmsformula(y ~ x + (1|C|school), sigma ~ 1 + (1|C|school))
Would include a fixed effect of x on the location of y, with a random intercept across schools; it would also include a random residual SD across schools; the location and scale intercepts would covary as well.
Is this what you are referring to?
4 Likes
So here is I’m trying to do in more details: So basically, I’m currently modelling a 1-1-1 multilevel mediation from bauer et al., (2006). That is, it is a stacked mediation model, meaning that the mediation and outcome equation are stacked into each other.
Here the equation:
In the first parenthesis lies the mediator equation and in the second one, the outcome equation. The S_Mij and SY_ij variables are specification variables. That is when Z is equal to the mediator then S_Mij is equal to 1, and the equation becomes
. The same thing happens when S_Yij = 1. The problem with this equation is with the residual. That is, the models needs a different residual when S_Mij = 1 or S_Yij = 1. Currently, I’m trying to model the heterogeneity of the residual variance in brms, but I’m not sure how to do it.
So there’s sort of like a dummy code per model?
Just model the scale (sigma) using similar dummy codes. You could even have it set up nearly identically:
brmsformula(yourTypicalFormulaHere, sigma ~ 0 + S_M + S_Y) – This means “Fit one intercept when S_M is 1, and one intercept when S_Y is 1, for the (log) residual SD”.
Ok, I’ll try that. Also. would you know a way to get the residuals in the original scale?
Thanks Stephen and Solomon!
So @Stephen_Martin, I used " hyp <- c(“exp(sigma_Sm) = 0”,
“exp(sigma_Sy) = 0”)
from equation: Stacked_Model_B <- brm(bf(Z ~ 0 + Sm + Smx + Sy + Sym + Syx + (0 + Sm + Smx + Sy + Sym + Syx | id) , sigma ~ 0 + Sm + Sy),data = data_stacked2, family = gaussian, iter = 100, control = list(adapt_delta=0.95))
To transform the residuals back to their original scale, and I think it outputs the correct residuals.
One more favour:
hyp1 <- c(“exp(sigma_Intercept) = 0”,
“exp(sigma_Intercept + sigma_variabley) = 0”)
from equation:
Stacked_Model_B <- brm(bf(Z ~ 0 + Sm + Smx + Sy + Sym + Syx + (0 + Sm + Smx + Sy + Sym + Syx | id) , sigma ~ variable),data = data_stacked2, family = gaussian, iter = 100, control = list(adapt_delta=0.95))
outputs the same thing as hyp, and was wondering if someone could give an intuitive expalanation please?
In one case, you use ~0 + Sm + Sy, which says “Don’t fit an intercept, but add these dummy codes as predictors”; in effect, this fits two intercepts - One for M and one for Y.
In the second case, you are fitting an intercept for when variable = 0; if variable = 1 when data are Y, then sigma_Intercept is the log SD when variable = 0 (When data are M, I assume). It’s two ways of estimating the same thing; it’s just a different parameterization. Param 1: Two intercepts, one for each model. Param 2: One intercept, and the ‘effect’ of variable is the difference in intercepts. Just like in regular, basic regression.
1 Like