functions.R 15.4 KB
Newer Older
Armin Rauschenberger's avatar
Armin Rauschenberger committed
1

Armin Rauschenberger's avatar
Armin Rauschenberger committed
2
#--- Main function -------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
3

Armin Rauschenberger's avatar
Armin Rauschenberger committed
4 5
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
6
#' Multivariate Elastic Net Regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
7 8
#' 
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
9 10 11 12 13 14 15
#' 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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
16 17
#' 
#' @param X
Armin Rauschenberger's avatar
Armin Rauschenberger committed
18
#' inputs\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
19
#' numeric matrix with \eqn{n} rows (samples)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
20
#' and \eqn{p} columns (variables)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
21 22 23 24 25 26 27 28 29
#'
#' @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
#'
Armin Rauschenberger's avatar
Armin Rauschenberger committed
30 31
#' @param foldid
#' fold identifiers\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
32
#' vector of length \eqn{n} with entries between \eqn{1} and \code{nfolds};
Armin Rauschenberger's avatar
Armin Rauschenberger committed
33
#' or \code{NULL} (balance)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
34
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
35
#' @param type.measure
Armin Rauschenberger's avatar
Armin Rauschenberger committed
36 37 38 39 40 41 42 43
#' 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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
44
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
45 46 47
#' @param alpha.meta
#' elastic net mixing parameter for meta learner\strong{:}
#' numeric between \eqn{0} (ridge) and \eqn{1} (lasso),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
48
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
49 50 51
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}}
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
52
#' @references 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
53 54 55 56 57 58 59
#' 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.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
60
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
61
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
62 63
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
64
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
65
#' object <- mixnet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
66
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
67
mixnet <- function(Y,X,family="gaussian",nfolds=10,foldid=NULL,type.measure="deviance",alpha.base=1,alpha.meta=0,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
68
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
69
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
70
  # family <- "gaussian"; nfolds <- 10; foldid <- NULL; type.measure <- "deviance"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
71
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
72
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
73 74 75 76 77 78 79 80 81 82 83 84
  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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
85 86
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
87 88 89 90 91 92 93 94 95 96
  #--- 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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
97
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
98
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
99 100 101 102
  #--- fold identifiers ---
  # provide foldid as matrix?
  if(is.null(foldid)){
    foldid <- palasso:::.folds(y=Y[,1],nfolds=nfolds) # temporary Y[,1]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
103
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
104
    nfolds <- length(unique(foldid))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
105 106
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
107 108 109 110 111 112 113 114 115 116
  #--- 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)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
117
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
118 119 120 121
  #--- predictions ---
  link <- list()
  for(i in seq_len(q)){
    link[[i]] <- matrix(data=NA,nrow=n,ncol=nlambda[i])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
122
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
123 124
  
  #--- base cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
125
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
126 127
    Y0 <- Y[foldid!=k,,drop=FALSE]
    Y1 <- Y[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
128 129
    X0 <- X[foldid!=k,,drop=FALSE]
    X1 <- X[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
130 131 132 133 134 135 136
    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
Armin Rauschenberger's avatar
Armin Rauschenberger committed
137
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
138 139
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
140 141 142 143 144 145 146 147
  #--- 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)]
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
148
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
149 150 151 152 153
  #--- 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]
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
154
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
155 156 157 158 159 160 161 162 163 164 165 166 167
  #--- 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
Armin Rauschenberger's avatar
Armin Rauschenberger committed
168 169
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
170
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
171 172 173 174 175
  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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
176 177
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
178 179 180 181 182 183 184 185 186 187 188
.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)
  }
}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
189

Armin Rauschenberger's avatar
Armin Rauschenberger committed
190 191 192 193 194 195 196 197 198 199
.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)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
200 201
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
202 203
#--- Methods for class "mixnet" -----------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
204 205
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
206
#' Makes Predictions
Armin Rauschenberger's avatar
Armin Rauschenberger committed
207 208
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
209
#' Predicts outcome from features with stacked model.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
210
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
211 212
#' @param object
#' \link[mixnet]{mixnet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
213
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
214 215 216 217
#' @param newx
#' covariates\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
#' and \eqn{p} columns (variables)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
218
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
219 220
#' @param type
#' character "link" or "response"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
221 222
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
223
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
224
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
225
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
226 227 228
#' 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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
229
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
230 231
#' object <- mixnet(Y=Y,X=X,family="binomial")
#' y_hat <- predict(object,newx=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
232
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
233 234
predict.mixnet <- function(object,newx,type="response",...){
  if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
235
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
236
  x <- object; rm(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
237
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
238
  cornet:::.check(x=newx,type="matrix")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
239
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
240 241 242
  q <- length(x$base)
  n <- nrow(newx)
  base <- meta <- matrix(NA,nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
243
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
244 245 246 247 248 249
  # base learners
  for(i in seq_len(q)){
    #base[,i] <- as.numeric(stats::predict(object=x$base[[i]]$glmnet.fit,newx=newx,s=x$base[[i]]$lambda.min,type="link"))
    base[,i] <- as.numeric(glmnet::predict.cv.glmnet(object=x$base[[i]],newx=newx,s="lambda.min",type="link"))
    # check whether fine for "binomial" family
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
250
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
251 252 253 254 255
  # meta learners
  for(i in seq_len(q)){
    meta[,i] <- as.numeric(stats::predict(object=x$meta[[i]],
                                          newx=base,s="lambda.min",type="link"))
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
256
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
257
  list <- list(base=base,meta=meta)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
258
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
259 260 261 262 263 264 265 266
  
  if(type=="response"){
    for(i in seq_len(q)){
      base[,i] <- .mean.function(x=base[,i],family=x$info$family[i])
      meta[,i] <- .mean.function(x=meta[,i],family=x$info$family[i])
    }
  } else if(type!="link"){
    stop("Invalid type.",call.=FALSE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
267 268
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
269 270 271
  list <- list(base=base,meta=meta)
  return(list)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
272 273
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
274 275
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
276
#' Extract Coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
277 278
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
279 280
#' Extracts pooled coefficients.
#' (The meta learners weights the coefficients from the base learners.)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
281
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
282
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
283
#' \link[mixnet]{mixnet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
284 285
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
286
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
287 288
#' 
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
289 290
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
291
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
#' object <- mixnet(Y=Y,X=X)
#' coef <- coef(object)
#' 
coef.mixnet <- function(object,...){
  if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
  
  # base coefficients
  base <- list()
  coef <- sapply(object$base,function(x) glmnet::coef.glmnet(object=x$glmnet.fit,s=x$lambda.min))
  base$alpha <- sapply(coef,function(x) x[1,])
  base$beta <- sapply(coef,function(x) x[-1,])
  names(base$alpha) <- colnames(base$beta) <- names(object$base)
  
  # meta coefficients
  meta <- list()
  weights <- weights.mixnet(object)
  meta$alpha <- weights[1,]
  meta$beta <- weights[-1,]
  
  # pooled coefficients
  pool <- list()
  pool$alpha <- meta$alpha + base$alpha %*% meta$beta
  pool$beta <- base$beta %*% meta$beta
  
  # q <- unique(object$info$q)
  # p <- unique(object$info$p)
  # pool$alpha <- rep(NA,times=q)
  # for(i in seq_len(q)){
  #   pool$alpha[i] <-  meta$alpha[i] + sum(meta$beta[,i] * base$alpha)
  # }
  # pool$beta <- matrix(NA,nrow=p,ncol=q)
  # for(i in seq_len(p)){
  #   for(j in seq_len(q)){
  #     pool$beta[i,j] <-  sum(meta$beta[,j] * base$beta[i,])
  #   }
  # }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
328
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
329
  return(pool)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
330 331
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
332 333
#' @export
#' @importFrom stats weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
334
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
335
#' Extract Weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
336 337
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
338 339
#' Extracts coefficients from the meta learner,
#' i.e. the weights for the base learners.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
340
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
341 342
#' @param object
#' \link[mixnet]{mixnet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
343
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
344 345
#' @param ...
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
346
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362
#' @examples
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
#' object <- mixnet(Y=Y,X=X)
#' weights(object)
#' 
# object <- mixnet(Y=Y,X=X)
weights.mixnet <- function(object,...){
  if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
  x <- object$meta
  coef <- lapply(object$meta,function(x) glmnet::coef.glmnet(object=x,s=x$lambda.min))
  coef <- do.call(what="cbind",args=coef)
  coef <- as.matrix(coef)
  colnames(coef) <- names(object$meta)
  return(coef)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
363 364
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
365 366 367 368 369 370 371 372 373 374 375 376 377 378 379
print.mixnet <- function(x,...){
  cat(paste0("mixnet object"),"\n")
}

#--- Manuscript functions ------------------------------------------------------

# @param spls
# sparse partial least squares (logical)
# @param mnorm
# multivariate normal regression (logical)
# @param mrce
# multivariate regression with covariance estimation (logical)



Armin Rauschenberger's avatar
Armin Rauschenberger committed
380
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
381
#' Model comparison
Armin Rauschenberger's avatar
Armin Rauschenberger committed
382 383
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
384
#' Compares univariate and multivariate regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
385
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
386
#' @inheritParams mixnet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
387
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
388 389
#' @param nfolds.ext
#' number of external folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
390
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391 392
#' @param nfolds.int
#' number of internal folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
393
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
394 395 396 397 398
#' @param foldid.ext
#' external fold identifiers\strong{:}
#' vector of length \eqn{n} with entries
#' between \eqn{1} and \code{nfolds.ext};
#' or \code{NULL}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
399
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
400 401 402 403 404
#' @param foldid.int
#' internal fold identifiers\strong{:}
#' vector of length \eqn{n} with entries
#' between \eqn{1} and \code{nfolds.int};
#' or \code{NULL}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
405
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
406 407
#' @param mnorm
#' multivariate normal regression\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
408 409
#' logical
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
410 411
#' @param spls
#' sparse partial least squares\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
412 413
#' logical
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
414 415 416
#' @param sier
#' high-dimensional multivariate regression\strong{:}
#' logical
Armin Rauschenberger's avatar
Armin Rauschenberger committed
417
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
418 419 420
#' @param mrce
#' multivariate regression with covariance estimation\strong{:}
#' logical
Armin Rauschenberger's avatar
Armin Rauschenberger committed
421
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
422 423
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}} and \code{\link[glmnet]{cv.glmnet}}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
424
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
425
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
426
#' NA
Armin Rauschenberger's avatar
Armin Rauschenberger committed
427
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
428
compare <- function(Y,X,family="gaussian",nfolds.ext=5,nfolds.int=10,foldid.ext=NULL,foldid.int=NULL,type.measure="deviance",alpha.base=1,alpha.meta=0,mnorm=FALSE,spls=FALSE,sier=FALSE,mrce=FALSE,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
429
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
430
  # family <- "gaussian"; nfolds <- 5; foldid <- NULL; type.measure <- "deviance"; alpha.base=1;alpha.meta=0;spls=FALSE;mnorm=FALSE;sier=FALSE;mrce=FALSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
431
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
432 433 434
  n <- nrow(Y)
  q <- ncol(Y)
  p <- ncol(X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
435
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
436 437 438 439 440
  #--- fold identifiers ---
  if(is.null(foldid.ext)){
    foldid.ext <- palasso:::.folds(y=Y[,1],nfolds=nfolds.ext) # temporary Y[,1]
  } else {
    nfolds.ext <- length(unique(foldid.ext))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
441 442
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
443 444 445 446 447
  #--- family ---
  if(length(family)==1){
    family <- rep(family,times=q)
  } else if(length(family)!=q){
    stop("Invalid argument family",call.=FALSE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
448 449
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
450
  #--- cross-validated predictions ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
451
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
452 453 454
  models <- c("base","meta","spls"[spls],"mnorm"[mnorm],"sier"[sier],"mrce"[mrce],"none")
  pred <- lapply(X=models,function(x) matrix(data=NA,nrow=n,ncol=q))
  names(pred) <- models
Armin Rauschenberger's avatar
Armin Rauschenberger committed
455
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503
  for(i in seq_len(nfolds.ext)){
    
    Y0 <- Y[foldid.ext!=i,,drop=FALSE]
    Y1 <- Y[foldid.ext==i,,drop=FALSE]
    X0 <- X[foldid.ext!=i,,drop=FALSE]
    X1 <- X[foldid.ext==i,,drop=FALSE]
    if(is.null(foldid.int)){
      foldid <- palasso:::.folds(y=Y0[,1],nfolds=nfolds.int) # temporary Y0[,1]
    } else {
      foldid <- foldid.int[foldid.ext!=i]
    }
    
    # base and meta learners
    fit <- mixnet(Y=Y0,X=X0,family=family,type.measure=type.measure,alpha.base=alpha.base,alpha.meta=alpha.meta,foldid=foldid) # add ,...
    temp <- predict.mixnet(fit,newx=X1)
    pred$base[foldid.ext==i,] <- temp$base
    pred$meta[foldid.ext==i,] <- temp$meta
    
    # other learners
    cond <- apply(X0,2,function(x) stats::sd(x)!=0)
    x0 <- X0[,cond]
    x1 <- X1[,cond]
    y0 <- apply(X=Y0,MARGIN=2,FUN=function(x) ifelse(is.na(x),sample(x[!is.na(x)],size=1),x))
    all(Y0==y0,na.rm=TRUE)
    if(mnorm){
      net <- glmnet::cv.glmnet(x=X0,y=y0,family="mgaussian",foldid=foldid,...) # ellipsis
      pred$mnorm[foldid.ext==i,] <- glmnet::predict.cv.glmnet(object=net,newx=X1,s="lambda.min",type="response")
    }
    if(spls){
      cv.spls <- spls::cv.spls(x=x0,y=y0,fold=nfolds.int,K=seq_len(10),
                               eta=seq(from=0.1,to=0.9,by=0.1),scale.x=FALSE,plot.it=FALSE)
      mspls <- spls::spls(x=x0,y=y0,K=cv.spls$K.opt,cv.spls$eta.opt,scale.x=FALSE)
      pred$spls[foldid.ext==i,] <- spls::predict.spls(object=mspls,newx=x1,type="fit")
    }
    if(sier){
      object <- SiER::cv.SiER(X=X0,Y=y0,K.cv=10)
      pred$sier[foldid.ext==i,] <- SiER::pred.SiER(cv.fit=object,X.new=X1)
    }
    if(mrce){
      lam1 <- rev(10^seq(from=-2,to=0,by=0.5))
      lam2 <- rev(10^seq(from=-2,to=0,by=0.5))
      object <- MRCE::mrce(X=x0,Y=y0,lam1.vec=lam1,lam2.vec=lam2,method="cv")
      pred$mrce[foldid.ext==i,] <- object$muhat + x1 %*% object$Bhat
    }
    
    pred$none[foldid.ext==i,] <- matrix(colMeans(Y0,na.rm=TRUE),nrow=sum(foldid.ext==i),ncol=ncol(Y0),byrow=TRUE)
    
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
504
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
505 506 507 508 509 510 511 512 513
  #--- cross-validated loss ---
  loss <- matrix(data=NA,nrow=length(models),ncol=ncol(Y),
                 dimnames=list(models,colnames(Y)))
  for(j in models){
    for(i in seq_len(q)){
      cond <- !is.na(Y[,i])
      loss[j,i] <- palasso:::.loss(y=Y[cond,i],fit=pred[[j]][cond,i],family=family[i],type.measure=type.measure)[[1]]
    }
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
514
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
515 516 517
  #--- model refit ---
  #fit <- mixnet(Y=Y,X=X,family=family,type.measure=type.measure,alpha.base=alpha.base,alpha.meta=alpha.meta) # add ,...
  #list <- list(loss=loss,fit=fit)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
518
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
519
  return(loss)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
520
}