# Model not recovering simulated parameters

**URL:** <https://discourse.mc-stan.org/t/model-not-recovering-simulated-parameters/32668>\
**Category:** Modeling\
**Tags:** rethinking, fitting-issues\
**Created:** [September 5, 2023, 11:27pm UTC](https://discourse.mc-stan.org/t/model-not-recovering-simulated-parameters/32668 "2023-09-05T23:27:03Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![e\_b](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/e_b/32/17071_2.png) [@e\_b](https://discourse.mc-stan.org/u/e_b)\
**Post date:** [September 5, 2023, 11:27pm UTC](https://discourse.mc-stan.org/t/model-not-recovering-simulated-parameters/32668/1 "2023-09-05T23:27:03Z")

</div>

I have been futzing around with a model for a while using simulated data, and I cannot recovered the exact parameters I set using the model. However, I am recovering the correct differences. I am not sure if I am messing up the model, or if I am doing everything correctly and my expectations for this model was incorrect.

Essentially, I want to use a model to determine how a behavior changes with the addition of a treatment. The behavioral response is a continuous measurement between 0 and 1 (or 0.99 to avoid certain math errors with a logit link), which I’ve simulated using the beta2 function in rethinking. In the absence of a treatment, the behavioral response should center around .4. With the addition of a treatment, the response should shift up by .1 (low quality treatments) or .3 (high quality treatments). This is how I’ve simulated the data in R:

```plaintext
#create df with 4 treatments, 6 replicates per treatment
treat <- c("A", "B", "C", "D")
ID <- c(1:6)
df <- expand_grid(treat, ID)

#expand each trial for 10 rows, representing 10 measurements per replicate
df <- df %>% 
  slice(rep(1:n(), 10)) 

#simulate behavioral response in absence of treatment
#simulating from beta to only get values between 0 and 1
df$absent <- rbeta2(n=240, .4 , 10)

#Simulate behavioral response with treatment
df <- df %>%
mutate(m_s = case_when(
  treat== "A" ~ .3,
  treat == "B" ~ .3,
  treat == "C" ~ .1,
  treat == "D" ~ .1)
)
df$present <- rbeta2(n=240, (.4 + df$m_s), 10)

#Replace "1" with ".99" in absent and present to avoid math errors
df$absent[df$absent == 1] <- .99
df$present[df$present == 1] <- .99

```

Then I tried to recover the m\_s values using this model:

```stan
parameters{
    real a;
    vector[4] t;
    real<lower=0> scale;
}
model{
    real mu;
    vector[240] mu2;
    scale ~ gamma( 15/1 , 1/1 );
    s ~ normal( 1 , 1 );
    a ~ normal( -0.4 , 0.5 );
    for ( i in 1:240 ) {
        mu2[i] = a + s[species[i]];
        mu2[i] = inv_logit(mu2[i]);
    }
    mu = a;
    mu = inv_logit(mu);
    present ~ beta( mu2*scale , (1-mu2)*scale );
    absent ~ beta( mu*scale , (1-mu)*scale );
}
generated quantities{
    vector[240] log_lik;
    real mu;
    vector[240] mu2;
    for ( i in 1:240 ) {
        mu2[i] = a + t[treat[i]];
        mu2[i] = inv_logit(mu2[i]);
    }
    mu = a;
    mu = inv_logit(mu);
    for ( i in 1:240 ) log_lik[i] = beta_lpdf( absent[i] | mu*scale , (1-mu)*scale );
}

```

However, every time I run this the output values for t, which should correspond to my m\_s values when I simulated the data, are off by ~.4 (after accounting for the logit link). The difference between treatments is the same, ~.2. I think my models estimates of t do not separate out the intercept, a, and I have no idea how to get the model to do so. Or even if it’s supposed to do so–I always assumed that models would estimate t given intercept, a, but maybe I’m wrong? I’m concerned if I can’t figure this out, I won’t understand my results when I apply the model to real data where I don’t know the true parameters. Does anyone know what’s going on here? Is the model the problem, or my expectations of the model?

---

<div class="post-metadata">

**Author:** ![desislava](https://avatars.discourse-cdn.com/v4/letter/d/5f8ce5/32.png) [@desislava](https://discourse.mc-stan.org/u/desislava)\
**Post date:** [September 9, 2023, 11:14am UTC](https://discourse.mc-stan.org/t/model-not-recovering-simulated-parameters/32668/2 "2023-09-09T11:14:22Z")

</div>

It may help to simplify the Stan code. You introduce the mean parameters `mu` and `mu2` in the model block. And then you recalculate those parameters in the generated quantities block. It’s easier to define `mu` and `mu2` in the transformed parameters block instead.

And I don’t think it’s true that

> [@e\_b](#):
>
> the output values for t, which should correspond to my m\_s values when I simulated the data

The `m_s` parameters are defined on the original scale in the interval (0, 1):

```r
df$present <- rbeta2(n=240, (.4 + df$m_s), 10) # The second param in rbeta2 is between 0 and 1.

```

While the `t` parameters are defined on the link (logit) scale:

```stan
for ( i in 1:240 ) {
    mu2[i] = a + t[treat[i]];
    mu2[i] = inv_logit(mu2[i]);
}

```

It’s easier to recover the `m_s` parameters on the original scale in the generated quantities block: `m_s = mu2 - mu`. I’m not sure you can simplify `m_s = inv_logit(a + t) - inv_logit(a)` into a function of `t` only.

```stan
data{
  int N;
  vector<lower=0, upper=1>[N] present;
  vector<lower=0, upper=1>[N] absent;
  array[N] int treat;
}
parameters{
    real a;
    vector[4] t;
    real<lower=0> scale;
}
transformed parameters{
    real<lower=0, upper=1> mu;
    vector<lower=0, upper=1>[4] mu2;
    mu = inv_logit( a );
    mu2 = inv_logit( a + t );
}
model{
    // priors
    scale ~ gamma( 15/1 , 1/1 );
    t ~ normal( 1 , 1 );
    a ~ normal( -0.4 , 0.5 );
    // likelihood
    present ~ beta( mu2[treat]*scale , (1-mu2[treat])*scale );
    absent ~ beta( mu*scale , (1-mu)*scale );
}
generated quantities{
  vector[4] m_s;
  m_s = mu2 - mu;
}

```