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

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

## ----data---------------------------------------------------------------------
data(randomREH3)   # event history
data(info3)        # actor attributes

head(randomREH3)
head(info3)

## ----tie-basic----------------------------------------------------------------
# Step 1: process the event history
reh <- remify(edgelist = randomREH3, model = "tie", directed = TRUE)
summary(reh)

# Step 2: define effects and compute statistics
effects <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
             outdegreeSender(scaling = "std") + isp(scaling = "std")
stats <- remstats(reh = reh, tie_effects = effects)
stats

# Step 3: estimate the model
fit <- remstimate(reh = reh, stats = stats)
summary(fit)

## ----tie-basic-diag, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white")----
diag <- diagnostics(object = fit, reh = reh, stats = stats)
diag

plot(x = fit, reh = reh, diagnostics = diag)

## ----tie-ordinal--------------------------------------------------------------
reh_ord <- remify(edgelist = randomREH3, model = "tie", ordinal = TRUE)

stats_ord <- remstats(reh = reh_ord,
                      tie_effects = ~ inertia(scaling = "std") +
                                      reciprocity(scaling = "std"))
fit_ord <- remstimate(reh = reh_ord, stats = stats_ord)
summary(fit_ord)

## ----tie-active---------------------------------------------------------------
reh_active <- remify(edgelist = randomREH3, model = "tie", riskset = "active")

cat("Full risk set:  ", reh$D, "dyads\n")
cat("Active risk set:", reh_active$activeD, "dyads\n")

stats_active <- remstats(reh = reh_active,
                         tie_effects = ~ inertia(scaling = "std") +
                                         reciprocity(scaling = "std"))
fit_active <- remstimate(reh = reh_active, stats = stats_active)
summary(fit_active)

## ----tie-active-saturated-----------------------------------------------------
reh_sat <- remify(edgelist = randomREH3, model = "tie",
                  riskset = "active_saturated")
cat("Active-saturated risk set:", reh_sat$activeD, "dyads\n")

## ----tie-manual---------------------------------------------------------------
# Define a manual risk set: a subset of directed dyads among actors 1-5
my_riskset <- data.frame(
  actor1 = c(1, 2, 3, 4, 5),
  actor2 = c(2, 3, 4, 5, 1)
)

reh_manual <- remify(
  edgelist       = randomREH3,
  model          = "tie",
  riskset        = "manual",
  manual_riskset = my_riskset
)
cat("Manual risk set:", reh_manual$activeD, "dyads\n")

## ----tie-typed----------------------------------------------------------------
reh_typed <- remify(
  edgelist   = randomREH3,
  model      = "tie",
  event_type = "setting"
)

## ----typed-ignore-------------------------------------------------------------
stats_ign <- remstats(
  reh = reh_typed,
  tie_effects = ~ inertia(consider_type = "ignore")
)
dimnames(stats_ign)[[3]]

## ----typed-separate-----------------------------------------------------------
stats_sep <- remstats(
  reh = reh_typed,
  tie_effects = ~ inertia(consider_type = "separate")
)
dimnames(stats_sep)[[3]]

## ----typed-interact-----------------------------------------------------------
stats_int <- remstats(
  reh = reh_typed,
  tie_effects = ~ inertia(consider_type = "interact")
)
dimnames(stats_int)[[3]]

## ----typed-fit----------------------------------------------------------------
stats_typed <- remstats(
  reh = reh_typed,
  tie_effects = ~ inertia(consider_type = "separate") +
                   reciprocity(consider_type = "separate")
)

fit_typed <- remstimate(reh = reh_typed, stats = stats_typed)
summary(fit_typed)

## ----tie-typed-ext------------------------------------------------------------
reh_typed_ext <- remify(
  edgelist               = randomREH3,
  model                  = "tie",
  event_type             = "setting",
  extend_riskset_by_type = TRUE
)
cat("Dyad-only risk set:      ", reh_typed$D, "\n")
cat("Type-extended risk set:  ", reh_typed_ext$D, "\n")

## ----exo-actor----------------------------------------------------------------
stats_exo <- remstats(
  reh = reh,
  tie_effects = ~ inertia(scaling = "std") +
                   reciprocity(scaling = "std") +
                   send("age", attr_actors = info3, scaling = "std") +
                   receive("age", attr_actors = info3, scaling = "std") +
                   difference("age", attr_actors = info3, scaling = "std")
)

fit_exo <- remstimate(reh = reh, stats = stats_exo)
summary(fit_exo)

## ----exo-interaction----------------------------------------------------------
stats_ix <- remstats(
  reh = reh,
  tie_effects = ~ inertia(scaling = "std") +
                   send("age", attr_actors = info3, scaling = "std") +
                   inertia(scaling = "std"):send("age", attr_actors = info3, scaling = "std")
)

fit_ix <- remstimate(reh = reh, stats = stats_ix)
summary(fit_ix)

## ----memory-------------------------------------------------------------------
# Full memory (default): entire event history
stats_full <- remstats(reh = reh, tie_effects = ~ inertia(), memory = "full")

# Window memory: only events within the last 2000 time units
stats_win <- remstats(reh = reh, tie_effects = ~ inertia(),
                      memory = "window", memory_value = 2000)

