diff --git a/.Rhistory b/.Rhistory
index 66c2813..a050964 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,512 +1,512 @@
-set.seed(1)
-X <- sample(1:2, N * 2, replace = TRUE) |> matrix(ncol = 2)
-y <- rnorm(N)
-for (i in 1:N) {
-if ((X[i, 1] == 2) & (X[i, 2] == 2)) {
-y[i] = y[i] + 1
+for (i in 1:length(y)) {
+# if (x[i] > 0) {
+# f[i] <- phi * dpois(x[i] - 1, lambda) + (1 - phi) * dpois(x[i], lambda)
+# } else {
+# f[i] <- dpois(x[i], lambda)
+# }
 }
+f
 }
-example_01 <- tibble::tibble(y = y, x1 = X[, 1], x2 = X[, 2])
-usethis::use_data(example_01, example_01, overwrite = TRUE)
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-library(ANOVABNPTestR)
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-pkgdown::build_site()
-?pkgdown::build_site
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-library(ANOVABNPTestR)
-options(pillar.print_max = 3)
-example_01
+dberpo(1L, 0.5, 0.5)
+if (x[i] > 0) {
+# f[i] <- phi * dpois(x[i] - 1, lambda) + (1 - phi) * dpois(x[i], lambda)
+# } else {
+# f[i] <- dpois(x[i], lambda)
+}
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+f <- x
+for (i in 1:length(y)) {
+if (x[i] > 0) {
+# f[i] <- phi * dpois(x[i] - 1, lambda) + (1 - phi) * dpois(x[i], lambda)
+# } else {
+# f[i] <- dpois(x[i], lambda)
+}
+}
+f
+}
+dberpo(1L, 0.5, 0.5)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+f <- x
+for (i in 1:length(y)) {
+print(x[i])
+# if (x[i] > 0) {
+# f[i] <- phi * dpois(x[i] - 1, lambda) + (1 - phi) * dpois(x[i], lambda)
+# } else {
+# f[i] <- dpois(x[i], lambda)
+# }
+}
+f
+}
+dberpo(1L, 0.5, 0.5)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+f <- x
+for (i in 1:length(x)) {
+print(x[i])
+# if (x[i] > 0) {
+# f[i] <- phi * dpois(x[i] - 1, lambda) + (1 - phi) * dpois(x[i], lambda)
+# } else {
+# f[i] <- dpois(x[i], lambda)
+# }
+}
+f
+}
+dberpo(1L, 0.5, 0.5)
 library(ANOVABNPTestR)
-options(pillar.print_max = 3)
-example_01
-head(example_01)
-head(example_01)
-options(pillar.print_max = 3)
-pkgdown::build_site(preview = FALSE)
-?ggplot2::facet_grid
+library(dplyr)
+library(ggplot2)
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+f <- x
+for (i in 1:length(x)) {
+if (x[i] > 0) {
+f[i] <- phi * dpois(x[i] - 1, lambda) + (1 - phi) * dpois(x[i], lambda)
+} else {
+f[i] <- dpois(x[i], lambda)
+}
+}
+f
+}
+dberpo(1L, 0.5, 0.5)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) #|>
 library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
