R/aomstats.R

In remstats: Computes Statistics for Relational Event History Data

#' aomstats
#'
#' Computes statistics for the sender activity rate step and receiver choice
#' step in actor-oriented relational event models (e.g., see Stadtfeld & Block,
#' 2017).
#'
#' @inheritParams remstats
#'
#' @section Effects:
#' The statistics to be computed are defined symbolically and should be
#' supplied to the \code{sender_effects} and/or \code{receiver_effects}
#' arguments in the form \code{~ effects}. The terms are separated by +
#' operators. For example: \code{receiver_effects = ~ inertia() + otp()}.
#' Interactions between two effects can be included with * or :
#' operators. For example: \code{receivereffects = ~ inertia():otp()}.  A list
#' of available effects can be obtained with \code{\link{actor_effects}()}.
#'
#' The majority of the statistics can be scaled in some way, see
#' the documentation of the \code{scaling} argument in the separate effect
#' functions for more information on this.
#'
#' @section Memory:
#' The default `memory` setting is `"full"`, which implies that at each time
#' point $t$ the entire event history before $t$ is included in the computation
#' of the statistics. Alternatively, when `memory` is set to `"window"`, only
#' the past event history within a given time window is considered (see
#' Mulders & Leenders, 2019). This length of this time window is set by the
#' `memory_value` parameter. For example, when `memory_value = 100` and `memory
#' = "window"`, at time point $t$ only the past events that happened at most
#' 100 time units ago are included in the computation of the statistics.
#' A third option is to set `memory` to `"interval"`. In this case, the past
#' event history within a given time interval is considered. For example, when
#' `"memory_value" = c(50, 100)` and `memory = "interval"`, at time point $t$
#' only the past events that happened between 50 and 100 time units ago are
#' included in the computation of the statistics. Finally, the fourth option is
#' to set `memory` to `"decay"`. In this case, the weight of the past event in
#' the computation of the statistics depend on the elapsed time between $t$ and
#' the past event. This weight is determined based on an exponential decay
#' function with half-life parameter `memory_value` (see Brandes et al., 2009).
#'
#' @section Event weights:
#' Note that if the relational event history contains a column that is named
#' ``weight'', it is assumed that these affect the endogenous statistics. These
#' affect the computation of all endogenous statistics with a few exceptions
#' that follow logically from their definition (e.g., the recenyContinue
#' statistic does depend on time since the event and not on event weights).
#'
#' @section Subset of the relational event history:
#' Optionally, statistics can be computed for a slice of the relational event
#' sequence - but based on the entire history. This is achieved by setting the
#' first and last values equal to the index of the first and last event for
#' which statistics are requested. For example, first = 5 and last = 5 computes
#' the statistics for only the 5th event in the relational event sequence,
#' based on the history that consists of events 1-4.
#'
#' @return An object of class 'aomstats'. List with in the first element the
#' statistics for the sender activity rate step and in the second element the
#' statistics for the receiver choice step. Statistics are computed once per 
#' unique time point (per-timepoint "pt" method), so the number of rows in the 
#' output equals \code{reh$M} (the number of  unique time points), which may be 
#' less than the total number of observed events when simultaneous events are 
#' present. The 'aomstats' object has the following attributes: 
#'   \describe{
#'     \item{\code{model}}{Type of model that is estimated.}
#'     \item{\code{formula}}{Model formula(s), obtained from the formula(s) 
#'     inputted to 'sender_effects' and/or 'receiver_effects'.}
#'     \item{\code{actors}}{The set of actors used to construct the statistics, 
#'     obtained from the remify object inputted to 'reh'.}
#'   }
#'
#' @examples
#' library(remstats)
#' 
#' # Load the data
#' data(history)
#' data(info)
#' 
#' # Prepare the data
#' reh <- remify::remify(edgelist = history, model = "actor")
#' 
#' # Define the sender effects
#' seff <- ~ send("extraversion")
#' 
#' # Define the receiver_effects
#' reff <- ~ receive("agreeableness") + inertia() + otp()
#' 
#' # Compute the statistics
#' aomstats(
#'   reh = reh, sender_effects = seff, receiver_effects = reff,
#'   attr_actors = info
#' )
#'
#' @references Stadtfeld, C., & Block, P. (2017). Interactions, actors, and
#' time: Dynamic network actor models for relational events. Sociological
#' Science, 4, 318–352. \doi{10.15195/v4.a14}
#'
#' Meijerink-Bosman, M., Back, M., Geukes, K., Leenders, R., & Mulder, J.
#' (2023). Discovering trends of social interaction behavior over time: An
#' introduction to relational event modeling: Trends of social interaction.
#' Behavior Research Methods, 55(3), 997-1023. \doi{10.3758/s13428-022-01821-8}
#'
#' @export
aomstats <- function(reh,
                     sender_effects = NULL,
                     receiver_effects = NULL,
                     memory = c("full", "window", "decay", "interval"),
                     memory_value = NA,
                     first = 2,
                     last = Inf,
                     display_progress = FALSE,
										 attr_actors = NULL,
										 attr_dyads = NULL) {

	start <- first
	stop <- last
	method <- "pt"
	
	if (!is.null(attr_actors)) {
		warning("'attr_actors' is deprecated. Supply attributes directly in the stat functions (e.g., send()). This argument is ignored.",
						call. = FALSE)
	}
	if (!is.null(attr_dyads)) {
		warning("'attr_dyads' is deprecated. Supply attributes directly in the stat functions (e.g., send()). This argument is ignored.",
						call. = FALSE)
	}
	
  # Validate remaining aomstats arguments
	attr_actors <- validate_aomstats_arguments(attr_actors, reh)	

  # Prepare the edgelist
  edgelist <- prepare_aomstats_edgelist(reh)

  # Prepare the network actors
  actors <- prepare_aomstats_actors(reh)

  # Prepare event weights
  weights <- prepare_aomstats_event_weights(reh)

  # Validate the memory argument
  memory <- match.arg(memory)
  memory_value <- validate_memory(memory, memory_value)
  if (memory == "window") {
    # Change memory to interval (window is a special case)
    memory <- "interval"
    memory_value <- c(0, memory_value)
  }

  # Validate the method
  method <- "pt"
  N_actors <- nrow(actors)

  # Types info — needed for consider_type="separate" per effect
  reh_norm  <- normalize_reh(reh)
  with_type <- isTRUE(reh_norm$meta$with_type)
  types_df  <- reh_norm$meta$dictionary$types  # data.frame: typeName, typeID
  C         <- if (with_type && !is.null(types_df)) nrow(types_df) else 1L
  # Per-event type IDs (1-based) — used when any effect has consider_type="separate"
  event_type_ids <- if (with_type && C > 1L) unlist(reh_norm$ids$type) else NULL

  # Initialize stats
  sender_stats <- NULL
  receiver_stats <- NULL

  # Sender model ------------------------------------------------------------
  sender_formula <- sender_effects
  if (!is.null(sender_formula)) {
    # Prepare sender_effects
    temp <- prepare_sender_effects(sender_formula, ordinal = reh$meta$ordinal)
    sender_effects <- temp$sender_effects
    sender_effects_names <- temp$sender_effects_names
    sender_interactions <- temp[[3]]

    # Prepare subset arguments
    subset <- prepare_subset(start, stop, edgelist, method, model = "sender")

    # Prepare sender covariate information
    sender_covar <- prepare_sender_covariates(sender_effects, attr_actors,
      actors, edgelist, reh)

    # Prepare sender_effects scaling
    sender_scaling <- prepare_aomstats_scaling(sender_effects,
      sender_interactions)

    # Extract consider_type per effect (default "ignore")
    s_consider_type <- sapply(sender_effects, function(x) {
      ct <- if ("consider_type" %in% names(x)) x$consider_type else "ignore"
      if (!ct %in% c("ignore", "separate"))
        stop("consider_type for sender effects must be 'ignore' or 'separate'")
      ct
    })
    s_consider_type <- c(s_consider_type,
      rep("ignore", sum(!sapply(sender_interactions, is.null))))

    sender_stats <- compute_aomstats_with_type(
      compute_fn     = function(eff_names, el, wts, s0, s1)
        compute_stats_sender(eff_names, el, actors[, 2], wts, sender_covar,
          sender_interactions, memory, memory_value, sender_scaling,
          s0, s1, method, display_progress),
      add_names_fn   = function(st) add_variable_names(st, sender_effects_names,
                         sender_effects, sender_interactions),
      edgelist       = edgelist,
      weights        = weights,
      event_type_ids = event_type_ids,
      types_df       = types_df,
      effects_names  = sender_effects_names,
      consider_type  = s_consider_type,
      with_type      = with_type,
      C              = C,
      subset_start   = subset$start,
      subset_stop    = subset$stop,
      N              = N_actors
    )
  }

  # Receiver model ----------------------------------------------------------
  receiver_formula <- receiver_effects
  if (!is.null(receiver_effects)) {
    # Prepare receiver_effects
    temp <- prepare_receiver_effects(receiver_formula)
    receiver_effects <- temp$receiver_effects
    receiver_effects_names <- temp$receiver_effects_names
    receiver_interactions <- temp$receiver_interactions

    # Prepare subset arguments
    subset <- prepare_subset(start, stop, edgelist, method, model = "sender") # if model is set to 'receiver' the dimension is based on number of events rather than time points.

    # Prepare receiver covariate information
    receiver_covar <- prepare_receiver_covariates(receiver_effects,
      attr_actors, attr_dyads, actors, edgelist, reh)

    # Prepare receiver_effects scaling
    receiver_scaling <- prepare_aomstats_scaling(receiver_effects,
      receiver_interactions)

    # Extract consider_type per effect (default "ignore")
    r_consider_type <- sapply(receiver_effects, function(x) {
      ct <- if ("consider_type" %in% names(x)) x$consider_type else "ignore"
      # backward compatibility
      if (isTRUE(ct))  ct <- "separate"
      if (isFALSE(ct)) ct <- "ignore"
      if (!ct %in% c("ignore", "separate"))
        stop("consider_type for receiver effects must be 'ignore' or 'separate'")
      ct
    })
    r_consider_type <- c(r_consider_type,
      rep("ignore", sum(!sapply(receiver_interactions, is.null))))

    receiver_stats <- compute_aomstats_with_type(
      compute_fn     = function(eff_names, el, wts, s0, s1)
        compute_stats_receiver(eff_names, el, actors[, 2], wts, receiver_covar,
          receiver_interactions, memory, memory_value, receiver_scaling,
          s0, s1, method, display_progress),
      add_names_fn   = function(st) add_variable_names(st, receiver_effects_names,
                         receiver_effects, receiver_interactions),
      edgelist       = edgelist,
      weights        = weights,
      event_type_ids = event_type_ids,
      types_df       = types_df,
      effects_names  = receiver_effects_names,
      consider_type  = r_consider_type,
      with_type      = with_type,
      C              = C,
      subset_start   = subset$start,
      subset_stop    = subset$stop,
      N              = N_actors
    )
  }

  # Subset output
  subset <- prepare_subset(start, stop, edgelist, method, "sender")
  subset$start <- subset$start + 1
  subset$stop <- subset$stop + 1
  subset <- as.data.frame(subset)

  # Output
  out <- list(sender_stats = sender_stats, receiver_stats = receiver_stats)
  class(out) <- c("aomstats", "remstats")
  attr(out, "model") <- "actor"
  attr(out, "formula") <- list(rate = sender_formula, choice = receiver_formula)
  attr(out, "actors") <- actors
  attr(out, "subset") <- subset
  attr(out, "method") <- method
  out
}