The Blume-Capel model is the full Spin distribution model. In addition to the linear node potentials (tau) and pairwise interactions (omega) of the Ising model, it estimates a quadratic node potential delta (the Blume-Capel single-ion / crystal-field term). The Ising model is the special case obtained by fixing all delta parameters to zero, so the two models are nested. By convention the Blume-Capel model is defined for the ternary c(-1, 0, 1) response coding (a warning is raised for any other response set), but any number of ordered numeric response options is allowed (the same set for every variable, and not necessarily integers).
BlumeCapel(data, omega = "full", tau, delta, beta, beta_model =
c("beta", "log_beta"), vars, groups, covs, means,
nobs, covtype = c("choose", "ML", "UB"), responses,
missing = "listwise", equal = "none",
baseline_saturated = TRUE, estimator = "default",
optimizer, storedata = FALSE, WLS.W, sampleStats,
identify = TRUE, verbose = FALSE, maxNodes = 20,
maxStates = 2^maxNodes, min_sum = -Inf,
bootstrap = FALSE, boot_sub,
boot_resample, penalty_lambda = NA,
penalty_alpha = 1, penalize_matrices)An object of the class psychonetrics
A data frame encoding the data used in the analysis. Can be missing if covs and nobs are supplied.
The network structure. Either "full" to estimate every element freely, "zero" to set all elements to zero, or a matrix of the dimensions nNode x nNode with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Optional vector encoding the threshold/intercept (linear node potential) structure. Set elements to 0 to indicate fixed to zero constrains, 1 to indicate free intercepts, and higher integers to indicate equality constraints. For multiple groups, this argument can be a list or array with each element/column encoding such a vector.
Optional vector encoding the quadratic node potential (Blume-Capel) structure. Set elements to 0 to indicate fixed to zero constrains, 1 to indicate free parameters, and higher integers to indicate equality constraints. For multiple groups, this argument can be a list or array with each element/column encoding such a vector. When missing, every delta parameter is freely estimated (fixing all of them to zero would instead yield the Ising model).
Optional scalar encoding the inverse temperature. 1 indicate free beta parameters, and higher integers to indicate equality constraints. For multiple groups, this argument can be a list or array with each element/column encoding such scalars.
How should beta be modeled? Set beta_model = "log_beta" to model the log of beta rather than beta directly.
An optional character vector encoding the variables used in the analysis. Must equal names of the dataset in data.
An optional string indicating the name of the group variable in the data.
A sample variance--covariance matrix, or a list/array of such matrices for multiple groups. Make sure covtype argument is set correctly to the type of covariances used.
A vector of sample means, or a list/matrix containing such vectors for multiple groups.
The number of observations used in covs and means, or a vector of such numbers of observations for multiple groups.
If 'covs' is used, this is the type of covariance (maximum likelihood or unbiased) the input covariance matrix represents. Set to "ML" for maximum likelihood estimates (denominator n) and "UB" to unbiased estimates (denominator n-1). The default will try to find the type used, by investigating which is most likely to result from integer valued datasets.
A vector of the response options used, encoded identically across all variables. Defaults to (and is conventionally) c(-1, 0, 1); a warning is raised for any other set. Automatically detected from the data when not supplied; required when covs is used. Any number of distinct response options is allowed, and the values need not be integers.
How should missingness be handled when data is used. Only "listwise" (listwise deletion) is currently supported for the Blume-Capel model; other options (e.g. "pairwise") are rejected at construction because they produce undefined sufficient statistics.
A character vector indicating which matrices should be constrained equal across groups.
A logical indicating if the baseline and saturated model should be included. Mostly used internally and NOT Recommended to be used manually.
The estimator to be used. Only "ML" (maximum likelihood) estimation is currently supported for the Blume-Capel model.
The optimizer to be used. Can be one of "nlminb" (the default R nlminb function), "ucminf" (from the optimr package), "nloptr_TNEWTON" (preconditioned truncated Newton via nloptr), and "LBFGS++" (pure C++ L-BFGS-B). Defaults to "nlminb".
Logical, should the raw data be stored? Needed for bootstrapping (see bootstrap).
Optional WLS weights matrix. CURRENTLY NOT USED.
An optional sample statistics object. Mostly used internally.
Logical, should the model be identified?
Logical, should messages be printed?
The maximum number of nodes allowed in the analysis. Used to set the default of maxStates (2^maxNodes); for binary data this reproduces the historical node limit. It is not recommended to set this higher.
The maximum number of response patterns the exact ML estimator may enumerate. Exact ML estimation sums over every possible response pattern when computing the partition function, expected values and (expected) Hessian, so the cost grows as length(responses)^nNode. The function stops with an error when length(responses)^nNode exceeds maxStates. The default 2^maxNodes reproduces the historical binary node limit while accounting for more than two response options. Raise this only if the computation is feasible (note that the number of states---and hence the run time---grows very quickly with both the number of nodes and the number of response options).
The minimum sum score that is artificially possible in the dataset. Defaults to -Inf. Set this only if you know a lower sum score is not possible in the data, for example due to selection bias.
Should the data be bootstrapped? If TRUE the data are resampled and a bootstrap sample is created. These must be aggregated using aggregate_bootstraps! Can be TRUE or FALSE. Can also be "nonparametric" (which sets boot_sub = 1 and boot_resample = TRUE) or "case" (which sets boot_sub = 0.75 and boot_resample = FALSE).
Proportion of cases to be subsampled (round(boot_sub * N)).
Logical, should the bootstrap be with replacement (TRUE) or without replacement (FALSE)
Numeric penalty strength for penalized ML estimation (PML/PFIML). NA (default) triggers automatic selection via EBIC-based grid search when a penalized estimator is used; set to a specific numeric value to use a fixed penalty strength (0 = no penalty). See find_penalized_lambda and penalize.
Elastic net mixing parameter: 1 = LASSO (default), 0 = ridge.
Character vector of matrix names to penalize. If missing, defaults are selected based on the model type (the network omega).
Sacha Epskamp <mail@sachaepskamp.com>
The Blume-Capel model is built on the Spin distribution:
\(\Pr(\boldsymbol{Y} = \boldsymbol{y}) = \frac{\exp\left( -\beta H\left(\boldsymbol{y}; \boldsymbol{\tau}, \boldsymbol{\delta}, \boldsymbol{\Omega}\right)\right)}{Z(\boldsymbol{\tau}, \boldsymbol{\delta}, \boldsymbol{\Omega})}\)
With Hamiltonian:
\(H\left(\boldsymbol{y}; \boldsymbol{\tau}, \boldsymbol{\delta}, \boldsymbol{\Omega}\right) = -\sum_{i=1}^{m} \tau_i y_{i} + \sum_{i=1}^{m} \delta_i y_{i}^2 - \sum_{i=2}^{m} \sum_{j=1}^{i-1} \omega_{ij} y_i y_j.\)
and Z representing the partition function or normalizing constant. Equivalently,
\(\Pr(\boldsymbol{Y} = \boldsymbol{y}) \propto \exp\left( \sum_i \tau_i y_i - \sum_i \delta_i y_i^2 + \sum_{i<j} \omega_{ij} y_i y_j \right).\)
For the ternary states c(-1, 0, 1) this is the usual Blume-Capel parameterization: tau controls the tendency toward the positive versus the negative state, delta controls the tendency toward the zero/middle state versus the active/extreme states (positive delta favours the middle category), and omega controls pairwise alignment. For more general ordered numeric states the same parameterization defines an ordinal spin model.
The Ising model is nested within the Blume-Capel model by fixing all delta parameters to zero. Exact ML estimation enumerates every response pattern, so the cost grows as length(responses)^nNode (see maxStates); the model is therefore limited to relatively small networks.
Epskamp, S., Maris, G., Waldorp, L. J., & Borsboom, D. (2018). Network Psychometrics. In: Irwing, P., Hughes, D., & Booth, T. (Eds.), The Wiley Handbook of Psychometric Testing, 2 Volume Set: A Multidisciplinary Reference on Survey, Scale and Test Development. New York: Wiley.
Ising
# \donttest{
library("dplyr")
# Simulate a small Blume-Capel data set (requires the development version of
# IsingSampler that supports the 'delta' argument):
if (FALSE) {
library("IsingSampler")
nNode <- 5
graph <- 0.4 * (matrix(1, nNode, nNode) - diag(nNode)) *
(matrix(c(0,1,0,0,1, 1,0,1,0,0, 0,1,0,1,0, 0,0,1,0,1, 1,0,0,1,0), 5, 5))
data <- IsingSampler(1000, graph = graph, thresholds = rep(0, nNode),
beta = 1, delta = 1, responses = c(-1, 0, 1))
# Fit the Blume-Capel model:
mod <- BlumeCapel(data, responses = c(-1, 0, 1)) %>% runmodel
# Estimated network:
getmatrix(mod, "omega")
# Estimated quadratic node potentials:
getmatrix(mod, "delta")
}
# }
Run the code above in your browser using DataLab