Performs fused ridge estimation of multiple precision matrices in cases where multiple classes of data is present for given a penalty matrix.

ridgeP.fused(
  Slist,
  ns,
  Tlist = default.target.fused(Slist, ns),
  lambda,
  Plist,
  maxit = 100L,
  verbose = TRUE,
  relative = TRUE,
  eps = sqrt(.Machine$double.eps)
)

Arguments

Slist

A list of length \(G\) of covariance matrices, i.e. square, symmetric numeric matrices of the same size. The \(g\)th matrix should correspond to the \(g\)th class.

ns

A numeric vector of sample sizes on which the matrices in Slist are based. I.e. ns[g] correspond to Slist[[g]].

Tlist

A list of length \(G\) of numeric p.d. target matrices corresponding to the matrices in Slist. If not supplied, the default is given by default.target.

lambda

The \(G\) by \(G\) penalty matrix. That is, a symmetric, non-negative numeric matrix of size \(G\) times \(G\) giving the class- and pair-specific penalties. The diagonal entries are the class specific ridge penalties. I.e. lambda[i, i] is the ridge penalty for class \(i\). The off-diagonal entries are the pair-specific fusion penalties. I.e. lambda[i, j] is the fusion penalty applied on the pair of classes \(i\) and \(j\).

Alternatively, can be supplied as a numeric of length 1 or 2. If a single number, a diagonal penalty with lambda in the diagonal is used If supplied as a numeric vector of two numbers, the first is used as a common ridge penalty and the second as a common fusion penalty.

Plist

An optional list of initial precision matrices for the fused ridge algorithm the same size as Slist. Can be omitted. Default is the nonfused ridge precision estimate using the pooled covariance matrix corresponding to setting all fusion penalties to zero.

maxit

A single integer giving the maximum number of allowed iterations. Can be set to Inf. If maxit is hit, a warning is given.

verbose

logical. Set to TRUE for extra output.

relative

logical indicating if the convergence criterion should be on a relative scale.

eps

A single positive numeric giving the convergence threshold.

Value

Returns a list as Slist with precision estimates of the corresponding classes.

Details

Performs a coordinate ascent to find the maximum likelihood of the fused likelihood problem for a given ridge penalty \(lambda\) and fused penalty matrix \(Lambda_f\).

Note

For extreme fusion penalties in lambda the algorithm is quite sensitive to the initial values given in Plist.

References

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2020). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. Journal of Machine Learning Research, 21(26): 1-52.

See also

default.penalty
ridgeP for the regular ridge estimate

Author

Anders Ellern Bilgrau, Carel F.W. Peeters <carel.peeters@wur.nl>, Wessel N. van Wieringen

Examples

# Create some (not at all high-dimensional) data on three classes
p <- 5  # Dimension
ns <- c(4, 6, 8)  # Sample sizes (K = 3 classes)
Slist <- createS(ns, p = p)
str(Slist, max.level = 2)  # The structure of Slist
#> List of 3
#>  $ class1: num [1:5, 1:5] 1.512 1.002 -0.812 0.276 0.261 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ class2: num [1:5, 1:5] 0.27 0.249 0.109 0.145 0.141 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ class3: num [1:5, 1:5] 1.956 1.057 -0.337 -0.692 0.902 ...
#>   ..- attr(*, "dimnames")=List of 2

#
# Estimate the precisions (using the complete penalty graph)
#

res1 <- ridgeP.fused(Slist, ns, lambda = c(1.3, 2.1))
#> i = 1   | max diff = 1.4002010950e-01
#> i = 2   | max diff = 1.0702011953e-01
#> i = 3   | max diff = 2.3734578607e-02
#> i = 4   | max diff = 4.0420482876e-03
#> i = 5   | max diff = 7.5721388746e-04
#> i = 6   | max diff = 1.4092417828e-04
#> i = 7   | max diff = 2.7060429387e-05
#> i = 8   | max diff = 5.4094720853e-06
#> i = 9   | max diff = 1.1344020719e-06
#> i = 10  | max diff = 2.5018649527e-07
#> i = 11  | max diff = 5.7933222898e-08
#> i = 12  | max diff = 1.4009241402e-08
#> Converged in 12 iterations, max diff < 1.49e-08.
print(res1)
#> $class1
#>             A          B           C           D           E
#> A  1.03597692 -0.7444659  0.10487143  0.06337964 -0.51367511
#> B -0.74446592  1.2427458  0.21726099  0.20546627  0.16747427
#> C  0.10487143  0.2172610  0.88468511 -0.22615653  0.07943883
#> D  0.06337964  0.2054663 -0.22615653  2.10615121 -0.11799362
#> E -0.51367511  0.1674743  0.07943883 -0.11799362  1.21975425
#> 
#> $class2
#>              A            B            C           D           E
#> A 794.42776683  -1.36095956   0.09142925   0.2561648  -1.0091060
#> B  -1.36095956 794.45413280  -0.03395939   0.3564319   0.1359921
#> C   0.09142925  -0.03395939 794.48929468  -0.3481282   0.7274477
#> D   0.25616483   0.35643188  -0.34812820 796.1015732  -0.1389002
#> E  -1.00910603   0.13599215   0.72744769  -0.1389002 794.8779662
#> 
#> $class3
#>            A          B          C          D          E
#> A  4.9327303 -2.1078911  0.4455113  1.0218273 -1.7272875
#> B -2.1078911  5.6947508  0.9728766  0.8425952 -0.6157010
#> C  0.4455113  0.9728766  5.5829703 -0.7185520  0.8790972
#> D  1.0218273  0.8425952 -0.7185520  8.0517058  0.4600716
#> E -1.7272875 -0.6157010  0.8790972  0.4600716  6.8532396
#> 

