knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
jmbr
(pronounced jimber) is an R package to facilitate analyses using Just Another Gibbs Sampler (JAGS
).
It is part of the mbr family of packages.
The first part of the model is where priors, random effects and the relationships of interest are set in JAGS.
Example model:
library(jmbr) library(embr)
model <- model("model { # Priors alpha ~ dnorm(0, 10^-2) T(0,) beta1 ~ dnorm(0, 10^-2) beta2 ~ dnorm(0, 10^-2) beta3 ~ dnorm(0, 10^-2) # Random Effect log_sAnnual ~ dnorm(0, 10^-2) log(sAnnual) <- log_sAnnual for(i in 1:nAnnual) { bAnnual[i] ~ dnorm(0, sAnnual^-2) } # Prediction of Interest for (i in 1:length(Pairs)) { log(ePairs[i]) <- alpha + beta1 * Year[i] + beta2 * Year[i]^2 + beta3 * Year[i]^3 + bAnnual[Annual[i]] Pairs[i] ~ dpois(ePairs[i]) } }")
The new expression is written in R Code and is used to calculate derived parameters.
new_expr = " for (i in 1:length(Pairs)) { log(prediction[i]) <- alpha + beta1 * Year[i] + beta2 * Year[i]^2 + beta3 * Year[i]^3 + bAnnual[Annual[i]] fit[i] <- prediction[i] residual[i] <- res_pois(Pairs[i], fit[i]) }"
This section modifies a data frame to the form it will be passed to the analysis code. The modified data is passed in list form.
modify_data = function(data) { data <- data |> select(-Eyasses) data }
Select data is a named list specifying the columns to select and their associated classes and values as well as transformations and scaling options.
Random effects gets the random effects definitions for an object as a named list, where bAnnual
refers to the column name Annual
in the data.
select_data = list("Pairs" = c(15L, 200L), "Year*" = 1L, Annual = factor()), random_effects = list(bAnnual = "Annual"),
All parameters in the data that are included in the model must be listed here.
- If there are values in the Pairs column outside of the specified range, including NA's, an error is thrown.
- "Year*" = 1L
indicates Year is of class integer.
Year-
= subtracts the minimum value (the first year) Year+
= subtracts the average value (centering)Year*
= subtracts the average value and divides by the SD (standardizing)Initial values of a parameter can be set prior to the analysis as a single argument function taking the modified data and returning a named list of initial values.
Unspecified initial values for each chain are drawn from the prior distributions.
gen_inits = function(data) { inits <- list() inits$ePairs <- data$Pairs + 1 inits },
At the end of the script is where the thinning rate is set, i.e. how much the MCMC chains should be thinned out before storing them.
Setting nthin = 1
corresponds to keeping all values.
Setting nthin = 100
would result in keeping every 100th value and discarding all other values.
model <- model("model { alpha ~ dnorm(0, 10^-2) beta1 ~ dnorm(0, 10^-2) beta2 ~ dnorm(0, 10^-2) beta3 ~ dnorm(0, 10^-2) log_sAnnual ~ dnorm(0, 10^-2) log(sAnnual) <- log_sAnnual for(i in 1:nAnnual) { bAnnual[i] ~ dnorm(0, sAnnual^-2) } for (i in 1:length(Pairs)) { log(ePairs[i]) <- alpha + beta1 * Year[i] + beta2 * Year[i]^2 + beta3 * Year[i]^3 + bAnnual[Annual[i]] Pairs[i] ~ dpois(ePairs[i]) } }", new_expr = " for (i in 1:length(Pairs)) { log(prediction[i]) <- alpha + beta1 * Year[i] + beta2 * Year[i]^2 + beta3 * Year[i]^3 + bAnnual[Annual[i]] fit[i] <- prediction[i] residual[i] <- res_pois(Pairs[i], fit[i]) }", modify_data = function(data) { data$nObs <- length(data$Annual) data }, select_data = list("Pairs" = c(15L, 200L), "Year*" = 1L, Annual = factor()), random_effects = list(bAnnual = "Annual"), nthin = 10L) data <- bauw::peregrine data$Annual <- factor(data$Year) set_analysis_mode("report")
Analysis mode can be set depending on the desired output.
set_analysis_mode("report")
Modes:
quick
: To quickly test code runs.report
: To produce results for a report. paper
: To produce results for a peer-reviewed paper.Analyse or reanalyse the model.
analysis <- analyse(model, data = data) analysis <- reanalyse(analysis)
Analysis Table:
n
autocorrelated samples.ess
.par(mar=c(1, 1, 1, 1)) plot(analysis)
Coefficient Table
Summary table of the posterior probability distribution.
coef(analysis)
The estimate is the median by default.
The zscore is $mean / sd$.
coef(analysis, simplify = TRUE)
The s-value is the suprisal value, which is a measure of directionality with respect to zero.
The s-value is zero (unsurprising) when p-value = 1.0 and increases exponentially as p approaches zero.
$$s = -log_2(p-value)$$ Example: How surprising it would be to throw 10 heads in 10 coin tosses.
A larger s-value provides more evidence against the null hypothesis and support that the data is in the direction of the posterior.
Example prediction:
Make predictions by varying Year
with other predictors, including the random effect of Annual
held constant.
year <- predict(analysis, new_data = "Year") library(ggplot2) ggplot(data = year, aes(x = Year, y = estimate)) + geom_point(data = bauw::peregrine, aes(y = Pairs)) + geom_line() + geom_line(aes(y = lower), linetype = "dotted") + geom_line(aes(y = upper), linetype = "dotted") + expand_limits(y = 0)
Predict()
queries the model and tells you what the expected number would be for that combination of values specified by [new data].
The example below would calculate the annual number of pairs for a typical number of fledged young of 50 (if Eyasses
was a parameter in the model).
year <- new_data(data, "Year", ref = list(Eyasses = 50L), obs_only = TRUE) %>% predict(analysis, new_data = ., ref_data = ref)
Arguments
new_data
: Creates a new data frame to calculate the predictions for.ref_data
: A data frame with 1 row indicating the reference values for calculating the effects size. ref = list(Eyasses = 50L)
.Predict can also take the form:
year <- predict(analysis, new_data = character(0), term = "ePairs")
Where term
calls the string of a term in the [new expression] of the model. By default it is the prediction[i]
.
Creates a new data frame to be passed to the [predict] function.
The idea is that most variables are held constant at a reference level while the variables of interest vary across their range.
year <- new_data( data, seq = "Year", ref = list(Eyasses = 50L), obs_only = TRUE) %>% predict(analysis, new_data = ., ref_data = ref)
Arguments
seq
: The name of columns to vary over. In this example, Year
. seq
then all levels of the factor are represented.ref
: A named list of reference values for variables not in seq
. Eyasses
constant at 50L.obs_only
: A list of character vectors indicating the sets of variables to only keep combinations for, i.e. combinations that were observed in the data. obs_only = TRUE
then obs_only
is set to be seq
.length_out
: A count indicating the length of numeric/integer sequences.seq
then all levels of the factor are represented (length_out
is ignored). obs_only
, then only observed factor levels are represented in sequences.Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.