Commit 5b89a40c authored by Armin Rauschenberger's avatar Armin Rauschenberger
Browse files

initialise

parent bcfdcc05
Package: cornet
Package: mixnet
Version: 0.0.1
Title: Elastic Net with Dichotomised Outcomes
Description: Implements lasso and ridge regression for dichotomised outcomes (Rauschenberger et al. 2019). Such outcomes are not naturally but artificially binary. They indicate whether an underlying measurement is greater than a threshold.
Title: Multivariate Regression
Description: Implements high-dimensional multivariate regression by stacking (Rauschenberger et al. 2019).
Depends: R (>= 3.0.0)
Imports: glmnet, palasso
Suggests: knitr, testthat
Enhances: RColorBrewer
Authors@R: person("Armin","Rauschenberger",email="a.rauschenberger@vumc.nl",role=c("aut","cre"))
Imports: glmnet, palasso, cornet
Suggests: knitr, testthat, MASS
Enhances: RColorBrewer, spls, SiER, MRCE
Authors@R: person("Armin","Rauschenberger",email="armin.rauschenberger@uni.lu",role=c("aut","cre"))
VignetteBuilder: knitr
License: GPL-3
LazyData: true
Language: en-GB
RoxygenNote: 6.1.1
URL: https://github.com/rauschenberger/cornet
BugReports: https://github.com/rauschenberger/cornet/issues
URL: https://github.com/rauschenberger/stanet
BugReports: https://github.com/rauschenberger/stanet/issues
# Generated by roxygen2: do not edit by hand
S3method(coef,cornet)
S3method(plot,cornet)
S3method(predict,cornet)
S3method(print,cornet)
export(cornet)
S3method(coef,mixnet)
S3method(predict,mixnet)
S3method(weights,mixnet)
export(mixnet)
importFrom(stats,weights)
## cornet 0.0.1 (2019-03-15)
## mixnet 0.0.1 (2019-06-01)
* first submission
\ No newline at end of file
#--- Workhorse function --------------------------------------------------------
#--- Main function -------------------------------------------------------------
#' @export
#' @aliases cornet-package
#' @title
#' Combined regression
#' Multivariate Elastic Net Regression
#'
#' @description
#' Implements lasso and ridge regression for dichotomised outcomes.
#' Such outcomes are not naturally but artificially binary.
#' They indicate whether an underlying measurement is greater than a threshold.
#'
#' @param y
#' continuous outcome\strong{:}
#' vector of length \eqn{n}
#'
#' @param cutoff
#' cut-off point for dichotomising outcome into classes\strong{:}
#' \emph{meaningful} value between \code{min(y)} and \code{max(y)}
#' Implements multivariate elastic net regression.
#'
#' @param Y
#' outputs\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
#' and \eqn{q} columns (variables),
#' with positive correlation (see details)
#'
#' @param X
#' features\strong{:}
#' inputs\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
#' and \eqn{p} columns (variables)
#'
#' @param alpha
#' elastic net mixing parameter\strong{:}
#' numeric between \eqn{0} (ridge) and \eqn{1} (lasso)
#'
#'
#' @param family
#' distribution\strong{:}
#' vector of length \eqn{1} or \eqn{q} with entries
#' \code{"gaussian"}, \code{"binomial"} or \code{"poisson"}
#'
#' @param nfolds
#' number of folds
#'
#' @param foldid
#' fold identifiers\strong{:}
#' vector with entries between \eqn{1} and \code{nfolds};
#' vector of length \eqn{n} with entries between \eqn{1} and \code{nfolds};
#' or \code{NULL} (balance)
#'
#' @param nfolds
#' number of folds\strong{:}
#' integer between \eqn{3} and \eqn{n}
#'
#' @param type.measure
#' loss function for binary classification\strong{:}
#' character \code{"deviance"}, \code{"mse"}, \code{"mae"},
#' or \code{"class"} (see \code{\link[glmnet]{cv.glmnet}})
#'
#' @param pi
#' pi sequence\strong{:}
#' vector of increasing values in the unit interval;
#' or \code{NULL} (default sequence)
#'
#' @param npi
#' number of \code{pi} values (weighting)
#'
#' @param sigma
#' sigma sequence\strong{:}
#' vector of increasing positive values;
#' or \code{NULL} (default sequence)
#' loss function\strong{:}
#' vector of length \eqn{1} or \eqn{q} with entries
#' \code{"deviance"}, \code{"class"}, \code{"mse"} or \code{"mae"}
#' (see \code{\link[glmnet]{cv.glmnet}})
#'
#' @param alpha.base
#' elastic net mixing parameter for base learners\strong{:}
#' numeric between \eqn{0} (ridge) and \eqn{1} (lasso)
#'
#' @param nsigma
#' number of \code{sigma} values (scaling)
#' @param alpha.meta
#' elastic net mixing parameter for meta learner\strong{:}
#' numeric between \eqn{0} (ridge) and \eqn{1} (lasso),
#'
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}}
#'
#' @details
#' The argument \code{family} is unavailable, because
#' this function fits a \emph{gaussian} model for the numeric response,
#' and a \emph{binomial} model for the binary response.
#'
#' Linear regression uses the loss function \code{"deviance"} (or \code{"mse"}),
#' but the loss is incomparable between linear and logistic regression.
#'
#' The loss function \code{"auc"} is unavailable for internal cross-validation.
#' If at all, use \code{"auc"} for external cross-validation only.
#'
#' @return
#' Returns an object of class \code{cornet}, a list with multiple slots:
#' \itemize{
#' \item \code{gaussian}: fitted linear model, class \code{glmnet}
#' \item \code{binomial}: fitted logistic model, class \code{glmnet}
#' \item \code{sigma}: scaling parameters \code{sigma},
#' vector of length \code{nsigma}
#' \item \code{pi}: weighting parameters \code{pi},
#' vector of length \code{npi}
#' \item \code{cvm}: evaluation loss,
#' matrix with \code{nsigma} rows and \code{npi} columns
#' \item \code{sigma.min}: optimal scaling parameter,
#' positive scalar
#' \item \code{pi.min}: optimal weighting parameter,
#' scalar in unit interval
#' \item \code{cutoff}: threshold for dichotomisation
#' }
#'
#' @seealso
#' Methods for objects of class \code{cornet} include
#' \code{\link[=coef.cornet]{coef}} and
#' \code{\link[=predict.cornet]{predict}}.
#'
#' @references
#' Armin Rauschenberger and Enrico Glaab (2019).
#' "Lasso and ridge regression for dichotomised outcomes".
#' \emph{Manuscript in preparation}.
#' A Rauschenberger, E Glaab (2019)
#' "Multivariate elastic net regression through stacked generalisation"
#' \emph{Manuscript in preparation.}
#'
#' @details
#' The \eqn{q} outcomes should be positively correlated.
#' Avoid negative correlations by changing the sign of the variable.
#'
#' @examples
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
#' net <- cornet(y=y,cutoff=0,X=X)
#' net
#' object <- mixnet(Y=Y,X=X)
#'
cornet <- function(y,cutoff,X,alpha=1,npi=101,pi=NULL,nsigma=99,sigma=NULL,nfolds=10,foldid=NULL,type.measure="deviance",...){
mixnet <- function(Y,X,family="gaussian",nfolds=10,foldid=NULL,type.measure="deviance",alpha.base=1,alpha.meta=0,...){
#--- temporary ---
# cutoff <- 0; npi <- 101; pi <- NULL; nsigma <- 99; sigma <- NULL; nfolds <- 10; foldid <- NULL; type.measure <- "deviance"; logistic <- TRUE
test <- list()
test$combined <- TRUE
# family <- "gaussian"; nfolds <- 10; foldid <- NULL; type.measure <- "deviance"
#--- checks ---
n <- length(y)
.check(x=y,type="vector")
if(all(y %in% c(0,1))){warning("Binary response.",call.=FALSE)}
.check(x=cutoff,type="scalar",min=min(y),max=max(y))
if(length(y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)}
.check(x=X,type="matrix",dim=c(n,NA))
.check(x=alpha,type="scalar",min=0,max=1)
.check(x=npi,type="scalar",min=1)
.check(x=pi,type="vector",min=0,max=1,null=TRUE)
.check(x=nsigma,type="scalar",min=1)
.check(x=sigma,type="vector",min=.Machine$double.eps,null=TRUE)
.check(x=nfolds,type="scalar",min=3,max=n)
.check(x=foldid,type="vector",dim=n,values=seq_len(nfolds),null=TRUE)
.check(x=type.measure,type="string",values=c("deviance","mse","mae","class"))
# never auc (min/max confusion)!
if(!is.null(list(...)$family)){stop("Reserved argument \"family\".",call.=FALSE)}
# binarisation
z <- 1*(y > cutoff)
# fold identifiers
if(is.null(foldid)){
foldid <- palasso:::.folds(y=z,nfolds=nfolds)
cornet:::.check(x=Y,type="matrix",miss=TRUE)
if(any(stats::cor(Y,use="pairwise.complete.obs")<0,na.rm=TRUE)){warning("Negative correlation!",call.=FALSE)}
cornet:::.check(x=X,type="matrix")
#cornet:::.check(x=family,type="string",values=c("gaussian","binomial","poisson"))
if(nrow(Y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)}
cornet:::.check(x=nfolds,type="scalar",min=3)
cornet:::.check(x=foldid,type="vector",values=seq_len(nfolds),null=TRUE)
cornet:::.check(x=type.measure,type="string",values=c("deviance","class","mse","mae")) # not auc (min/max confusion)
cornet:::.check(x=alpha.base,type="scalar",min=0,max=1)
cornet:::.check(x=alpha.meta,type="scalar",min=0,max=1)
if(!is.null(c(list(...)$lower.limits,list(...)$upper.limits))){
stop("Reserved arguments \"lower.limits\" and \"upper.limits\".",call.=FALSE)
}
#--- model fitting ---
fit <- list()
fit$gaussian <- glmnet::glmnet(y=y,x=X,family="gaussian",alpha=alpha,...)
fit$binomial <- glmnet::glmnet(y=z,x=X,family="binomial",alpha=alpha,...)
#--- sigma sequence ---
if(is.null(sigma)){
fit$sigma <- exp(seq(from=log(0.05*stats::sd(y)),
to=log(10*stats::sd(y)),length.out=nsigma))
} else {
fit$sigma <- sigma
nsigma <- length(sigma)
#--- dimensionality ---
n <- nrow(Y)
q <- ncol(Y)
p <- ncol(X)
#--- family ---
if(length(family)==1){
family <- rep(family,times=q)
} else if(length(family)!=q){
stop("Invalid argument family",call.=FALSE)
}
#--- pi sequence ---
if(is.null(pi)){
fit$pi <- seq(from=0,to=1,length.out=npi)
#--- fold identifiers ---
# provide foldid as matrix?
if(is.null(foldid)){
foldid <- palasso:::.folds(y=Y[,1],nfolds=nfolds) # temporary Y[,1]
} else {
fit$pi <- pi
npi <- length(pi)
nfolds <- length(unique(foldid))
}
#--- tuning parameters ---
lab.pi <- paste0("pi",seq_len(npi))
lab.sigma <- paste0("si",seq_len(nsigma))
names(fit$sigma) <- lab.sigma
names(fit$pi) <- lab.pi
#--- full fit ---
nlambda <- numeric()
base <- lapply(seq_len(q),function(x) list())
for(i in seq_len(q)){
cond <- !is.na(Y[,i])
#if(sum(cond)==0){nlambda[i] <- 0; next}
base[[i]]$glmnet.fit <- glmnet::glmnet(y=Y[cond,i],x=X[cond,],family=family[i],alpha=alpha.base,...) # ellipsis
base[[i]]$lambda <- base[[i]]$glmnet.fit$lambda
nlambda[i] <- length(base[[i]]$glmnet.fit$lambda)
}
#--- cross-validation ---
pred <- list()
pred$y <- matrix(data=NA,nrow=n,ncol=length(fit$gaussian$lambda))
pred$z <- matrix(data=NA,nrow=n,ncol=length(fit$binomial$lambda))
if(test$combined){
dimnames <- list(NULL,lab.sigma,lab.pi)
pred$combined <- array(data=NA,dim=c(n,nsigma,npi),dimnames=dimnames)
#--- predictions ---
link <- list()
for(i in seq_len(q)){
link[[i]] <- matrix(data=NA,nrow=n,ncol=nlambda[i])
}
#--- base cross-validation ---
for(k in seq_len(nfolds)){
y0 <- y[foldid!=k]
y1 <- y[foldid==k]
z0 <- z[foldid!=k]
z1 <- z[foldid==k]
Y0 <- Y[foldid!=k,,drop=FALSE]
Y1 <- Y[foldid==k,,drop=FALSE]
X0 <- X[foldid!=k,,drop=FALSE]
X1 <- X[foldid==k,,drop=FALSE]
# linear regression
net <- glmnet::glmnet(y=y0,x=X0,family="gaussian",alpha=alpha,...)
temp_y <- stats::predict(object=net,newx=X1,type="response",s=fit$gaussian$lambda)
pred$y[foldid==k,seq_len(ncol(temp_y))] <- temp_y
cvm <- palasso:::.loss(y=y1,fit=temp_y,family="gaussian",type.measure="deviance")[[1]]
y_hat <- temp_y[,which.min(cvm)]
# logistic regression
net <- glmnet::glmnet(y=z0,x=X0,family="binomial",alpha=alpha,...)
temp_z <- stats::predict(object=net,newx=X1,type="response",s=fit$binomial$lambda)
pred$z[foldid==k,seq_len(ncol(temp_z))] <- temp_z
cvm <- palasso:::.loss(y=z1,fit=temp_z,family="binomial",type.measure=type.measure)[[1]]
z_hat <- temp_z[,which.min(cvm)]
# combined regression
if(test$combined){
for(i in seq_along(fit$sigma)){
for(j in seq_along(fit$pi)){
cont <- stats::pnorm(q=y_hat,mean=cutoff,sd=fit$sigma[i])
pred$combined[foldid==k,i,j] <- fit$pi[j]*cont + (1-fit$pi[j])*z_hat
}
}
for(i in seq_len(q)){
cond <- !is.na(Y0[,i])
#if(sum(cond)==0){next}
object <- glmnet::glmnet(y=Y0[cond,i],x=X0[cond,],family=family[i],alpha=alpha.base,...) # ellipsis
temp <- stats::predict(object=object,newx=X1,type="link",
s=base[[i]]$glmnet.fit$lambda)
link[[i]][foldid==k,seq_len(ncol(temp))] <- temp
}
}
#--- evaluation ---
#--- tune base lambdas ---
for(i in seq_len(q)){
fit <- .mean.function(link[[i]],family=family[i])
cond <- !is.na(Y[,i])
base[[i]]$cvm <- palasso:::.loss(y=Y[cond,i],fit=fit[cond,],
family=family[i],type.measure=type.measure)[[1]]
base[[i]]$lambda.min <- base[[i]]$lambda[which.min(base[[i]]$cvm)]
}
# linear loss
fit$gaussian$cvm <- palasso:::.loss(y=y,fit=pred$y,family="gaussian",type.measure="deviance")[[1]]
fit$gaussian$lambda.min <- fit$gaussian$lambda[which.min(fit$gaussian$cvm)]
#--- predictions ---
hat <- matrix(NA,nrow=n,ncol=q)
for(i in seq_len(q)){
hat[,i] <- link[[i]][,base[[i]]$lambda==base[[i]]$lambda.min]
}
# logistic loss
fit$binomial$cvm <- palasso:::.loss(y=z,fit=pred$z,family="binomial",type.measure=type.measure)[[1]]
fit$binomial$lambda.min <- fit$binomial$lambda[which.min(fit$binomial$cvm)]
# combined loss
if(test$combined){
dimnames <- list(lab.sigma,lab.pi)
fit$cvm <- matrix(data=NA,nrow=nsigma,ncol=npi,dimnames=dimnames)
for(i in seq_len(nsigma)){
for(j in seq_len(npi)){
fit$cvm[i,j] <- palasso:::.loss(y=z,fit=pred$combined[,i,j],family="binomial",type.measure=type.measure)[[1]]
}
}
temp <- which(fit$cvm==min(fit$cvm),arr.ind=TRUE,useNames=TRUE)
if(nrow(temp)>1){warning("Multiple!",call.=FALSE);temp <- temp[1,,drop=FALSE]}
fit$sigma.min <- fit$sigma[temp[1]]
fit$pi.min <- fit$pi[temp[2]]
if(fit$cvm[names(fit$sigma.min),names(fit$pi.min)]!=min(fit$cvm)){stop("Internal error.")}
#--- meta cross-validation ---
meta <- list()
for(i in seq_len(q)){
cond <- !is.na(Y[,i])
meta[[i]] <- glmnet::cv.glmnet(y=Y[cond,i],x=hat[cond,],
lower.limits=0, # important: 0
upper.limits=Inf, # important: Inf
foldid=foldid[cond],
family=family[i],
type.measure=type.measure,
intercept=TRUE, # with intercept
alpha=alpha.meta,...) # ellipsis
# consider trying different alpha.meta and selecting best one
}
#--- return ---
fit$cutoff <- cutoff
fit$info <- list(type.measure=type.measure,
sd.y=stats::sd(y),
"+"=sum(z==1),
"-"=sum(z==0),
table=table(z),
n=n,p=ncol(X),
test=as.data.frame(test))
class(fit) <- "cornet"
return(fit)
names(base) <- names(meta) <- paste0("y",seq_len(q))
info <- data.frame(q=q,p=p,family=family,type.measure=type.measure)
list <- list(base=base,meta=meta,info=info)
class(list) <- "mixnet"
return(list)
}
#--- Methods -------------------------------------------------------------------
.mean.function <- function(x,family){
if(family=="gaussian"){
return(x)
} else if(family=="binomial"){
return(1/(1+exp(-x)))
} else if(family=="poisson"){
return(exp(x))
} else {
stop("Family not implemented.",call.=FALSE)
}
}
#' @export
#' @title
#' Combined regression
#'
#' @description
#' Prints summary of cornet object.
#'
#' @param x
#' \link[cornet]{cornet} object
#'
#' @param ...
#' further arguments (not applicable)
#'
#' @return
#' Returns sample size \eqn{n},
#' number of covariates \eqn{p},
#' information on dichotomisation,
#' tuned scaling parameter (sigma),
#' tuned weighting parameter (pi),
#' and corresponding loss.
#'
#' @examples
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
#' net <- cornet(y=y,cutoff=0,X=X)
#' print(net)
#'
print.cornet <- function(x,...){
cat("cornet object:\n")
cat(paste0("n = ",x$info$n,","," p = ",x$info$p),"\n")
cat(paste0("z = I(y > ",signif(x$cutoff,digits=2),"): "))
cat(paste0(x$info$"+","+"," vs ",x$info$"-","-"),"\n")
cat(paste0("sigma.min = ",signif(x$sigma.min,digits=1)),"\n")
cat(paste0("pi.min = ",round(x$pi.min,digits=2)),"\n")
type <- x$info$type.measure
value <- signif(x$cvm[names(x$sigma.min),names(x$pi.min)],digits=2)
cat(paste0(type," = ",value))
return(invisible(NULL))
.link.function <- function(x,family){
if(family=="gaussian"){
return(x)
} else if(family=="binomial"){
return(log(1/(1-x)))
} else if(family=="poisson"){
return(log(x))
} else {
stop("Family not implemented.",call.=FALSE)
}
}
#--- Methods for class "mixnet" -----------------------------------------------
#' @export
#' @title
#' Extract estimated coefficients
#' Makes Predictions
#'
#' @description
#' Extracts estimated coefficients from linear and logistic regression,
#' under the penalty parameter that minimises the cross-validated loss.
#'
#' @param object
#' \link[cornet]{cornet} object
#' Predicts outcome from features with stacked model.
#'
#' @param ...
#' further arguments (not applicable)
#' @param object
#' \link[mixnet]{mixnet} object
#'
#' @return
#' This function returns a matrix with \eqn{n} rows and two columns,
#' where \eqn{n} is the sample size. It includes the estimated coefficients
#' from linear regression (1st column: \code{"beta"})
#' and logistic regression (2nd column: \code{"gamma"}).
#'
#' @examples
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
#' net <- cornet(y=y,cutoff=0,X=X)
#' coef(net)
#' @param newx
#' covariates\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
#' and \eqn{p} columns (variables)
#'
coef.cornet <- function(object,...){
if(length(list(...))!=0){warning("Ignoring arguments.")}
s <- object$gaussian$lambda.min
beta <- glmnet::coef.glmnet(object=object$gaussian,s=s)
s <- object$binomial$lambda.min
gamma <- glmnet::coef.glmnet(object=object$binomial,s=s)
coef <- cbind(beta,gamma)
colnames(coef) <- c("beta","gamma")
return(coef)
}
#' @export
#' @title
#' Plot loss matrix
#'
#' @description
#' Plots the loss for different combinations of
#' scaling (sigma) and weighting (pi) parameters.
#'
#' @param x
#' \link[cornet]{cornet} object
#' @param type
#' character "link" or "response"
#'
#' @param ...
#' further arguments (not applicable)
#'
#' @return
#' This function plots the evaluation loss (\code{cvm}).
#' Whereas the matrix has sigma in the rows, and pi in the columns,
#' the plot has sigma on the \eqn{x}-axis, and pi on the \eqn{y}-axis.
#' For all combinations of sigma and pi, the colour indicates the loss.
#' If the R package \code{RColorBrewer} is installed,
#' blue represents low. Otherwise, red represents low.
#' White always represents high.
#'
#' @examples
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
#' Y <- matrix(rbinom(n=n*q,size=1,prob=0.5),nrow=n,ncol=q)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
#' net <- cornet(y=y,cutoff=0,X=X)
#' plot(net)
#' object <- mixnet(Y=Y,X=X,family="binomial")
#' y_hat <- predict(object,newx=X)
#'
plot.cornet <- function(x,...){
predict.mixnet <- function(object,newx,type="response",...){
if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
if(length(list(...))!=0){warning("Ignoring arguments.")}
k <- 100
levels <- stats::quantile(x$cvm,probs=seq(from=0,to=1,length.out=k+1))
x <- object; rm(object)