In jchiquet/spring: Fit the SPRING model

First we load the library:

library(spring)
t0 <- proc.time()

We specify a model (basically, a set of parameter): this stands for the regression coefficient and the covariance matrix.

 ## problem dimension
p <- 20 # number of regressors
q <- 3  # number of responses
k <- 10 # number of non null direct effects

## direct coefficients
Omega.xy <- Matrix(0,p,q)
Omega.xy[sample(1:(p*q),k)] <- sample(c(-1,1),k,rep=T)

## covariance of the noise (residual)
cor <- 0.5
R <- toeplitz(cor^(1:q-1))

## turn this into a "spring" model
true.model <- model$new(
     coef.direct = Omega.xy,
     coef.regres = Matrix(- Omega.xy %*% R),
     covariance  = Matrix(R),
     intercept   = Matrix(rep(0,q),q,1))

Now, produce an arbitrary design matrix and use the drawResponse method to get some data:

nb.samples <- 100
X <- matrix(rnorm(nb.samples * p), nb.samples, p)

## draw some data with a given design matrix
Y <- true.model$drawResponse(X)

out <- spring(X, Y)
BIC <- pen.crit(out)$BIC
plot(BIC$criterion)

BIC$model$plot()

Stability selection:

library(stabs)
set.seed(5678)
fitfunSpring <- function(x, y, q, ...){
  fit <- spring(x=x, y=y, comp.df = FALSE, verbose = FALSE, ...)[[1]]
  id_selected <- lapply([e-mail], function(x) which(rowSums(x) != 0))
  qvals <- sapply(lapply([e-mail], function(x) rowSums(x) != 0), sum)
  id_selected <- id_selected[qvals <= q]
  selected <- rep(FALSE, ncol(x))
  selected[id_selected[[length(id_selected)]]] <- TRUE

  path <- matrix(FALSE, ncol(x), length(id_selected))
  for(idx in seq_along(id_selected)) path[id_selected[[idx]], idx] <- TRUE
  rownames(path) <- colnames(as.data.frame(x))
  res <- list(selected = selected, path = path)
  res
}

stab.spring <- stabsel(x=X, y=as.matrix(Y),
                       PFER=1, B=4, q=2, mc.cores=1, 
                       fitfun=fitfunSpring,
                       args.fitfun=list(lambda2 = 0, min.ratio = 0.5, nlambda1 = 20,
                                        struct = sparseMatrix(i=1:p, j=1:p,x=rep(1, p))))

plot(stab.spring, "spring stability path")

t1 <- proc.time()
t1 - t0
print(sessionInfo(), locale=FALSE)

jchiquet/spring documentation built on May 18, 2019, 10:22 p.m.