The generic function ‘mle’ fits the graphical model selected by glasso, mglasso or cglasso.
mle(object, …)# S3 method for glasso
mle(object, …, maxit = 1.0e+04, thr = 1.0e-04, trace = 0L)
# S3 method for mglasso
mle(object, …, maxit_em = 1.0e+03, thr_em = 1.0e-4, maxit_bcd = 1.0e+4,
thr_bcd = 1.0e-4, trace = 0L)
# S3 method for cglasso
mle(object, …, maxit_em = 1.0e+03, thr_em = 1.0e-4, maxit_bcd = 1.0e+4,
thr_bcd = 1.0e-4, trace = 0L)
an object of class ‘glasso’, ‘mglasso’ or ‘cglasso’.
maximum number of iterations of the glasso algorithm. Default is 1.0E+4.
threshold for the convergence of the glasso algorithm. Default is 1.0E-4.
maximum number of iterations of the EM algorithm. Default is 1.0E+03.
threshold for the convergence of the EM algorithm. Default is 1.0E-4.
maximum number of iterations of the glasso algorithm. Default is 1.0E+4.
threshold for the convergence of the glasso algorithm. Default is 1.0E-4.
integer for printing out information as iterations proceed: trace = 0 no information is printed out on video; trace = 1 basic information is printed out on video; trace = 2 detailed information is printed out on video.
additional argument added for backward compatibility with the generic function ‘mle’.
If ‘object’ has class ‘glasso’, then ‘mle’ returns and object with S3 class ‘ggm’, which inherits the class ‘glasso’. An object with class ‘ggm’ is a list containing the following components:
the call that produced this object.
the matrix used to compute the covariance matrix.
the covariance matrix used to fit the ggm model.
the number of fitted ggm model.
the \(p\)-dimensional vector reporting the values of the tuning parameter used to fit the glasso model.
the maximum number of iterations of the glasso algorithm.
the threshold for the convergence of the glasso algorithm.
an array of dimension \((p\times p\times\texttt{nrho})\). Sgm[, , k] is the estimate of the covariance matrix of the \(k\)th ggm model.
an array of dimension \((p\times p\times\texttt{nrho})\). Tht[, , k] is the estimate of the precision matrix of the \(k\)th ggm model.
an array of dimension \((p\times p\times\texttt{nrho})\). Adj[, , k] is the adjacency matrix associated to Tht[, , k], i.e. Adj[i, j, k] \(= 1\) iff Tht[i, j, k] \(\neq 0\) and 0 otherwise.
the \(p\)-dimensional vector reporting the number of non-zero partial correlation coefficients.
the \(p\)-dimensional vector reporting the values of the measure \(R^2\). See section Details, in glasso.
the \(p\)-dimensional vector reporting the number of connected components (for internal purposes only).
the \((p\times\texttt{nrho})\)-dimensional matrix encoding the connected components (for internal purposes only).
the \((p\times\texttt{nrho})\)-dimensional matrix reporting the number of vertices per connected component (for internal purposes only).
the \(p\)-dimensional vector reporting the number of iterations.
a description of the error that has occurred.
the integer used for printing out information.
If object has class mglasso, then mle returns and object with S3 class mggm, which inherits the class mglasso. An object with class mggm is a list containing the following components:
the call that produced this object.
the object with S3 class ‘datacggm’ used to fit the cggm model.
the number of fitted cggm model.
the \(p\)-dimensional vector reporting the values of the tuning parameter used to fit the cglasso model.
maximum number of iterations of the EM algorithm.
threshold for the convergence of the EM algorithm.
maximum number of iterations of the glasso algorithm.
threshold for the convergence of the glasso algorithm.
an array of dimension \(n\times p\times\texttt{nrho}\). Xipt[, , k] is the matrix where the censored vaules are replaced using the conditional expected vaules computed in the E-step of the algorithm describted in section Details.
an array of dimension \(p\times p\times\texttt{nrho}\). S[, , k] is the matrix \(S'\) used to fit the cggm model (see the section Details).
a matrix of dimension \(p\times\texttt{nrho}\). The \(k\)th column is the estimate of the expected values of the \(k\)th cggm model.
an array of dimension \(p\times p\times\texttt{nrho}\). Sgm[, , k] is the estimate of the covariance matrix of the \(k\)th cggm model.
an array of dimension \(p\times p\times\texttt{nrho}\). Tht[, , k] is the estimate of the precision matrix of the \(k\)th cggm model.
an array of dimension \(p\times p\times\texttt{nrho}\). Adj[, , k] is the adjacency matrix associated to Tht[, , k], i.e. Adj[i, j, k] \(= 1\) iff Tht[i, j, k] \(\neq 0\) and 0 otherwise.
the \(p\)-dimensional vector reporting the number of non-zero partial correlation coefficients.
the \(p\)-dimensional vector reporting the values of the measure \(R^2\). See section Details, in cglasso.
the \(p\)-dimensional vector reporting the number of connected components (for internal purposes only).
the \((p\times\texttt{nrho})\)-dimensional matrix encoding the connected components (for internal purposes only).
the \((p\times\texttt{nrho})\)-dimensional matrix reporting the number of vertices per connected component (for internal purposes only).
the \((\texttt{nrho}\times 2)\)-dimensional matrix reporting the number of iterations.
a description of the error that has occurred.
the name of the Fortran subroutine where the error has occurred (for internal debug only).
the integer used for printing out information.
If object has class cglasso, then mle returns and object with S3 class cggm, which inherits the class cglasso. An object with class cggm is a list containing the following components:
the call that produced this object.
the object with S3 class ‘datacggm’ used to fit the cggm model.
the number of fitted cggm model.
the \(p\)-dimensional vector reporting the values of the tuning parameter used to fit the cglasso model.
maximum number of iterations of the EM algorithm.
threshold for the convergence of the EM algorithm.
maximum number of iterations of the glasso algorithm.
threshold for the convergence of the glasso algorithm.
an array of dimension \(n\times p\times\texttt{nrho}\). Xipt[, , k] is the matrix where the censored vaules are replaced using the conditional expected vaules computed in the E-step of the algorithm describted in section Details.
an array of dimension \(p\times p\times\texttt{nrho}\). S[, , k] is the matrix \(S'\) used to fit the cggm model (see the section Details).
a matrix of dimension \(p\times\texttt{nrho}\). The \(k\)th column is the estimate of the expected values of the \(k\)th cggm model.
an array of dimension \(p\times p\times\texttt{nrho}\). Sgm[, , k] is the estimate of the covariance matrix of the \(k\)th cggm model.
an array of dimension \(p\times p\times\texttt{nrho}\). Tht[, , k] is the estimate of the precision matrix of the \(k\)th cggm model.
an array of dimension \(p\times p\times\texttt{nrho}\). Adj[, , k] is the adjacency matrix associated to Tht[, , k], i.e. Adj[i, j, k] \(= 1\) iff Tht[i, j, k] \(\neq 0\) and 0 otherwise.
the \(p\)-dimensional vector reporting the number of non-zero partial correlation coefficients.
the \(p\)-dimensional vector reporting the values of the measure \(R^2\). See section Details, in cglasso.
the \(p\)-dimensional vector reporting the number of connected components (for internal purposes only).
the \((p\times\texttt{nrho})\)-dimensional matrix encoding the connected components (for internal purposes only).
the \((p\times\texttt{nrho})\)-dimensional matrix reporting the number of vertices per connected component (for internal purposes only).
the \((\texttt{nrho}\times 2)\)-dimensional matrix reporting the number of iterations.
a description of the error that has occurred.
the name of the Fortran subroutine where the error has occurred (for internal debug only).
the integer used for printing out information.
The generic function ‘mle’ computes the maximum likelihood estimates of the graphical model selected by the function glasso, mglasso or cglasso.
If ‘object’ has class ‘glasso’, the method function ‘mle.glasso’ computes the maximum likelihood estimates of the parameters of the Gaussian graphical models (ggm) associated to the sequence of glasso estimates. Formally, for a given value of the tuning parameter let \(\hat{\Theta}^\rho\) be the glasso estimate of the precision matrix, then ‘mle.glasso’ solves the following maximization problem: $$\max_{\bar{\Theta}}\log det\bar{\Theta} - tr\{S\bar{\Theta}\},$$ where \(\bar{\Theta} = \{\bar{\theta}_{hk}\}\) and \(\bar{\theta}_{hk} = 0\) if \(\hat{\theta}^{\rho}_{hk} = 0\) otherwise it is estimated.
If ‘object’ has class ‘mglasso’, the method function ‘mle.mglasso’ computes the maximum likelihood estimates of the parameters of the Gaussian graphical models with missing-at-random data (mggm) associated to the sequence of mglasso estimates. Formally, for a given value of the tuning parameter let \(\hat{\Theta}^\rho\) be the mglasso estimate of the precision matrix, then ‘mle.mglasso’ computes the maximum likelihood estimate by the following EM-like algorithm:
| Step | Description |
| 1. | let \(\hat{\Theta}^\rho\) be the mglasso estimate of the precision matrix; |
| 2. | E-step |
| use the moments of the conditional normal distribution to impute the missing values; | |
| 3. | M-step |
| let \(X'\) the completed matrix and \(S'\) the corresponding empirical variance matrix, then: | |
| let \(\hat{\mu}_h = \sum_{i = 1}^n x'_{ih} / n\) | |
| estimate \(\bar{\Theta}\) maximixing the function: \(\log det\bar{\Theta} - tr\{S'\bar{\Theta}\},\) where \(\bar{\theta}_{hk} = 0\) if \(\hat{\theta}^{\rho}_{hk} = 0\) | |
| otherwise it is estimated; |
If ‘object’ has class ‘cglasso’, the method function ‘mle.cglasso’ computes the maximum likelihood estimates of the parameters of the censored Gaussian graphical models (cggm) associated to the sequence of cglasso estimates. Formally, for a given value of the tuning parameter let \(\hat{\Theta}^\rho\) be the cglasso estimate of the precision matrix, then ‘mle.cglasso’ computes the maximum likelihood estimate by the following EM-like algorithm:
| Step | Description |
| 1. | let \(\hat{\Theta}^\rho\) be the cglasso estimate of the precision matrix; |
| 2. | E-step |
| use the moments of the truncated normal distribution to compute the current estimates of the | |
| marginal means, denoted by \(\bar{x}'\), and to complete the empirical covariance matrix \(S'\); | |
| 3. | M-step |
| let \(\hat{\mu} = \bar{x}'\); | |
| estimate \(\bar{\Theta}\) maximixing the \(Q\)-function: \(\log det\bar{\Theta} - tr\{S'\bar{\Theta}\},\) where \(\bar{\theta}_{hk} = 0\) if \(\hat{\theta}^{\rho}_{hk} = 0\) | |
| otherwise it is estimated; |
Augugliaro, L., Abbruzzo, A., and Vinciotti, V. (2018). \(\ell_1\)-Penalized gaussian graphical model. Biostatistics (to appear).
Stadler N., and Buhlmann P. (2012) <DOI:10.1007/s11222-010-9219-7>. Missing values: sparse inverse covariance estimation and an extension to sparse regression. Statistics and Computing 22, 219--235.
glasso, mglasso, cglasso, to_graph, and the method functions summary, coef, plot, aic, bic and ebic.
# NOT RUN {
library("cglasso")
set.seed(123)
#################
# cglasso model #
#################
n <- 100L
p <- 5L
mu <- rep.int(0L, times = p)
X <- rdatacggm(n = n, mu = mu, probr = 0.05)
out <- cglasso(X = X)
out_mle <- mle(out)
out_mle
inherits(out_mle, "cglasso")
class(out_mle)
# inheriting method functions from 'cglasso': some examples
coef(out_mle, nrho = 10L, print.info = TRUE)
to_graph(out_mle, nrho = 10L, weighted = TRUE)
out_aic <- aic(out_mle)
out_aic
plot(out_mle, typeplot = "graph", gof = out_aic)
#################
# mglasso model #
#################
R <- event(X)
X <- as.matrix(X)
X[R == 1L] <- NA
out <- mglasso(X = X)
out_mle <- mle(out)
out_mle
inherits(out_mle, "mglasso")
class(out_mle)
# inheriting method functions from 'mglasso': some examples
coef(out_mle, nrho = 10L, print.info = TRUE)
to_graph(out_mle, nrho = 10L, weighted = TRUE)
out_aic <- aic(out_mle)
out_aic
plot(out_mle, typeplot = "graph", gof = out_aic)
################
# glasso model #
################
X <- MASS::mvrnorm(n = n, mu = mu, Sigma = diag(p))
out <- glasso(X)
out_mle <- mle(out)
out_mle
inherits(out_mle, "glasso")
class(out_mle)
# inheriting method functions from 'glasso': some examples
coef(out_mle, nrho = 10L, print.info = TRUE)
to_graph(out_mle, nrho = 10L, weighted = TRUE)
out_aic <- aic(out_mle)
out_aic
plot(out_mle, typeplot = "graph", gof = out_aic)
# }
Run the code above in your browser using DataLab