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()
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(fit@coef.regres, function(x) which(rowSums(x) != 0)) qvals <- sapply(lapply(fit@coef.regres, 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.