## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(dpi = 72, collapse = TRUE, comment = "#>")
# Mixture fits need flexmix (declared in Suggests). If it is unavailable - e.g.
# a CRAN check run with _R_CHECK_FORCE_SUGGESTS_=false - the modelling chunks
# render code-only instead of erroring the vignette build.
has_flexmix <- requireNamespace("flexmix", quietly = TRUE)

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

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

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

stats <- remstats(
  reh = reh,
  tie_effects = ~ inertia(scaling = "std") + reciprocity(scaling = "std"),
  first = 100, memory = "decay", memory_value = 2000
)

## ----fit-k2, message=FALSE, eval=has_flexmix----------------------------------
fit_k2 <- dlcrem(reh = reh, stats = stats, k = 2)
fit_k2

## ----coefs, eval=has_flexmix--------------------------------------------------
coef(fit_k2)

## ----select-k, message=FALSE, eval=has_flexmix--------------------------------
fit_multi <- remixture(
  reh    = reh,
  stats  = stats,
  random = ~ (1 | dyad),
  k      = 1:4,
  nrep   = 5
)

bic_table(fit_multi)

## ----plot-k, out.width="60%", fig.align="center", dev=c("jpeg"), dev.args = list(bg = "white"), eval=has_flexmix----
plot(fit_multi)

## ----best-model, eval=has_flexmix---------------------------------------------
bics   <- bic_table(fit_multi)
best_k <- bics$k[which.min(bics$BIC)]
cat("Best K:", best_k, "\n")

fit_best <- fit_multi[[paste0("k", best_k)]]
fit_best

## ----diagnostics, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white"), eval=has_flexmix----
diag_k2 <- diagnostics(fit_k2, reh, stats)
diag_k2

plot(diag_k2)

## ----compare-loglik, message=FALSE, eval=has_flexmix--------------------------
fit_mle <- remstimate(reh = reh, stats = stats)

cat("MLE logLik:    ", round(as.numeric(logLik(fit_mle)), 2), "\n")
cat("MIXREM logLik: ", round(as.numeric(logLik(fit_k2)),  2), "\n")