Get posterior draws of implicit latent variables in nonlinear model

Modeling
,
1

I fitted a model using the following code in R (data are attached).

Is it possible to get posterior distributions for :

Y1 = inv_logit(Araw) *
exp(-exp(lnBraw) * exp(-c0 * inv_logit(Mraw) * inv_logit(Sraw) * adj_age)

for the original observations, and for:

Y2 = inv_logit(Graw) * inv_logit(Craw)

using new data with a constant gdd?

newdata <- ppari_data |> mutate(gdd = 2500)

I want all fixed and random effects in the posterior distributions.

library(tidyverse)
library(brms)
library(rstan)
library(cmdstanr)


ppari_data <- read_csv("ppari_data.csv") |>
  mutate(
    srv = factor(srv)
  )

bf_lgt_ppari_e <- bf(
  lgt_ppari ~ logit(
    inv_logit(Araw) * inv_logit(Graw) * inv_logit(Craw) *
      exp(-exp(lnBraw) * exp(-c0 * inv_logit(Mraw) * inv_logit(Sraw) * adj_age))
  ),
  Araw ~ 1 + (1|srv),
  Graw ~ 1 + s(gdd, k = 5) + (1|srv),
  Craw ~ 1 + s(chill28Feb, k = 5) + (1|srv),
  lnBraw ~ 1 + s(m2_tree, k = 5) + (1|srv),
  Sraw ~ 1 + s(shape, k = 5) + (1|srv),
  Mraw ~ 1 + s(m2_tree, k = 5) + (1|srv),
  c0 ~ 1,
  nl = TRUE
)

priors_lgt_ppari_e <- c(
  prior(normal(logit(0.99), 1), nlpar = "Araw"),
  prior(normal(logit(0.98), 1), nlpar = "Graw"),
  prior(normal(logit(0.98), 1), nlpar = "Craw"),
  prior(normal(log(54),      1), nlpar = "lnBraw"),
  prior(normal(logit(0.5),  1), nlpar = "Sraw"),
  prior(normal(logit(0.5),  1), nlpar = "Mraw"),
  prior(normal(1.35,      0.01), nlpar = "c0")
)

fit_ppari <- brm(
  formula = bf_lgt_ppari_e,
  data    = ppari_data,
  family  = gaussian(),
  prior   = priors_lgt_ppari_e,
  init    = "random",
  control = list(adapt_delta = 0.99,
                 max_treedepth = 15),
  chains  = 4,
  cores   = 4,
  iter    = 4000,
  warmup  = 500,
  backend = "cmdstanr",
  threads = threading(5),
  file    = "fit_ppari.rds"
  )

I am trying to use the data and restrictions imposed by the model and priors to “elicit” values for the nonlinear parameters, which (in proper functions) represent latent variables. Loosely, Y1 represents the actual proportion of canopy size for a given age achieved with the realized gdd, etc., whereas Y2 represents the maximum canopy size achievable based on the rest of the conditions, when gdd and chill are not limiting. In addition to exploring those two variables directly, I want to use then as predictors in another model. I have not been able to formulate and fit the complete DAG with explicit latent variables and all components yet. I am working on a step-by-step approach.
Any comments on anything that comes to mind will be greatly appreciated.

I have tried many things, including getting posterior_linpred for each parameter, but that apparently breaks the correlations among parameters. I am at the point of operating directly on

fit_lgt_ppari_e_parm_draws <- as_draws(fit_ppari) |>
 map(as_tibble) |>
 bind_rows()

which seems extremely cumbersome, and would not give me the values for the Y1 I need.

Thank you!

ppari_data.csv (360.1 KB)

2

I think I found the solution. My earlier attempts with posterior_linpred were thwarted by an unrelated error. Correlations are not broken by obtaining each linpred separately. Using the “transformed parameters” approach helped me double check results, but I do not recommend it in this specific case because it creates unnecessary large objects.

Comments are encouraged.

# Solution ========

## New fit with same data but explicitly getting draws of nonlinear
## parameters plus residual variance.

stan2get_la_scaf <- "
  vector[N] nlp_Araw;
  vector[N] nlp_Graw;
  vector[N] nlp_Craw;
  vector[N] nlp_lnBraw;
  vector[N] nlp_Sraw;
  vector[N] nlp_Mraw;
  vector[N] nlp_c0;

  // (you already have s_*/r_* computed above)
  for (n in 1:N) {
    // fixed–effects + group–effects for Araw
    nlp_Araw[n]  = X_Araw[n] * b_Araw
                 + r_1_Araw_1[J_1[n]] * Z_1_Araw_1[n];

    // Graw: population + spline + group
    nlp_Graw[n]  = X_Graw[n] * b_Graw
                 + Xs_Graw[n] * bs_Graw
                 + Zs_Graw_1_1[n] * s_Graw_1_1
                 + r_2_Graw_1[J_2[n]] * Z_2_Graw_1[n];

    // Craw
    nlp_Craw[n]  = X_Craw[n] * b_Craw
                 + Xs_Craw[n] * bs_Craw
                 + Zs_Craw_1_1[n] * s_Craw_1_1
                 + r_3_Craw_1[J_3[n]] * Z_3_Craw_1[n];

    // lnBraw
    nlp_lnBraw[n] = X_lnBraw[n] * b_lnBraw
                  + Xs_lnBraw[n] * bs_lnBraw
                  + Zs_lnBraw_1_1[n] * s_lnBraw_1_1
                  + r_4_lnBraw_1[J_4[n]] * Z_4_lnBraw_1[n];

    // Sraw
    nlp_Sraw[n]   = X_Sraw[n] * b_Sraw
                  + Xs_Sraw[n] * bs_Sraw
                  + Zs_Sraw_1_1[n] * s_Sraw_1_1
                  + r_5_Sraw_1[J_5[n]] * Z_5_Sraw_1[n];

    // Mraw
    nlp_Mraw[n]   = X_Mraw[n] * b_Mraw
                  + Xs_Mraw[n] * bs_Mraw
                  + Zs_Mraw_1_1[n] * s_Mraw_1_1
                  + r_6_Mraw_1[J_6[n]] * Z_6_Mraw_1[n];

    // c0
    nlp_c0[n]     = X_c0[n] * b_c0;
  }
  "

bf_lgt_ppari_f <- bf(
  lgt_ppari ~ logit(
    inv_logit(Araw) * inv_logit(Graw) * inv_logit(Craw) *
      exp(-exp(lnBraw) * exp(-c0 * inv_logit(Mraw) * inv_logit(Sraw) * adj_age))
  ),
  Araw ~ 1 + (1|srv),
  Graw ~ 1 + s(gdd, k = 5) + (1|srv),
  Craw ~ 1 + s(chill28Feb, k = 5) + (1|srv),
  lnBraw ~ 1 + s(m2_tree, k = 5) + (1|srv),
  Sraw ~ 1 + s(shape, k = 5) + (1|srv),
  Mraw ~ 1 + s(m2_tree, k = 5) + (1|srv),
  c0 ~ 1,
  nl = TRUE
)

priors_lgt_ppari_f <- c(
  prior(normal(logit(0.99), 1), nlpar = "Araw"),
  prior(normal(logit(0.98), 1), nlpar = "Graw"),
  prior(normal(logit(0.98), 1), nlpar = "Craw"),
  prior(normal(log(54),      1), nlpar = "lnBraw"),
  prior(normal(logit(0.5),  1), nlpar = "Sraw"),
  prior(normal(logit(0.5),  1), nlpar = "Mraw"),
  prior(normal(1.35,      0.001), nlpar = "c0")
)


fit_lgt_ppari_f <- brm(
  formula = bf_lgt_ppari_f,
  data    = ppari_data,
  stanvars = stanvar(
    scode = stan2get_la_scaf,
    block = "tparameters",
    position = "end"
  ),
  family  = gaussian(),
  prior   = priors_lgt_ppari_f,
  init    = "random",
  control = list(adapt_delta = 0.99,
                 max_treedepth = 15),
  chains  = 4,
  cores   = 4,
  iter    = 4000,
  warmup  = 500,
  backend = "cmdstanr",
  threads = threading(5),
  file    = paste0(
    proj_path,
    "/models/fit_lgt_ppari_f.rds"
  )
)


## Get posterior_predict ======

post_pred_lgt_ppari <- posterior_predict(
  fit_lgt_ppari_f
) |>
  unname()

## Reconstruct draws using posterior_linpred =====
## Note that each is obtained "independently" in different
## calls to posterior_linpred()
A_lin   <- posterior_linpred(fit_lgt_ppari_f,
                             nlpar      = "Araw",
                             newdata    = ppari_data,
                             re_formula = NULL,
                             transform  = FALSE)

G_lin   <- posterior_linpred(fit_lgt_ppari_f,
                             nlpar      = "Graw",
                             newdata    = ppari_data,
                             re_formula = NULL,
                             transform  = FALSE)

C_lin   <- posterior_linpred(fit_lgt_ppari_f,
                             nlpar      = "Craw",
                             newdata    = ppari_data,
                             re_formula = NULL,
                             transform  = FALSE)

B_lin   <- posterior_linpred(fit_lgt_ppari_f,
                             nlpar      = "lnBraw",
                             newdata    = ppari_data,
                             re_formula = NULL,
                             transform  = FALSE)

S_lin   <- posterior_linpred(fit_lgt_ppari_f,
                             nlpar      = "Sraw",
                             newdata    = ppari_data,
                             re_formula = NULL,
                             transform  = FALSE)

M_lin   <- posterior_linpred(fit_lgt_ppari_f,
                             nlpar      = "Mraw",
                             newdata    = ppari_data,
                             re_formula = NULL,
                             transform  = FALSE)

c0_lin  <- posterior_linpred(fit_lgt_ppari_f,
                             nlpar      = "c0",
                             newdata    = ppari_data,
                             re_formula = NULL,
                             transform  = FALSE)

A_draws  <- plogis(A_lin)
G_draws  <- plogis(G_lin)
C_draws  <- plogis(C_lin)
b_draws  <- exp(B_lin)
S_draws  <- plogis(S_lin)
M_draws  <- plogis(M_lin)
c0_draws <-      c0_lin

adj_age_mat <- matrix(
  rep(ppari_data$adj_age, nrow(A_draws)),
  nrow = nrow(A_draws),
  byrow = TRUE
)

## A. Reconstruct draws of response using _lin and covariate from data
reconst_lgt_ppari_A <- qlogis(
  A_draws *
    G_draws *
    C_draws *
    exp(
      -b_draws *
        exp(-c0_draws * M_draws * S_draws * adj_age_mat)
    )
)
## Summarize response draws for each observation
mu_lin_mean <- colMeans(reconst_lgt_ppari_A)

## Reconstruct draws using transformed parameter draws =====
## Note that these were obtained as vectors, explicitly
## preserving their joint distribution.

Araw <- fit_lgt_ppari_f |>
  spread_draws(nlp_Araw[obs]) |>
  pivot_wider(
    names_from = obs,
    names_prefix = "obs",
    values_from = nlp_Araw
  ) |>
  dplyr::select(starts_with("obs")) |>
  as.matrix() |>
  unname()

Craw <- fit_lgt_ppari_f |>
  spread_draws(nlp_Craw[obs]) |>
  pivot_wider(
    names_from = obs,
    names_prefix = "obs",
    values_from = nlp_Craw
  ) |>
  dplyr::select(starts_with("obs")) |>
  as.matrix() |>
  unname()

Graw <- fit_lgt_ppari_f |>
  spread_draws(nlp_Graw[obs]) |>
  pivot_wider(
    names_from = obs,
    names_prefix = "obs",
    values_from = nlp_Graw
  ) |>
  dplyr::select(starts_with("obs")) |>
  as.matrix() |>
  unname()

lnBraw <- fit_lgt_ppari_f |>
  spread_draws(nlp_lnBraw[obs]) |>
  pivot_wider(
    names_from = obs,
    names_prefix = "obs",
    values_from = nlp_lnBraw
  ) |>
  dplyr::select(starts_with("obs")) |>
  as.matrix() |>
  unname()

Sraw <- fit_lgt_ppari_f |>
  spread_draws(nlp_Sraw[obs]) |>
  pivot_wider(
    names_from = obs,
    names_prefix = "obs",
    values_from = nlp_Sraw
  ) |>
  dplyr::select(starts_with("obs")) |>
  as.matrix() |>
  unname()

Mraw <- fit_lgt_ppari_f |>
  spread_draws(nlp_Mraw[obs]) |>
  pivot_wider(
    names_from = obs,
    names_prefix = "obs",
    values_from = nlp_Mraw
  ) |>
  dplyr::select(starts_with("obs")) |>
  as.matrix() |>
  unname()

c0 <- fit_lgt_ppari_f |>
  spread_draws(nlp_c0[obs]) |>
  pivot_wider(
    names_from = obs,
    names_prefix = "obs",
    values_from = nlp_c0
  ) |>
  dplyr::select(starts_with("obs")) |>
  as.matrix() |>
  unname()

## B. Reconstruct draws of response using transformed parameter draws
## and covariate from data
reconst_lgt_ppari_B <- qlogis(
  plogis(Araw) *
    plogis(Graw) *
    plogis(Craw) *
    exp(
      -exp(lnBraw) *
        exp(
          -c0 *
            plogis(Mraw) *
            plogis(Sraw) *
            adj_age_mat
        )
    )
)

## Compare results =====

hist(post_pred_lgt_ppari - reconst_lgt_ppari_A, breaks = 200)
hist(post_pred_lgt_ppari - reconst_lgt_ppari_B, breaks = 200)
hist(reconst_lgt_ppari_A - reconst_lgt_ppari_B, breaks = 200)

# Reconstructions A and B seem to be the same down to "rounding error."
# Reconstructions differ from posterior predict due to normal residual with sigma = 0.39
# By using posterior_linpred for any of the parameters with the
# appropriate newdata I can get posterior_draws for any values of
# gdd, chill, m2tree and shape desired.

Reconstructions A and B areseem to be the same down to “rounding error.”
Reconstructions differ from posterior predict due to normal residual with sigma = 0.39
By using posterior_linpred for any of the parameters with the appropriate newdata I can get posterior_draws for any values of gdd, chill, m2tree and shape desired.

3

Thanks for following up and sharing the solution you found.