Compute the maximum likelihood estimate of the precision matrix, given a known graphical structure (i.e., an adjacency matrix). This approach was originally described in "The Elements of Statistical Learning" (see pg. 631, Hastie et al. 2009) .
constrained(Sigma, adj) mle_known_graph(Sigma, adj)
| Sigma | Covariance matrix |
|---|---|
| adj | Adjacency matrix that encodes the constraints, where a zero indicates that element should be zero. |
A list containing the following:
Theta: Inverse of the covariance matrix (precision matrix)
Sigma: Covariance matrix.
wadj: Weighted adjacency matrix, corresponding to the partial correlation network.
The algorithm is written in c++, and should scale to high dimensions
nicely.
Note there are a variety of algorithms for this purpose. Simulation studies indicated that this approach is both accurate and computationally efficient (HFT therein, Emmert-Streib et al. 2019)
Emmert-Streib F, Tripathi S, Dehmer M (2019).
“Constrained covariance matrices with a biologically realistic structure: Comparison of methods for generating high-dimensional Gaussian graphical models.”
Frontiers in Applied Mathematics and Statistics, 5, 17.
Hastie T, Tibshirani R, Friedman J (2009).
The elements of statistical learning: data mining, inference, and prediction.
Springer Science \& Business Media.
# \donttest{ # data y <- ptsd # fit model fit <- ggmncv(cor(y), n = nrow(y), penalty = "lasso", progress = FALSE) # set negatives to zero (sign restriction) adj_new <- ifelse( fit$P <= 0, 0, 1) check_zeros <- TRUE # track trys iter <- 0 # iterate until all positive while(check_zeros){ iter <- iter + 1 fit_new <- constrained(cor(y), adj = adj_new) check_zeros <- any(fit_new$wadj < 0) adj_new <- ifelse( fit_new$wadj <= 0, 0, 1) } # } # alias # data y <- ptsd # nonreg (lambda = 0) fit <- ggmncv(cor(y), n = nrow(y), lambda = 0, progress = FALSE) # set values less than |0.1| to zero adj_new <- ifelse( abs(fit$P) <= 0.1, 0, 1) # mle given the graph mle_known_graph(cor(y), adj_new) #> $Theta #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] #> [1,] 2.040 -0.469 0.000 -0.625 -0.301 0.000 0.000 0.000 0.000 0.000 #> [2,] -0.469 2.308 -1.000 0.000 0.000 0.000 -0.157 0.000 0.000 0.000 #> [3,] 0.000 -1.000 2.421 -0.597 -0.333 0.000 0.000 0.000 0.335 0.000 #> [4,] -0.625 0.000 -0.597 2.237 -0.570 -0.141 -0.282 0.000 -0.311 0.000 #> [5,] -0.301 0.000 -0.333 -0.570 2.199 -0.269 -0.218 0.000 0.000 0.000 #> [6,] 0.000 0.000 0.000 -0.141 -0.269 1.338 -0.331 -0.162 0.205 0.000 #> [7,] 0.000 -0.157 0.000 -0.282 -0.218 -0.331 1.636 0.000 0.000 0.000 #> [8,] 0.000 0.000 0.000 0.000 0.000 -0.162 0.000 1.134 0.000 0.000 #> [9,] 0.000 0.000 0.335 -0.311 0.000 0.205 0.000 0.000 1.743 0.000 #> [10,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.729 #> [11,] -0.609 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.552 -1.058 #> [12,] 0.000 0.000 -0.425 0.223 0.000 0.000 -0.368 0.000 -0.458 0.153 #> [13,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 #> [14,] 0.000 0.000 -0.165 0.372 0.000 0.000 -0.127 0.000 -0.332 0.000 #> [15,] 0.000 0.000 0.000 0.000 -0.327 0.000 -0.236 0.000 -0.174 0.000 #> [16,] 0.000 -0.209 0.000 0.000 0.000 0.000 0.000 -0.193 0.000 -0.138 #> [17,] 0.000 0.000 -0.212 0.000 0.000 -0.135 0.000 0.000 -0.199 0.000 #> [18,] 0.141 -0.244 0.000 0.000 -0.358 0.000 0.000 0.000 0.000 0.059 #> [19,] 0.000 0.000 0.000 0.000 0.000 0.000 0.272 -0.219 0.000 -0.222 #> [20,] 0.000 -0.369 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 #> [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] #> [1,] -0.609 0.000 0.000 0.000 0.000 0.000 0.000 0.141 0.000 0.000 #> [2,] 0.000 0.000 0.000 0.000 0.000 -0.209 0.000 -0.244 0.000 -0.369 #> [3,] 0.000 -0.425 0.000 -0.165 0.000 0.000 -0.212 0.000 0.000 0.000 #> [4,] 0.000 0.223 0.000 0.372 0.000 0.000 0.000 0.000 0.000 0.000 #> [5,] 0.000 0.000 0.000 0.000 -0.327 0.000 0.000 -0.358 0.000 0.000 #> [6,] 0.000 0.000 0.000 0.000 0.000 0.000 -0.135 0.000 0.000 0.000 #> [7,] 0.000 -0.368 0.000 -0.127 -0.236 0.000 0.000 0.000 0.272 0.000 #> [8,] 0.000 0.000 0.000 0.000 0.000 -0.193 0.000 0.000 -0.219 0.000 #> [9,] -0.552 -0.458 0.000 -0.332 -0.174 0.000 -0.199 0.000 0.000 0.000 #> [10,] -1.058 0.153 0.000 0.000 0.000 -0.138 0.000 0.059 -0.222 0.000 #> [11,] 2.460 0.000 0.000 0.000 -0.391 0.000 0.000 0.000 0.000 0.000 #> [12,] 0.000 2.009 -0.383 0.000 0.000 0.000 0.000 0.000 -0.668 0.000 #> [13,] 0.000 -0.383 2.352 -1.084 0.000 0.000 0.000 -0.373 -0.333 -0.282 #> [14,] 0.000 0.000 -1.084 1.985 -0.331 0.000 0.000 0.147 0.000 0.000 #> [15,] -0.391 0.000 0.000 -0.331 1.812 -0.399 0.000 0.000 0.000 0.000 #> [16,] 0.000 0.000 0.000 0.000 -0.399 1.423 0.000 -0.241 0.000 0.000 #> [17,] 0.000 0.000 0.000 0.000 0.000 0.000 1.450 -0.571 0.000 0.000 #> [18,] 0.000 0.000 -0.373 0.147 0.000 -0.241 -0.571 1.829 -0.278 0.000 #> [19,] 0.000 -0.668 -0.333 0.000 0.000 0.000 0.000 -0.278 1.892 -0.370 #> [20,] 0.000 0.000 -0.282 0.000 0.000 0.000 0.000 0.000 -0.370 1.431 #> #> $Sigma #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] #> [1,] 1.012 0.525 0.471 0.595 0.530 0.250 0.365 0.125 0.322 0.343 0.531 0.263 #> [2,] 0.525 1.007 0.679 0.478 0.496 0.264 0.403 0.164 0.257 0.236 0.348 0.379 #> [3,] 0.471 0.679 1.002 0.555 0.555 0.295 0.408 0.155 0.203 0.200 0.304 0.423 #> [4,] 0.595 0.478 0.555 1.003 0.610 0.337 0.443 0.123 0.280 0.236 0.366 0.217 #> [5,] 0.530 0.496 0.555 0.610 1.005 0.395 0.479 0.161 0.283 0.243 0.373 0.323 #> [6,] 0.250 0.264 0.295 0.337 0.395 1.000 0.392 0.193 0.036 0.094 0.144 0.160 #> [7,] 0.365 0.403 0.408 0.443 0.479 0.392 1.005 0.127 0.256 0.172 0.287 0.364 #> [8,] 0.125 0.164 0.155 0.123 0.161 0.193 0.127 1.000 0.108 0.109 0.127 0.174 #> [9,] 0.322 0.257 0.203 0.280 0.283 0.036 0.256 0.108 1.009 0.327 0.513 0.440 #> [10,] 0.343 0.236 0.200 0.236 0.243 0.094 0.172 0.109 0.327 1.006 0.644 0.169 #> [11,] 0.531 0.348 0.304 0.366 0.373 0.144 0.287 0.127 0.513 0.644 1.008 0.291 #> [12,] 0.263 0.379 0.423 0.217 0.323 0.160 0.364 0.174 0.440 0.169 0.291 1.004 #> [13,] 0.222 0.343 0.330 0.153 0.285 0.129 0.261 0.163 0.370 0.197 0.280 0.533 #> [14,] 0.170 0.252 0.249 0.052 0.210 0.082 0.248 0.119 0.403 0.188 0.278 0.412 #> [15,] 0.381 0.369 0.351 0.339 0.461 0.209 0.405 0.159 0.420 0.330 0.491 0.342 #> [16,] 0.279 0.370 0.307 0.256 0.322 0.170 0.254 0.244 0.247 0.268 0.310 0.254 #> [17,] 0.239 0.330 0.361 0.269 0.330 0.228 0.240 0.124 0.276 0.148 0.225 0.283 #> [18,] 0.261 0.435 0.404 0.299 0.440 0.221 0.277 0.173 0.264 0.165 0.246 0.370 #> [19,] 0.218 0.328 0.312 0.164 0.256 0.113 0.143 0.254 0.317 0.257 0.280 0.558 #> [20,] 0.236 0.412 0.321 0.196 0.251 0.123 0.192 0.140 0.221 0.166 0.217 0.347 #> [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] #> [1,] 0.222 0.170 0.381 0.279 0.239 0.261 0.218 0.236 #> [2,] 0.343 0.252 0.369 0.370 0.330 0.435 0.328 0.412 #> [3,] 0.330 0.249 0.351 0.307 0.361 0.404 0.312 0.321 #> [4,] 0.153 0.052 0.339 0.256 0.269 0.299 0.164 0.196 #> [5,] 0.285 0.210 0.461 0.322 0.330 0.440 0.256 0.251 #> [6,] 0.129 0.082 0.209 0.170 0.228 0.221 0.113 0.123 #> [7,] 0.261 0.248 0.405 0.254 0.240 0.277 0.143 0.192 #> [8,] 0.163 0.119 0.159 0.244 0.124 0.173 0.254 0.140 #> [9,] 0.370 0.403 0.420 0.247 0.276 0.264 0.317 0.221 #> [10,] 0.197 0.188 0.330 0.268 0.148 0.165 0.257 0.166 #> [11,] 0.280 0.278 0.491 0.310 0.225 0.246 0.280 0.217 #> [12,] 0.533 0.412 0.342 0.254 0.283 0.370 0.558 0.347 #> [13,] 1.005 0.655 0.359 0.265 0.281 0.431 0.515 0.420 #> [14,] 0.655 1.004 0.402 0.227 0.199 0.254 0.354 0.286 #> [15,] 0.359 0.402 1.005 0.444 0.255 0.320 0.277 0.238 #> [16,] 0.265 0.227 0.444 1.002 0.237 0.363 0.255 0.214 #> [17,] 0.281 0.199 0.255 0.237 1.002 0.511 0.262 0.208 #> [18,] 0.431 0.254 0.320 0.363 0.511 1.003 0.413 0.304 #> [19,] 0.515 0.354 0.277 0.255 0.262 0.413 1.003 0.445 #> [20,] 0.420 0.286 0.238 0.214 0.208 0.304 0.445 1.003 #> #> $wadj #> [,1] [,2] [,3] [,4] [,5] [,6] #> [1,] 0.00000000 0.21614234 0.00000000 0.29257145 0.1421146 0.00000000 #> [2,] 0.21614234 0.00000000 0.42304318 0.00000000 0.0000000 0.00000000 #> [3,] 0.00000000 0.42304318 0.00000000 0.25653342 0.1443226 0.00000000 #> [4,] 0.29257145 0.00000000 0.25653342 0.00000000 0.2569977 0.08150009 #> [5,] 0.14211463 0.00000000 0.14432256 0.25699771 0.0000000 0.15682367 #> [6,] 0.00000000 0.00000000 0.00000000 0.08150009 0.1568237 0.00000000 #> [7,] 0.00000000 0.08079604 0.00000000 0.14740916 0.1149350 0.22372179 #> [8,] 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000 0.13151656 #> [9,] 0.00000000 0.00000000 -0.16307925 0.15749933 0.0000000 -0.13423851 #> [10,] 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [11,] 0.27185333 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [12,] 0.00000000 0.00000000 0.19270881 -0.10519180 0.0000000 0.00000000 #> [13,] 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [14,] 0.00000000 0.00000000 0.07526729 -0.17653454 0.0000000 0.00000000 #> [15,] 0.00000000 0.00000000 0.00000000 0.00000000 0.1638159 0.00000000 #> [16,] 0.00000000 0.11532566 0.00000000 0.00000000 0.0000000 0.00000000 #> [17,] 0.00000000 0.00000000 0.11314998 0.00000000 0.0000000 0.09692185 #> [18,] -0.07299568 0.11875860 0.00000000 0.00000000 0.1785104 0.00000000 #> [19,] 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [20,] 0.00000000 0.20304329 0.00000000 0.00000000 0.0000000 0.00000000 #> [,7] [,8] [,9] [,10] [,11] [,12] #> [1,] 0.00000000 0.0000000 0.00000000 0.00000000 0.2718533 0.00000000 #> [2,] 0.08079604 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [3,] 0.00000000 0.0000000 -0.16307925 0.00000000 0.0000000 0.19270881 #> [4,] 0.14740916 0.0000000 0.15749933 0.00000000 0.0000000 -0.10519180 #> [5,] 0.11493498 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [6,] 0.22372179 0.1315166 -0.13423851 0.00000000 0.0000000 0.00000000 #> [7,] 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000 0.20298605 #> [8,] 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [9,] 0.00000000 0.0000000 0.00000000 0.00000000 0.2665771 0.24475232 #> [10,] 0.00000000 0.0000000 0.00000000 0.00000000 0.5130038 -0.08209259 #> [11,] 0.00000000 0.0000000 0.26657709 0.51300384 0.0000000 0.00000000 #> [12,] 0.20298605 0.0000000 0.24475232 -0.08209259 0.0000000 0.00000000 #> [13,] 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000 0.17619361 #> [14,] 0.07047447 0.0000000 0.17848805 0.00000000 0.0000000 0.00000000 #> [15,] 0.13706964 0.0000000 0.09790874 0.00000000 0.1851954 0.00000000 #> [16,] 0.00000000 0.1519317 0.00000000 0.08797902 0.0000000 0.00000000 #> [17,] 0.00000000 0.0000000 0.12517580 0.00000000 0.0000000 0.00000000 #> [18,] 0.00000000 0.0000000 0.00000000 -0.03317781 0.0000000 0.00000000 #> [19,] -0.15460256 0.1495123 0.00000000 0.12274249 0.0000000 0.34263043 #> [20,] 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 #> [,13] [,14] [,15] [,16] [,17] [,18] #> [1,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 -0.07299568 #> [2,] 0.0000000 0.00000000 0.00000000 0.11532566 0.00000000 0.11875860 #> [3,] 0.0000000 0.07526729 0.00000000 0.00000000 0.11314998 0.00000000 #> [4,] 0.0000000 -0.17653454 0.00000000 0.00000000 0.00000000 0.00000000 #> [5,] 0.0000000 0.00000000 0.16381590 0.00000000 0.00000000 0.17851041 #> [6,] 0.0000000 0.00000000 0.00000000 0.00000000 0.09692185 0.00000000 #> [7,] 0.0000000 0.07047447 0.13706964 0.00000000 0.00000000 0.00000000 #> [8,] 0.0000000 0.00000000 0.00000000 0.15193168 0.00000000 0.00000000 #> [9,] 0.0000000 0.17848805 0.09790874 0.00000000 0.12517580 0.00000000 #> [10,] 0.0000000 0.00000000 0.00000000 0.08797902 0.00000000 -0.03317781 #> [11,] 0.0000000 0.00000000 0.18519536 0.00000000 0.00000000 0.00000000 #> [12,] 0.1761936 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 #> [13,] 0.0000000 0.50168414 0.00000000 0.00000000 0.00000000 0.17983877 #> [14,] 0.5016841 0.00000000 0.17452942 0.00000000 0.00000000 -0.07714899 #> [15,] 0.0000000 0.17452942 0.00000000 0.24847995 0.00000000 0.00000000 #> [16,] 0.0000000 0.00000000 0.24847995 0.00000000 0.00000000 0.14938526 #> [17,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.35062692 #> [18,] 0.1798388 -0.07714899 0.00000000 0.14938526 0.35062692 0.00000000 #> [19,] 0.1578574 0.00000000 0.00000000 0.00000000 0.00000000 0.14944362 #> [20,] 0.1537130 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 #> [,19] [,20] #> [1,] 0.0000000 0.0000000 #> [2,] 0.0000000 0.2030433 #> [3,] 0.0000000 0.0000000 #> [4,] 0.0000000 0.0000000 #> [5,] 0.0000000 0.0000000 #> [6,] 0.0000000 0.0000000 #> [7,] -0.1546026 0.0000000 #> [8,] 0.1495123 0.0000000 #> [9,] 0.0000000 0.0000000 #> [10,] 0.1227425 0.0000000 #> [11,] 0.0000000 0.0000000 #> [12,] 0.3426304 0.0000000 #> [13,] 0.1578574 0.1537130 #> [14,] 0.0000000 0.0000000 #> [15,] 0.0000000 0.0000000 #> [16,] 0.0000000 0.0000000 #> [17,] 0.0000000 0.0000000 #> [18,] 0.1494436 0.0000000 #> [19,] 0.0000000 0.2248647 #> [20,] 0.2248647 0.0000000 #>