The main function provided by the package.
gergm(formula, covariate_data = NULL, beta_correlation_model = FALSE,
distribution_estimator = c("none", "rowwise-marginal", "joint"),
network_is_directed = TRUE, include_diagonal = FALSE,
number_of_networks_to_simulate = 500, MCMC_burnin = 100, thin = 1,
proposal_variance = 0.1, target_accept_rate = 0.25, seed = 123,
hyperparameter_optimization = FALSE, convex_hull_proportion = 0.9,
convex_hull_convergence_proportion = 0.9, sample_edges_at_a_time = 0,
parallel = FALSE, parallel_statistic_calculation = FALSE, cores = 1,
use_stochastic_MH = FALSE, stochastic_MH_proportion = 0.25,
slackr_integration_list = NULL, convergence_tolerance = 0.5,
MPLE_gain_factor = 0, acceptable_fit_p_value_threshold = 0.05,
normalization_type = c("log", "division"),
transformation_type = c("Cauchy", "LogCauchy", "Gaussian", "LogNormal"),
estimation_method = c("Metropolis", "Gibbs"),
downweight_statistics_together = TRUE, force_x_theta_updates = 1,
force_x_lambda_updates = 1, output_directory = NULL, output_name = NULL,
generate_plots = TRUE, verbose = TRUE, use_previous_thetas = TRUE,
fine_grained_pv_optimization = FALSE, use_MPLE_only = c(FALSE, TRUE),
start_with_zeros = FALSE, stop_for_degeneracy = FALSE,
maximum_number_of_lambda_updates = 10,
maximum_number_of_theta_updates = 10,
user_specified_initial_thetas = NULL, integration_intervals = 150,
distribution_mple_regularization_weight = 0.05,
theta_grid_optimization_list = NULL, weighted_MPLE = FALSE,
estimate_model = TRUE, optimization_method = c("BFGS", "L-BFGS-B",
"Nelder-Mead", "CG", "SANN", "Brent"), ...)A formula object that specifies the relationship between statistics and the observed network. Currently, the user may specify a model using any combination of the following statistics: `out2stars(alpha = 1)`, `in2stars(alpha = 1)`, `ctriads(alpha = 1)`, `mutual(alpha = 1)`, `ttriads(alpha = 1)`, `absdiff(covariate = "MyCov")`, `sender(covariate = "MyCov")`, `reciever(covariate = "MyCov")`, `nodematch(covariate)`, `nodemix(covariate, base = "MyBase")`, `netcov(network)` and `edges(alpha = 1, method = c("regression","endogenous"))`. If the user specifies `nodemix(covariate, base = NULL)`, then all levels of the covariate will be matched on. Note that the `edges` term must be specified if the user wishes to include an intercept (strongly recommended). The user may select the "regression" method (default) to include an intercept in the lambda transformation of the network, or "endogenous" to include the intercept as in a traditional ERGM model. To use exponential down-weighting for any of the network level terms, simply specify a value for alpha less than 1. The `(alpha = 1)` term may be omitted from the structural terms if no exponential down weighting is required. In this case, the terms may be provided as: `out2star`, `in2star`, `ctriads`, `mutual`, `ttriads`. If the network is undirected the user may only specify the following terms: `twostars(alpha = 1)`, `ttriads(alpha = 1)`, `absdiff(covariate = "MyCov")`, `sender(covariate = "MyCov")`, `nodematch(covariate)`, `nodemix(covariate, base = "MyBase")`, `netcov(network)` and `edges(alpha = 1, method = c("regression","endogenous"))`. In some cases, the user may only wish to calculate endogenous statistics for edges between some subset of the nodes in the network. For each of the endogenous statistics, the user may optionally specify a `covariate` and `base` field such as in `in2stars(covariate = "Type", base = "C")`. This will add an in-2star statistic for the subnetwork defined by actors who match each level of the categorical variable "Type" (in this example), and exclude the subnetwork for type "C", if the `base` argument is provided. If the `base` argument is excluded, then terms will be added to the model for all levels of the statistic. This can be a useful option if the user believes that a network property varies with some property of nodes.
A data frame containing node level covariates the user wished to transform into sender or receiver effects. It must have row names that match every entry in colnames(raw_network), should have descriptive column names. If left NULL, then no sender or receiver effects will be added.
Defaults to FALSE. If TRUE, then the beta correlation model is estimated. A correlation network must be provided, but all covariates and undirected statistics may be supplied as normal.
Provides an option to estimate the structure of row-wise marginal and joint distribtuions using a uniform-dirichlet proposal distribution. THIS FEATURE IS EXPERIMENTAL. Defaults to "none", in which case a normal GERGM is estimated, but can be set to one of "rowwise-marginal" and "joint" to propose either row-wsie marginal distribtuions or joint distributions. If an option other than "none" is selected, the beta correlation model will be turned off, estimation will automatically be set to Metropolis, no covariate data will be allowed, and the network will be set to directed. Furthermore, a "diagonal" statistic will be added to the model which simply records the sum of the diagonal of the network. The "mutual" statistic will also be adapted to include the diagonal elements. In the future, more statistics which take account of the network diagonal will be included.
Logical specifying whether or not the observed network is directed. Default is TRUE.
Logical indicating whether the diagonal should be included in the estimation proceedure. If TRUE, then a "diagonal" statistic is added to the model. Defaults to FALSE.
Number of simulations generated for estimation via MCMC. Default is 500.
Number of samples from the MCMC simulation procedure that will be discarded before drawing the samples used for estimation. Default is 100.
The proportion of samples that are kept from each simulation. For example, thin = 1/200 will keep every 200th network in the overall simulated sample. Default is 1.
The variance specified for the Metropolis Hastings simulation method. This parameter is inversely proportional to the average acceptance rate of the M-H sampler and should be adjusted so that the average acceptance rate is approximately 0.25. Default is 0.1.
The target Metropolis Hastings acceptance rate. Defaults to 0.25
Seed used for reproducibility. Default is 123.
Logical indicating whether automatic hyperparameter optimization should be used. Defaults to FALSE. If TRUE, then the algorithm will automatically seek to find an optimal burnin and number of networks to simulate, and if using Metropolis Hasings, will attempt to select a proposal variance that leads to a acceptance rate within +-0.05 of target_accept_rate. Furthermore, if degeneracy is detected, the algorithm will attempt to adress the issue automatically. WARNING: This feature is experimental, and may greatly increase runtime. Please monitor console output!
Defaults to 0.9. Must be a number between 0 and 1, with values between 0.5 and 0.9 preferred. Setting a value for this parameter makes use of the initialization method introduced by Hummel, Hunter and Handcock (2012), which can provide a much more stable initialization method than MPLE. Selecting this method will not supercede other initialization methods, as it will be run after MPLE or grid search. However, in practice, it should often preclude the need for other initialization methods. If NULL, then standard MPLE an grid search methods may be used. Not recommended
Defaults to 0.9. This must be a number between 0 and 1, with values between 0.7 and 0.95 preferred. This parameter controls the stopping behavior of the convex hull initialization proceedure, with a higher value requiring a better fit before moving to standard estimation.
Defaults to 0. If greater than zero, then this is the number of edges to be updated at once during MCMCMLE. The lower this number is set, the higher the Metropolis Hastings acceptance rate should be. This option will primarily be relevant for large networks.
Logical indicating whether the weighted MPLE objective and any other operations that can be easily parallelized should be calculated in parallel. Defaults to FALSE. If TRUE, a significant speedup in computation may be possible.
Logical indicating whether network statistics should be calculated in parallel. This will tend to be slower for networks with les than ~30 nodes but may provide a substantial speedup for larger networks.
Numeric value defaulting to 1. Can be set to any number up to the number of threads/cores available on your machine. Will be used to speed up computations if parallel = TRUE.
A logical indicating whether a stochastic approximation to the h statistics should be used under Metropolis Hastings in-between thinned samples. This may dramatically speed up estimation. Defaults to FALSE. HIGHLY EXPERIMENTAL!
Percentage of dyads/triads to use for approximation, defaults to 0.25.
An optional list object that contains information necessary to provide updates about model fitting progress to a Slack channel (https://slack.com/). This can be useful if models take a long time to run, and you wish to receive updates on their progress (or if they become degenerate). The list object must be of the following form: list(model_name = "descriptive model name", channel = "#yourchannelname", incoming_webhook_url = "https://hooks.slack.com/services/XX/YY/ZZ"). You will need to set up incoming webhook integration for your slack channel and then paste in the URL you get from slack into the incoming_webhook_url field. If all goes well, and the computer you are running the GERGM estimation on has internet access, your slack channel will receive updates when you start estimation, after each lambda/theta parameter update, if the model becomes degenerate, and when it completes running.
Threshold designated for stopping criterion. If the difference of parameter estimates from one iteration to the next all have a p -value (under a paired t-test) greater than this value, the parameter estimates are declared to have converged. Default is 0.5, which is quite conservative.
Multiplicative constant between 0 and 1 that controls how far away the initial theta estimates will be from the standard MPLEs via a one step Fisher update. In the case of strongly dependent data, it is suggested to use a value of 0.10. Default is 0.
A p-value threshold for how closely statistics of observed network conform to statistics of networks simulated from GERGM parameterized by converged final parameter estimates. Default value is 0.05.
If only a raw_network is provided the function will automatically check to determine if all edges fall in the [0,1] interval. If edges are determined to fall outside of this interval, then a trasformation onto the interval may be specified. If "division" is selected, then the data will have a value added to them such that the minimum value is at least zero (if necessary) and then all edge values will be divided by the maximum to ensure that the maximum value is in [0,1]. If "log" is selected, then the data will have a value added to them such that the minimum value is at least zero (if necessary), then 1 will be added to all edge values before they are logged and then divided by the largest value, again ensuring that the resulting network is on [0,1]. Defaults to "log" and need not be set to NULL if providing covariates as it will be ignored.
Specifies how covariates are transformed onto the raw network. When working with heavy tailed data that are not strictly positive, select "Cauchy" to transform the data using a Cauchy distribution. If data are strictly positive and heavy tailed (such as financial data) it is suggested the user select "LogCauchy" to perform a Log-Cauchy transformation of the data. For a tranformation of the data using a Gaussian distribution, select "Gaussian" and for strictly positive raw networks, select "LogNormal". The Default value is "Cauchy".
Simulation method for MCMC estimation. Default is "Metropolis", which allows for the most flexible model specifications, but may also be set to "Gibbs", if the user wishes to use Gibbs sampling.
Logical specifying whether or not the weights should be applied inside or outside the sum. Default is TRUE and user should not select FALSE under normal circumstances.
Defaults to 1 where theta estimation is not allowed to converge until thetas have updated for x iterations . Useful when model is not degenerate but simulated statistics do not match observed network well when algorithm stops after first y updates.
Defaults to 1 where lambda estimation is not allowed to converge until lambdas have updated for x iterations . Useful when model is not degenerate but simulated statistics do not match observed network well when algorithm stops after first y updates.
The directory where you would like output generated by the GERGM estimation procedure to be saved (if output_name is specified). This includes, GOF, trace, and parameter estimate plots, as well as a summary of the estimation procedure and an .Rdata file containing the GERGM object returned by this function. May be left as NULL if the user would prefer all plots be printed to the graphics device.
The common name stem you would like to assign to all objects output by the gergm() function. Default value of NULL will not save any output directly to .pdf files, it will be printed to the console instead. Must be a character string or NULL. For example, if "Test" is supplied as the output_name, then 4 files will be output: "Test_GOF.pdf", "Test_Parameter_Estim ates.pdf", "Test_GERGM_Object.Rdata", "Test_Estimation_Log.txt", and "Test_Trace_Plot.pdf"
Defaults to TRUE, if FALSE, then no diagnostic or parameter plots are generated.
Defaults to TRUE (providing lots of output while model is running). Can be set to FALSE if the user wishes to see less output.
Logical, defaults to TRUE. IF TRUE, then at each iteration of covariate parameter updates beyond the first, the MPLE initialization for theta parameters will be skipped, and the previous iteration thetas will be used as a starting point. This option will only work with convex hull initialization. The intuition here is that if MPLE is poor, the previous theta estimates may be a better starting point for convex hull initialization. In some cases this can lead to a very large speedup in estimation. However, this can lead to arbitrarily poor performance in some situations.
Logical indicating whether fine grained proposal variance optimization should be used. This will often slow down proposal variance optimization, but may provide better results. Highly recommended if running a correlation model.
Logical specifying whether or not only the maximum pseudo likelihood estimates should be obtained. In this case, no simulations will be performed. Default is FALSE.
Defaults to FALSE. IF TRUE, then no MPLE is used and allendogenous parameters are initialized to zero. This method will still work with convex_hull_proportion.
When TRUE, automatically stops estimation when degeneracy is detected, even when hyperparameter_optimization is set to TRUE. Defaults to FALSE.
Maximum number of iterations of outer MCMC loop which alternately estimates transform parameters and ERGM parameters. In the case that data_transformation = NULL, this argument is ignored. Default is 10.
Maximum number of iterations within the MCMC inner loop which estimates the ERGM parameters. Default is 100.
Optional numeric vector of user specified theta values to be used for initialization (instead of MPLE). This option should not be used exept when MPLE performance is poor. Defaults to NULL.
The number of intervals to be used for numerical integration in weighted MPLE and in MPLE for the distribution estimator. Defaults to 150 but may be increased if more accuracy is required.
L2 regularization of the MPLE theta estimates may be necessary to provide an adequate starting position when using the distribtuion estimator. We suggest a value of 0.05, but the optimal L2 penalty will be application specific. Setting to zero removes all L2 regularization.
Defaults to NULL. This highly experimental feature may allow the user to address model degeneracy arising from a suboptimal theta initialization. It performs a grid search around the theta values calculated via MPLE to select a potentially improved initialization. The runtime complexity of this feature grows exponentially in the size of the grid and number of parameters -- use with great care. This feature may only be used if hyperparameter_optimization = TRUE, and if a list object of the following form is provided: list(grid_steps = 2, step_size = 0.5, cores = 2, iteration_fraction = 0.5). grid_steps indicates the number of steps out the grid search will perform, step_size indicates the fraction of the MPLE theta estimate that each grid search step will change by, cores indicates the number of cores to be used for parallel optimization, and iteration_fraction indicates the fraction of the number of MCMC iterations that will be used for each grid point (should be set less than 1 to speed up optimization). In general grid_steps should be smaller the more structural parameters the user wishes to specify. For example, with 5 structural parameters (mutual, ttriads, etc.), grid_steps = 3 will result in a (2*3+1)^5 = 16807 parameter grid search. Again this feature is highly experimental and should only be used as a last resort (after playing with exponential down weighting and the MPLE_gain_factor).
Defaults to FALSE. Should be used whenever the user is specifying statistics with alpha down weighting. Tends to provide better initialization when downweight_statistics_together = FALSE.
Logical indicating whether a model should be estimated. Defaults to TRUE, but can be set to FALSE if the user simply wishes to return a GERGM object containing the model specification. Useful for debugging.
The optimization method used by the 'optim' function in estimating theta and lambda parameter estimates. Defualts to "BFGS", but can also be any one of "L-BFGS-B", "Nelder-Mead", "CG", "SANN", or "Brent". "L-BFGS-B" is preferred for fitting beta correlation models.
Optional arguments, currently unsupported.
A gergm object containing parameter estimates.
# NOT RUN {
set.seed(12345)
net <- matrix(rnorm(100,0,20),10,10)
colnames(net) <- rownames(net) <- letters[1:10]
formula <- net ~ mutual + ttriads
test <- gergm(formula,
normalization_type = "division",
network_is_directed = TRUE,
use_MPLE_only = FALSE,
estimation_method = "Metropolis",
number_of_networks_to_simulate = 40000,
thin = 1/10,
proposal_variance = 0.5,
downweight_statistics_together = TRUE,
MCMC_burnin = 10000,
seed = 456,
convergence_tolerance = 0.01,
MPLE_gain_factor = 0,
force_x_theta_updates = 4)
# }
Run the code above in your browser using DataLab