# Decay memory: exponential decay with half-life 1000
stats_dec <- remstats(reh = reh, tie_effects = ~ inertia(),
                      memory = "decay", memory_value = 1000)

# Compare the inertia statistic at the last time point for the first dyad
M <- dim(stats_full)[1]
cat("Full memory:   ", stats_full[M, 1, "inertia"], "\n")
cat("Window (2000): ", stats_win[M, 1, "inertia"], "\n")
cat("Decay (1000):  ", stats_dec[M, 1, "inertia"], "\n")

## ----sampling-----------------------------------------------------------------
stats_sampled <- remstats(
  reh = reh,
  tie_effects = ~ inertia(scaling = "std") + reciprocity(scaling = "std"),
  sampling = TRUE,
  samp_num = 5,
  seed = 42
)

fit_sampled <- remstimate(reh = reh, stats = stats_sampled)
summary(fit_sampled)

## ----assess-helper------------------------------------------------------------
assess <- function(formula, memory_value = 2000) {
  s <- remstats(reh, tie_effects = formula, first = 200,
                memory = "decay", memory_value = memory_value)
  m <- remstimate(reh, s)
  d <- diagnostics(m, reh, s)
  data.frame(
    npar          = length(coef(m)),
    AIC           = round(m$AIC, 1),
    BIC           = round(m$BIC, 1),
    median_recall = round(d$recall$summary$median_rel_rank, 3)
  )
}

## ----assess-models, message=FALSE, warning=FALSE------------------------------
f_empty <- ~ 1
f_noage <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
             outdegreeSender(scaling = "std") + isp(scaling = "std")
f_wage  <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
             outdegreeSender(scaling = "std") +
             send("age", attr_actors = info3, scaling = "std") +
             isp(scaling = "std")
f_rich  <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
             outdegreeSender(scaling = "std") + indegreeReceiver(scaling = "std") +
             send("age", attr_actors = info3, scaling = "std") +
             receive("age", attr_actors = info3, scaling = "std") +
             isp(scaling = "std") + osp(scaling = "std")

comparison <- rbind(
  empty  = assess(f_empty),
  no_age = assess(f_noage),
  wage   = assess(f_wage),
  rich   = assess(f_rich)
)
comparison

## ----assess-summary-----------------------------------------------------------
stats_wage <- remstats(reh, tie_effects = f_wage, first = 200,
                       memory = "decay", memory_value = 2000)
summary(remstimate(reh, stats_wage))

## ----assess-memory, message=FALSE, warning=FALSE, out.width="60%", fig.align="center", dev=c("jpeg"), dev.args = list(bg = "white")----
half_lives <- c(500, 1000, 2000, 3000, 4000)
mem_path <- do.call(rbind, lapply(half_lives, function(h)
  cbind(half_life = h, assess(f_noage, memory_value = h))))
mem_path[, c("half_life", "BIC", "median_recall")]

plot(mem_path$half_life, mem_path$BIC, type = "b", pch = 19,
     xlab = "decay half-life", ylab = "BIC", main = "BIC vs. memory half-life")

## ----actor-basic--------------------------------------------------------------
reh_ao <- remify(edgelist = randomREH3, model = "actor", directed = TRUE)

rate_effects   <- ~ outdegreeSender(scaling = "std") + indegreeSender(scaling = "std")
choice_effects <- ~ inertia(scaling = "std") + reciprocity(scaling = "std")

stats_ao <- remstats(
  reh = reh_ao,
  sender_effects   = rate_effects,
  receiver_effects = choice_effects
)
stats_ao

fit_ao <- remstimate(reh = reh_ao, stats = stats_ao)
summary(fit_ao)

## ----actor-exo----------------------------------------------------------------
rate_effects_exo   <- ~ outdegreeSender(scaling = "std") +
                        send("age", attr_actors = info3, scaling = "std")
choice_effects_exo <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") +
                        receive("age", attr_actors = info3, scaling = "std")

stats_ao_exo <- remstats(
  reh = reh_ao,
  sender_effects   = rate_effects_exo,
  receiver_effects = choice_effects_exo
)

fit_ao_exo <- remstimate(reh = reh_ao, stats = stats_ao_exo)
summary(fit_ao_exo)

## ----actor-diag, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white")----
diag_ao <- diagnostics(object = fit_ao, reh = reh_ao, stats = stats_ao)
plot(x = fit_ao, reh = reh_ao, diagnostics = diag_ao)

## ----hmc, message=FALSE-------------------------------------------------------
fit_hmc <- remstimate(
  reh      = reh,
  stats    = stats,
  approach = "Bayesian",
  bayes    = list(
    nsim    = 300L,
    nchains = 2L,
    burnin  = 300L,
    L       = 100L,
    epsilon = 0.001,
    thin    = 2L
  ),
  seed = 42
)

summary(fit_hmc)

## ----hmc-diag, out.width="50%", dev=c("jpeg"), dev.args = list(bg = "white")----
diag_hmc <- diagnostics(object = fit_hmc, reh = reh, stats = stats)
plot(x = fit_hmc, reh = reh, diagnostics = diag_hmc)