## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## -----------------------------------------------------------------------------
library(gipsDA)

## -----------------------------------------------------------------------------
set.seed(42)

train_id <- unlist(
  lapply(split(seq_len(nrow(iris)), iris$Species), sample, size = 35),
  use.names = FALSE
)

train <- iris[train_id, ]
test <- iris[-train_id, ]

table(train$Species)
table(test$Species)

## -----------------------------------------------------------------------------
str(train)

## -----------------------------------------------------------------------------
is.factor(train$Species)

## -----------------------------------------------------------------------------
vapply(train[, 1:4], is.numeric, logical(1))

## -----------------------------------------------------------------------------
x_train <- train[, 1:4]
x_test <- test[, 1:4]

train_center <- vapply(x_train, mean, numeric(1))
train_scale <- vapply(x_train, sd, numeric(1))

train_scaled <- train
test_scaled <- test

train_scaled[, 1:4] <- scale(
  x_train,
  center = train_center,
  scale = train_scale
)

test_scaled[, 1:4] <- scale(
  x_test,
  center = train_center,
  scale = train_scale
)

## -----------------------------------------------------------------------------
fit_scaled <- gipslda(Species ~ ., data = train_scaled)
pred_scaled <- predict(fit_scaled, test_scaled)

mean(pred_scaled$class == test_scaled$Species)

## -----------------------------------------------------------------------------
lda_fit <- gipslda(Species ~ ., data = train)
qda_fit <- gipsqda(Species ~ ., data = train)
joint_qda_fit <- gipsmultqda(Species ~ ., data = train)

## -----------------------------------------------------------------------------
lda_pred <- predict(lda_fit, test)
qda_pred <- predict(qda_fit, test)
joint_qda_pred <- predict(joint_qda_fit, test)

c(
  gipslda = mean(lda_pred$class == test$Species),
  gipsqda = mean(qda_pred$class == test$Species),
  gipsmultqda = mean(joint_qda_pred$class == test$Species)
)

## ----model-hierarchy, echo = FALSE, fig.cap = "The diagram illustrates the hierarchical relationships between the models.", out.width = "95%"----
knitr::include_graphics("figures/models_hierarchy.png")

## -----------------------------------------------------------------------------
lda_map <- gipslda(
  Species ~ .,
  data = train,
  MAP = TRUE
)

lda_map

## -----------------------------------------------------------------------------
lda_avg <- gipslda(
  Species ~ .,
  data = train,
  MAP = FALSE
)

lda_avg

## -----------------------------------------------------------------------------
fit_bf <- gipsqda(
  Species ~ .,
  data = train,
  optimizer = "BF"
)

fit_bf

## ----eval = FALSE-------------------------------------------------------------
# fit_mh <- gipsqda(
#   Species ~ .,
#   data = train,
#   optimizer = "MH",
#   max_iter = 1000
# )

## ----eval = FALSE-------------------------------------------------------------
# fit_mh_100 <- gipsqda(
#   Species ~ .,
#   data = train,
#   optimizer = "MH",
#   max_iter = 100
# )
# 
# fit_mh_1000 <- gipsqda(
#   Species ~ .,
#   data = train,
#   optimizer = "MH",
#   max_iter = 1000
# )

## -----------------------------------------------------------------------------
lda_fit$prior

## -----------------------------------------------------------------------------
equal_prior <- rep(1 / length(levels(train$Species)), length(levels(train$Species)))
names(equal_prior) <- levels(train$Species)

lda_equal_prior <- gipslda(
  Species ~ .,
  data = train,
  prior = equal_prior
)

lda_equal_prior$prior

## -----------------------------------------------------------------------------
sum(equal_prior)

## -----------------------------------------------------------------------------
lda_classic <- gipslda(
  Species ~ .,
  data = train,
  weighted_avg = FALSE
)

lda_weighted <- gipslda(
  Species ~ .,
  data = train,
  weighted_avg = TRUE
)

## -----------------------------------------------------------------------------
pred_classic <- predict(lda_classic, test)
pred_weighted <- predict(lda_weighted, test)

c(
  classic = mean(pred_classic$class == test$Species),
  weighted = mean(pred_weighted$class == test$Species)
)

## -----------------------------------------------------------------------------
fit_formula <- gipslda(
  Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
  data = train
)

fit_formula

## -----------------------------------------------------------------------------
fit_formula_short <- gipslda(
  Species ~ .,
  data = train
)

fit_formula_short

## -----------------------------------------------------------------------------
fit_subset <- gipslda(
  Species ~ .,
  data = iris,
  subset = Species != "setosa"
)

fit_subset

## -----------------------------------------------------------------------------
x <- as.matrix(iris[, 1:4])
grouping <- iris$Species

