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

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

## ----data---------------------------------------------------------------------
data(edgelist0)
data(edgelist0_actors)

head(edgelist0)
head(edgelist0_actors)

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

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

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

## ----actor-basic--------------------------------------------------------------
# Step 1: process the event history
reh_ao <- remify(edgelist = edgelist0, model = "actor", directed = TRUE, riskset = "active")
summary(reh_ao)

# Step 2: specify statistics for sender and receiver choice models separately
sender_effects <- ~ outdegreeSender() + indegreeSender()
choice_effects <- ~ inertia() + reciprocity()

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

# Step 3: estimate the model
fit_ao <- remstimate(reh = reh_ao, stats = stats_ao, method = "MLE")
summary(fit_ao)

## ----ordinal-tie--------------------------------------------------------------
# Tie model with ordinal timing
reh_ord <- remify(edgelist = edgelist0, model = "tie", ordinal = TRUE)
stats_ord <- remstats(reh = reh_ord, tie_effects = ~ inertia() + reciprocity())
fit_ord <- remstimate(reh = reh_ord, stats = stats_ord, method = "MLE")
summary(fit_ord)

## ----ordinal-actor------------------------------------------------------------
# Actor model with ordinal timing
reh_ao_ord <- remify(edgelist = edgelist0, model = "actor", ordinal = TRUE, riskset = "active")
stats_ao_ord <- remstats(
  reh = reh_ao_ord,
  sender_effects = ~ outdegreeSender(),
  receiver_effects = ~ inertia() + reciprocity()
)
fit_ao_ord <- remstimate(reh = reh_ao_ord, stats = stats_ao_ord, method = "MLE")
summary(fit_ao_ord)

## ----undirected---------------------------------------------------------------
reh_undir <- remify(edgelist = edgelist0, model = "tie", directed = FALSE)
cat("Directed dyads:", reh$D, "\n")
cat("Undirected pairs:", reh_undir$D, "\n")

stats_undir <- remstats(
  reh = reh_undir,
  tie_effects = ~ inertia() + sp()
)

fit_undir <- remstimate(reh = reh_undir, stats = stats_undir, method = "MLE")
summary(fit_undir)

## ----active-------------------------------------------------------------------
reh_active <- remify(edgelist = edgelist0, 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() + reciprocity())
fit_active <- remstimate(reh = reh_active, stats = stats_active, method = "MLE")
summary(fit_active)

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

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

## ----typed-setup--------------------------------------------------------------
reh_typed <- remify(
  edgelist   = edgelist0,
  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-----------------------------------------------------------
reh_typed_ext <- remify(
  edgelist   = edgelist0,
  model      = "tie",
  event_type = "setting",
  extend_riskset_by_type = TRUE
)
stats_int <- remstats(
  reh = reh_typed_ext,
  tie_effects = ~ inertia(consider_type = "interact")
)
dimnames(stats_int)[[3]]

## ----typed-fit----------------------------------------------------------------
stats_typed <- remstats(
  reh = reh_typed_ext,
  tie_effects = ~ inertia(consider_type = "interact") +
                   reciprocity(consider_type = "ignore")
)
fit_typed <- remstimate(reh = reh_typed_ext, stats = stats_typed, method = "MLE")
summary(fit_typed)

## ----typed-actor--------------------------------------------------------------
reh_ao_typed <- remify(
  edgelist   = edgelist0,
  model      = "actor",
  event_type = "setting"
)
stats_ao_typed <- remstats(
  reh = reh_ao_typed,
  sender_effects = ~ outdegreeSender(),
  receiver_effects = ~ inertia(consider_type = "separate") + reciprocity()
)
fit_ao_typed <- remstimate(reh = reh_ao_typed, stats = stats_ao_typed, method = "MLE")
summary(fit_ao_typed)

## ----scaling------------------------------------------------------------------
reh <- remify(edgelist = edgelist0, model = "tie", directed = TRUE)

stats_none <- remstats(reh = reh, tie_effects = ~ inertia(scaling = "none"))
stats_prop <- remstats(reh = reh, tie_effects = ~ inertia(scaling = "prop"))
stats_std  <- remstats(reh = reh, tie_effects = ~ inertia(scaling = "std"))

# Compare values at the last time point for the first dyad
cat("Raw count:      ", stats_none[114, 1, "inertia"], "\n")
cat("Proportional:   ", stats_prop[114, 1, "inertia"], "\n")
cat("Standardized:   ", stats_std[114, 1, "inertia"], "\n")

remst_none <- remstimate(reh, stats = stats_none)
coef(remst_none)
remst_prop <- remstimate(reh, stats = stats_prop)
coef(remst_prop)
remst_std <- remstimate(reh, stats = stats_std)
coef(remst_std)

## ----exo-tie------------------------------------------------------------------
exo_effects <- ~ inertia(scaling = "std") +
                   reciprocity(scaling = "std") +
                   send("job", edgelist0_actors) +
                   same("job", edgelist0_actors)
stats_exo <- remstats(
  reh = reh,
  tie_effects = exo_effects
)
fit_exo <- remstimate(reh = reh, stats = stats_exo, method = "MLE")
summary(fit_exo)

## ----exo-actor----------------------------------------------------------------
rate_effects_exo <- ~ outdegreeSender(scaling = "std") + send("job", edgelist0_actors)
choice_effects_exo <- ~ inertia(scaling = "std") + reciprocity(scaling = "std") + same("job", edgelist0_actors)

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, method = "MLE")
summary(fit_ao_exo)

## ----exo-interaction----------------------------------------------------------
stats_ix <- remstats(
  reh = reh,
  tie_effects = ~ inertia() +
                   send("job", edgelist0_actors) +
                   inertia():send("job", edgelist0_actors)
)

fit_ix <- remstimate(reh = reh, stats = stats_ix, method = "MLE")
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 500 time units
stats_win <- remstats(reh = reh, tie_effects = ~ inertia(),
                      memory = "window", memory_value = 500)

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

# Compare the inertia statistic at the last time point for a single dyad
cat("Full memory:  ", stats_full[900, 2, "inertia"], "\n")
cat("Window (500): ", stats_win[900, 2, "inertia"], "\n")
cat("Decay (200):  ", stats_dec[900, 2, "inertia"], "\n")

## ----sampling-----------------------------------------------------------------
reh_samp <- remify(edgelist = edgelist0, model = "tie", directed = TRUE, event_type = "setting")
stats_sampled <- remstats(
  reh = reh_samp,
  tie_effects = exo_effects,
  sampling = TRUE,
  samp_num = 20,
  seed = 42
)
fit_sampled <- remstimate(reh = reh_samp, stats = stats_sampled, method = "MLE")
summary(fit_sampled)

## ----hmc, message=FALSE-------------------------------------------------------
P <- length(dimnames(stats_exo)[[3]])
my_prior <- list(
  mean = rep(0, P),
  vcov = diag(c(100, rep(1, P - 1)))  # wide prior on intercept, N(0,1) on effects
)

fit_exo_hmc <- remstimate(
  reh    = reh,
  stats  = stats_exo,
  method = "HMC",
  prior  = my_prior,
  burnin = 200,
  nsim   = 500,
  thin   = 2,
  L      = 100,
  seed   = 42
)

summary(fit_exo_hmc)

## ----diag-tie, out.width="80%", fig.width=8, fig.height=4, dev=c("jpeg"), dev.args = list(bg = "white")----
diag <- diagnostics(object = fit_exo, reh = reh, stats = stats_exo)
diag
plot(x = fit_exo, reh = reh, diagnostics = diag)

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

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