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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
2 3 4 5
#.loss <- get(".loss",envir=asNamespace("palasso"))
#.folds <- get(".folds",envir=asNamespace("palasso"))
#.check <- get(".check",envir=asNamespace("cornet"))

Armin Rauschenberger's avatar
Armin Rauschenberger committed
6
#--- Main function -------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
7

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

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
204 205 206 207
.link.function <- function(x,family){
  if(family=="gaussian"){
    return(x)
  } else if(family=="binomial"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
208 209
    if(any(x<0|x>1)){stop("Invalid!",call.=FALSE)}
    return(log(x/(1-x)))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
210
  } else if(family=="poisson"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
211
    if(any(x<0)){stop("Invalid!",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
212 213 214 215
    return(log(x))
  } else {
    stop("Family not implemented.",call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
216 217
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
218
#--- Methods for class "joinet" -----------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
219

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
291 292
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
293
#' Extract Coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
294 295
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
296
#' Extracts pooled coefficients.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
297 298
#' (The meta learners linearly combines
#' the coefficients from the base learners.)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
299
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
300
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
301
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
302 303
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
304
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
305 306
#' 
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
307 308
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
309
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
310
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
311 312
#' coef <- coef(object)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
313
coef.joinet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
314 315 316 317 318 319 320 321 322 323 324
  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()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
325
  weights <- weights.joinet(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
  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
346
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
347
  return(pool)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
348 349
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
350 351
#' @export
#' @importFrom stats weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
352
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
353
#' Extract Weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
354 355
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
356 357
#' Extracts coefficients from the meta learner,
#' i.e. the weights for the base learners.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
358
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
359
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
360
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
361
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
362 363
#' @param ...
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
364
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
365 366 367 368
#' @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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
369
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
370 371
#' weights(object)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
372
weights.joinet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
373 374 375 376 377 378 379
  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
380 381
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
382 383
print.joinet <- function(x,...){
  cat(paste0("joinet object"),"\n")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
384 385 386 387
}

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
388
#' @export
Armin Rauschenberger's avatar
Armin Rauschenberger committed
389
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
390
#' Model comparison
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391 392
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
393
#' Compares univariate and multivariate regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
394
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
395
#' @inheritParams joinet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
396
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
397 398
#' @param nfolds.ext
#' number of external folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
399
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
400 401
#' @param nfolds.int
#' number of internal folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
402
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
403 404 405 406 407
#' @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
408
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
409 410 411 412 413
#' @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
414
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
415 416 417
#' @param mnorm,spls,sier,mrce
#' experimental arguments\strong{:}
#' logical (install packages \code{spls}, \code{SiER}, or \code{MRCE})
Armin Rauschenberger's avatar
Armin Rauschenberger committed
418
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
419 420
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}} and \code{\link[glmnet]{cv.glmnet}}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
421
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
422
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
423 424 425
#' n <- 40; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
426
#' cv.joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
427
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
428
cv.joinet <- 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 431 432
  n <- nrow(Y)
  q <- ncol(Y)
  p <- ncol(X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
433
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
434 435 436 437 438
  #--- 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
439 440
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
441 442 443 444 445
  #--- 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
446 447
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
448
  #--- cross-validated predictions ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
449
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
450 451 452
  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
453
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
454 455 456 457 458 459 460 461 462 463 464 465 466
  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
Armin Rauschenberger's avatar
Armin Rauschenberger committed
467 468
    fit <- joinet(Y=Y0,X=X0,family=family,type.measure=type.measure,alpha.base=alpha.base,alpha.meta=alpha.meta,foldid=foldid) # add ,...
    temp <- predict.joinet(fit,newx=X1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
469 470 471 472 473 474 475 476 477
    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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
478
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
    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
503
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
504 505 506 507 508 509 510 511 512
  #--- 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
513
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
514
  #--- model refit ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
515
  #fit <- joinet(Y=Y,X=X,family=family,type.measure=type.measure,alpha.base=alpha.base,alpha.meta=alpha.meta) # add ,...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
516
  #list <- list(loss=loss,fit=fit)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
517
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
518
  return(loss)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
519
}