Fits a Discrete Exponential-Family Model using Maximum Likelihood.
An object of class DEFM.
The parameters of the model.
The row and column of the array to turn on for the log odds.
An object of class DEFM.
Double vector or named list. Starting point for the MLE (defaults to zero). When named, values are matched to the model terms.
Lower and upper limits for the optimization (only used
with method = "L-BFGS-B"). Default to -20 and 20.
Character scalar. Either "L-BFGS-B" (default) or "BFGS",
passed to stats::optim(), or "newton" for Newton-Raphson with step
halving (see details).
Named list or named numeric vector. Parameters held fixed at the given values.
Integer scalar. Number of threads used to compute the log-likelihood and its derivatives (only relevant if the package was compiled with OpenMP).
List of control parameters. Passed to stats::optim();
with method = "newton", maxit (default 100), reltol (default
1e-10), and grtol (default 1e-8, maximum absolute gradient) are used.
Further arguments passed to summary_table (with the exception of as_texreg, which is set to TRUE).
Either a single or a list of defm fit objects.
Function to be called from the texreg package, e.g., texreg::screenreg.
Whether or not to skip the intercept (logit) terms when printing the table
logodds returns a numeric vector with the log-odds for each observation in the data.
An object of class defm_mle, which extends stats4::mle
(all the methods for mle objects, e.g., coef(), vcov(), and
summary(), work with it). It also stores the names of the outcomes and
covariates, which are used by summary_table() to build a compact table
of the estimates.
An object of class texreg with additional attributes: custom.coef.map,
reorder.coef, and groups.
The likelihood function of the DEFM is closely-related to the Exponential-Family Random Graph Model [ERGM]. Furthermore, the DEFM can be treated as a generalization of the ERGM. The model implemented here can be viewed as an ERGM for a bipartite network, where the actors are individuals and the events are the binary outputs.
If the model features no markov terms, i.e., terms that depend on more than one output, then the model is equivalent to a logistic regression. The example below shows this equivalence.
The optimization uses the exact gradient of the log-likelihood
(loglike_defm_grad()), and the variance-covariance matrix is the inverse
of minus the exact Hessian (loglike_defm_hess()) at the estimates. Since
the log-likelihood is concave, method = "newton" usually converges in a
few iterations, but, unlike "L-BFGS-B", it does not bound the
parameters, so it may drift when the MLE does not exist (e.g., under
separation).
Vega Yon, G. G., Pugh, M. J., & Valente, T. W. (2022). Discrete Exponential-Family Models for Multivariate Binary Outcomes (arXiv:2211.00627). arXiv. https://arxiv.org/abs/2211.00627
DEFM for objects of class DEFM and loglike_defm() for the
log-likelihood function of DEFMs.
#' Using Valente's SNS data
data(valentesnsList)
# Creating the DEFM object
logit_0 <- new_defm(
id = valentesnsList$id,
X = valentesnsList$X,
Y = valentesnsList$Y[,1,drop=FALSE],
order = 0
)
# Building the model
td_logit_intercept(logit_0)
td_logit_intercept(logit_0, covar = "Hispanic")
td_logit_intercept(
logit_0,
covar = "exposure_smoke"
)
td_logit_intercept(logit_0, covar = "Grades")
init_defm(logit_0) # Needs to be initialized
# Fitting the model
res_0 <- defm_mle(logit_0)
# Refitting the model using GLM
res_glm <- with(
valentesnsList,
glm(Y[,1] ~ X[,1] + X[,3] + X[,7], family = binomial())
)
# Comparing results
summary_table(res_0)
#> ====================================================================
#> Outcome-specific (logit) terms
#> --------------------------------------------------------------------
#> alcohol
#> --------------------------------------------------------------------
#> (Intercept) -0.273
#> (0.316)
#> Hispanic 1.335***
#> (0.133)
#> exposure_smoke 0.681***
#> (0.104)
#> Grades -0.284***
#> (0.070)
#> ====================================================================
#> N = 1730; AIC = 2058.186; BIC = 2080.010; Log-likelihood = -1025.093
#> *** p < 0.001; ** p < 0.01; * p < 0.05
summary(res_0)
#> Maximum likelihood estimation
#>
#> Call:
#> defm_mle(object = logit_0)
#>
#> Coefficients:
#> Estimate Std. Error
#> Logit intercept alcohol -0.2734291 0.31600926
#> Logit intercept alcohol x Hispanic 1.3351276 0.13263076
#> Logit intercept alcohol x exposure_smoke 0.6812039 0.10420829
#> Logit intercept alcohol x Grades -0.2840048 0.06962068
#>
#> -2 log L: 2050.186
#>
#> Tip: call summary_table() on the fitted model for a compact table by outcome (console, HTML, or Markdown).
summary(res_glm)
#>
#> Call:
#> glm(formula = Y[, 1] ~ X[, 1] + X[, 3] + X[, 7], family = binomial())
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.27338 0.31600 -0.865 0.387
#> X[, 1] 1.33511 0.13262 10.067 < 2e-16 ***
#> X[, 3] 0.68121 0.10421 6.537 6.27e-11 ***
#> X[, 7] -0.28401 0.06962 -4.080 4.51e-05 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 2294.3 on 1729 degrees of freedom
#> Residual deviance: 2050.2 on 1726 degrees of freedom
#> AIC: 2058.2
#>
#> Number of Fisher Scoring iterations: 3
#>
# Comparing the logodds
head(logodds(logit_0, par = coef(res_0), i = 0, j = 0))
#> [1] -2.3746571 -1.4094484 -1.5514508 -0.3295174 1.3168951 -0.3295174