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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
214
#--- Methods for class "joinet" -----------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
215

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

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

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

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

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

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