fit_matrix_lda <- gipslda(x, grouping)
fit_matrix_qda <- gipsqda(x, grouping)
fit_matrix_joint <- gipsmultqda(x, grouping)

## -----------------------------------------------------------------------------
predict(fit_matrix_lda, x[1:5, ])$class
predict(fit_matrix_qda, x[1:5, ])$class
predict(fit_matrix_joint, x[1:5, ])$class

## -----------------------------------------------------------------------------
x_df <- iris[, 1:4]
y <- iris$Species

fit_df_lda <- gipslda(x_df, y)
fit_df_qda <- gipsqda(x_df, y)
fit_df_joint <- gipsmultqda(x_df, y)

## -----------------------------------------------------------------------------
predict(fit_df_lda, x_df[1:5, ])$class
predict(fit_df_qda, x_df[1:5, ])$class
predict(fit_df_joint, x_df[1:5, ])$class

## -----------------------------------------------------------------------------
pred <- predict(lda_fit, test)

names(pred)

## -----------------------------------------------------------------------------
head(pred$class)
head(pred$posterior)

## -----------------------------------------------------------------------------
names(qda_pred)
names(joint_qda_pred)

## -----------------------------------------------------------------------------
pred_plugin <- predict(lda_fit, test, method = "plug-in")
pred_predictive <- predict(lda_fit, test, method = "predictive")
pred_debiased <- predict(lda_fit, test, method = "debiased")

## -----------------------------------------------------------------------------
c(
  plugin = mean(pred_plugin$class == test$Species),
  predictive = mean(pred_predictive$class == test$Species),
  debiased = mean(pred_debiased$class == test$Species)
)

## -----------------------------------------------------------------------------
head(pred_plugin$posterior)
head(pred_predictive$posterior)
head(pred_debiased$posterior)

## -----------------------------------------------------------------------------
qda_loo <- predict(qda_fit, method = "looCV")
joint_qda_loo <- predict(joint_qda_fit, method = "looCV")

c(
  gipsqda_loo_accuracy = mean(qda_loo$class == train$Species),
  gipsmultqda_loo_accuracy = mean(joint_qda_loo$class == train$Species)
)

## -----------------------------------------------------------------------------
print(lda_fit)

## -----------------------------------------------------------------------------
print(qda_fit)

## -----------------------------------------------------------------------------
print(joint_qda_fit)

## -----------------------------------------------------------------------------
summary(lda_fit)
summary(qda_fit)
summary(joint_qda_fit)

## -----------------------------------------------------------------------------
names(lda_fit)
names(qda_fit)
names(joint_qda_fit)

## -----------------------------------------------------------------------------
lda_fit$optimization_info

## -----------------------------------------------------------------------------
qda_fit$optimization_info

## -----------------------------------------------------------------------------
joint_qda_fit$optimization_info

## -----------------------------------------------------------------------------
inspect_model <- function(object) {
  data.frame(
    component = names(object),
    class = vapply(
      object,
      function(x) paste(class(x), collapse = ", "),
      character(1)
    ),
    length = vapply(object, length, integer(1)),
    dim = vapply(
      object,
      function(x) {
        d <- dim(x)
        if (is.null(d)) "" else paste(d, collapse = " x ")
      },
      character(1)
    ),
    row.names = NULL
  )
}

## -----------------------------------------------------------------------------
inspect_model(lda_fit)

## -----------------------------------------------------------------------------
inspect_model(qda_fit)

## -----------------------------------------------------------------------------
inspect_model(joint_qda_fit)

## -----------------------------------------------------------------------------
coef(lda_fit)

## ----fig.width = 6, fig.height = 5, fig.alt = "Plot of the fitted gipslda model in discriminant space."----
plot(lda_fit)

## ----fig.width = 6, fig.height = 5, fig.alt = "Pairs plot of discriminant coordinates for the fitted gipslda model."----
pairs(lda_fit, type = "std")

## -----------------------------------------------------------------------------
fit_lda <- gipslda(Species ~ ., data = train)
fit_qda <- gipsqda(Species ~ ., data = train)
fit_joint <- gipsmultqda(Species ~ ., data = train)

pred_lda <- predict(fit_lda, test)
pred_qda <- predict(fit_qda, test)
pred_joint <- predict(fit_joint, test)

c(
  gipslda = mean(pred_lda$class == test$Species),
  gipsqda = mean(pred_qda$class == test$Species),
  gipsmultqda = mean(pred_joint$class == test$Species)
)

## ----eval = FALSE-------------------------------------------------------------
# fit <- gipsqda(
#   Species ~ .,
#   data = train,
#   optimizer = "MH",
#   max_iter = 1000
# )

## -----------------------------------------------------------------------------
head(pred_lda$posterior)
head(pred_qda$posterior)
head(pred_joint$posterior)