-ggplot2::ggplot(example1, aes(x = y)) +
-ggplot2::geom_histogram() +
-ggplot2::facet_grid(vars(x1), vars(x2))
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_histogram() +
-ggplot2::facet_grid(vars(x1), vars(x2))
-?ggplot2::facet_wrap
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2))
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), ncol = 1)
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), nrow = 1)
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), ncol = 1)
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), ncol = 1, strip.position = "right")
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), ncol = 1, strip.position = "right", labeller = "label_both")
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), now = 1, labeller = "label_both")
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), nrow = 1, labeller = "label_both")
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), nrow = 1, labeller = "label_both", dir = "v")
-ggplot2::ggplot(example_01, aes(x = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), nrow = 1, labeller = "label_both", dir = "h")
-ggplot2::ggplot(example_01, aes(y = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(vars(x1, x2), nrow = 1, labeller = "label_both")
-ggplot2::ggplot(example_01, aes(y = y)) +
-ggplot2::geom_boxplot() +
-ggplot2::facet_wrap(
-nrow = 1,
-vars(x1, x2),
-labeller = "label_both"
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
+n <- 1000
+phi <- 0.5
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) #|>
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
 ) +
-theme_classic()
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-pkgdown::build_site(preview = FALSE)
-ggplot(example_01, aes(y = y)) +
-geom_boxplot() +
-facet_wrap(
-nrow = 1,
-vars(x1, x2),
-labeller = "label_both"
+ggplot2::geom_line()
+library(ANOVABNPTestR)
+library(dplyr)
+library(ggplot2)
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
+n <- 1000
+phi <- 0.5
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
 ) +
-theme_classic() +
-theme(panel.spacing = unit(2, "lines"))
-ggplot(example_01, aes(y = y)) +
-geom_boxplot() +
-facet_wrap(
-nrow = 1,
-vars(x1, x2),
-labeller = "label_both"
+ggplot2::geom_line()
+n <- 1000
+phi <- 0.2
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
 ) +
-theme_classic() +
-theme(panel.spacing = unit(1, "lines"))
-theme(panel.spacing = unit(1, "lines"))
-pkgdown::build_site(preview = FALSE)
-library(ANOVABNPTestR)
-ANOVABNPTestR::setup()
-setup()
-library(ANOVABNPTestR)
-ANOVABNPTestR::setup()
-library(ANOVABNPTestR)
-library(ANOVABNPTestR)
-library(ANOVABNPTestR)
-ANOVABNPTestR::setup()
-library(ANOVABNPTestR)
-library(ANOVABNPTestR)
-devtools::load_all()
-install.packages("devtools")
-devtools::load_all()
-library(ANOVABNPTestR)
-library(ANOVABNPTestR)
+ggplot2::geom_line()
 library(ANOVABNPTestR)
+library(dplyr)
+library(ggplot2)
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
 library(ANOVABNPTestR)
-library(ANOVABNPTestR)
-D <- 2
-N <- 1000
-X <- sample(1:2, N * D, replace = TRUE) |> matrix(nrow = N, ncol = D)
-y <- runif(N) < 0.75
-library(ANOVABNPTestR)
-anova_bnp_normal_bernoulli(y, X)
-library(ANOVABNPTestR)
-anova_bnp_normal_bernoulli(y, X)
-anova_bnp_bernoulli(y, X)
-pkgdown::build_site(preview = FALSE)
-library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
-remotes::install_github("igutierrezm/ANOVABNPTestR")
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
 library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
-foodspoilage <- read_excel("FoodSpoilage.xlsx")
-foodspoilage <- readxl::read_excel("FoodSpoilage.xlsx")
-# Transforming covariate to integer type
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula1"] <- 1
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula2"] <- 2
-foodspoilage$Preservative <- as.integer(foodspoilage$Preservative)
-# Model
-yvec <- foodspoilage$Spoilage
-Xmat <- foodspoilage[,1:4] |> as.matrix()
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-ANOVABNPTestR::setup()
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-View(foodspoilage)
-# Model
-yvec <- foodspoilage$Spoilage == 1
-Xmat <- foodspoilage[,1:4] |> as.matrix()
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-# Transforming covariate to integer type
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula1"] <- 1L
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula2"] <- 2L
-foodspoilage$Preservative <- as.integer(foodspoilage$Preservative)
-# Model
-yvec <- foodspoilage$Spoilage == 1
-Xmat <- foodspoilage[,1:4] |> as.matrix()
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-Xmat
-Xmat <- foodspoilage[,1:4] |> as.matrix() |> as.integer()
-Xmat
-Xmat <- foodspoilage[,1:4] |> as.matrix()
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-# Model
-yvec <- foodspoilage$Spoilage == 1
-Xmat <- foodspoilage[,1:4] |> as.matrix()
-mode(Xmat) <- "integer"
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-View(foodspoilage)
-foodspoilage <- as.integer(foodspoilage)
-foodspoilage <- as_integer(foodspoilage)
-foodspoilage <- lapply(as.integer, foodspoilage)
-foodspoilage <- lapply(foodspoilage, as.integer)
-foodspoilage <- lapply(foodspoilage, as.integer) |> as.data.frame()
-View(foodspoilage)
-# Model
-yvec <- foodspoilage$Spoilage == 1
-Xmat <- foodspoilage[,1:4] |> as.matrix()
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-View(Xmat)
-unique(yvec)
-unique(Xmat)
-fit <- anova_bnp_bernoulli(yvec, Xmat)
-yvec <- rnorm(length(yvec))
-Xmat <- foodspoilage[, 1:4] |> as.matrix()
-fit <- anova_bnp_normal(yvec, Xmat)
-View(Xmat)
-View(Xmat)
-View(Xmat)
-write.csv(Xmat, file = "Xmat.csv", row.names = FALSE)
-write.csv(yvec, file = "yvec.csv", row.names = FALSE)
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
+n <- 1000
+phi <- 0.2
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+devtools::load_all(".")
+library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
+n <- 1000
+phi <- 0.2
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
+) +
+ggplot2::geom_line()
+devtools::load_all(".")
 library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
+library(tidyr)
+ANOVABNPTestR::setup()
+options(download.file.method = "wininet")
+ANOVABNPTestR::setup()
+1+1
 library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
+library(tidyr)
+options(download.file.method = "wininet")
 ANOVABNPTestR::setup()
-foodspoilage <- readxl::read_excel("FoodSpoilage.xlsx")
-# Transforming covariate to integer type
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula1"] <- 1L
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula2"] <- 2L
-foodspoilage <- lapply(foodspoilage, as.integer) |> as.data.frame()
-# Model
-yvec <- foodspoilage$Spoilage == 1
-Xmat <- foodspoilage[, 1:4] |> as.matrix()
-fit <- anova_bnp_normal(yvec, Xmat)
-yvec <- rnorm(length(yvec))
-Xmat <- foodspoilage[, 1:4] |> as.matrix()
-fit <- anova_bnp_normal(yvec, Xmat)
-hypothesis_post_simple(fit)
-hypothesis_post_interaction(fit)
-predictive_plot_simple(fit, 1) + theme_classic()
+devtools::load_all(".")
 library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
+library(tidyr)
+options(download.file.method = "wininet")
 ANOVABNPTestR::setup()
+devtools::load_all(".")
 library(ANOVABNPTestR)
+library(dplyr)
 library(ggplot2)
+library(tidyr)
 ANOVABNPTestR::setup()
-foodspoilage <- readxl::read_excel("FoodSpoilage.xlsx")
-# Transforming covariate to integer type
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula1"] <- 1L
-foodspoilage$Preservative[foodspoilage$Preservative=="Formula2"] <- 2L
-foodspoilage <- lapply(foodspoilage, as.integer) |> as.data.frame()
-# Model
-yvec <- foodspoilage$Spoilage == 1
-yvec <- rnorm(length(yvec))
-Xmat <- foodspoilage[, 1:4] |> as.matrix()
-fit <- anova_bnp_normal(yvec, Xmat)
-hypothesis_post_simple(fit)
-hypothesis_post_interaction(fit)
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) + theme_classic()
-fit
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1) + theme_classic()
+library(ANOVABNPTestR)
 library(dplyr)
