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

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
200 201 202 203 204 205 206 207 208 209
.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
210 211
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
212
#--- Methods for class "joinet" -----------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
213

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

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
344 345
#' @export
#' @importFrom stats weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
346
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
347
#' Extract Weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
348 349
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
350 351
#' Extracts coefficients from the meta learner,
#' i.e. the weights for the base learners.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
352
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
353
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
354
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
355
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
356 357
#' @param ...
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
358
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
359 360 361 362
#' @examples
#' n <- 30; q <- 2; p <- 20
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
363
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
364 365
#' weights(object)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
366
weights.joinet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
367 368 369 370 371 372 373
  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
374 375
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
376 377
print.joinet <- function(x,...){
  cat(paste0("joinet object"),"\n")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
378 379 380 381
}

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

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