functions.R 18.3 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
#' Armin Rauschenberger, Enrico Glaab (2019)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
59 60
#' "joinet: predicting correlated outcomes jointly
#' to improve clinical prognosis"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
61
#' \emph{Manuscript in preparation}.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
62 63
#' 
#' @details
Armin Rauschenberger's avatar
Armin Rauschenberger committed
64
#' \strong{correlation:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
65 66
#' 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
67
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
68
#' \strong{elastic net:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
69 70
#' \code{alpha.base} controls input-output effects,
#' \code{alpha.meta} controls output-output effects;
Armin Rauschenberger's avatar
Armin Rauschenberger committed
71 72 73 74 75 76 77 78 79 80 81
#' lasso renders sparse models (\code{alpha}\eqn{=1}),
#' ridge renders dense models (\code{alpha}\eqn{=0})
#' 
#' @return
#' This function returns an object of class \code{joinet}.
#' Available methods include
#' \code{\link[=predict.joinet]{predict}},
#' \code{\link[=coef.joinet]{coef}},
#' and \code{\link[=weights.joinet]{weights}}.
#' The slots \code{base} and \code{meta} each contain
#' \eqn{q} \code{\link[glmnet]{cv.glmnet}}-like objects.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
82
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
83
#' @seealso
Armin Rauschenberger's avatar
Armin Rauschenberger committed
84
#' \code{\link{cv.joinet}}, vignette
Armin Rauschenberger's avatar
Armin Rauschenberger committed
85
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
86
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
87
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
88
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
89
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
90
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
91
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
92 93 94
#' \dontrun{
#' browseVignettes("joinet") # further examples}
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
95
joinet <- 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
96
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
97
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
98
  # family <- "gaussian"; nfolds <- 10; foldid <- NULL; type.measure <- "deviance"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
99
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
100
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
101 102 103
  Y <- as.matrix(Y)
  X <- as.matrix(X)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
104 105 106
  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")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
107
  #cornet:::.check(x=family,type="vector",values=c("gaussian","binomial","poisson"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
108 109 110 111 112 113 114 115
  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
116 117
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
118 119 120 121 122 123 124 125 126 127
  #--- 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
128
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
129
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
130 131 132 133
  #--- 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
134
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
135
    nfolds <- length(unique(foldid))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
136 137
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
138 139 140 141 142 143 144 145 146 147
  #--- 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
148
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
149 150 151 152
  #--- 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
153
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
154 155
  
  #--- base cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
