remstimate - optimization of tie-oriented and actor-oriented likelihood

A function for the optimization of tie-oriented and actor-oriented likelihood. There are four optimization algorithms: two Frequentists, Maximum Likelihood Estimation (MLE) and Adaptive Gradient Descent (GDADAMAX), and two Bayesian, Bayesian Sampling Importance Resampling (BSIR) and Hamiltonian Monte Carlo (HMC).

Usage

remstimate(
  reh,
  stats,
  method = c("MLE", "GDADAMAX", "BSIR", "HMC"),
  ncores = attr(reh, "ncores"),
  prior = NULL,
  nsim = 1000L,
  nchains = 1L,
  burnin = 500L,
  thin = 10L,
  init = NULL,
  epochs = 1000L,
  L = 50L,
  epsilon = ifelse(method == "GDADAMAX", 0.001, 0.002),
  seed = NULL,
  WAIC = FALSE,
  silent = TRUE,
  ...
)

Arguments

reh

a remify object of the processed relational event history. Output object of the function remify::remify().

stats

a remstats object: when `attr(reh,"model")` is `"tie"`, stats is an array of statistics with dimensions [M x D x P]: where M is the number of events, D is the number of possible dyads (full riskset), P is the number of statistics; if `attr(reh,"model")` is `"actor"`, stats is a list that can contain up to two arrays named "sender_stats" and "receiver_stats" with dimensions [M x N x P], where N are the actors (senders in the array "sender_stats", receivers in the array "receiver_stats"). Furthermore, it is possible to only estimate the sender rate model or only the receiver choice model, by using the correct naming of the arrays.

method

the optimization method to estimate model parameters. Methods available are: Maximum Likelihood Estimation ("MLE", and also the default method), Adaptive Gradient Descent ("GDADAMAX") based on the work of Diederik P. Kingma and Jimmy Ba, 2014 (<doi:10.48550/arXiv.1412.6980>), Bayesian Sampling Importance Resampling ("BSIR"), Hamiltonian Monte Carlo ("HMC"). (default method is "MLE").

ncores

[optional] number of threads for the parallelization. (default value is 1, which means no parallelization).

prior

[optional] prior distribution when method is "BSIR". Default value is NULL, which means that no prior is assumed. For the tie-oriented modeling, the argument prior is the name of the function in the format name_package::name_density_function. The parameters of the prior distribution can be supplied as inputs to the remstimate function (e.g., remstimate::remstimate(reh=reh,stats=stats,method="BSIR",ncores=5,prior=mvnfast::dmvn,mu=rep(0,3),sigma=diag(3)*2,log=TRUE) ). For actor-oriented modeling the argument prior is a named list of two objects "sender_model", which calls the prior function for the sender rate model, and, "receiver_model", which calls the prior function for the receiver choice model. For the specification of the prior parameters, the user must define an optional argument called prior_args, which is also a named list (with names "sender_model" and "receiver_model"): each list is a list of objects named after the prior arguments and with value of the prior argument (e.g., prior_args$sender_model = list(mu = rep(1.5,3), sigma = diag(3)*0.5, log = TRUE)). Finally, both in tie-oriented and actor-oriented modeling prior functions must have an argument that returns the value of the density on a logarithmic scale (i.e., log=TRUE). log=TRUE is already set up internally by remstimate().

nsim

[optional] when method is "HMC", nsim is the number of simulations (iterations) in each chain, when method is "BSIR", then nsim is the number of samples from the proposal distribution. Default value is 1000.

nchains

[optional] number of chains to generate in the case of method = "HMC". Default value is 1.

burnin

[optional] number of initial iterations to be added as burnin for method = "HMC". Default value is 500.

thin

[optional] number of steps to skip in the posterior draws for method = "HMC". Default value is 10. If nsim<100, thin is set to 1.

init

[optional] vector of initial values if tie-oriented model, or a named list of two vectors ('sender_model' and 'receiver_model') if both models of the actor-oriented framework are specified. init can also be a list of only one vector (named 'sender_model' or 'receiver_model'), if the interest is to estimate one specific model of the actor-oriented framework. init is used for the methods "GDADAMAX" and "HMC". If init is NULL, then it will be assigned internally.

epochs

[optional] It is the number of iteration used in the method "GDADAMAX". Default value is 1000.

L

[optional] number of leap-frog steps to use in the method "HMC". Default value is 50.

epsilon

[optional] It is a parameter used in two methods: if method is "GDADAMAX", it represents the inter-iteration difference of the loss function and it is used as stop-rule within the algorithm (default value is 0.001), if method is "HMC" (default value is 0.002), it is a parameter used in the leap-frog algorithm and it is proportional to the step size.

seed

[optional] seed value for reproducibility. If NULL, seed will be assigned by the machine and saved in the output object.

WAIC

