Yes. I will first explain the reasoning to better understand my code. I start by storing the hurdlefit object. Afterwards, I store the simulation attributes in the object summary_data. With this, I essentially take the mean of every parameter prediction over all 4 chains of simulations. I do this for all 4 parameters, I call them beta0, beta1, gamma0, gamma1. Since the gamma parameters come from a lognormal regression, I take the exp() of them.

```
hurdlefit <- stan('backup.stan', data = list(N=N, Y_obs=Y_obs, K=K, X=X, K_hu=K_hu, X_hu=X_hu, prior_only=prior_only))
beta0 <- median((c(summary_data[[1]][[1]]$b_hu_Intercept, summary_data[[1]][[2]]$b_hu_Intercept, summary_data[[1]][[3]]$b_hu_Intercept, summary_data[[1]][[4]]$b_hu_Intercept)))
beta1 <- median((c(summary_data[[1]][[1]]$`b_hu[1]`, summary_data[[1]][[2]]$`b_hu[1]`, summary_data[[1]][[3]]$`b_hu[1]`, summary_data[[1]][[4]]$`b_hu[1]`)))
gamma0 <- median(exp(c(summary_data[[1]][[1]]$b_Intercept, summary_data[[1]][[2]]$b_Intercept, summary_data[[1]][[3]]$b_Intercept, summary_data[[1]][[4]]$b_Intercept)))
gamma1 <- median(exp(c(summary_data[[1]][[1]]$`b[1]`, summary_data[[1]][[2]]$`b[1]`, summary_data[[1]][[3]]$`b[1]`, summary_data[[1]][[4]]$`b[1]`)))
```

Afterwards, it is a little strange, but I prefer plotting in python so I simply save these values to a file and load them up in python. After which I plot them using the theoretical equation for a lognormal hurdle model, which if Iâ€™m not mistaken is:

\frac{1}{1+e^-(\beta_0+\beta_1 x)}*(\gamma_0 + \gamma_1 x)

Where the first part of the equation is the logistic portion and the second is the linear portion.

Here is the code for plotting:

```
# logistic function
def logit(x, beta0, beta1):
return 1/(1+np.exp(-(beta0 + beta1*x)))
# log_linear function
def linear(x, gamma0, gamma1):
return gamma0 + gamma1*x
# hurdle model functional form
def hurdle(x, beta0, beta1, gamma0, gamma1):
return logit(x, beta0, beta1)*linear(x, gamma0, gamma1)
sigma_plot_values = np.linspace(1, 2.8, num=100)
fitted_line = hurdle(sigma_plot_values, beta0, beta1, gamma0, gamma1)
fig = plt.figure(figsize=(10,8), layout='constrained')
plt.xlim(1.0, 2.8)
plt.ylim(-0.5, 11)
plt.plot(sigma_plot_values, fitted_line, color='black')
plt.show()
```

And for some reason, this is never the same curve as the one obtained from doing conditional_effects() with the brms object, which should be the same.