rags2ridges (version 2.2.2)

ridgeP.fused: Fused ridge estimation

Description

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

Usage

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.

relative

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

verbose

logical. Set to TRUE for extra output.

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\).

References

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2015). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes, arXiv:1509.07982v1 [stat.ME].

See Also

default.penalty ridgeP for the regular ridge estimate

Examples

# NOT RUN {
# 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

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

res1 <- ridgeP.fused(Slist, ns, lambda = c(1.3, 2.1))
print(res1)

# 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)
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)
res3 <- ridgeP.fused(Slist, ns, lambda = mylambda)
# 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))
print(res4)

# 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))
print(res5)
# }