functions.R 15.6 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
#' @param alpha.meta
#' elastic net mixing parameter for meta learner\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
47
#' 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 62 63 64 65 66
#' elastic net mixing parameters:
#' \code{alpha.base} controls input-output effects,
#' \code{alpha.meta} controls output-output effects;
#' ridge (\eqn{0}) renders dense models,
#' lasso (\eqn{1}) renders sparse models
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
67
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
68 69
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
70
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
71
#' object <- mixnet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
72
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
73
mixnet <- function(Y,X,family="gaussian",nfolds=10,foldid=NULL,type.measure="deviance",alpha.base=0,alpha.meta=0,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
74
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
75
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
76
  # family <- "gaussian"; nfolds <- 10; foldid <- NULL; type.measure <- "deviance"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
77
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
78
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
79 80 81 82 83 84 85 86 87 88 89 90
  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
91 92
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
93 94 95 96 97 98 99 100 101 102
  #--- 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
103
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
104
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
105 106 107 108
  #--- 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
109
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
110
    nfolds <- length(unique(foldid))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
111 112
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
113 114 115 116 117 118 119 120 121 122
  #--- 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
123
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
124 125 126 127
  #--- 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
128
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
129 130
  
  #--- base cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
131
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
132 133
    Y0 <- Y[foldid!=k,,drop=FALSE]
    Y1 <- Y[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
134 135
    X0 <- X[foldid!=k,,drop=FALSE]
    X1 <- X[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
136 137 138 139 140 141 142
    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
143
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
144 145
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
146 147 148 149 150 151 152 153
  #--- 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
154
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
155 156 157 158 159
  #--- 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
160
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
161 162 163 164 165 166 167 168 169 170 171 172 173
  #--- 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
174 175
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
176
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
177 178 179 180 181
  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
182 183
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
184 185 186 187 188 189 190 191 192 193 194
.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
195

Armin Rauschenberger's avatar
Armin Rauschenberger committed
196 197 198 199 200 201 202 203 204 205
.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
206 207
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
208 209
#--- Methods for class "mixnet" -----------------------------------------------

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
279 280
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
281
#' Extract Coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
282 283
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
284 285
#' Extracts pooled coefficients.
#' (The meta learners weights the coefficients from the base learners.)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
286
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
287
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
288
#' \link[mixnet]{mixnet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
289 290
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
291
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
292 293
#' 
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
294 295
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
296
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
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 328 329 330 331 332
#' 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
333
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
334
  return(pool)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
335 336
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
337 338
#' @export
#' @importFrom stats weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
339
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
340
#' Extract Weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
341 342
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
343 344
#' Extracts coefficients from the meta learner,
#' i.e. the weights for the base learners.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
345
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
346 347
#' @param object
#' \link[mixnet]{mixnet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
348
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
349 350
#' @param ...
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
351
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367
#' @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
368 369
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384
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
385
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
386
#' Model comparison
Armin Rauschenberger's avatar
Armin Rauschenberger committed
387 388
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
389
#' Compares univariate and multivariate regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
390
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391
#' @inheritParams mixnet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
392
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
393 394
#' @param nfolds.ext
#' number of external folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
395
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
396 397
#' @param nfolds.int
#' number of internal folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
398
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
399 400 401 402 403
#' @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
404
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
405 406 407 408 409
#' @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
410
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
411 412
#' @param mnorm
#' multivariate normal regression\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
413 414
#' logical
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
415 416
#' @param spls
#' sparse partial least squares\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
417 418
#' logical
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
419 420 421
#' @param sier
#' high-dimensional multivariate regression\strong{:}
#' logical
Armin Rauschenberger's avatar
Armin Rauschenberger committed
422
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
423 424 425
#' @param mrce
#' multivariate regression with covariance estimation\strong{:}
#' logical
Armin Rauschenberger's avatar
Armin Rauschenberger committed
426
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
427 428
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}} and \code{\link[glmnet]{cv.glmnet}}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
429
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
430
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
431
#' NA
Armin Rauschenberger's avatar
Armin Rauschenberger committed
432
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
433
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
434
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
435
  # 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
436
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
437 438 439
  n <- nrow(Y)
  q <- ncol(Y)
  p <- ncol(X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
440
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
441 442 443 444 445
  #--- 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
446 447
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
448 449 450 451 452
  #--- 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
453 454
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
455
  #--- cross-validated predictions ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
456
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
457 458 459
  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
460
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
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 504 505 506 507 508
  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
509
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
510 511 512 513 514 515 516 517 518
  #--- 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
519
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
520 521 522
  #--- 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
523
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
524
  return(loss)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
525
}