-predictive_plot_simple(fit, 1) + theme_classic()
-group_codes(fit)
-?c_accros
-?c_accross
-?c_across
-group_codes(fit) |>
-mutate(touse = pmax(c_across(-c("group", var1))))
-group_codes(fit) |>
-mutate(touse = pmax(c_across(-c("group", "x1"))))
-group_codes(fit)
-pmax(c_across(-c("group", "x1")))
-dense_rank(c(2, 3, 4, 2, 1))
-a <- group_codes(fit)
-a
-a %>%
-mutate(touse = pmax(c_across(-c("group", "x1"))))
-a %>%
-mutate(touse = sum(c_across(-c("group", "x1"))))
-a %>%
-rowwise() %>%
-mutate(touse = sum(c_across(-c("group", "x1"))))
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-a <- group_codes(fit)
-a %>%
-rowwise() %>%
-mutate(touse = sum(c_across(-c("group", "x1"))))
-a %>%
-rowwise() %>%
-mutate(touse = pmax(c_across(-c("group", "x1"))))
-a %>%
-rowwise() %>%
-mutate(touse = sd(c_across(-c("group", "x1"))))
-pmax(c(1, 2, 1))
-max(c(1, 2, 1))
-a %>%
-rowwise() %>%
-mutate(touse = max(c_across(-c("group", "x1"))))
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) + theme_classic()
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-?mutate_at
-min(c(1,2, 4,-1))
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-predictive_plot_simple(fit, 1)
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1)
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) + theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) + theme_classic()
-group_codes(fit)
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-View(foodspoilage)
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-a <- predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-View(a)
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-a <- predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-View(a)
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-a <- predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-a <- predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-a
-mutate(
-a,
-touse = max(c_across(-c("group", "x1", "x1"))))
-mutate(
-a,
-touse = max(c_across(-c("group", "x1", "x2"))))
-mutate(
-a,
-touse = pmax(c_across(-c("group", "x1", "x2"))))
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-a <- predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-a
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-a <- predictive_plot_interaction(fit, 1, 2) #+ theme_classic()
-a
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) + theme_classic()
-source("R/predictive_plot_simple.R")
-source("R/predictive_plot_interaction.R")
-predictive_plot_simple(fit, 1) + theme_classic()
-predictive_plot_interaction(fit, 1, 2) + theme_classic()
-ggplot(example_01, aes(y = y)) +
-geom_boxplot() +
-facet_wrap(
-nrow = 1,
-vars(x1, x2),
-labeller = "label_both"
+library(ggplot2)
+library(tidyr)
+devtools::load_all(".")
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
+n <- 1000
+phi <- 0.2
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
+) +
+ggplot2::geom_line()
+n <- 1000
+phi <- 0.8
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
 ) +
-theme_classic() +
-theme(panel.spacing = unit(1, "lines"))
-fit <- anova_bnp_normal(yvec, Xmat)
-yvec <- example_01[[1]]
-Xmat <- example_01[, 2:3] |> as.matrix()
-fit <- anova_bnp_normal(yvec, Xmat)
-hypothesis_post_simple(fit)
-hypothesis_post_interaction(fit)
-predictive_plot_simple(fit, 1) +
-theme_classic()
-predictive_plot_simple(fit, 1) +
-theme_classic()
-predictive_plot_interaction(fit, 1, 2) +
-theme_classic()
+ggplot2::geom_line()
+n <- 2000
+phi <- 0.8
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
+) +
+ggplot2::geom_line()
+n <- 10000
+phi <- 0.8
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <- anova_bnp_berpoi(yvec, Xmat, iter = 20000L, warmup = 10000L)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
+) +
+ggplot2::geom_line()
 library(ANOVABNPTestR)
-data("example_01", package = "ANOVABNPTestR")
-?pull
-a <- data.frame(a = 1)
-pull(a, "a")
 library(dplyr)