156
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
157 158
    Y0 <- Y[foldid!=k,,drop=FALSE]
    Y1 <- Y[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
159 160
    X0 <- X[foldid!=k,,drop=FALSE]
    X1 <- X[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
161 162 163 164 165 166 167
    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
168
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
169 170
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
171 172 173 174 175 176 177 178
  #--- 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
179
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
180 181 182 183 184
  #--- 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
185
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
186 187 188 189 190 191 192 193 194 195 196 197 198
  #--- 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
199 200
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
201
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
202 203 204
  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
205
  class(list) <- "joinet"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
206
  return(list)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
207 208
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
209 210 211 212 213 214 215 216 217 218 219
.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
220

Armin Rauschenberger's avatar
Armin Rauschenberger committed
221 222 223 224
.link.function <- function(x,family){
  if(family=="gaussian"){
    return(x)
  } else if(family=="binomial"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
225 226
    if(any(x<0|x>1)){stop("Invalid!",call.=FALSE)}
    return(log(x/(1-x)))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
227
  } else if(family=="poisson"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
228
    if(any(x<0)){stop("Invalid!",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
229 230 231 232
    return(log(x))
  } else {
    stop("Family not implemented.",call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
233 234
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
235
#--- Methods for class "joinet" -----------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
236

Armin Rauschenberger's avatar
Armin Rauschenberger committed
237 238
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
239
#' Make Predictions
Armin Rauschenberger's avatar
Armin Rauschenberger committed
240 241
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
242
#' Predicts outcome from features with stacked model.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
243
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
244
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
245
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
246
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
247 248 249 250
#' @param newx
#' covariates\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
#' and \eqn{p} columns (variables)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
251
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
252 253
#' @param type
#' character "link" or "response"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
254 255
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
256
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
257
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
258 259 260 261 262
#' @return 
#' This function returns predictions from base and meta learners.
#' The slots \code{base} and \code{meta} each contain a matrix
#' with \eqn{n} rows (samples) and \eqn{q} columns (variables).
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
263
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
264
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
265
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
266 267 268
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
#' object <- joinet(Y=Y,X=X)
#' predict(object,newx=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
269
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
270
predict.joinet <- function(object,newx,type="response",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
271
  if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
272
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
273
  x <- object; rm(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
274
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
275
  newx <- as.matrix(newx)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
276
  cornet:::.check(x=newx,type="matrix")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
277
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
278 279 280
  q <- length(x$base)
  n <- nrow(newx)
  base <- meta <- matrix(NA,nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
281
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
282 283 284 285 286 287
  # 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
288
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
289 290 291 292 293
  # 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
294
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
295
  list <- list(base=base,meta=meta)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
296
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
297 298 299 300 301 302 303
  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
304 305
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
306 307 308
  list <- list(base=base,meta=meta)
  return(list)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
309 310
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
311 312
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
313
#' Extract Coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
314 315
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
316
#' Extracts pooled coefficients.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
317 318
#' (The meta learners linearly combines
#' the coefficients from the base learners.)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
319
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
320
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
321
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
322 323
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
324
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
325
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
326 327 328 329 330 331 332
#' @return
#' This function returns the pooled coefficients.
#' The slot \code{alpha} contains the intercepts
#' in a vector of length \eqn{q},
#' and the slot \code{beta} contains the slopes
#' in a matrix with \eqn{p} rows (inputs) and \eqn{q} columns.
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
333
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
334
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
335
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
336
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
337
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
338 339
#' coef <- coef(object)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
340
coef.joinet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
341 342 343 344 345 346 347 348 349 350 351
  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
352
  weights <- weights.joinet(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372
  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
373
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
374
  return(pool)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
375 376
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
377 378
#' @export
#' @importFrom stats weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
379
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
380
#' Extract Weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
381 382
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
383 384
#' Extracts coefficients from the meta learner,
#' i.e. the weights for the base learners.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
385
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
386
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
387
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
388
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
389 390
#' @param ...
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
392 393 394 395 396 397 398 399
#' @return
#' This function returns a matrix with
#' \eqn{1+q} rows and \eqn{q} columns.
#' The first row contains the intercepts,
#' and the other rows contain the slopes,
#' which are the effects of the outcomes
#' in the row on the outcomes in the column.
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
400
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
401
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
402
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
403
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
404
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
405 406
#' weights(object)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
407
weights.joinet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
408 409 410 411 412 413 414
  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
415 416
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
417 418
print.joinet <- function(x,...){
  cat(paste0("joinet object"),"\n")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
419 420 421 422
}

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
423
#' @export
Armin Rauschenberger's avatar
Armin Rauschenberger committed
424
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
425
#' Model comparison
Armin Rauschenberger's avatar
Armin Rauschenberger committed
426 427
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
428
#' Compares univariate and multivariate regression.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
429
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
430
#' @inheritParams joinet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
431
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
432 433
#' @param nfolds.ext
#' number of external folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
434
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
435 436
#' @param nfolds.int
#' number of internal folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
437
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
438 439 440 441 442
#' @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
443
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
444 445 446 447 448
#' @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
449
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
450 451
#' @param mnorm,spls,sier,mrce
#' experimental arguments\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
452
#' logical (requires packages \code{spls}, \code{SiER}, or \code{MRCE})
Armin Rauschenberger's avatar
Armin Rauschenberger committed
453
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
454
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
455 456 457 458 459
#' further arguments passed to \code{\link[glmnet]{glmnet}}
#' and \code{\link[glmnet]{cv.glmnet}}
#' 
#' @return 
#' This function returns a matrix with \eqn{q} columns,
Armin Rauschenberger's avatar
Armin Rauschenberger committed
460 461 462
#' including the cross-validated loss from the univariate models
#' (\code{base}), the multivariate models (\code{meta}),
#' and the intercept-only models (\code{none}).
Armin Rauschenberger's avatar
Armin Rauschenberger committed
463
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
464
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
465
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
466
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
467 468 469 470 471
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
#' cv.joinet(Y=Y,X=X)
#' 
#' \dontrun{
#' # correlated features
Armin Rauschenberger's avatar
Armin Rauschenberger committed
472
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
473 474 475 476
#' mu <- rep(0,times=p)
#' Sigma <- 0.90^abs(col(diag(p))-row(diag(p)))
#' X <- MASS::mvrnorm(n=n,mu=mu,Sigma=Sigma)
#' mu <- rowSums(X[,sample(seq_len(p),size=5)])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
477 478 479
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=mu))
#' #Y <- t(MASS::mvrnorm(n=q,mu=mu,Sigma=diag(n)))
#' cv.joinet(Y=Y,X=X)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
480 481 482
#' 
#' \dontrun{
#' # other distributions
Armin Rauschenberger's avatar
Armin Rauschenberger committed
483
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
484
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
485
#' eta <- rowSums(X[,1:5])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
486 487 488
#' Y <- replicate(n=q,expr=rbinom(n=n,size=1,prob=1/(1+exp(-eta))))
#' cv.joinet(Y=Y,X=X,family="binomial")
#' Y <- replicate(n=q,expr=rpois(n=n,lambda=exp(scale(eta))))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507
#' cv.joinet(Y=Y,X=X,family="poisson")}
#' 
#' \dontrun{
#' # uncorrelated outcomes
#' n <- 50; p <- 100; q <- 3
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
#' y <- rnorm(n=n,mean=rowSums(X[,1:5]))
#' Y <- cbind(y,matrix(rnorm(n*(q-1)),nrow=n,ncol=q-1))
#' cv.joinet(Y=Y,X=X)}
#' 
#' \dontrun{
#' # sparse and dense models
#' n <- 50; p <- 100; q <- 3
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
#' set.seed(1) # fix folds
#' cv.joinet(Y=Y,X=X,alpha.base=1) # lasso
#' set.seed(1)
#' cv.joinet(Y=Y,X=X,alpha.base=0) # ridge}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
508
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
509
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
510
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
511 512 513
  n <- nrow(Y)
  q <- ncol(Y)
  p <- ncol(X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
514
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
515 516 517 518 519
  #--- 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
520 521
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
522 523 524 525 526
  #--- 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
527 528
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
529
  #--- cross-validated predictions ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
530
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
531
  models <- c("base","meta","mnorm"[mnorm],"spls"[spls],"sier"[sier],"mrce"[mrce],"none")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
532 533
  pred <- lapply(X=models,function(x) matrix(data=NA,nrow=n,ncol=q))
  names(pred) <- models
Armin Rauschenberger's avatar
Armin Rauschenberger committed
534
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
535 536 537 538 539 540 541 542 543 544 545 546 547
  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
548 549
    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
550 551 552 553 554 555 556 557 558
    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
559
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
560 561 562 563 564 565 566 567 568 569 570 571 572 573 574
    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){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
575
      # bug?
Armin Rauschenberger's avatar
Armin Rauschenberger committed
576 577 578 579 580 581 582 583 584
      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
585
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
586 587 588 589 590 591 592 593 594
  #--- 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
595
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
596
  #--- model refit ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
597
  #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
598
  #list <- list(loss=loss,fit=fit)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
599
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
600
  return(loss)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
601
}