Lognormal family in brms

Interfaces brms
1

I simulated lognormal data, where the variation on the log scale is constant which I believe represents the lognormal family assumptions in brms.

When I simulated this data I recovered all parameters almost exactly, but when plotting the y obs by average yobs - yrep I see some trend. Is this normal?

I then used my model on my actual data and got the same trend in this plot equivalent to the “error_scatter_avg” in the pp_check. Should it be like this or is something wrong?

image
image697×510 44 KB

Or is this plot more informative:

image
image693×507 42.7 KB

ppc_error_scatter_avg(
  y = log(data2$y), 
  yrep = post_pred,
  x = colMeans(log_epred)  # Use the expected values (no noise) here
)
)

I know that posterior_predredict does the values on the raw scale so I used posterior_linpred which is the log_y and added the noise manually

# residual plot: log scale posterior predictive draws
posterior_draws <- as_draws_df(m1) |> as.data.frame() # posterior draws
log_epred <- posterior_linpred(m1, re.form = NULL)
post_pred <- log_epred
post_pred[] <- NA

# add noise
for( i in 1:nrow(post_pred)) {
  
  sd_group1 <- exp(posterior_draws$b_sigma_group1[i])
  names(sd_group1) <- "group1"
  sd_group2 <- exp(posterior_draws$b_sigma_group2[i])
  names(sd_group2) <- "group2"
  indx <- data2$group
  
  sd <- c(sd_group1 , sd_group2)[indx ]
  
  mu_ <- log_epred[i,]
  
  # alpha & beta
  post_pred[i,] <- rnorm(n = length(mu_), mean = mu_, sd = sd)
  
}

ppc_error_scatter_avg(
     y = log(data2$y), 
     yrep = post_pred
   )

Shouldn’t a lognormal model in brms produce homoscedastic residuals on the log scale, given its assumption of Gaussian errors, or is some heteroscedastic trend expected?

@jsocolar

2

For questions like this it would be informative to see a full reproducible example, including the data simulation, model fitting, processing, and plotting.

It looks to me like you’ve got a lot of data points but a limited number of unique predictor combinations in your simulation, which is why the points in your second plot line up on a limited number of lines.

The trends in your first figure are totally expected–given some set of covariate values, when the observed value of y is higher the residual is bigger. That’s all this figure is showing; it’s a tautology.

The pattern in the residuals is as much a feature of the dataset as a feature of the model. The question you should be asking is: shouldn’t my data simulation produce homoskedastic residuals? If you’re convinced that the residuals you’re seeing are not homoskedastic on the log scale, I would begin by interrogating your data simulation.

3

@jsocolar

Here is my model

fit.brms.log <- brm(bf(formula = y ~ X1 * X2, sigma ~ 1), 
                    family = brms::lognormal(),
                    #prior = prior(normal(0,0.05), class = "sigma", lb = 0),
                    data = data,
                    chains = 4,
                    iter = 5000,
                    warmup = 2000,
                    cores = 4,
                    threads = 4,
                    control = list(adapt_delta = 0.9, max_treedepth = 10),
                    backend = "cmdstanr"
)

Here is my simulation:
I set sd as constant on the log scale

{
  # seed
  set.seed(1995)
  
  # samples
  n <- 450
  
  # generate predictor variables (X1, X2)
  X1 <- rnorm(n)  # X1 ~ Normal(0, 1)
  X2 <- rnorm(n)  # X2 ~ Normal(0, 1)
  
  # set coefficients/model parameters: change in the log of the expected count
  beta_0 <- log(5)  # Intercept 
  beta_1 <- log(1.3)  # predictor 1
  beta_2 <- log(2.4)  # predictor 2
  beta_3 <- log(1.2)
  
  # compute the linear predictor (eta)
  eta <- beta_0 + (beta_1 * X1) + (beta_2 * X2) + (beta_3 * X1 * X2)
  
  mu <- exp(eta)
  sd <- .75
  
  sdlog2 <- log(1 + (sd^2))
  log_mu <- log(mu) - 0.5*sdlog2
  log_sd <- sqrt(sdlog2)
  
  y <- rlnorm(n, meanlog = log_mu, sdlog = log_sd)
  
  # dataset
  data <- data.frame(
    X1 = X1,
    X2 = X2,
    y = y
  )
}

and you can see that this plot shows the residuals are homoscedastic

image
image693×511 75.4 KB

Then here is the version of the simulation where I made sd a function of the mean on the log scale.

{
  # seed
  set.seed(1995)
  
  # samples
  n <- 450
  
  # generate predictor variables (X1, X2)
  X1 <- rnorm(n)  # X1 ~ Normal(0, 1)
  X2 <- rnorm(n)  # X2 ~ Normal(0, 1)
  
  # set coefficients/model parameters: change in the log of the expected count
  beta_0 <- log(5)  # Intercept 
  beta_1 <- log(1.3)  # predictor 1
  beta_2 <- log(2.4)  # predictor 2
  beta_3 <- log(1.2)
  
  # compute the linear predictor (eta)
  eta <- beta_0 + (beta_1 * X1) + (beta_2 * X2) + (beta_3 * X1 * X2)
  
  mu <- exp(eta)
  sd <- .75
  
  sdlog2 <- log(1 + (sd^2/mu^2))
  log_mu <- log(mu) - 0.5*sdlog2
  log_sd <- sqrt(sdlog2)
  
  y <- rlnorm(n, meanlog = log_mu, sdlog = log_sd)
  
  # dataset
  data <- data.frame(
    X1 = X1,
    X2 = X2,
    y = y
  )
}

With this simulation I got this plot

image
image693×514 46.7 KB

So to me it seems like this posterior predictive check is sensitive to detecting homoskedasticity assumption on the log scale.

Do you think this is valid? Do you recommend other diagnostic plots for lognormal glmm?

Thanks.

4

This all makes sense, but it doesn’t reproduce your plots from above? Do you consider your original question to be resolved?