I would think that you would simulate the proportion that survived and then for that proportion simulate the heights and the rest would be zero. This would then provide a simulation of overall distribution of plant heights.

Say you plant 100 seeds at time zero, and all seeds germinate. Let’s suppose you sent out some undergrad student to measure the plants at time T and they took a ruler with markings to every cm and so all you have is integer data to the nearest cm of plants that are alive. Those alive get a non-zero height in cm and those dead get a 0 cm height. I might imagine your generated quantities in a Stan model to look something like:

```
generated quantities {
vector[N] survived_preds;
vector[N] y_preds;
for (n in 1:N) {
survived_preds[n] = bernoulli_rng(alpha_s + beta1_s*Fire[n] + beta2_s*Drought[n] + beta3_s*Herbivory[n]);
if (survived_preds[n]==1) {
y_preds[n] = poisson_rng( alpha_h + beta1_h*Fire[n] + beta2_h*Drought[n] + beta3_h*Herbivory[n] );
}
else if (survived_preds[n]==0) {
y_preds[n] = 0.0;
}
}
}
```

In the generated quantities code above, to simulate the plant heights from your data generation process you simulate dead or alive and then if alive you simulate the height and if dead the height is just zero. (I did not test this code to make sure that else if will work this way in GQ block, but I think this would be the concept).

I’m not sure how you would do this in brms, but it seems to me that you could simply extract the MCMC samples and then use R’s rng functions to simulate predictions with similar code as shown above for the GQ block in a Stan model. For example, see below, where I simulate parameters that would be your MCMC samples extracted from the brms fit. Then plotting quantiles with some heavily modified code from Michael Betancourt’s excellent histograms (hopefully without botching the binning, but you should check for yourself).

```
M <- 1000 #number of MCMC samples from brms run
N <- 100 #number of observations in data
#data
Fire <- rbinom(N, 1, 0.1)
Drought <- rbinom(N, 1, 0.3)
Herbivory <- rbinom(N, 1, 0.4)
#MCMC samples from posterior probability distributions for model parameters
alpha_s <- rnorm(M, 1.75, 0.1)
beta1_s <- rnorm(M, -1, 0.1)
beta2_s <- rnorm(M, -0.5, 0.1)
beta3_s <- rnorm(M, -0.25, 0.1)
alpha_h <- rnorm(M, 30, 0.1)
beta1_h <- rnorm(M, -10, 0.1)
beta2_h <- rnorm(M, -5, 0.1)
beta3_h <- rnorm(M, -3, 0.1)
#simulate survivor status
p <- matrix(NA, nrow=M, ncol=N)
survived_preds <- matrix(NA, nrow=M, ncol=N)
for (n in 1:N){
for (m in 1:M) {
p[m,n] <- plogis(alpha_s[m] + beta1_s[m]*Fire[n] + beta2_s[m]*Drought[n] + beta3_s[m]*Herbivory[n])
survived_preds[m,n] <- rbinom(1, 1, p[m,n])
}
}
#simulate predictions of plant heights
lambda <- matrix(NA, nrow=M, ncol=N)
y_preds <- matrix(NA, nrow=M, ncol=N)
for (n in 1:N){
for (m in 1:M) {
if (survived_preds[m,n]==1) {
lambda[m,n] <- alpha_h[m] + beta1_h[m]*Fire[n] + beta2_h[m]*Drought[n] + beta3_h[m]*Herbivory[n]
y_preds[m,n] <- rpois(1, lambda=lambda[m,n])
} else {
y_preds[m,n] <- 0
}
}
}
#color settings
library(scales)
c_light_m1 <- c("#edf8e9")
c_light_highlight_m1 <- c("#c7e9c0")
c_mid_m1 <- c("#a1d99b")
c_mid_highlight_m1 <- c("#74c476")
c_dark_m1 <- c("#31a354")
c_dark_highlight_m1 <- c("#006d2c")
#plot quantiles of y_preds for each integer
B <- max(y_preds) + 1
idx <- rep(1:B, each=2)
x <- sapply(1:length(idx), function(b) if(b %% 2 == 0) idx[b] + 0.5 else idx[b] - 0.5)
x <- (x - 1)
pred_counts <- sapply(1:1000, function(n)
hist(y_preds[n,], breaks=(0:B) - 0.5, plot=FALSE)$counts)
probs = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
cred<- sapply(1:B, function(b) quantile(pred_counts[b,], probs=probs))
pad_cred <- do.call(cbind, lapply(idx, function(n) cred[1:9,n]))
plot(1, type="n", main="",
xlim=c(-0.5, B + 0.5), xlab="height(cm)",
ylim=c(0, 30), ylab="", cex.lab=2 , cex.axis=1.9)
title("Posterior predictions for plant height", line = 0.5, font.main=1, cex.main=2.2)
polygon(c(x, rev(x)), c(pad_cred[1,], rev(pad_cred[9,])),
col = alpha(c_light_m1, 0.7), border = NA)
polygon(c(x, rev(x)), c(pad_cred[2,], rev(pad_cred[8,])),
col = alpha(c_light_highlight_m1, 0.7), border = NA)
polygon(c(x, rev(x)), c(pad_cred[3,], rev(pad_cred[7,])),
col = alpha(c_mid_m1, 0.7), border = NA)
polygon(c(x, rev(x)), c(pad_cred[4,], rev(pad_cred[6,])),
col = alpha(c_mid_highlight_m1, 0.7), border = NA)
lines(x, pad_cred[5,], col=alpha(c_dark_m1, 0.7), lwd=2)
```

Then you would just overlay the histogram of your observations onto this plot to get your posterior retrodictive check.