# The same using the penalty matrix (the diagnal is ignored)
mylambda  <- matrix(c(1.3, 2.1, 2.1,
                      2.1, 1.3, 2.1,
                      2.1, 2.1, 1.3), 3, 3, byrow = TRUE)
res2 <- ridgeP.fused(Slist, ns, lambda = mylambda)
#> i = 1   | max diff = 1.4002010950e-01
#> i = 2   | max diff = 1.0702011953e-01
#> i = 3   | max diff = 2.3734578607e-02
#> i = 4   | max diff = 4.0420482876e-03
#> i = 5   | max diff = 7.5721388746e-04
#> i = 6   | max diff = 1.4092417828e-04
#> i = 7   | max diff = 2.7060429387e-05
#> i = 8   | max diff = 5.4094720853e-06
#> i = 9   | max diff = 1.1344020719e-06
#> i = 10  | max diff = 2.5018649527e-07
#> i = 11  | max diff = 5.7933222898e-08
#> i = 12  | max diff = 1.4009241402e-08
#> Converged in 12 iterations, max diff < 1.49e-08.
stopifnot(all.equal(res1, res2))


#
# Estimate the precisions (using a non-complete penalty graph)
#

# Say we only want to shrink pairs (1,2) and (2,3) and not (1,3)
mylambda[1,3] <- mylambda[3,1] <- 0
print(mylambda)
#>      [,1] [,2] [,3]
#> [1,]  1.3  2.1  0.0
#> [2,]  2.1  1.3  2.1
#> [3,]  0.0  2.1  1.3
res3 <- ridgeP.fused(Slist, ns, lambda = mylambda)
#> i = 1   | max diff = 2.8142966298e-01
#> i = 2   | max diff = 2.3624155485e-02
#> i = 3   | max diff = 3.8485663292e-02
#> i = 4   | max diff = 3.5156955750e-03
#> i = 5   | max diff = 4.5185832791e-04
#> i = 6   | max diff = 6.1386422646e-05
#> i = 7   | max diff = 8.5434453489e-06
#> i = 8   | max diff = 1.2139880452e-06
#> i = 9   | max diff = 1.7613303766e-07
#> i = 10  | max diff = 2.6104072928e-08
#> i = 11  | max diff = 3.9537010608e-09
#> Converged in 11 iterations, max diff < 1.49e-08.
# which similar to, but not the same as res1 and res2.


#
# Using other custom target matrices
#

