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

## ----setup--------------------------------------------------------------------
library(glmnet)
library(rACMEMEEV)
set.seed(987654)

## -----------------------------------------------------------------------------
simulate_gamma_mix <- function(n, a1 = 3, a2 = 4, a3 = 5) {
  shape1 <- 2
  rate1 <- 2
  shape2 <- 3
  rate2 <- 3
  shape3 <- 5
  rate3 <- 5

  x <- rgamma(n, shape1, rate1)
  y <- rgamma(n, shape2, rate2)
  z <- rgamma(n, shape3, rate3)
  output <- a1 * x + a2 * y + a3 * z + rnorm(n)

  frame <- data.frame(
    list(
      output = output,
      x = x,
      y = y,
      z = z
    )
  )
  return(frame)
}

## -----------------------------------------------------------------------------
df <- simulate_gamma_mix(100)
head(df, 10)

## -----------------------------------------------------------------------------
model <- glm("output ~ x + y + z", data = df)
summary(model)

## -----------------------------------------------------------------------------
tainted_df <- df
tainted_df$x <- df$x + rnorm(100, mean = 2)
tainted_df$y <- df$y + rnorm(100, mean = 2)
tainted_df$z <- df$z + rnorm(100, mean = 2)
head(tainted_df, 10)

## -----------------------------------------------------------------------------
tainted_model <- glm("output ~ x + y + z", data = tainted_df)
summary(tainted_model)

## -----------------------------------------------------------------------------
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)

cat(sprintf(
  "X validity coefficient: %s\nY validity coefficient: %s\nZ validity coefficient: %s", x_v_coef, y_v_coef, z_v_coef
))

## -----------------------------------------------------------------------------
output <- acme_model(tainted_df, c("x", "y", "z"))

## -----------------------------------------------------------------------------
output$covariance_matrix

## -----------------------------------------------------------------------------
output$model

## -----------------------------------------------------------------------------
validity_coefficients <- c(x_v_coef, y_v_coef, z_v_coef)
lambda <- attenuation_matrix(output, c("x", "y", "z"), validity_coefficients)

## -----------------------------------------------------------------------------
lambda$matrix

## -----------------------------------------------------------------------------
model_output <- multivariate_model(
  "output ~ x + y + z",
  data = tainted_df,
  columns = c("x", "y", "z"),
  a_c_matrix = lambda$matrix,
  sds = lambda$sds,
  variances = lambda$variances,
  univariate = TRUE
)

## -----------------------------------------------------------------------------
summary(tainted_model)
summary(model)

## -----------------------------------------------------------------------------
summary(model_output$multivariate)

## -----------------------------------------------------------------------------
quantile(model_output$multivariate$x, 0.975)
quantile(model_output$multivariate$y, 0.975)
quantile(model_output$multivariate$z, 0.975)

## -----------------------------------------------------------------------------
plots <- plot_covariates(model_output, c("x", "y", "z"))
plots$x
plots$y
plots$z

## -----------------------------------------------------------------------------
df <- simulate_gamma_mix(100, a3 = 2)
head(df, 10)

## -----------------------------------------------------------------------------
tainted_df_lower_assoc <- df
tainted_df_lower_assoc$x <- df$x + rnorm(100, mean = 2)
tainted_df_lower_assoc$y <- df$y + rnorm(100, mean = 2)
tainted_df_lower_assoc$z <- df$z + rnorm(100, mean = 2)
head(tainted_df_lower_assoc, 10)

## -----------------------------------------------------------------------------
output <- acme_model(tainted_df_lower_assoc, c("x", "y", "z"))
lambda <- attenuation_matrix(output, c("x", "y", "z"), validity_coefficients)
model_output <- multivariate_model(
  "output ~ x + y + z",
  data = tainted_df_lower_assoc,
  columns = c("x", "y", "z"),
  a_c_matrix = lambda$matrix,
  sds = lambda$sds,
  variances = lambda$variances,
  univariate = TRUE
)

## -----------------------------------------------------------------------------
summary(model_output$multivariate)

## -----------------------------------------------------------------------------
traceplots(model_output$naive, c("x", "y", "z"))

## -----------------------------------------------------------------------------
acf_plots(model_output$naive, c("x", "y", "z"))