This document shows how to:
brms,posterior_predict() andposterior_linpred(transform = TRUE).The Michaelis-Menten mean function is
\[ \mu_i = \frac{V_{max} \cdot x_i}{Half_{max} + x_i} \]
Because we fit this model with a Gamma family and a log link:
\[ \mu_i = \exp(\eta_i) \]
library(brms)Loading required package: RcppLoading 'brms' package (version 2.23.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'The following object is masked from 'package:stats':
arset.seed(123)
n <- 200
light <- runif(n, 0, 100)
Vmax <- 12
Halfmax <- 25
mu <- (Vmax * light) / (Halfmax + light)
shape <- 8
rate <- shape / mu
biomass <- rgamma(n, shape = shape, rate = rate)
dat <- data.frame(
biomass = biomass,
light = light
)
head(dat) biomass light
1 4.558295 28.75775
2 8.659189 78.83051
3 10.302748 40.89769
4 7.070707 88.30174
5 8.741957 94.04673
6 2.032083 4.55565plot(dat$light, dat$biomass,
pch = 16,
col = adjustcolor("black", 0.5),
xlab = "Light",
ylab = "Biomass")
The bf(...) function here is specifying that we are including a non-linear equation. The Vmax+Halfmax~1 part is saying that both of those parameters are “fixed effects,” meaning they do not vary by group. The nl=TRUE is specifying that this is a non-linear model. We set priors with lb=0, because it doesn’t make sense for either of these parameters to be negative.
fit <- brm(
bf(
biomass ~ (Vmax * light) / (Halfmax + light),
Vmax + Halfmax ~ 1,
nl = TRUE
),
data = dat,
family = Gamma(link = "log"),
prior = c(
prior(normal(10, 5), nlpar = "Vmax", lb = 0),
prior(normal(5, 2), nlpar = "Halfmax", lb = 0)
),
chains = 4,
iter = 2000,
seed = 123
)Compiling Stan program...WARNING: Rtools is required to build R packages, but is not currently installed.
Please download and install the appropriate version of Rtools for 4.5.2 from
https://cran.r-project.org/bin/windows/Rtools/.Trying to compile a simple C fileRunning "C:/PROGRA~1/R/R-45~1.2/bin/x64/Rcmd.exe" SHLIB foo.c
using C compiler: 'gcc.exe (GCC) 14.3.0'
gcc -I"C:/PROGRA~1/R/R-45~1.2/include" -DNDEBUG -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/Rcpp/include/" -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/RcppEigen/include/" -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported" -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/StanHeaders/include/src/" -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/StanHeaders/include/" -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include "C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp" -std=c++1y -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -std=gnu2x -mfpmath=sse -msse2 -mstackrealign -c foo.c -o foo.o
cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
In file included from C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
from C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
from C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
from <command-line>:
C:/Users/trevorcaughlin/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
679 | #include <cmath>
| ^~~~~~~
compilation terminated.
make: *** [C:/PROGRA~1/R/R-45~1.2/etc/x64/Makeconf:289: foo.o] Error 1WARNING: Rtools is required to build R packages, but is not currently installed.
Please download and install the appropriate version of Rtools for 4.5.2 from
https://cran.r-project.org/bin/windows/Rtools/.Start sampling
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000244 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.44 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.98 seconds (Warm-up)
Chain 1: 0.975 seconds (Sampling)
Chain 1: 1.955 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000142 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.42 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1.009 seconds (Warm-up)
Chain 2: 0.905 seconds (Sampling)
Chain 2: 1.914 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000111 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.11 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 1.036 seconds (Warm-up)
Chain 3: 0.833 seconds (Sampling)
Chain 3: 1.869 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 7.5e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.959 seconds (Warm-up)
Chain 4: 0.924 seconds (Sampling)
Chain 4: 1.883 seconds (Total)
Chain 4: This is setting up a grid for values that we are going to predict across.
pred_seq <- seq(min(dat$light), max(dat$light), length.out = 200)
nd <- data.frame(light = pred_seq)Note that the posterior_predict is drawing posterior values from the gamma distribution, representing full uncertainty: sampling uncertainty and parameter uncertainty. These are the values you would expect to sample if you collected data again. In contrast, posterior_linpred is sampling the posterior draws for the mean, including parameter uncertainty online. Note that posterior_linpred is the same as what we have been doing manually, writing something like:
\[\mu_i = \exp\left(\frac{V_{max} \cdot x_i}{Half_{max} + x_i}\right)\] This is what we have been doing manually all along. The built-in function makes it easier, but still good to understand what the link functions are doing! The transform=T in posterior_linpred is specifying that we want output displayed on the scale of the response.
ppred <- posterior_predict(fit, newdata = nd)
lpred <- posterior_linpred(fit, newdata = nd, transform = TRUE)ppred_med <- apply(ppred, 2, median)
ppred_lwr <- apply(ppred, 2, quantile, probs = 0.025)
ppred_upr <- apply(ppred, 2, quantile, probs = 0.975)
lpred_med <- apply(lpred, 2, median)
lpred_lwr <- apply(lpred, 2, quantile, probs = 0.025)
lpred_upr <- apply(lpred, 2, quantile, probs = 0.975)plot(dat$light, dat$biomass,
pch = 16,
col = adjustcolor("black", 0.5),
xlab = "Light",
ylab = "Biomass")
polygon(
c(pred_seq, rev(pred_seq)),
c(ppred_lwr, rev(ppred_upr)),
col = adjustcolor("skyblue", 0.25),
border = NA
)
polygon(
c(pred_seq, rev(pred_seq)),
c(lpred_lwr, rev(lpred_upr)),
col = adjustcolor("tomato", 0.25),
border = NA
)
lines(pred_seq, lpred_med, lwd = 2, col = "red")
points(dat$light, dat$biomass,
pch = 16,
col = adjustcolor("black", 0.5))