## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(dpi = 72, collapse = TRUE, comment = "#>")
# Penalized fits need glmnet (frequentist) / shrinkem (Bayesian), both declared
# in Suggests. If a backend is unavailable - e.g. a CRAN check run with
# _R_CHECK_FORCE_SUGGESTS_=false - the relevant chunks render code-only instead
# of erroring the vignette build.
has_glmnet   <- requireNamespace("glmnet",   quietly = TRUE)
has_shrinkem <- requireNamespace("shrinkem", quietly = TRUE)

## ----load---------------------------------------------------------------------
library(remverse)

## ----data---------------------------------------------------------------------
data(randomREH3)
data(info3)

reh <- remify(edgelist = randomREH3, model = "tie", directed = TRUE)

## ----kitchen-sink-------------------------------------------------------------
effects_rich <- ~ inertia(scaling = "std") +
                   reciprocity(scaling = "std") +
                   indegreeSender(scaling = "std") +
                   indegreeReceiver(scaling = "std") +
                   outdegreeSender(scaling = "std") +
                   outdegreeReceiver(scaling = "std") +
                   totaldegreeDyad(scaling = "std") +
                   isp(scaling = "std") +
                   osp(scaling = "std") +
                   itp(scaling = "std") +
                   otp(scaling = "std") +
                   send("age", attr_actors = info3, scaling = "std") +
                   receive("age", attr_actors = info3, scaling = "std") +
                   difference("age", attr_actors = info3, scaling = "std")

stats_rich <- remstats(reh = reh, tie_effects = effects_rich, first = 500,
                       memory = "decay", memory_value = 2000)
dimnames(stats_rich)[[3]]

## ----mle-rich-----------------------------------------------------------------
fit_mle <- remstimate(reh = reh, stats = stats_rich)
summary(fit_mle)

## ----shrinkem, message=FALSE, eval=has_shrinkem-------------------------------
fit_bayes <- rempenalty(reh = reh, stats = stats_rich, approach = "Bayesian")
summary(fit_bayes)
round(coef(fit_bayes), 3)

## ----glmnet-default, message=FALSE, eval=has_glmnet---------------------------
fit_glmnet <- rempenalty(reh = reh, stats = stats_rich, approach = "frequentist")
summary(fit_glmnet)
coef(fit_glmnet)

## ----compare-coefs, eval=has_glmnet-------------------------------------------
coefs_mle    <- coef(fit_mle)
coefs_bayes <- coef(fit_bayes)
coefs_glmnet <- coef(fit_glmnet)

shared <- intersect(intersect(names(coefs_mle), names(coefs_glmnet)), names(coefs_bayes))

comparison <- data.frame(
  statistic = shared,
  MLE       = round(coefs_mle[shared], 3),
  BAYES     = round(coefs_bayes[shared], 3),
  GLMNET    = round(coefs_glmnet[shared], 3),
  row.names = NULL
)
comparison

## ----alpha, message=FALSE, eval=has_glmnet------------------------------------
# Pure lasso (default): maximum sparsity
fit_lasso <- rempenalty(reh, stats_rich, approach = "frequentist", alpha = 1)

# Pure ridge: shrinkage without variable selection
fit_ridge <- rempenalty(reh, stats_rich, approach = "frequentist", alpha = 0)

cat("Lasso non-zero:", sum(coef(fit_lasso) != 0), "of", length(coef(fit_lasso)), "\n")
cat("Ridge non-zero:", sum(coef(fit_ridge) != 0), "of", length(coef(fit_ridge)), "\n")

## ----lambda, message=FALSE, eval=has_glmnet-----------------------------------
fit_min <- rempenalty(reh, stats_rich, approach = "frequentist",
                      lambda_select = "min")

cat("lambda.1se non-zero:", sum(coef(fit_glmnet) != 0), "\n")
cat("lambda.min non-zero:", sum(coef(fit_min)    != 0), "\n")

## ----diagnostics, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white"), eval=has_glmnet----
diag_mle    <- diagnostics(fit_mle,    reh, stats_rich)
diag_bayes <- diagnostics(fit_bayes,   reh, stats_rich)
diag_glmnet <- diagnostics(fit_glmnet, reh, stats_rich)

cat("MLE recall:    ", round(diag_mle$recall$summary$mean_rel_rank, 3),    "\n")
cat("BAYES recall: ", round(diag_bayes$recall$summary$mean_rel_rank, 3), "\n")
cat("GLMNET recall: ", round(diag_glmnet$recall$summary$mean_rel_rank, 3), "\n")

plot(diag_mle)
plot(diag_bayes)
plot(diag_glmnet)

## ----show-comparison, eval=has_glmnet-----------------------------------------
comparison