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)

Arguments

Sigma

Covariance matrix

adj

Adjacency matrix that encodes the constraints, where a zero indicates that element should be zero.

Value

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.

Note

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)

References

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.

Examples

# \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
#>