## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(BH)
library(rACMEMEEV)

## -----------------------------------------------------------------------------
x <- rgamma(100, shape = 1, rate = 0.01)
y <- rgamma(100, shape = 3, rate = 0.02)
z <- rgamma(100, shape = 3, rate = 0.3)

## -----------------------------------------------------------------------------
output <- rlnorm(100, meanlog = 3.5, sdlog = 0.2)

## -----------------------------------------------------------------------------
df <- data.frame(
  list(x = x, y = y, z = z, output = output)
)

head(df, 10)

x_v_coef <- generate_coefficient(1000, 0.4, 0.7, 0.95)
y_v_coef <- generate_coefficient(1000, 0.5, 0.7, 0.95)
z_v_coef <- generate_coefficient(1000, 0.3, 0.6, 0.95)

## -----------------------------------------------------------------------------
jags_output <- acme_model(df, c("x", "y", "z"))

## -----------------------------------------------------------------------------
stan_output <- acme_model(df, c("x", "y", "z"), stan = TRUE)

## -----------------------------------------------------------------------------
jags_output$model

## -----------------------------------------------------------------------------
stan_output$model

## -----------------------------------------------------------------------------
jags_output$covariance_matrix

## -----------------------------------------------------------------------------
stan_output$covariance_matrix

## -----------------------------------------------------------------------------
validity_coefficients <- c(x_v_coef, y_v_coef, z_v_coef)
jags_lambda <- attenuation_matrix(
  jags_output,
  c("x", "y", "z"),
  validity_coefficients
)
jags_model_output <- multivariate_model(
  "output ~ x + y + z",
  data = df,
  columns = c("x", "y", "z"),
  a_c_matrix = jags_lambda$matrix,
  sds = jags_lambda$sds,
  variances = jags_lambda$variances,
  univariate = TRUE
)

## -----------------------------------------------------------------------------
validity_coefficients <- c(x_v_coef, y_v_coef, z_v_coef)
stan_lambda <- attenuation_matrix(
  stan_output,
  c("x", "y", "z"),
  validity_coefficients,
  stan = TRUE
)
stan_model_output <- multivariate_model(
  "output ~ x + y + z",
  data = df,
  columns = c("x", "y", "z"),
  a_c_matrix = stan_lambda$matrix,
  sds = stan_lambda$sds,
  variances = stan_lambda$variances,
  univariate = TRUE
)

## -----------------------------------------------------------------------------
jags_plots <- plot_covariates(jags_model_output, c("x", "y", "z"))
jags_plots$x
jags_plots$y
jags_plots$z

## -----------------------------------------------------------------------------
stan_plots <- plot_covariates(stan_model_output, c("x", "y", "z"))
stan_plots$x
stan_plots$y
stan_plots$z

## -----------------------------------------------------------------------------
traceplots(stan_output$samples, c("x", "y", "z"), pre_model = TRUE, stan = TRUE)

## -----------------------------------------------------------------------------
traceplots(jags_output$samples, c("x", "y", "z"), pre_model = TRUE)