-pull(a, "a")
-library(ANOVABNPTestR)
-devtools::load_all(".")
-load("~/r-projects/ANOVABNPTestR/data/example_01.rda")
-View(example_01)
-N <- 1000
-y <- rpoi(1, 1000)
-y <- rpois(1, 1000)
-y <- rpois(1000, 1)
-N <- 1000
-x <- 1 + (runif(N) < 0.5, 2)
-N <- 1000
-x <- 1 + (runif(N, 2) < 0.5)
-x <- 1 + (runif(2 * N) < 0.5)
-N <- 1000
-x <- 1 + (runif(2 * N) < 0.5)
-X <- matrix(x, nrow = N, ncol = 2)
-y <- rpois(N, 1) + (X[, 1] * X[, 2])
-N <- 1000
-x <- 1 + (runif(2 * N) < 0.5)
-Xmat <- matrix(x, nrow = N, ncol = 2)
-yvec <- rpois(N, 1) + (Xmat[, 1] * Xmat[, 2])
-fit <- anova_bnp_berpoi(yvec, Xmat)
-N <- 1000
-x <- as.integer(1 + (runif(2 * N) < 0.5))
-typeof(x)
-typeof(y)
-N <- 1000
-x <- 1 + (runif(2 * N) < 0.5)
-Xmat <- matrix(x, nrow = N, ncol = 2)
-yvec <- rpois(N, 1) + (Xmat[, 1] * Xmat[, 2])
-fit <- anova_bnp_berpoi(as.integer(yvec), as.integer(Xmat))
-as.integer(Xmat)
-N <- 1000
-x <- as.integer(1 + (runif(2 * N) < 0.5))
-Xmat <- matrix(x, nrow = N, ncol = 2)
-yvec <- rpois(N, 1) + (Xmat[, 1] * Xmat[, 2])
-N <- 1000
-x <- as.integer(1 + (runif(2 * N) < 0.5))
-Xmat <- matrix(x, nrow = N, ncol = 2)
-yvec <- rpois(N, 1) + (Xmat[, 1] * Xmat[, 2])
-fit <- anova_bnp_berpoi(yvec, Xmat)
-hypothesis_post_simple(fit)
-N <- 1000
-x <- as.integer(1 + (runif(2 * N) < 0.5))
-Xmat <- matrix(x, nrow = N, ncol = 2)
-yvec <- rpois(N, 1) + (Xmat[, 1] == 2) * (Xmat[, 2] == 2)
-fit <- anova_bnp_berpoi(yvec, Xmat)
-hypothesis_post_simple(fit)
-hypothesis_post_interaction(fit)
+library(ggplot2)
+library(tidyr)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
+n <- 10000
+phi <- 0.8
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 0.1,
+b = 0.1,
+a1 = mean(yvec),
+b1 = 1,
+alpha0 = 1.0,
+beta0 = 1.0,
+lb = 0,
+ub = 10
+)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 0.1,
+b = 0.1,
+a1 = mean(yvec),
+b1 = 1,
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+n <- 10000
+phi <- 0.8
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 0.1,
+b = 0.1,
+a1 = numeric(mean(yvec)),
+b1 = numeric(1),
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 0.1,
+b = 0.1,
+a1 = numeric(mean(yvec)) + 1,
+b1 = numeric(1),
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 0.1,
+b = 0.1,
+a1 = numeric(mean(yvec)) + 1,
+b1 = numeric(1),
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+?anova_bnp_berpoi
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 1,
+b = 1,
+a1 = numeric(mean(yvec)) + 1,
+b1 = numeric(1),
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 1,
+b = 1,
+a1 = round(mean(yvec)) + 1,
+b1 = 1L,
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 1,
+b = 1,
+a1 = mean(yvec),
+b1 = 1,
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 1.0,
+b = 1.0,
+a1 = 1.0,
+b1 = 1.0,
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
+) +
+ggplot2::geom_line()
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 0.1,
+b = 0.1,
+a1 = 1.0,
+b1 = 1.0,
+alpha0 = 1.0,
+beta0 = 1.0,
+)
+fit$f_post |>
+dplyr::rename(fh = f) |>
+dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+tidyr::pivot_longer(cols = fh:f0) |>
+ggplot2::ggplot(
+ggplot2::aes(x = y, color = name, y = value)
+) +
+ggplot2::geom_line()
+fit <-
+anova_bnp_berpoi(
+yvec,
+Xmat,
+iter = 2000L,
+warmup = 1000L,
+rho = 1.0,
+a = 0.1,
+b = 0.1,
+a1 = mean(yvec),
+b1 = 1.0,
+alpha0 = 1.0,
+beta0 = 1.0,
+)
diff --git a/DESCRIPTION b/DESCRIPTION
index bdf6f66..54ef2c9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -19,7 +19,7 @@ Imports:
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.0
+RoxygenNote: 7.2.3
 Suggests:
     rmarkdown,
     knitr
