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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
2
#--- import unexported functions:
Armin Rauschenberger's avatar
Armin Rauschenberger committed
3
# FUNCTION <- get("FUNCTION",envir=asNamespace("PACKAGE"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
4

Armin Rauschenberger's avatar
Armin Rauschenberger committed
5
6
7
#--- deactivate on Solaris:
# if(!grepl('SunOS',Sys.info()['sysname'])){}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
8
#--- Main function -------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
9

Armin Rauschenberger's avatar
Armin Rauschenberger committed
10
11
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
12
#' Multivariate Elastic Net Regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
13
14
#' 
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
15
16
17
18
19
#' Implements multivariate elastic net regression.
#'  
#' @param Y
#' outputs\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
20
#' and \eqn{q} columns (variables)
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
#' @param alpha.meta
Armin Rauschenberger's avatar
Armin Rauschenberger committed
51
#' elastic net mixing parameter for meta learners\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
57
58
59
60
61
62
#' @param constraint
#' non-negativity constraints\strong{:}
#' logical (see details)
#' 
#' @param weight
#' inclusion/exclusion of variables\strong{:}
#' logical matrix with \eqn{q} rows and \eqn{p} columns,
#' or \code{NULL} (see details)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
63
64
65
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}}
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
66
#' @references 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
67
#' Armin Rauschenberger, Enrico Glaab (2021)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
68
#' "Predicting correlated outcomes from molecular data"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
69
#' \emph{Manuscript in preparation}.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
70
71
#' 
#' @details
Armin Rauschenberger's avatar
Armin Rauschenberger committed
72
73
74
75
76
77
78
79
80
81
82
83
84
#' \strong{non-negativity constraints:}
#' If it is reasonable to assume that the outcomes
#' are \emph{positively} correlated
#' (potentially after changing the sign of some outcomes)
#' we recommend to set \code{constraint=TRUE}.
#' Then non-negativity constraints are imposed on the meta learner.
#' 
#' \strong{inclusion/exclusion of variables:}
#' The entry in the \eqn{j}th column and the \eqn{k}th row
#' indicates whether the \eqn{j}th feature may be used for 
#' modelling the \eqn{k}th outcome
#' (where \eqn{0} means \code{FALSE} and
#' \eqn{1} means \code{TRUE}).
Armin Rauschenberger's avatar
Armin Rauschenberger committed
85
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
86
#' \strong{elastic net:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
87
88
#' \code{alpha.base} controls input-output effects,
#' \code{alpha.meta} controls output-output effects;
Armin Rauschenberger's avatar
Armin Rauschenberger committed
89
90
91
92
93
94
95
96
97
98
99
#' 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
100
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
101
#' @seealso
Armin Rauschenberger's avatar
Armin Rauschenberger committed
102
#' \code{\link{cv.joinet}}, vignette
Armin Rauschenberger's avatar
Armin Rauschenberger committed
103
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
104
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
105
#' \dontshow{
Armin Rauschenberger's avatar
Armin Rauschenberger committed
106
#' if(!grepl('SunOS',Sys.info()['sysname'])){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
107
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
108
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
109
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
110
#' object <- joinet(Y=Y,X=X)}}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
111
112
113
114
#' \dontrun{
#' 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])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
115
#' object <- joinet(Y=Y,X=X)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
116
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
117
118
119
#' \dontrun{
#' browseVignettes("joinet") # further examples}
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
120
121
joinet <- function(Y,X,family="gaussian",nfolds=10,foldid=NULL,type.measure="deviance",alpha.base=1,alpha.meta=1,constraint=TRUE,weight=NULL,...){
  # IMPLEMENT CODE FOR CONSTRAINT AND WEIGHT!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
122
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
123
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
124
  # family <- "gaussian"; nfolds <- 10; foldid <- NULL; type.measure <- "deviance"; alpha.base <- alpha.meta <- 1; constraint <- TRUE; weight <- NULL
Armin Rauschenberger's avatar
Armin Rauschenberger committed
125
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
126
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
127
128
129
  Y <- as.matrix(Y)
  X <- as.matrix(X)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
130
  cornet:::.check(x=Y,type="matrix",miss=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
131
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
132
133
134
135
136
137
138
  if(constraint){
    for(i in 1:(ncol(Y)-1)){
      for(j in i:ncol(Y)){
        cor <- stats::cor.test(Y[,i],Y[,j],use="pairwise.complete.obs")
        if(cor$statistic<0 & cor$p.value<0.05){
          warning(paste("Columns",i,"and",j,"are negatively correlated. Consider using constraint=FALSE."),call.=FALSE)
        }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
139
140
141
142
143
      }
    }
  }
  
  #if(any(stats::cor(Y,use="pairwise.complete.obs")<0,na.rm=TRUE)){warning("Negative correlation!",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
144
  cornet:::.check(x=X,type="matrix")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
145
  #cornet:::.check(x=family,type="vector",values=c("gaussian","binomial","poisson"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
146
147
148
149
150
151
  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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
152
  cornet:::.check(x=weight,type="matrix",min=0,max=1,null=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
153
154
  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
155
156
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
157
158
159
160
161
  #--- dimensionality ---
  n <- nrow(Y)
  q <- ncol(Y)
  p <- ncol(X)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
162
163
164
165
166
167
  if(is.null(weight)){
    pf <- matrix(1,nrow=q,ncol=p)
  } else {
    pf <- 1/weight
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
168
169
170
171
172
  #--- 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
173
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
174
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
175
176
177
178
  #--- 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
179
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
180
    nfolds <- length(unique(foldid))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
181
182
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
183
184
185
186
187
188
  #--- 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}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
189
    base[[i]]$glmnet.fit <- glmnet::glmnet(y=Y[cond,i],x=X[cond,],family=family[i],alpha=alpha.base,penalty.factor=pf[i,],...) # ellipsis
Armin Rauschenberger's avatar
Armin Rauschenberger committed
190
191
192
    base[[i]]$lambda <- base[[i]]$glmnet.fit$lambda
    nlambda[i] <- length(base[[i]]$glmnet.fit$lambda)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
193
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
194
195
196
197
  #--- 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
198
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
199
200
  
  #--- base cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
201
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
202
203
    Y0 <- Y[foldid!=k,,drop=FALSE]
    Y1 <- Y[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
204
205
    X0 <- X[foldid!=k,,drop=FALSE]
    X1 <- X[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
206
207
208
    for(i in seq_len(q)){
      cond <- !is.na(Y0[,i])
      #if(sum(cond)==0){next}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
209
      object <- glmnet::glmnet(y=Y0[cond,i],x=X0[cond,],family=family[i],alpha=alpha.base,penalty.factor=pf[i,],...) # ellipsis
Armin Rauschenberger's avatar
Armin Rauschenberger committed
210
211
212
      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
213
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
214
215
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
216
217
218
219
220
221
222
  #--- 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
223
    class(base[[i]]) <- "cv.glmnet" # trial 2020-01-10
Armin Rauschenberger's avatar
Armin Rauschenberger committed
224
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
225
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
226
227
228
229
230
  #--- 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
231
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
232
233
234
235
236
  #--- 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,],
Armin Rauschenberger's avatar
Armin Rauschenberger committed
237
                                   lower.limits=ifelse(constraint,0,-Inf), # important: 0 (was lower.limits=0)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
238
239
240
241
242
243
244
                                   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
245
246
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
247
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
248
249
250
  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
251
  class(list) <- "joinet"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
252
  return(list)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
253
254
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
255
.mean.function <- function(x,family){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
256
  if(family %in% c("gaussian","cox")){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
257
258
259
260
261
262
263
264
265
    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
266

Armin Rauschenberger's avatar
Armin Rauschenberger committed
267
.link.function <- function(x,family){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
268
  if(family %in% c("gaussian","cox")){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
269
270
    return(x)
  } else if(family=="binomial"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
271
272
    if(any(x<0|x>1)){stop("Invalid!",call.=FALSE)}
    return(log(x/(1-x)))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
273
  } else if(family=="poisson"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
274
    if(any(x<0)){stop("Invalid!",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
275
276
277
278
    return(log(x))
  } else {
    stop("Family not implemented.",call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
279
280
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
281
#--- Methods for class "joinet" -----------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
282

Armin Rauschenberger's avatar
Armin Rauschenberger committed
283
284
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
285
#' Make Predictions
Armin Rauschenberger's avatar
Armin Rauschenberger committed
286
287
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
288
#' Predicts outcome from features with stacked model.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
289
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
290
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
291
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
292
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
293
294
295
296
#' @param newx
#' covariates\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
#' and \eqn{p} columns (variables)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
297
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
298
299
#' @param type
#' character "link" or "response"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
300
301
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
302
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
303
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
304
305
306
307
308
#' @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
309
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
310
#' \dontshow{
Armin Rauschenberger's avatar
Armin Rauschenberger committed
311
#' if(!grepl('SunOS',Sys.info()['sysname'])){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
312
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
313
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
314
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
315
316
#' Y[,1] <- 1*(Y[,1]>median(Y[,1]))
#' object <- joinet(Y=Y,X=X,family=c("binomial","gaussian","gaussian"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
317
318
319
320
321
322
323
324
#' predict(object,newx=X)}}
#' \dontrun{
#' 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])))
#' Y[,1] <- 1*(Y[,1]>median(Y[,1]))
#' object <- joinet(Y=Y,X=X,family=c("binomial","gaussian","gaussian"))
#' predict(object,newx=X)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
325
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
326
predict.joinet <- function(object,newx,type="response",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
327
  if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
328
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
329
  x <- object; rm(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
330
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
331
  newx <- as.matrix(newx)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
332
  cornet:::.check(x=newx,type="matrix")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
333
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
334
335
336
  q <- length(x$base)
  n <- nrow(newx)
  base <- meta <- matrix(NA,nrow=n,ncol=q)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
337
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
338
339
  # base learners
  for(i in seq_len(q)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
340
341
    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"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
342
343
    # check whether fine for "binomial" family
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
344
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
345
346
347
348
349
  # 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
350
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
351
  list <- list(base=base,meta=meta)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
352
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
353
354
355
356
357
358
359
  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
360
361
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
362
363
364
  list <- list(base=base,meta=meta)
  return(list)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
365
366
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
367
368
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
369
#' Extract Coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
370
371
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
372
#' Extracts pooled coefficients.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
373
374
#' (The meta learners linearly combines
#' the coefficients from the base learners.)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
375
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
376
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
377
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
378
379
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
380
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
381
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
382
383
384
385
386
387
388
#' @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
389
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
390
#' \dontshow{
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391
#' if(!grepl('SunOS',Sys.info()['sysname'])){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
392
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
393
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
394
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
395
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
396
397
398
399
400
401
402
#' coef <- coef(object)}}
#' \dontrun{
#' 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])))
#' object <- joinet(Y=Y,X=X)
#' coef <- coef(object)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
403
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
404
coef.joinet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
405
406
407
408
  if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
  
  # base coefficients
  base <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
409
  coef <- sapply(object$base,function(x) stats::coef(object=x$glmnet.fit,s=x$lambda.min))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
410
411
412
413
414
415
  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
416
  weights <- weights.joinet(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
  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
437
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
438
  return(pool)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
439
440
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
441
442
#' @export
#' @importFrom stats weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
443
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
444
#' Extract Weights
Armin Rauschenberger's avatar
Armin Rauschenberger committed
445
446
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
447
448
#' Extracts coefficients from the meta learner,
#' i.e. the weights for the base learners.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
449
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
450
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
451
#' \link[joinet]{joinet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
452
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
453
454
#' @param ...
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
455
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
456
457
458
459
460
461
462
463
#' @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
464
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
465
#' \dontshow{
Armin Rauschenberger's avatar
Armin Rauschenberger committed
466
#' if(!grepl('SunOS',Sys.info()['sysname'])){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
467
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
468
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
469
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
470
#' object <- joinet(Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
471
472
473
474
475
476
477
#' weights(object)}}
#' \dontrun{
#' 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])))
#' object <- joinet(Y=Y,X=X)
#' weights(object)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
478
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
479
weights.joinet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
480
481
  if(length(list(...))!=0){warning("Ignoring argument.",call.=FALSE)}
  x <- object$meta
Armin Rauschenberger's avatar
Armin Rauschenberger committed
482
  coef <- lapply(object$meta,function(x) stats::coef(object=x,s=x$lambda.min))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
483
484
485
486
  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
487
488
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
489
490
print.joinet <- function(x,...){
  cat(paste0("joinet object"),"\n")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
491
492
493
494
}

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
495
#' @export
Armin Rauschenberger's avatar
Armin Rauschenberger committed
496
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
497
#' Model comparison
Armin Rauschenberger's avatar
Armin Rauschenberger committed
498
499
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
500
#' Compares univariate and multivariate regression.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
501
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
502
#' @inheritParams joinet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
503
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
504
505
#' @param nfolds.ext
#' number of external folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
506
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
507
508
#' @param nfolds.int
#' number of internal folds
Armin Rauschenberger's avatar
Armin Rauschenberger committed
509
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
510
511
512
513
514
#' @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
515
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
516
517
518
519
520
#' @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
521
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
522
#' @param compare
Armin Rauschenberger's avatar
Armin Rauschenberger committed
523
#' experimental arguments\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
524
525
526
#' character vector with entries "mnorm", "spls", "mrce",
#' "sier", "mtps", "rmtl", "gpm" and others
#' (requires packages \code{spls}, \code{MRCE}, \code{SiER}, \code{MTPS}, \code{RMTL} or \code{GPM})
Armin Rauschenberger's avatar
Armin Rauschenberger committed
527
528
529
530
#' 
#' @param mice
#' missing data imputation\strong{:}
#' logical (\code{mice=TRUE} requires package \code{mice})
Armin Rauschenberger's avatar
Armin Rauschenberger committed
531
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
532
#' @param cvpred
Armin Rauschenberger's avatar
Armin Rauschenberger committed
533
#' return cross-validated predictions: logical
Armin Rauschenberger's avatar
Armin Rauschenberger committed
534
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
535
536
537
538
#' @param times
#' measure computation time\strong{:}
#' logical
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
539
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
540
541
542
543
544
#' 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
545
546
547
#' 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
548
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
549
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
550
#' \dontshow{
Armin Rauschenberger's avatar
Armin Rauschenberger committed
551
#' if(!grepl('SunOS',Sys.info()['sysname'])){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
552
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
553
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
554
#' Y <- replicate(n=q,expr=rnorm(n=n,mean=rowSums(X[,1:5])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
555
556
557
558
559
560
#' cv.joinet(Y=Y,X=X)}}
#' \dontrun{
#' 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])))
#' cv.joinet(Y=Y,X=X)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
561
562
563
#' 
#' \dontrun{
#' # correlated features
Armin Rauschenberger's avatar
Armin Rauschenberger committed
564
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
565
566
567
568
#' 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
569
570
571
#' 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
572
573
574
#' 
#' \dontrun{
#' # other distributions
Armin Rauschenberger's avatar
Armin Rauschenberger committed
575
#' n <- 50; p <- 100; q <- 3
Armin Rauschenberger's avatar
Armin Rauschenberger committed
576
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
577
#' eta <- rowSums(X[,1:5])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
578
579
580
#' 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
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
#' 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
600
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
601
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=1,compare=FALSE,mice=FALSE,cvpred=FALSE,times=FALSE,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
602
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
603
  if(FALSE){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
604
  fold <- foldid.ext
Armin Rauschenberger's avatar
Armin Rauschenberger committed
605
  family <- "gaussian"; nfolds.ext <- 5; nfolds.int <- 10; foldid.ext <- foldid.int <- NULL; type.measure <- "deviance"; alpha.base <- alpha.meta <- 1; mice <- cvpred <- times <- FALSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
606
607
  foldid.ext <- fold; nfolds.ext <- 1
  #nfolds.ext <- 1; nfolds.int <- 10; foldid.int <- NULL; compare <- TRUE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
608
609
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
610
  if(length(compare)==1 && compare==TRUE){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
611
612
613
614
615
616
617
    if(all(family=="gaussian")){
      compare <- c("mnorm","mars","spls","mrce","map","mrf","sier","mcen","gpm","rmtl","mtps")
    } else if(all(family=="binomial")){
      compare <- c("mars","mcen","rmtl","mtps")
    } else {
      stop("Comparison not implemented for mixed families.",call.=FALSE)
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
618
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
619
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
620
621
622
  n <- nrow(Y)
  q <- ncol(Y)
  p <- ncol(X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
623
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
624
625
626
627
  #--- fold identifiers ---
  if(is.null(foldid.ext)){
    foldid.ext <- palasso:::.folds(y=Y[,1],nfolds=nfolds.ext) # temporary Y[,1]
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
628
629
    #nfolds.ext <- length(unique(foldid.ext))
    nfolds.ext <- max(foldid.ext)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
630
631
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
632
633
634
635
636
  #--- 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
637
638
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
639
640
  # check packages
  pkgs <- .packages(all.available=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
641
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
642
  if(is.character(compare)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
643
644
    for(i in seq_along(compare)){
      pkg <- switch(compare[i],mnorm="glmnet",mars="earth",spls="spls",
Armin Rauschenberger's avatar
Armin Rauschenberger committed
645
646
647
                  mrce="MRCE",map="remMap",mrf="MultivariateRandomForest",
                  sier="SiER",mcen="mcen",gpm="GPM",rmtl="RMTL",mtps="MTPS",
                  stop("Invalid method.",call.=FALSE))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
648
649
650
      if(!pkg %in% pkgs){
        stop("Method \"",compare[i],"\" requires package \"",pkg,"\".",call.=FALSE)
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
651
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
652
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
653
654
655
656
657
658
659
660

  #--- checks ---
  #if(any( & any(family!="gaussian")){
  #  stop("\"mnorm\", \"spls\", \"mrce\" and \"sier\" require family \"gaussian\".",call.=FALSE)
  #}
  #if(any(mtps,rmtl) & any(!family %in% c("gaussian","binomial"))){
  #  stop("\"mtps\" and \"rmtl\" require family \"gaussian\" or \"binomial\".",call.=FALSE)
  #}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
661
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
662
  #--- cross-validated predictions ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
663
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
664
  models <- c("base","meta",compare,"none")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
665
  pred <- lapply(X=models,function(x) matrix(data=NA,nrow=n,ncol=q))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
666
667
  time <- lapply(X=models,function(x) NA)
  names(pred) <- names(time) <- models
Armin Rauschenberger's avatar
Armin Rauschenberger committed
668
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
669
670
671
  for(i in seq_len(nfolds.ext)){
    
    Y0 <- Y[foldid.ext!=i,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
672
    #Y1 <- Y[foldid.ext==i,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
673
674
    X0 <- X[foldid.ext!=i,,drop=FALSE]
    X1 <- X[foldid.ext==i,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
675
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
676
677
678
679
680
681
682
    # standardise features (trial)
    #mu <- apply(X=X0,MARGIN=2,FUN=function(x) mean(x))
    #sd <- apply(X=X0,MARGIN=2,FUN=function(x) stats::sd(x))
    #X0 <- t((t(X0)-mu)/sd)[,sd!=0]
    #X1 <- t((t(X1)-mu)/sd)[,sd!=0]
    # or standardise once before cv?
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
683
684
685
686
    # remove constant features
    cond <- apply(X=X0,MARGIN=2,FUN=function(x) stats::sd(x)!=0)
    X0 <- X0[,cond]; X1 <- X1[,cond]
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
687
688
689
690
691
692
693
    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
694
    start <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
695
    fit <- joinet(Y=Y0,X=X0,family=family,type.measure=type.measure,alpha.base=alpha.base,alpha.meta=alpha.meta,foldid=foldid) # add ellipsis (...)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
696
    # also do not standardise!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
697
    temp <- predict.joinet(fit,newx=X1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
698
699
    pred$base[foldid.ext==i,] <- temp$base
    pred$meta[foldid.ext==i,] <- temp$meta
Armin Rauschenberger's avatar
Armin Rauschenberger committed
700
701
    end <- Sys.time()
    time$meta <- as.numeric(difftime(end,start,units="secs"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
702
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
703
    # missing values
Armin Rauschenberger's avatar
Armin Rauschenberger committed
704
705
706
707
    if(mice & any(is.na(Y0))){
      if(requireNamespace("mice",quietly=TRUE)){
        y0 <- as.matrix(mice::complete(data=mice::mice(Y0,m=1,method="pmm",seed=1,printFlag=FALSE),action="all")[[1]])
      } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
708
        stop("Imputation by PMM requires install.packages(\"mice\").",call.=FALSE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
709
710
711
712
      }
    } else {
      y0 <- apply(X=Y0,MARGIN=2,FUN=function(x) ifelse(is.na(x),stats::median(x[!is.na(x)]),x))
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
713
    all(Y0==y0,na.rm=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
714
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
715
    # other learners
Armin Rauschenberger's avatar
Armin Rauschenberger committed
716
717
    
    if("mnorm" %in% compare){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
718
719
720
721
722
      cat("mnorm"," ")
      start <- Sys.time()
      if(any(family!="gaussian")){
        stop("mnorm requires \"gaussian\" family.",call.=FALSE)
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
723
      net <- glmnet::cv.glmnet(x=X0,y=y0,family="mgaussian",foldid=foldid,alpha=alpha.base) # add ellipsis (...)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
724
      pred$mnorm[foldid.ext==i,] <- stats::predict(object=net,newx=X1,s="lambda.min",type="response")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
725
726
727
728
      end <- Sys.time()
      time$mnorm <- as.numeric(difftime(end,start,units="secs"))
    } else {
      net <- glmnet::glmnet(x=X0,y=y0,family="mgaussian")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
729
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
730
731
    
    if("mars" %in% compare){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
732
733
      cat("mars"," ")
      start <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
734
      if(all(family=="gaussian")){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
735
        object <- earth::earth(x=X0,y=y0)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
736
        # equivalent: object <- mda::mars(x=X0,y=y0)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
737
      } else if(all(family=="binomial")){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
738
        object <- earth::earth(x=X0,y=y0,glm=list(family=stats::binomial))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
739
      } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
740
        stop("MARS requires either \"gaussian\" or \"binomial\" family.",call.=FALSE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
741
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
742
      # pmethod="cv" not available for multivariate outputs
Armin Rauschenberger's avatar
Armin Rauschenberger committed
743
      
Armin Rauschenberger's avatar
Armin Rauschenberger committed
744
745
746
747
748
749
750
751
752
      ### start trial ###
      if(FALSE){
      #nk <- min(200, max(20, 2 * ncol(X0))) + 1
      #nprune <- round(seq(from=2,to=nk,length.out=10))
      #object <- list()
      #for(j in seq_along(nprune)){
      #  object[[j]] <- earth::earth(x=X0,y=y0,nprune=nprune[j],pmethod="cv",nfold=nfolds.int)
      #}
      #sapply(object,function(x) x$gcv)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
753
754
755
      ## i.e. run earth/mars with tryCatch for each nprune
      ## and select run with best cvm (here gcv)
      # tune nprune (use default nk)!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
756
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
757
      
Armin Rauschenberger's avatar
Armin Rauschenberger committed
758
759
      #pred$mars[foldid.ext==i,] <- earth:::predict.earth(object=object,newdata=X1,type="response") # original
      pred$mars[foldid.ext==i,] <- stats::predict(object=object,newdata=X1,type="response") # trial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
760
761
      end <- Sys.time()
      time$mars <- as.numeric(difftime(end,start,units="secs"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
762
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
763
764
    
    if("spls" %in% compare){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
765
766
767
768
769
770
      cat("spls"," ")
      start <- Sys.time()
      if(any(family!="gaussian")){
        stop("spls requires \"gaussian\" family.")
      }
      invisible(utils::capture.output(cv.spls <- spls::cv.spls(x=X0,y=y0,fold=nfolds.int,K=seq_len(min(ncol(X0),10)),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
771
                               eta=seq(from=0.0,to=0.9,by=0.1),plot.it=FALSE))) #scale.x=FALSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
772
      object <- spls::spls(x=X0,y=y0,K=cv.spls$K.opt,eta=cv.spls$eta.opt) #scale.x=FALSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
773
      pred$spls[foldid.ext==i,] <- spls::predict.spls(object=object,newx=X1,type="fit")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
774
775
      end <- Sys.time()
      time$spls <- as.numeric(difftime(end,start,units="secs"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
776
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
777
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
778
779
    if("mrce" %in% compare){
      cat("mrce"," ")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
780
      start <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
781
782
      if(any(family!="gaussian")){
        stop("MRCE requires \"gaussian\" family.",call.=FALSE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
783
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
784
      lam1 <- lam2 <- 10^seq(from=1,to=-4,length.out=11)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
785
      invisible(utils::capture.output(trials <- lapply(lam2,function(x) tryCatch(expr=MRCE::mrce(X=X0,Y=y0,lam1.vec=lam1,lam2.vec=x,method="cv",kfold=nfolds.int),error=function(x) NULL))))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
786
787
788
      cv.err <- sapply(trials,function(x) ifelse(is.null(x),Inf,min(x$cv.err)))
      object <- trials[[which.min(cv.err)]]
      pred$mrce[foldid.ext==i,] <- matrix(object$muhat,nrow=nrow(X1),ncol=q,byrow=TRUE) + X1 %*% object$Bhat
Armin Rauschenberger's avatar
Armin Rauschenberger committed
789
      end <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
790
      time$mrce <- as.numeric(difftime(end,start,units="secs"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
791
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
792
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
793
794
    if("map" %in% compare){
      cat("map"," ")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
795
      start <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
796
797
      if(any(family!="gaussian")){
        stop("map requires \"gaussian\" family.")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
798
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
799
800
      mean <- colMeans(y0)
      y0s <- y0-matrix(data=mean,nrow=nrow(X0),ncol=ncol(y0),byrow=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
801
      lamL1.v <- lamL2.v <- exp(seq(from=0,to=5,length.out=11))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
802
      cv <- remMap::remMap.CV(X=X0,Y=y0s,lamL1.v=lamL1.v,lamL2.v=lamL2.v,fold=nfolds.int)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
803
      #graphics::plot(x=lamL1.v,y=log(as.numeric(cv$ols.cv[,3])))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
804
805
806
807
808
809
810
811
      pick <- which.min(as.vector(cv$ols.cv))
      lamL1 <- cv$l.index[1,pick]
      lamL2 <- cv$l.index[2,pick]
      # index <- which(cv$ols.cv==min(cv$ols.cv),arr.ind=TRUE)[1,]
      # rev(lamL1.v)[index[1]]
      # rev(lamL2.v)[index[2]]
      ##cat("lam1:",lamL1,", lam2:",lamL2)
      object <- remMap::remMap(X.m=X0,Y.m=y0s,lamL1=lamL1,lamL2=lamL2)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
812
      pred$map[foldid.ext==i,] <- matrix(data=mean,nrow=nrow(X1),ncol=ncol(y0),byrow=TRUE) + X1 %*% object$phi
Armin Rauschenberger's avatar
Armin Rauschenberger committed
813
      end <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
814
      time$map <- as.numeric(difftime(end,start,units="secs"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
815
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
816
817
    
    if("mrf" %in% compare){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
818
819
820
821
822
823
      cat("mrf"," ")
      start <- Sys.time()
      if(any(family!="gaussian")){
        stop("mrf requires \"gaussian\" family.")
      }
      pred$mrf[foldid.ext==i,] <- MultivariateRandomForest::build_forest_predict(trainX=X0,trainY=y0,
Armin Rauschenberger's avatar
Armin Rauschenberger committed
824
                                   n_tree=100,m_feature=min(ncol(X0)-1,5),min_leaf=min(nrow(X0),5),testX=X1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
825
      # use n_tree=500, m_feature=floor(ncol(X0)/3)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
826
      # alternative: IntegratedMRF
Armin Rauschenberger's avatar
Armin Rauschenberger committed
827
828
829
      # Check why this does not work well!
      end <- Sys.time()
      time$mrf <- as.numeric(difftime(end,start,units="secs"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
830
831
    }
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
832
833
834
835
836
837
    if("sier" %in% compare){
      cat("sier"," ")
      start <- Sys.time()
      if(any(family!="gaussian")){
        stop("SiER requires \"gaussian\" family.",call.=FALSE)
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
838
      invisible(utils::capture.output(object <- SiER::cv.SiER(X=X0,Y=y0,K.cv=3)))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
839
840
841
842
843
844
845
      # trial with K.cv=3 (for spped-up)
      # use upper.comp=10 and thres=0.01  (changed for speed-up)
      pred$sier[foldid.ext==i,] <- SiER::pred.SiER(cv.fit=object,X.new=X1)
      end <- Sys.time()
      time$sier <- as.numeric(difftime(end,start,units="secs"))
    }
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
846
    if("mcen" %in% compare){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
847
848
849
850
851
852
853
854
855
856
      cat("mcen"," ")
      start <- Sys.time()
      if(all(family=="gaussian")){
        type <- "mgaussian"
      } else if(all(family=="binomial")){
        type <- "mbinomial"
      } else {
        stop("mcen requires either \"gaussian\" or \"binomial\".",call.=FALSE)
      }
      object <- mcen::cv.mcen(x=X0,y=y0,family=type,folds=foldid,ky=1,
Armin Rauschenberger's avatar
Armin Rauschenberger committed
857
858
                              gamma_y=seq(from=0.1,to=5.1,by=1),ndelta=5)
      # TEMPORARY gamma_y=seq(from=0.1,to=5.1,length.out=3) and ndelta=3 (for speed-up)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
859
860
      #temp <- mcen:::predict.cv.mcen(object=object,newx=X1) # original
      temp <- stats::predict(object=object,newx=X1) # trial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
861
      pred$mcen[foldid.ext==i,] <- as.matrix(temp)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
862
863
864
      # single cluster (ky=1) due to setting and error
      end <- Sys.time()
      time$mcen <- as.numeric(difftime(end,start,units="secs"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
865
866
    }
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
867
868
869
    if("gpm" %in% compare){
      cat("gpm"," ")
      start <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
870
871
872
      if(any(family!="gaussian")){
        stop("GPM requires \"gaussian\" family.",call.=FALSE)
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
873
      object <- GPM::Fit(X=X0,Y=y0)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
874
      pred$gpm[foldid.ext==i,] <- GPM::Predict(XF=X1,Model=object)$YF
Armin Rauschenberger's avatar
Armin Rauschenberger committed
875
876
877
878
      end <- Sys.time()
      time$gpm <- as.numeric(difftime(end,start,units="secs"))
    }
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
879
880
    if("rmtl" %in% compare){
      cat("rmtl"," ")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
881
      start <- Sys.time()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
882
883
884
885
886
887
888
889
      if(all(family=="gaussian")){
        type <- "Regression"
        y0l <- lapply(seq_len(ncol(y0)),function(i) y0[,i,drop=FALSE])
      } else if(all(family=="binomial")){
        type <- "Classification"
        y0l <- lapply(seq_len(ncol(y0)),function(i) 2*y0[,i,drop=FALSE]-1)
      } else {
        stop("RMTL requires either \"gaussian\" or \"binomial\".",call.=FALSE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
890
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
891
892
893
894
895
896
897
898
899
900
      X0l <- lapply(seq_len(ncol(y0)),function(i) X0)
      X1l <- lapply(seq_len(ncol(y0)),function(i) X1)
      #---------------------------
      #--- manual tuning start ---
      Lam1_seq <- c(10^seq(from=1,to=-4,length.out=11),0)
      Lam2_seq <- c(10^seq(from=1,to=-4,length.out=11),0)
      cvMTL <- list()
      seed <- .Random.seed
      for(j in seq_along(Lam2_seq)){
        .Random.seed <- seed
Armin Rauschenberger's avatar
Armin Rauschenberger committed
901
        cvMTL[[j]] <- RMTL::cvMTL(X=X0l,Y=y0l,type=type,Lam1_seq=Lam1_seq,Lam2=Lam2_seq[j],nfolds=nfolds.int)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
902
903
904
905
906
907
908
909
910
      }
      cvm <- vapply(X=cvMTL,FUN=function(x) min(x$cvm),FUN.VALUE=numeric(1))
      Lam1 <- cvMTL[[which.min(cvm)]]$Lam1.min
      Lam2 <- Lam2_seq[which.min(cvm)]
      #graphics::plot(x=Lam2_seq,y=cvm)
      #cat(Lam1," ",Lam2,"\n")
      #--- manual tuning end ----
      #--------------------------
      MTL <- RMTL::MTL(X=X0l,Y=y0l,type=type,Lam1=Lam1,Lam2=Lam2)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
911
912
      #temp <- RMTL:::predict.MTL(object=MTL,newX=X1l) # original
      temp <- stats::predict(object=MTL,newX=X1l)
Armin Rauschenberger's avatar
Armin Rauschenberger committed