#############################################################
# _compare_helpers.R -- cross-family coefficient extractors.
#
# The take-home deliverable of the capstone's payoff section: every
# model family in this course prints its coefficient table in its own
# dialect (summary() matrices, posterior draws, tidy() frames). Each
# row_*() function below speaks one dialect and returns the SAME
# five-column data.frame --
#
#     model | term | estimate | lo | hi
#
# -- so you can rbind() the lot and draw one forest plot. Each needs
# only the fitted object (plus, for row_dcr, the dcr output): point
# them at your own fits and the plot below works unchanged.
#
# THE ONE RULE OF USE: these helpers align *verdicts*, not magnitudes.
# glm and ergm coefficients live on the logit scale, lame's on the
# probit scale (roughly 1.6x smaller for the same effect), and the
# ergm's are conditional on the rest of the graph besides. Read sign
# and interval-vs-zero across rows; never read size across families.
#############################################################

# strip each family's naming decoration so the same covariate lines up
# across models: lame appends _dyad/_node, ergm prepends edgecov./nodecov.
clean_term <- function(x) {
  x <- sub("^(edgecov|nodecov|absdiff)\\.", "", x)
  x <- sub("_(dyad|node|row|col)$", "", x)
  x[x == "averageId_sum"] <- "averageId_actor1"
  x[x == "size_sum"] <- "size_actor1"
  x[x %in% c("(Intercept)", "intercept", "edges")] <- "(baseline)"
  x
}

# glm / lm: point estimate +/- 1.96 * naive SE
row_glm <- function(fit, model = "glm (naive)") {
  s <- summary(fit)$coefficients
  data.frame(model = model, term = clean_term(rownames(s)),
             estimate = s[, 1],
             lo = s[, 1] - 1.96 * s[, 2],
             hi = s[, 1] + 1.96 * s[, 2],
             row.names = NULL)
}

# same glm, dcr's dyadic cluster-robust variance: same estimate,
# different uncertainty calculation
row_dcr <- function(fit, dcr_out, model = "glm + dcr SEs") {
  est <- coef(fit); se <- dcr_out$dcrse
  data.frame(model = model, term = clean_term(names(est)),
             estimate = est,
             lo = est - 1.96 * se, hi = est + 1.96 * se,
             row.names = NULL)
}

# lame / amen fits (srm, ame): posterior median and 95% credible
# interval straight from the BETA draws
row_ame <- function(fit, model = "ame") {
  q <- t(apply(fit$BETA, 2, quantile, c(0.5, 0.025, 0.975)))
  data.frame(model = model, term = clean_term(rownames(q)),
             estimate = q[, 1], lo = q[, 2], hi = q[, 3],
             row.names = NULL)
}

# ergm: MLE +/- 1.96 * SE from the asymptotic covariance
row_ergm <- function(fit, model = "ergm") {
  est <- coef(fit); se <- sqrt(diag(vcov(fit)))
  data.frame(model = model, term = clean_term(names(est)),
             estimate = est,
             lo = est - 1.96 * se, hi = est + 1.96 * se,
             row.names = NULL)
}

# saom (RSiena sienaFit): estimate +/- 1.96 * SE from the covariance.
# THE LOUD WARNING: saom parameters are rate and evaluation-function
# quantities -- statements about *changes between waves*, not about
# tie probabilities in a snapshot. never facet these beside the
# tie-model rows above: rbind() them into their OWN rows object and
# draw a separate forest_plot(), comparing verdicts only within the
# family. (this is the estimand lesson of the payoff table -- "how do
# ties change?" is a different question, so its coefficients answer
# a different question too.)
row_saom <- function(sfit, model = "saom") {
  est <- sfit$theta
  se  <- sqrt(diag(sfit$covtheta))
  data.frame(model = model, term = sfit$requestedEffects$effectName,
             estimate = est, lo = est - 1.96 * se, hi = est + 1.96 * se,
             row.names = NULL)
}

# one forest plot from any rbind() of the rows above. baselines and
# structural terms (gwesp/gwdegree/gwdsp at any decay, sociality, ...)
# are dropped by default -- they have no cross-family analogue.
forest_plot <- function(rows, drop = "(baseline)") {
  rows <- rows[!rows$term %in% drop &
               !grepl("^gw(esp|degree|dsp)|sociality", rows$term), ]
  rows$model <- factor(rows$model, levels = rev(unique(rows$model)))
  msu_palette <- c("#18453B", "#008934", "#7BBD00",
                   "#0B9A6D", "#535054", "#A7A8AA")
  ggplot2::ggplot(rows,
                  ggplot2::aes(x = estimate, y = model, color = model)) +
    ggplot2::geom_vline(xintercept = 0, linetype = "dashed",
                        color = "grey55") +
    ggplot2::geom_pointrange(ggplot2::aes(xmin = lo, xmax = hi),
                             linewidth = 0.6, show.legend = FALSE) +
    ggplot2::scale_color_manual(values = msu_palette) +
    ggplot2::facet_wrap(~ term, scales = "free_x") +
    ggplot2::labs(x = "estimate (each family's own scale)", y = NULL)
}
