## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 6,
  fig.height = 4,
  dpi = 96
)

## -----------------------------------------------------------------------------
library(int3ract)

set.seed(1402)
dat <- data.frame(x = rnorm(100), z = rnorm(100), w = rnorm(100))
dat$y <- dat$x + 0.5 * dat$z - 0.5 * dat$w +
  0.5 * dat$x * dat$z * dat$w + rnorm(100, sd = 4)

## -----------------------------------------------------------------------------
fit2 <- lm(y ~ x * z, data = dat)
jn2  <- JN(fit2, theta_1 = "x", theta_2 = "z")
jn2

## -----------------------------------------------------------------------------
summary(jn2)

## ----fig.alt = "Conditional effect of x across the range of z, with a confidence band shaded by significance and a histogram of observed z values below."----
plot(jn2, which = "x")

## -----------------------------------------------------------------------------
jn_regions(jn2)

## -----------------------------------------------------------------------------
head(as.data.frame(jn2), 3)

## -----------------------------------------------------------------------------
fit3 <- lm(y ~ x * z * w, data = dat)
jn3  <- JN(fit3, theta_1 = "x", theta_2 = "z", theta_3 = "w",
           range_size = 20)
jn3

## ----fig.width = 6, fig.height = 5, fig.alt = "Heatmap of the conditional effect of x over the two-dimensional grid of z and w, crosshatched where the effect is not significant, with marginal histograms of observed values along both axes."----
plot(jn3, which = "x")

## -----------------------------------------------------------------------------
jn3_fdr <- JN(fit3, theta_1 = "x", theta_2 = "z", theta_3 = "w",
              range_size = 20, control_fdr = TRUE)
head(jn_regions(jn3_fdr), 3)

## ----eval = requireNamespace("lme4", quietly = TRUE)--------------------------
d <- dat
d$g <- rep(letters[1:5], each = 20)
fit_mer <- lme4::lmer(y ~ x * z + (1 | g), data = d)
JN(fit_mer, theta_1 = "x", theta_2 = "z")

## ----bayes, eval = requireNamespace("MCMCpack", quietly = TRUE), message = FALSE----
post <- MCMCpack::MCMCregress(y ~ x * z, data = dat,
                              burnin = 500, mcmc = 2000, verbose = 0)

jnb <- JN(post, theta_1 = "x", theta_2 = "z",
          theta_1_vals = seq(-3, 3, 0.5),
          theta_2_vals = seq(-3, 3, 0.5))
jnb

## ----eval = requireNamespace("MCMCpack", quietly = TRUE), fig.alt = "Overlaid conditional posterior densities of the effect of x at a sequence of values of z."----
plot(jnb, which = "x")

## ----eval = requireNamespace("MCMCpack", quietly = TRUE), fig.alt = "Posterior mean effect of x across the range of z with a shaded credible band."----
plot(jnb, which = "x", type = "band")

## ----fig.alt = "The same conditional effect plot for x, drawn in different colours and without the data-density panel."----
plot(jn2, which = "x",
     style = jn_style(sig_color = "steelblue", non_sig_color = "grey70",
                      show_density = FALSE))

## -----------------------------------------------------------------------------
figs <- jn_plots(jn2)
names(figs)

## ----eval = FALSE-------------------------------------------------------------
# jn_save(jn2, folder = "figures", device = "pdf")

## ----eval = FALSE-------------------------------------------------------------
# # siena07() results carry estimates and a covariance matrix -> Wald tests
# JN(saom_fit, theta_1 = 4, theta_2 = 7, theta_int_12 = 12,
#    theta_1_vals = c(0, 6), theta_2_vals = c(-2, 2))
# 
# # sienaBayes() results carry draws -> conditional posteriors
# JN(bayes_fit, theta_1 = 4, theta_2 = 7, theta_int_12 = 12,
#    theta_1_vals = seq(0, 6, 1), theta_2_vals = seq(-2, 2, 1))

## ----eval = FALSE-------------------------------------------------------------
# jn_input.myfit <- function(object, theta_1, theta_2, theta_3 = NULL, ...) {
#   idx <- c(theta_1, theta_2, paste(theta_1, theta_2, sep = ":"))
#   jn_wald(coefficients = object$estimates[idx],
#           vcov         = object$covariance[idx, idx],
#           labels       = c(theta_1, theta_2))
# }