# Construct a custom target list
myTlist <- list(diag(p), matrix(1, p, p), matrix(0, p, p))
res4 <- ridgeP.fused(Slist, ns, Tlist = myTlist, lambda = c(1.3, 2.1))
#> i = 1   | max diff = 3.1257889996e+01
#> i = 2   | max diff = 3.4409827451e+00
#> i = 3   | max diff = 6.8012210972e-01
#> i = 4   | max diff = 6.7082874983e-01
#> i = 5   | max diff = 7.6327877176e-01
#> i = 6   | max diff = 7.8392918969e-01
#> i = 7   | max diff = 7.3694397431e-01
#> i = 8   | max diff = 4.9799334039e-01
#> i = 9   | max diff = 2.5046310639e-01
#> i = 10  | max diff = 6.7723607315e-02
#> i = 11  | max diff = 1.2022395832e-02
#> i = 12  | max diff = 1.7562206922e-03
#> i = 13  | max diff = 2.3338896503e-04
#> i = 14  | max diff = 2.9862421611e-05
#> i = 15  | max diff = 3.7716355751e-06
#> i = 16  | max diff = 4.7469274413e-07
#> i = 17  | max diff = 5.9738296332e-08
#> i = 18  | max diff = 7.5260405998e-09
#> Converged in 18 iterations, max diff < 1.49e-08.
print(res4)
#> $class1
#>            A            B          C           D            E
#> A  1.4001801 -0.768203923  0.2943821 -0.16611466 -0.378855549
#> B -0.7682039  1.747537878  0.2492091 -0.02477384  0.005000194
#> C  0.2943821  0.249209095  1.6251428 -0.16973280 -0.196420943
#> D -0.1661147 -0.024773838 -0.1697328  2.36227912 -0.253810347
#> E -0.3788555  0.005000194 -0.1964209 -0.25381035  1.752663100
#> 
#> $class2
#>           A         B         C         D         E
#> A 1.9040215 0.3079378 0.7597049 0.7138363 0.4966479
#> B 0.3079378 1.8671581 0.6353090 0.8075746 0.9546634
#> C 0.7597049 0.6353090 1.9737974 0.7540297 1.0124113
#> D 0.7138363 0.8075746 0.7540297 2.5020341 0.6444428
#> E 0.4966479 0.9546634 1.0124113 0.6444428 2.0137004
#> 
#> $class3
#>             A           B           C           D           E
#> A  0.86617971 -0.42853113 -0.05883382  0.20106627 -0.42673730
#> B -0.42853113  0.99650400  0.20034523  0.11603158 -0.07879847
#> C -0.05883382  0.20034523  0.80337872 -0.25067769  0.12606867
#> D  0.20106627  0.11603158 -0.25067769  1.66211100  0.05613467
#> E -0.42673730 -0.07879847  0.12606867  0.05613467  1.31782482
#> 

# Alternative, see ?default.target.fused
myTlist2 <- default.target.fused(Slist, ns, type = "Null")  # For the null target
res5 <- ridgeP.fused(Slist, ns, Tlist = myTlist2, lambda = c(1.3, 2.1))
#> i = 1   | max diff = 3.1789287061e+01
#> i = 2   | max diff = 3.5058518047e+00
#> i = 3   | max diff = 7.0527998738e-01
#> i = 4   | max diff = 6.7057088293e-01
#> i = 5   | max diff = 7.6258026156e-01
#> i = 6   | max diff = 7.9625679777e-01
#> i = 7   | max diff = 7.9418677424e-01
#> i = 8   | max diff = 7.0796931357e-01
#> i = 9   | max diff = 4.1580490657e-01
#> i = 10  | max diff = 1.1030966909e-01
#> i = 11  | max diff = 1.7946878411e-02
#> i = 12  | max diff = 2.2562875966e-03
#> i = 13  | max diff = 2.5298644600e-04
#> i = 14  | max diff = 2.7225692389e-05
#> i = 15  | max diff = 2.8922309980e-06
#> i = 16  | max diff = 3.0623759458e-07
#> i = 17  | max diff = 3.2422824794e-08
#> i = 18  | max diff = 3.4360847420e-09
#> Converged in 18 iterations, max diff < 1.49e-08.
print(res5)
#> $class1
#>             A           B           C           D           E
#> A  0.94773211 -0.49249803  0.22725182 -0.07529617 -0.24908314
#> B -0.49249803  1.18969215  0.21109864  0.04493370  0.06829704
#> C  0.22725182  0.21109864  1.06996036 -0.07654167 -0.09757846
#> D -0.07529617  0.04493370 -0.07654167  1.64429502 -0.15573996
#> E -0.24908314  0.06829704 -0.09757846 -0.15573996  1.13107810
#> 
#> $class2
#>             A           B           C           D          E
#> A  1.33681377 -0.36667832 -0.02448179 -0.04813310 -0.2613155
#> B -0.36667832  1.22917966 -0.16039971  0.02744357  0.1858189
#> C -0.02448179 -0.16039971  1.25169252 -0.07359242  0.2211603
#> D -0.04813310  0.02744357 -0.07359242  1.74066264 -0.1858643
#> E -0.26131552  0.18581895  0.22116025 -0.18586428  1.2771818
#> 
#> $class3
#>             A           B           C          D           E
#> A  0.92406642 -0.41872351 -0.04795204  0.2621847 -0.43141317
#> B -0.41872351  1.07733367  0.24409918  0.1818829 -0.06404697
#> C -0.04795204  0.24409918  0.89920929 -0.1981665  0.18679648
#> D  0.26218466  0.18188290 -0.19816645  1.8174454  0.11296282
#> E -0.43141317 -0.06404697  0.18679648  0.1129628  1.42893361
#>