[optional] logical value. The Watanabe Akaike's Information Criterion (WAIC) will be calculated is WAIC = TRUE. The default number of simulations used to calculate the WAIC is 500. In order to specify a different number of simulations, the user must supply an additional argument nsimWAIC to the function.

silent

[optional-not-yet-implemented] a TRUE/FALSE value. If FALSE, progress of optimization status will be printed out.

...

additional parameters. They can be parameters of other functions defined as input in some of the arguments above. (e.g., arguments of the prior distribution)

Value

'remstimate' S3 object.

Examples

# ------------------------------------ #
#       tie-oriented model: "MLE"      #
# ------------------------------------ #

# loading data
data(tie_data)

# processing event sequence with remify
tie_reh <- remify::remify(edgelist = tie_data$edgelist, model = "tie")
  
# specifying linear predictor
tie_model <- ~ 1 + 
               remstats::indegreeSender()+
               remstats::inertia()+
               remstats::reciprocity() 

# calculating statistics
tie_reh_stats <- remstats::remstats(reh = tie_reh, 
                                    tie_effects = tie_model)

# running estimation
tie_mle <- remstimate::remstimate(reh = tie_reh,
                                  stats = tie_reh_stats,
                                  method = "MLE",
                                  ncores = 1)
# summary
summary(tie_mle)
#> Relational Event Model (tie oriented) 
#> 
#> Call:
#> ~1 + remstats::indegreeSender() + remstats::inertia() + remstats::reciprocity()
#> 
#> 
#> Coefficients (MLE with interval likelihood):
#> 
#>                  Estimate   Std. Err    z value Pr(>|z|) Pr(=0)
#> baseline        -4.910454   0.187555 -26.181372   0.0000 <2e-16
#> indegreeSender   0.043490   0.036449   1.193170   0.2328 0.8307
#> inertia         -0.201506   0.088154  -2.285831   0.0223 0.4231
#> reciprocity     -0.052137   0.098237  -0.530728   0.5956 0.8968
#> Null deviance: 1216.739 on 100 degrees of freedom
#> Residual deviance: 1210.625 on 96 degrees of freedom
#> Chi-square: 6.11449 on 4 degrees of freedom, asymptotic p-value 0.1907597 
#> AIC: 1218.625 AICC: 1219.046 BIC: 1229.045

# ------------------------------------ #
#      actor-oriented model: "MLE"     #
# ------------------------------------ #

# loading data
data(ao_data)

# processing event sequence with remify
ao_reh <- remify::remify(edgelist = ao_data$edgelist, model = "actor")
  
# specifying linear predictor (for sender rate and receiver choice model)
rate_model <- ~ 1 + remstats::indegreeSender()
choice_model <- ~ remstats::inertia() + remstats::reciprocity()

# calculating statistics
ao_reh_stats <- remstats::remstats(reh = ao_reh, 
                                   sender_effects = rate_model, 
                                   receiver_effects = choice_model)

# running estimation
ao_mle <- remstimate::remstimate(reh = ao_reh,
                                 stats = ao_reh_stats,
                                 method = "MLE",
                                 ncores = 1)
# summary
summary(ao_mle)
#> Relational Event Model (actor oriented) 
#> 
#> Call rate model **for sender**:
#> 
#> 	~1 + remstats::indegreeSender()
#> 
#> 
#> Coefficients rate model (MLE with interval likelihood):
#> 
#>                   Estimate    Std. Err     z value Pr(>|z|) Pr(=0)
#> baseline        -4.8062015   0.1832563 -26.2266709      0.0 <2e-16
#> indegreeSender  -0.0083633   0.0159467  -0.5244534      0.6 0.8971
#> Null deviance: 1177.625 on 100 degrees of freedom
#> Residual deviance: 1177.348 on 98 degrees of freedom
#> Chi-square: 0.2770484 on 2 degrees of freedom, asymptotic p-value 0.8706422 
#> AIC: 1181.348 AICC: 1181.472 BIC: 1186.558
#> -------------------------------------------------------------------------------- 
#> 
#> Call choice model **for receiver**:
#> 
#> 	~remstats::inertia() + remstats::reciprocity()
#> 
#> 
#> Coefficients choice model (MLE with interval likelihood):
#> 
#>              Estimate  Std. Err   z value Pr(>|z|) Pr(=0)
#> inertia     -0.032508  0.077037 -0.421983   0.6730 0.9015
#> reciprocity  0.018287  0.071630  0.255305   0.7985 0.9064
#> Null deviance: 2.772589 on 100 degrees of freedom
#> Residual deviance: 277.0551 on 98 degrees of freedom
#> Chi-square: -274.2826 on 2 degrees of freedom, asymptotic p-value 1 
#> AIC: 281.0551 AICC: 281.1789 BIC: 286.2655

# ------------------------------------ #
#   for more examples check vignettes  #
# ------------------------------------ #