K_sigma; // number of population-level effects
matrix[N, K_sigma] X_sigma; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
int Kc_sigma = K_sigma - 1;
matrix[N, Kc_sigma] Xc_sigma; // centered version of X_sigma without an intercept
vector[Kc_sigma] means_X_sigma; // column means of X_sigma before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
for (i in 2:K_sigma) {
means_X_sigma[i - 1] = mean(X_sigma[, i]);
Xc_sigma[, i - 1] = X_sigma[, i] - means_X_sigma[i - 1];
}
}
parameters {
vector[Kc] b; // population-level effects
real Intercept; // temporary intercept for centered predictors
vector[Kc_sigma] b_sigma; // population-level effects
real Intercept_sigma; // temporary intercept for centered predictors
}
transformed parameters {
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = Intercept + Xc * b;
// initialize linear predictor term
vector[N] sigma = Intercept_sigma + Xc_sigma * b_sigma;
for (n in 1:N) {
// inverse function
sigma[n] = sqrt(sigma[n])^2;
}
target += normal_lpdf(Y | mu, sigma);
}
// priors including constants
target += student_t_lpdf(Intercept | 3, 0, 2.5);
target += student_t_lpdf(Intercept_sigma | 3, 0, 2.5);
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
// actual population-level intercept
real b_sigma_Intercept = Intercept_sigma - dot_product(means_X_sigma, b_sigma);
}
```
---
### Author:  [@kholsinger](https://discourse.mc-stan.org/u/kholsinger)
#### Post date: [April 16, 2021, 5:39pm UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/3 "2021-04-16T17:39:02Z")
I’ve only skimmed your Stan code, but it appears to me that you want
```
sigma[n] = sqrt(sigma[n]^2);
```
instead of
```
sigma[n] = sqrt(sigma[n])^2
```
Even if I’ve got that right, it may not fix your problem, but it might be a start.
Kent
---
### Author:  [@sdaza](https://discourse.mc-stan.org/u/sdaza)
#### Post date: [April 17, 2021, 10:23am UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/4 "2021-04-17T10:23:43Z")
Thanks for your reply, @kholsinger !
The model below using brms gives me a good approximation of the coefficients that define sigma (vest):
```
f1 = bf(y ~ g + E + g * E) +
nlf(sigma ~ (sqrt(vest))^2) +
lf(vest ~ 1 + E)
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.04 0.02 -0.00 0.08 1.00 7694 6459
vest_Intercept 0.60 0.01 0.58 0.63 1.00 7736 5898
vest_E 0.15 0.01 0.13 0.17 1.00 6958 6162
g 0.49 0.02 0.45 0.53 1.00 7650 5821
E 0.53 0.02 0.50 0.56 1.00 7592 5897
g:E 0.13 0.02 0.09 0.16 1.00 7499 6472
```
But I get these messages:
```
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: normal_lpdf: Scale parameter[1] is nan, but must be > 0! (in 'model65de1eb63ec1_a2b8f3797cffb15953b3bf0a1f419319' at line 42)
```
If I use:
```
f2 = bf(y ~ g + E + g * E) +
nlf(sigma ~ sqrt(vest^2)) +
lf(vest ~ 1 + E)
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.74 0.34 -1.12 -0.39 3.04 2 11
vest_Intercept 1.24 2.00 -0.34 6.25 3.04 2 11
vest_E -3.25 3.82 -13.38 -0.49 3.04 2 11
g 0.11 1.55 -1.44 1.66 3.04 2 11
E -0.15 0.52 -0.67 0.37 3.03 2 11
g:E 1.69 0.26 1.43 1.96 3.04 2 11
```
The coefficients are way off. As far as I understand, the sigma function should include the inverse function.
Thank you again!
Sebastian
---
### Author:  [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [April 20, 2021, 3:57pm UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/5 "2021-04-20T15:57:30Z")
Hi,
I think the main problem is probably this:
> [@sdaza](#):
>
> `sigma ~ (sqrt(vest))^2)`
This is equivalent to taking the absolute value of `vest`, hower, absolute value is a problem for Stan as it makes the posterior not smooth (the derivative is not continuous).
I understand you use it to be able to work with the “identity” link, but what is the reason you don’t want to use the default “log” link (which is used precisely because it avoids those issue)? Maybe understanding the motivation will help us find a better solution for you.
Alternatively you can use the “softplus” link ([Rectifier (neural networks) - Wikipedia](https://en.wikipedia.org/wiki/Rectifier_(neural_networks)#Softplus)), which behaves almost like identity for positive values of the linear predictor.
Best of luck with your model!
---
### Author:  [@sdaza](https://discourse.mc-stan.org/u/sdaza)
#### Post date: [April 21, 2021, 6:52am UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/6 "2021-04-21T06:52:30Z")
Thanks, @martinmodrak. I am trying to verify and reproduce a test of the proportional contribution of `h` and `e` in the variability of an unobserved `y`. The purpose of the model is to test the proportional contribution of `h` and `e` using this hypothesis:
`hyp <- "g * vest_E = g:E * vest_Intercept"`
Using the data generated below, the difference between the `g * vest_E = g:E * vest_Intercept` should be zero.
```
b0 = 0.8
b1 = 0.2
a = 0.5
h2 = 0.4
h = sqrt(h2)
e = sqrt(1-h^2)
E = rnorm(1000, 0, 1)
N = length(E)
g = rnorm(N,0,1)
sigma = sqrt((b0*e + b1*e*E)^2)
y = a*E + b0*h*g + b1*h*E*g + rnorm(N,0,sigma)
dat = data.frame(y, g, E)
"b0*e = 0.62 and b1*e = 0.155"
"G = b0*h = 0.506 and GxE = b1*h = 0.126"
"diff1 = 0.816 and diff2 = 0.816"
```
If I use the default log link, the scale of the coefficients for sigma is multiplicative (not additive), and I cannot recover the coefficients `b0*e` and `b1*e` that define sigma and then perform the test. I can use the `exp` in sigma to generate the data, and the log link to recover the coefficients, but that changes a bit the data generation I am trying to test (btw, it’s not my model, I am just testing it and the scale of the model).
---
### Author:  [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [April 21, 2021, 11:57am UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/7 "2021-04-21T11:57:26Z")
So I think the model as you present it is _very_ weird. For example for fixed `b1` increase in `E` leads to increase in `sigma` for some observations and decrease in `sigma` for others, just to demonstrate this with the values you have:
```r
b0 = 0.8
b1 = 0.2
h2 = 0.4
h = sqrt(h2)
e = sqrt(1-h^2)
E = c(1, - 4.5)
sqrt((b0*e + b1*e*E)^2)
#> [1] 0.77459667 0.07745967
sqrt((b0*e + b1*e*(E + 0.5))^2)
#> [1] 0.8520563 0.0000000
```
So with the identity link + absolute value, the coefficients for `sigma` are very hard to interpret unless you can make sure that the linear predictor is never smaller than 0 (e.g. by constraining all coefficients AND predictors to be \> 0). The `softplus` link also provides a similar effect where it makes the result always positive, but “softly”.
You also assume the coefficients for `sigma` and `y` are shared, which is currently not supported by `brms` and also makes IMHO little sense - why would something have _exactly_ the same effect on mean and variance? (if there is a reason for this, then maybe better understanding the context would help to see how one can handle this better)
So even if we find a way to make the model work as you’ve written it, I think any result from such a model would be very likely meaningless. If you are 100% sure you want to use the model _exactly_ like that, we can think of ways to make it work, but I would really first stop and ask why are we using such a model and think about alternatives that would make more sense and likely be also easier to fit.
---
### Author:  [@sdaza](https://discourse.mc-stan.org/u/sdaza)
#### Post date: [April 21, 2021, 12:57pm UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/8 "2021-04-21T12:57:34Z")
Thanks for taking the time to review this, @martinmodrak!
Let me go into details. [This is a scaling model for gene x environment (GxE)](https://www.biorxiv.org/content/10.1101/2020.09.08.287888v2) interactions and a test to check if the GxE interaction you get from your model comes from a scaling model where the proportion variation of the phenotype remains constant across E.
Here the generative scaling model (assuming all the variables have variance 1 and mean 0). We have an unobserved phenotype (G = genetic info):
Y^\*\_i = h G\_i + e \epsilon\_i
Where e = \sqrt{1 - h^2}, and \epsilon\_i is an unobserved error term.
The observed Y takes the form (E = environment):
Y\_i = a\_0 + a\_1 E\_i + (b\_0 + b\_1 E\_i) Y\_i^\*
Y\_i = a\_0 + a\_1 E\_i + b\_0hG\_i + b\_1hE\_iG\_i + b\_0e\epsilon\_i + b\_1e E\_i\epsilon\_i
The model to run with real data would be (this is what I am trying to estimate using brms):
Y\_i = \tau\_0 + \tau\_1 E\_i + \pi\_0G\_i + \pi\_1E\_iG\_i + \lambda\_0\epsilon\_i + \lambda\_1 E\_i\epsilon\_i
Using the coefficients above, I can run a test to check if my _empirical_ model comes from the _scaling_ where the contribution to the phenotype remains constant across E:
\frac{\pi\_0}{\lambda\_0} = \frac{\pi\_1}{\lambda\_1} = \frac{b\_0h}{b\_0e} = \frac{b\_1h}{b\_1e}
H\_0 : \pi\_0\lambda\_1 - \pi\_1\lambda\_0 = 0
---
### Author:  [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [April 27, 2021, 2:13pm UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/9 "2021-04-27T14:13:55Z")
Sorry, I dropped the ball a bit on this. So rewriting the model a bit just for my own clarity:
Y\_i \sim N(\mu\_i, \sigma\_i) \\ \mu\_i = \tau\_0 + \tau\_1 E\_i + \pi\_0G\_i + \pi\_1E\_iG\_i \\ \sigma\_i = \lambda\_0+ \lambda\_1 E\_i
Some things to note:
- if E\_i is a discrete predictor (as appears to be the case in the paper you linked), you just end up having a different fixed variance for each level of E. So changing link function to `log` still allows each level to have its own variance. The only thing that changed is the structure of the prior, but that can IMHO be easily adjusted to not be very different than when assuming identity link. If E is to be a smooth predictor than you can’t match it exactly - and as I said above you can get some very weird behaviour from such a model.
- If you allow \lambda\_0, \lambda\_1 \< 0 the model is implicitly multimodal - flipping signs on both lambdas gives you IMHO _exactly_ the same likelihood. Weirdly enough the paper reports fitted coefficients for some lambdas less than zero, so they are probably not worried about this - possible precisely because they only have discrete (categorical) E.
The hypothesis test as you envision it will have its own challenges, but let’s keep that for later.
So the big question is: do you need continuous E? If not, just use the `log` link. If yes, I think the `softplus` link is going to give you behaviour quite close to the paper in the model while hopefully not breaki ng `brms`.
Hope that helps at least a bit.
---
### Author:  [@sdaza](https://discourse.mc-stan.org/u/sdaza)
#### Post date: [April 28, 2021, 7:32am UTC](https://discourse.mc-stan.org/t/retrieving-sigma-coefficients-when-sigma-comes-from-a-function/21759/10 "2021-04-28T07:32:42Z")
Thanks again, @martinmodrak. In principle, E\_i can be continuous or discrete. In their simulation, though, E\_i is continuous. If I use the `softplus` link (based on my analysis), I can’t reproduce the test they propose because of the scale of the coefficients is different (multiplicative vs additive).
The key question for me would be: let’s assume I get some data and estimate the model using a `softplus` link. Then I try to use the test they propose to discard the scaling model is generating the GxE coefficient, but because the scale of the coefficients is different, the test seems useless. This would be independent whether the scaling model is or not a good model to test (that’s is my actual research question).
For now, I just changed the generation data definition so that the log link works and I can reproduce the scaling test:
```
sigma = exp(b0*e + b1*e*E)
y = rnorm(N, a0 + a1*E + b0*h*G + b1*h*E*G, sigma)
```
Thank you so much!