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)
)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.
A numeric vector of sample sizes on which the matrices in
Slist are based. I.e. ns[g] correspond to Slist[[g]].
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.
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.
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.
A single integer giving the maximum number of allowed
iterations. Can be set to Inf. If maxit is hit, a warning
is given.
logical. Set to TRUE for extra output.
logical indicating if the convergence criterion should
be on a relative scale.
A single positive numeric giving the convergence threshold.
Returns a list as Slist with precision estimates of the
corresponding classes.
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\).
For extreme fusion penalties in lambda the algorithm is quite
sensitive to the initial values given in Plist.
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.
default.penalty ridgeP for the
regular ridge estimate
# 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.546 1.072 -0.897 0.163 0.161 ...
#> ..- attr(*, "dimnames")=List of 2
#> $ class2: num [1:5, 1:5] 0.508 0.277 0.265 0.156 -0.413 ...
#> ..- attr(*, "dimnames")=List of 2
#> $ class3: num [1:5, 1:5] 1.566 0.469 -0.456 -0.125 0.71 ...
#> ..- 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 = 3.0367895493e-01
#> i = 2 | max diff = 1.8727083322e-02
#> i = 3 | max diff = 3.4766604330e-03
#> i = 4 | max diff = 6.3247273882e-04
#> i = 5 | max diff = 1.1531201870e-04
#> i = 6 | max diff = 2.1394821563e-05
#> i = 7 | max diff = 4.0747115949e-06
#> i = 8 | max diff = 8.0199482189e-07
#> i = 9 | max diff = 1.6404990566e-07
#> i = 10 | max diff = 3.5027458680e-08
#> i = 11 | max diff = 7.8261410914e-09
#> Converged in 11 iterations, max diff < 1.49e-08.
print(res1)
#> $class1
#> A B C D E
#> A 9.4975989 -1.49128356 0.6957991 -0.31747217 -0.1956141
#> B -1.4912836 10.63983272 0.6973545 0.08555426 0.4562045
#> C 0.6957991 0.69735450 10.2926366 -0.32508052 0.4043885
#> D -0.3174722 0.08555426 -0.3250805 12.20187351 -0.1563601
#> E -0.1956141 0.45620446 0.4043885 -0.15636009 10.2099139
#>
#> $class2
#> A B C D E
#> A 582.29863021 -1.153538605 0.018545143 -0.35633520 0.21367572
#> B -1.15353860 583.027411239 -0.008986424 0.17191072 0.94727569
#> C 0.01854514 -0.008986424 582.750267173 -0.30114793 1.19332408
#> D -0.35633520 0.171910724 -0.301147929 584.55616917 -0.03149909
#> E 0.21367572 0.947275691 1.193324078 -0.03149909 582.11233401
#>
#> $class3
#> A B C D E
#> A 1.1803062 -0.7396866 0.1106980 -0.17035973 -0.42346203
#> B -0.7396866 1.9505614 0.5917630 0.25788342 0.35790463
#> C 0.1106980 0.5917630 1.5194527 -0.49385549 0.68022457
#> D -0.1703597 0.2578834 -0.4938555 3.29199107 0.02844718
#> E -0.4234620 0.3579046 0.6802246 0.02844718 1.85210630
#>
# 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 = 3.0367895493e-01
#> i = 2 | max diff = 1.8727083322e-02
#> i = 3 | max diff = 3.4766604330e-03
#> i = 4 | max diff = 6.3247273882e-04
#> i = 5 | max diff = 1.1531201870e-04
#> i = 6 | max diff = 2.1394821563e-05
#> i = 7 | max diff = 4.0747115949e-06
#> i = 8 | max diff = 8.0199482189e-07
#> i = 9 | max diff = 1.6404990566e-07
#> i = 10 | max diff = 3.5027458680e-08
#> i = 11 | max diff = 7.8261410914e-09
#> Converged in 11 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 = 4.0122942069e-01
#> i = 2 | max diff = 7.3627641176e-03
#> i = 3 | max diff = 2.3634364338e-03
#> i = 4 | max diff = 2.3185281854e-04
#> i = 5 | max diff = 2.3785358600e-05
#> i = 6 | max diff = 2.5863349741e-06
#> i = 7 | max diff = 2.9566282889e-07
#> i = 8 | max diff = 3.5137581769e-08
#> i = 9 | max diff = 4.3064653738e-09
#> Converged in 9 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.2665407850e+01
#> i = 2 | max diff = 3.5457820795e+00
#> i = 3 | max diff = 7.0068257338e-01
#> i = 4 | max diff = 7.2429885604e-01
#> i = 5 | max diff = 8.4709534890e-01
#> i = 6 | max diff = 8.4070776325e-01
#> i = 7 | max diff = 6.3960573506e-01
#> i = 8 | max diff = 3.3744898320e-01
#> i = 9 | max diff = 1.1180494878e-01
#> i = 10 | max diff = 2.4130378665e-02
#> i = 11 | max diff = 3.9601545661e-03
#> i = 12 | max diff = 5.6656994838e-04
#> i = 13 | max diff = 7.6836375497e-05
#> i = 14 | max diff = 1.0261461462e-05
#> i = 15 | max diff = 1.3693289378e-06
#> i = 16 | max diff = 1.8349573770e-07
#> i = 17 | max diff = 2.4726744801e-08
#> i = 18 | max diff = 3.3509183752e-09
#> Converged in 18 iterations, max diff < 1.49e-08.
print(res4)
#> $class1
#> A B C D E
#> A 1.3041441 -0.73754139 0.3363622 -0.22340612 -0.20361615
#> B -0.7375414 1.84856542 0.3104566 -0.05351622 0.04772841
#> C 0.3363622 0.31045660 1.6249836 -0.21055935 -0.10838058
#> D -0.2234061 -0.05351622 -0.2105593 2.49321804 -0.15728225
#> E -0.2036162 0.04772841 -0.1083806 -0.15728225 1.68636401
#>
#> $class2
#> A B C D E
#> A 1.7236118 0.4078962 0.7482677 0.6214714 0.8891218
#> B 0.4078962 2.0453675 0.7459626 0.8983823 1.1544750
#> C 0.7482677 0.7459626 1.9964389 0.7522931 1.2619198
#> D 0.6214714 0.8983823 0.7522931 2.6127972 0.8901378
#> E 0.8891218 1.1544750 1.2619198 0.8901378 1.9165903
#>
#> $class3
#> A B C D E
#> A 0.76339094 -0.30021508 0.04307929 -0.081956811 -0.337791557
#> B -0.30021508 1.09419399 0.29119415 0.080841973 0.018322805
#> C 0.04307929 0.29119415 0.87222565 -0.297370748 0.227757917
#> D -0.08195681 0.08084197 -0.29737075 1.684798441 0.009372881
#> E -0.33779156 0.01832280 0.22775792 0.009372881 1.260161592
#>
# 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.3444160465e+01
#> i = 2 | max diff = 3.6383661303e+00
#> i = 3 | max diff = 7.3820857508e-01
#> i = 4 | max diff = 7.2360444799e-01
#> i = 5 | max diff = 8.4500858618e-01
#> i = 6 | max diff = 8.9090121627e-01
#> i = 7 | max diff = 8.0079289102e-01
#> i = 8 | max diff = 5.3907965324e-01
#> i = 9 | max diff = 1.8809562490e-01
#> i = 10 | max diff = 3.8975794046e-02
#> i = 11 | max diff = 5.7107767311e-03
#> i = 12 | max diff = 7.0973068633e-04
#> i = 13 | max diff = 8.3297683956e-05
#> i = 14 | max diff = 9.6442563456e-06
#> i = 15 | max diff = 1.1180544917e-06
#> i = 16 | max diff = 1.3033117814e-07
#> i = 17 | max diff = 1.5285158604e-08
#> i = 18 | max diff = 1.8024837296e-09
#> Converged in 18 iterations, max diff < 1.49e-08.
print(res5)
#> $class1
#> A B C D E
#> A 0.8792906 -0.46163350 0.25483586 -0.12774353 -0.12080558
#> B -0.4616335 1.26409038 0.25980669 0.01482165 0.07774399
#> C 0.2548359 0.25980669 1.08596253 -0.12174527 -0.03575802
#> D -0.1277435 0.01482165 -0.12174527 1.73976750 -0.10028210
#> E -0.1208056 0.07774399 -0.03575802 -0.10028210 1.03373825
#>
#> $class2
#> A B C D E
#> A 1.14827426 -0.29913416 -0.04178625 -0.15389463 0.07155385
#> B -0.29913416 1.35014726 -0.07895991 0.08356129 0.30036714
#> C -0.04178625 -0.07895991 1.27622462 -0.08679001 0.40516656
#> D -0.15389463 0.08356129 -0.08679001 1.82208663 -0.01714969
#> E 0.07155385 0.30036714 0.40516656 -0.01714969 1.02022752
#>
#> $class3
#> A B C D E
#> A 0.82756800 -0.26551333 0.07501802 -0.03878723 -0.33842203
#> B -0.26551333 1.19646672 0.35289732 0.13709729 0.03520421
#> C 0.07501802 0.35289732 0.98140935 -0.25602667 0.28531030
#> D -0.03878723 0.13709729 -0.25602667 1.79410383 0.03387104
#> E -0.33842203 0.03520421 0.28531030 0.03387104 1.34417265
#>