diff --git a/Manifest.toml b/Manifest.toml
new file mode 100644
index 0000000..929954b
--- /dev/null
+++ b/Manifest.toml
@@ -0,0 +1,443 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.7.3"
+manifest_format = "2.0"
+
+[[deps.ANOVADDPTest]]
+deps = ["DPMNeal3", "DataFrames", "Distributions", "ExtractMacro", "Random", "SpecialFunctions", "StatsBase", "StatsModels"]
+path = "C:\\Users\\Iegutierrezm\\.julia\\dev\\ANOVADDPTest"
+uuid = "39310853-6f77-48f5-bb9a-ec22dd77aa3b"
+version = "0.1.0"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.Calculus]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
+uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
+version = "0.5.1"
+
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "e7ff6cadf743c098e08fca25c91103ee4303c9bb"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.15.6"
+
+[[deps.ChangesOfVariables]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "38f7a08f19d8810338d4f5085211c7dfa5d5bdd8"
+uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+version = "0.1.4"
+
+[[deps.Compat]]
+deps = ["Dates", "LinearAlgebra", "UUIDs"]
+git-tree-sha1 = "3ca828fe1b75fa84b021a7860bd039eaea84d2f2"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "4.3.0"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[deps.Crayons]]
+git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
+uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+version = "4.1.1"
+
+[[deps.DPMNeal3]]
+deps = ["Distributions", "Parameters", "Random"]
+git-tree-sha1 = "0fc231ebece26a068ab19226412cb5b2f927ef99"
+uuid = "af0b21b6-ea01-48a9-a6c6-8ff6351a99f2"
+version = "0.1.0"
+
+[[deps.DataAPI]]
+git-tree-sha1 = "46d2680e618f8abd007bce0c3026cb0c4a8f2032"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.12.0"
+
+[[deps.DataFrames]]
+deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SnoopPrecompile", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
+git-tree-sha1 = "5b93f1b47eec9b7194814e40542752418546679f"
+uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+version = "1.4.2"
+
+[[deps.DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.13"
+
+[[deps.DataValueInterfaces]]
+git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
+uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
+version = "1.0.0"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[deps.DensityInterface]]
+deps = ["InverseFunctions", "Test"]
+git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
+uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
+version = "0.4.0"
+
+[[deps.Distributions]]
+deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
+git-tree-sha1 = "04db820ebcfc1e053bd8cbb8d8bccf0ff3ead3f7"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.25.76"
+
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "c36550cb29cbe373e95b3f40486b9a4148f89ffd"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.9.2"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[deps.DualNumbers]]
+deps = ["Calculus", "NaNMath", "SpecialFunctions"]
+git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
+uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
+version = "0.6.8"
+
+[[deps.ExtractMacro]]
+git-tree-sha1 = "35616bd2fb03730d2700ebe71870af7169d9e6a2"
+uuid = "f86d3d12-fd5b-522c-99e9-61577282a1e9"
+version = "1.0.0"
+
+[[deps.FileWatching]]
+uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
+
+[[deps.FillArrays]]
+deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
+git-tree-sha1 = "802bfc139833d2ba893dd9e62ba1767c88d708ae"
+uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
+version = "0.13.5"
+
+[[deps.Formatting]]
+deps = ["Printf"]
+git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
+uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
+version = "0.4.2"
+
+[[deps.Future]]
+deps = ["Random"]
+uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
+
+[[deps.HypergeometricFunctions]]
+deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions", "Test"]
+git-tree-sha1 = "709d864e3ed6e3545230601f94e11ebc65994641"
+uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
+version = "0.3.11"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.8"
+
+[[deps.InvertedIndices]]
+git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
+uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
+version = "1.1.0"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.1"
+
+[[deps.IteratorInterfaceExtensions]]
+git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
+uuid = "82899510-4779-5014-852e-03e436cf321d"
+version = "1.0.0"
+
+[[deps.JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.4.1"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.LogExpFunctions]]
+deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "94d9c52ca447e23eac0c0f074effbcd38830deb5"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.18"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[deps.Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "1.0.2"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[deps.NaNMath]]
+deps = ["OpenLibm_jll"]
+git-tree-sha1 = "a7c3d1da1189a1c2fe843a3bfa04d18d20eb3211"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "1.0.1"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[deps.OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[deps.PDMats]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "cf494dca75a69712a72b80bc48f59dcf3dea63ec"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.11.16"
+
+[[deps.Parameters]]
+deps = ["OrderedCollections", "UnPack"]
+git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
+uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+version = "0.12.3"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[deps.PooledArrays]]
+deps = ["DataAPI", "Future"]
+git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7"
+uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
+version = "1.4.2"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
+
+[[deps.PrettyTables]]
+deps = ["Crayons", "Formatting", "Markdown", "Reexport", "StringManipulation", "Tables"]
+git-tree-sha1 = "460d9e154365e058c4d886f6f7d6df5ffa1ea80e"
+uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+version = "2.1.2"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[deps.QuadGK]]
+deps = ["DataStructures", "LinearAlgebra"]
+git-tree-sha1 = "97aa253e65b784fd13e83774cadc95b38011d734"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.6.0"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.Reexport]]
+git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.2.2"
+
+[[deps.Rmath]]
+deps = ["Random", "Rmath_jll"]
+git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.7.0"
+
+[[deps.Rmath_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
+uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
+version = "0.3.0+0"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.ShiftedArrays]]
+git-tree-sha1 = "503688b59397b3307443af35cd953a13e8005c16"
+uuid = "1277b4bf-5013-50f5-be3d-901d8477a67a"
+version = "2.0.0"
+
+[[deps.SnoopPrecompile]]
+git-tree-sha1 = "f604441450a3c0569830946e5b33b78c928e1a85"
+uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+version = "1.0.1"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[deps.SortingAlgorithms]]
+deps = ["DataStructures"]
+git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "1.0.1"
+
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "d75bda01f8c31ebb72df80a46c88b25d1c79c56d"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.1.7"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.StatsAPI]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "f9af7f195fb13589dd2e2d57fdb401717d2eb1f6"
+uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+version = "1.5.0"
+
+[[deps.StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
+git-tree-sha1 = "d1bf48bfcc554a3761a133fe3a9bb01488e06916"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.33.21"
+
+[[deps.StatsFuns]]
+deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
+git-tree-sha1 = "5783b877201a82fc0014cbf381e7e6eb130473a4"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "1.0.1"
+
+[[deps.StatsModels]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Printf", "REPL", "ShiftedArrays", "SparseArrays", "StatsBase", "StatsFuns", "Tables"]
+git-tree-sha1 = "a5e15f27abd2692ccb61a99e0854dfb7d48017db"
+uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
+version = "0.6.33"
+
+[[deps.StringManipulation]]
+git-tree-sha1 = "46da2434b41f41ac3594ee9816ce5541c6096123"
+uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e"
+version = "0.3.0"
+
+[[deps.SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[deps.TableTraits]]
+deps = ["IteratorInterfaceExtensions"]
+git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
+uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
+version = "1.0.1"
+
+[[deps.Tables]]
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
+git-tree-sha1 = "c79322d36826aa2f4fd8ecfa96ddb47b174ac78d"
+uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+version = "1.10.0"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[deps.UnPack]]
+git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
+uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+version = "1.0.2"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
diff --git a/NAMESPACE b/NAMESPACE
index 9d95016..c83c53b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,8 @@ export(hypothesis_post_simple)
 export(predictive_plot_interaction)
 export(predictive_plot_simple)
 export(setup)
+export(shift_plot)
+export(shift_post)
 importFrom(JuliaConnectoR,juliaEval)
 importFrom(JuliaConnectoR,juliaImport)
 importFrom(dplyr,across)
@@ -23,6 +25,7 @@ importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,pull)
 importFrom(dplyr,rowwise)
+importFrom(ggplot2,aes)
 importFrom(ggplot2,aes_string)
 importFrom(ggplot2,geom_line)
 importFrom(ggplot2,ggplot)
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..e88973f
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+ANOVADDPTest = "39310853-6f77-48f5-bb9a-ec22dd77aa3b"
diff --git a/R/anova_bnp_normal.R b/R/anova_bnp_normal.R
index 8a0093d..f228f0a 100644
--- a/R/anova_bnp_normal.R
+++ b/R/anova_bnp_normal.R
@@ -38,6 +38,7 @@ anova_bnp_normal <- function(
   out <-
     list(
       f_post = fit$fpost,
+      shift_post = fit$shiftpost,
       gamma_post = fit$group_probs,
       group_codes = fit$group_codes,
       hypothesis_post_simple = fit$effects1,
diff --git a/R/f_post.R b/R/f_post.R
index 0361c94..824fe48 100644
--- a/R/f_post.R
+++ b/R/f_post.R
@@ -11,3 +11,17 @@
 f_post <- function(fit) {
   fit$f_post
 }
+
+#' Return the shift functions
+#'
+#' @param fit An object of class anova_bnp_model.
+#' @return A tibble with the following variables:
+#' \itemize{
+#'   \item group the group id.
+#'   \item y the value of the response variable.
+#'   \item the shift function evaluated at y, given the group.
+#' }
+#' @export
+shift_post <- function(fit) {
+  fit$shift_post
+}
diff --git a/R/shift_plot.R b/R/shift_plot.R
new file mode 100644
index 0000000..bcb2a5f
--- /dev/null
+++ b/R/shift_plot.R
@@ -0,0 +1,19 @@
+#' Plot a shift function
+#'
+#' Plot the shift function for one particular group.
+#' @param fit An object of class anova_bnp_model.
+#' @param group The group id.
+#' @return A ggplot2 plot.
+#' @importFrom dplyr filter
+#' @importFrom ggplot2 aes geom_line ggplot
+#' @importFrom rlang := .data
+#' @export
+shift_plot <- function(fit, group) {
+  out <-
+    fit |>
+    shift_post() |>
+    dplyr::filter(.data$group == group) |>
+    ggplot2::ggplot(ggplot2::aes(x = y, y = shift)) +
+    ggplot2::geom_line()
+  return(out)
+}
diff --git a/extra/dberpo.R b/extra/dberpo.R
new file mode 100644
index 0000000..d9b1a5a
--- /dev/null
+++ b/extra/dberpo.R
@@ -0,0 +1,16 @@
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+  ngrid <- length(x)
+  pmf <- rep(0, ngrid)
+  for (i in 1:ngrid) {
+    if (x[i] == 0) {
+      pmf[i] = dpois(x = x[i], lambda = lambda)
+    } else {
+      pmf[i] = (
+        (0 + phi) * dpois(x = x[i] - 1, lambda = lambda) +
+        (1 - phi) * dpois(x = x[i] - 0, lambda = lambda)
+      )
+    }
+  }
+  return(pmf)
+}
diff --git a/extra/demo_berpoi.R b/extra/demo_berpoi.R
index bce0aa4..fc3deca 100644
--- a/extra/demo_berpoi.R
+++ b/extra/demo_berpoi.R
@@ -1,9 +1,54 @@
 library(ANOVABNPTestR)
+library(dplyr)
+library(ggplot2)
+library(tidyr)
 
-N <- 1000
-x <- as.integer(1 + (runif(2 * N) < 0.5))
-Xmat <- matrix(x, nrow = N, ncol = 2)
-yvec <- rpois(N, 1) + (Xmat[, 1] == 2) * (Xmat[, 2] == 2)
-fit <- anova_bnp_berpoi(yvec, Xmat)
-hypothesis_post_simple(fit)
-hypothesis_post_interaction(fit)
+# Probability mass function Bernoulli Poisson distribution
+dberpo <- function(x, phi, lambda) {
+  phi * dpois(x - 1, lambda) + (1 - phi) * dpois(x, lambda)
+}
+
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo <- function(n, phi, lambda) {
+  rpois(n, lambda) + rbinom(n, size = 1, prob = phi)
+}
+
+n <- 10000
+phi <- 0.8
+lambda <- 0.5
+yvec <- rberpo(n, phi, lambda)
+Xmat <- matrix(1L, n, 1)
+hist(yvec)
+fit <-
+  anova_bnp_berpoi(
+    yvec,
+    Xmat,
+    iter = 2000L,
+    warmup = 1000L,
+    rho = 1.0,
+    a = 0.1,
+    b = 0.1,
+    a1 = mean(yvec),
+    b1 = 1.0,
+    alpha0 = 1.0,
+    beta0 = 1.0,
+  )
+
+fit$f_post |>
+  dplyr::rename(fh = f) |>
+  dplyr::mutate(f0 = dberpo(y, phi, lambda)) |>
+  tidyr::pivot_longer(cols = fh:f0) |>
+  ggplot2::ggplot(
+    ggplot2::aes(x = y, color = name, y = value)
+  ) +
+  ggplot2::geom_line()
+
+
+
+# N <- 1000
+# x <- as.integer(1 + (runif(2 * N) < 0.5))
+# Xmat <- matrix(x, nrow = N, ncol = 2)
+# yvec <- rpois(N, 1) + (Xmat[, 1] == 2) * (Xmat[, 2] == 2)
+# fit <- anova_bnp_berpoi(yvec, Xmat)
+# hypothesis_post_simple(fit)
+# hypothesis_post_interaction(fit)
diff --git a/extra/demo_normal.R b/extra/demo_normal.R
new file mode 100644
index 0000000..7e37f0b
--- /dev/null
+++ b/extra/demo_normal.R
@@ -0,0 +1,11 @@
+N <- 1000
+x <- as.integer(1 + (runif(2 * N) < 0.5))
+Xmat <- matrix(x, nrow = N, ncol = 2)
+yvec <- rnorm(N)
+fit <- anova_bnp_normal(yvec, Xmat)
+
+tbl_shift <- shift_post(fit)
+hypothesis_post_simple(fit)
+hypothesis_post_interaction(fit)
+
+shift_plot(fit, group = 2)
diff --git a/extra/mc_count_data.R b/extra/mc_count_data.R
new file mode 100644
index 0000000..16c8b44
--- /dev/null
+++ b/extra/mc_count_data.R
@@ -0,0 +1,156 @@
+# Monte Carlo - Coun-data
+library(ANOVABNPTestR)
+library(ggplot2)
+library(MASS) # negative binomial regression
+source("extra/rberpo.R")
+source("extra/dberpo.R")
+
+# Design of the experiment
+# Two factor with two levels
+r1 <- 48  # factor 1 level 1
+r2 <- 48 # factor 1 level 2
+r3 <- 48  # factor 2 level 1
+r4 <- 48  # factor 2 level 2
+n <- r1 + r2 + r3 + r4
+
+f1 <- as.factor(c(rep(1, r1), rep(1, r2), rep(2, r3), rep(2, r4)))
+f2 <- as.factor(c(rep(1, r1), rep(2, r2), rep(1, r3), rep(2, r4)))
+Xmat <- cbind(as.integer(f1), as.integer(f2))
+group <-
+  as.factor(
+    (f1 == 1 & f2 == 1) * 1 +
+    (f1 == 1 & f2 == 2) * 2 +
+    (f1 == 2 & f2 == 1) * 3 +
+    (f1 == 2 & f2 == 2) * 4
+  )
+
+# Simulation of the data
+comp <- seq(1, 2)
+nm <- 100
+alpha <- 0.05
+count_simple_f1 <- 0
+count_simple_f2 <- 0
+count_inter <- 0
+f_sum <- 0
+count_f1_pois <- 0
+count_f2_pois <- 0
+count_f1f2_pois <- 0
+count_f1_nb <- 0
+count_f2_nb <- 0
+count_f1f2_nb <- 0
+
+
+for( k in 1: nm) {
+  # Experiment (1,1)_1
+  y1 <- rberpo(n = r1, phi = 0.5, lambda = 0.5)
+
+  # Experiment (1,2)_2
+  y2 <- rberpo(n = r1, phi = 0.5, lambda = 1.0)
+
+  # Experiment (2,1)_3
+  y3 <- rberpo(n = r1, phi = 0.5, lambda = 1.5)
+
+  # Experiment (2,2)_4
+  y4 <- rberpo(n = r1, phi = 0.5, lambda = 2.0)
+
+  # Fitting our model
+  Y <- as.integer(c(y1, y2, y3, y4))
+  data <- data.frame(Y, f1, f2, group)
+  lb <- as.integer(0)
+  ub <- as.integer(25)
+  a1 <- as.numeric(round(mean(Y),0))  # solo acepta valores enteros, debería aceptar Float 64
+  b1 <- as.numeric(1.0) # solo acepta valores enteros, debería aceptar Float 64
+
+  fit <-
+    anova_bnp_berpoi(
+      Y,
+      Xmat,
+      iter = 10000L,
+      warmup = 2000L,
+      rho = 1.0,
+      a = 0.1,
+      b = 0.1,
+      a1 = a1,
+      b1 = b1,
+      alpha0 = 1.0,
+      beta0 = 1.0,
+      lb = lb,
+      ub = ub
+    )
+  f_pred <- fit$f_post
+  simple <- hypothesis_post_simple(fit)
+  inter <- hypothesis_post_interaction(fit)
+  count_simple_f1 <- count_simple_f1 + (simple[1,2] > 0.5) * 1
+  count_simple_f2 <- count_simple_f2 + (simple[2,2] > 0.5) * 1
+  count_inter <- count_inter + (inter[1,3] > 0.5) * 1
+  f_sum <- f_sum + f_pred$f
+
+  # classical models
+  # Poisson regression
+  fit_glm1 <- summary(glm(Y ~ f1 * f2, data = data, family = poisson()))$coeff
+  count_f1_pois <- count_f1_pois + (fit_glm1[2,4] <= alpha) * 1
+  count_f2_pois <- count_f2_pois + (fit_glm1[3,4] <= alpha) * 1
+  count_f1f2_pois <- count_f1f2_pois + (fit_glm1[4,4] <= alpha) * 1
+
+  # Negative binomial regression
+  fit_glm2 <- summary(glm.nb(Y ~ f1 + f2 + f1 * f2, data = data))$coeff
+  count_f1_nb <- count_f1_nb + (fit_glm2[2, 4] <= alpha) * 1
+  count_f2_nb <- count_f2_nb + (fit_glm2[3, 4] <= alpha) * 1
+  count_f1f2_nb <- count_f1f2_nb + (fit_glm2[4, 4] <= alpha) * 1
+  print(k)
+}
+
+f_mean <- f_sum / nm
+pred_mc <- data.frame(group = f_pred$group, y = f_pred$y, fy = f_mean)
+#pdf(file=paste0("pred_size_",r1,".pdf"))
+#plot(pred_mc$y[pred_mc$group == 1], pred_mc$fy[pred_mc$group==1], pch = 20, col=1, ylim=c(0,0.2), xlab = "y", ylab="pmf")
+#lines(pred_mc$y[pred_mc$group == 1], pred_mc$fy[pred_mc$group==1], lty=1)
+#points(pred_mc$y[pred_mc$group == 2], pred_mc$fy[pred_mc$group==2], pch = 20, col=2)
+#lines(pred_mc$y[pred_mc$group == 2], pred_mc$fy[pred_mc$group==2], lty=1, col=2)
+#points(pred_mc$y[pred_mc$group == 3], pred_mc$fy[pred_mc$group==3], pch = 20, col=3)
+#lines(pred_mc$y[pred_mc$group == 3], pred_mc$fy[pred_mc$group==3], lty=1, col=3)
+#points(pred_mc$y[pred_mc$group == 4], pred_mc$fy[pred_mc$group==4], pch = 20, col=4)
+#lines(pred_mc$y[pred_mc$group == 4], pred_mc$fy[pred_mc$group==4], lty=1, col=4)
+#legend(17,0.2, legend=c("E = (1,1)", "E = (1,2)", "E = (2,1)", "E = (2,2)"), lty = rep(1,4), col=c(1,2,3,4))
+#dev.off()
+
+names <- c("Q1", "Q2", "Q1Q2")
+freq_bnp <- c(count_simple_f1, count_simple_f2, count_inter)
+freq_pois <- c(count_f1_pois, count_f2_pois, count_f1f2_pois)
+freq_nb <- c(count_f1_nb, count_f2_nb, count_f1f2_nb)
+
+
+mc_result <- data.frame(names, freq_bnp, freq_pois, freq_nb)
+
+write.table(x = mc_result, file=paste0("mc_size_", r1, ".csv"), row.names = FALSE,  sep=",")
+
+# True pmf
+grid <- seq(0, max(Y) + 3, 1)
+pmf1 <- dberpo(x = grid, phi = 0.5, lambda = 0.5)
+pmf2 <- dberpo(x = grid, phi = 0.5, lambda = 1.0)
+pmf3 <- dberpo(x = grid, phi = 0.5, lambda = 1.5)
+pmf4 <- dberpo(x = grid, phi = 0.5, lambda = 2.0)
+
+
+pdf(file=paste0("pred_size_",r1,".pdf"))
+plot(grid, pmf1, pch = 20, col = 1, ylim = c(0,0.8), xlab = "y", ylab="pmf")
+lines(grid, pmf1, lty=1, col=1, lwd = 1.5)
+points(grid, pmf2, pch = 20, col = 2)
+lines(grid, pmf2, lty=1, col=2, lwd = 1.5)
+points(grid, pmf3, pch = 20, col = 3)
+lines(grid, pmf3, lty=1, col=3, lwd = 1.5)
+points(grid, pmf4, pch = 20, col = 4)
+lines(grid, pmf4, lty=1, col=4, lwd = 1.5)
+points(pred_mc$y[pred_mc$group == 1], pred_mc$fy[pred_mc$group==1], pch = 20, col=1)
+lines(pred_mc$y[pred_mc$group == 1], pred_mc$fy[pred_mc$group==1], lty=2, col=1, lwd=2)
+points(pred_mc$y[pred_mc$group == 2], pred_mc$fy[pred_mc$group==2], pch = 20, col=2)
+lines(pred_mc$y[pred_mc$group == 2], pred_mc$fy[pred_mc$group==2], lty=2, col=2, lwd=2)
+points(pred_mc$y[pred_mc$group == 3], pred_mc$fy[pred_mc$group==3], pch = 20, col=3)
+lines(pred_mc$y[pred_mc$group == 3], pred_mc$fy[pred_mc$group==3], lty=2, col=3, lwd=2)
+points(pred_mc$y[pred_mc$group == 4], pred_mc$fy[pred_mc$group==4], pch = 20, col=4)
+lines(pred_mc$y[pred_mc$group == 4], pred_mc$fy[pred_mc$group==4], lty=2, col=4, lwd=2)
+legend(6,0.7, legend=c("E = (1,1) --> REF", "E = (1,2) --> Q2", "E = (2,1) --> Q1", "E = (2,2) --> Q1Q2"), lty = rep(2,4), lwd=rep(2,4), col=c(1,2,3,4))
+text(7.0,0.19, labels=c(paste0("Replicates = ",r1)))
+dev.off()
+
+
diff --git a/extra/pkgdown.R b/extra/pkgdown.R
new file mode 100644
index 0000000..43f75dc
--- /dev/null
+++ b/extra/pkgdown.R
@@ -0,0 +1,13 @@
+# Deploy locally
+devtools::document()
+pkgdown::build_site()
+
+# # Delete to gh-pages branch
+# shell("git branch -d gh-pages")
+# shell("git push origin --delete gh-pages")
+
+# Deploy to gh-pages (locally)
+pkgdown::deploy_to_branch()
+
+# Deploy to gh-pages (remotely)
+# TODO
diff --git a/extra/rberpo.R b/extra/rberpo.R
new file mode 100644
index 0000000..2f164c4
--- /dev/null
+++ b/extra/rberpo.R
@@ -0,0 +1,16 @@
+# Generates random draws from a Bernoulli Poisson distribution
+rberpo = function(n, phi, lambda) {
+  y <- rep(0, n)
+  index <- 1:2
+  weights <- c(phi, 1 - phi)
+  for (i in 1:n) {
+    j <- sample(x = index, size = 1, prob = weights)
+    if (j == 1) {
+      y[i] <- rpois(1, lambda = lambda) + 1
+    } else {
+      y[i] <- rpois(n = 1, lambda = lambda)
+    }
+  }
+  return(y)
+}
+
diff --git a/man/shift_plot.Rd b/man/shift_plot.Rd
new file mode 100644
index 0000000..7a8efeb
--- /dev/null
+++ b/man/shift_plot.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/shift_plot.R
+\name{shift_plot}
+\alias{shift_plot}
+\title{Plot a shift function}
+\usage{
+shift_plot(fit, group)
+}
+\arguments{
+\item{fit}{An object of class anova_bnp_model.}
+
+\item{group}{The group id.}
+}
+\value{
+A ggplot2 plot.
+}
+\description{
+Plot the shift function for one particular group.
+}
diff --git a/man/shift_post.Rd b/man/shift_post.Rd
new file mode 100644
index 0000000..ae2e100
--- /dev/null
+++ b/man/shift_post.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/f_post.R
+\name{shift_post}
+\alias{shift_post}
+\title{Return the shift functions}
+\usage{
+shift_post(fit)
+}
+\arguments{
+\item{fit}{An object of class anova_bnp_model.}
+}
+\value{
+A tibble with the following variables:
+\itemize{
+\item group the group id.
+\item y the value of the response variable.
+\item the shift function evaluated at y, given the group.
+}
+}
+\description{
+Return the shift functions
+}