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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
2
#--- Workhorse function --------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
3

Armin Rauschenberger's avatar
Armin Rauschenberger committed
4
5
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
6
#' Logistic regression with a continuous response
Armin Rauschenberger's avatar
Armin Rauschenberger committed
7
8
#' 
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
9
#' Implements logistic regression with a continuous response.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
10
11
12
13
14
15
#'  
#' @param y
#' continuous response\strong{:}
#' vector of length \eqn{n}
#' 
#' @param cutoff
Armin Rauschenberger's avatar
Armin Rauschenberger committed
16
#' cutoff point for dichotomising response into classes\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
17
18
19
20
#' value between \code{min(y)} and \code{max(y)}
#' 
#' @param X
#' covariates\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
21
#' numeric matrix with \eqn{n} rows (samples)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
22
#' and \eqn{p} columns (variables)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
23
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
24
25
26
27
#' @param alpha
#' elastic net mixing parameter\strong{:}
#' numeric between \eqn{0} (ridge) and \eqn{1} (lasso)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
28
29
30
31
#' @param foldid
#' fold identifiers\strong{:}
#' vector with entries between \eqn{1} and \code{nfolds};
#' or \code{NULL} (balance)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
32
33
34
35
#' 
#' @param nfolds
#' number of folds
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
36
#' @param type.measure
Armin Rauschenberger's avatar
Armin Rauschenberger committed
37
38
39
#' loss function for binary classification
#' (linear regression uses the deviance)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
40
41
42
43
44
45
46
47
#' @param pi
#' pi sequence\strong{:}
#' vector of values in the unit interval;
#' or \code{NULL} (default sequence)
#' 
#' @param npi
#' number of \code{pi} values
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
48
49
50
51
52
53
54
55
#' @param sigma
#' sigma sequence\strong{:}
#' vector of increasing positive values;
#' or \code{NULL} (default sequence)
#' 
#' @param nsigma
#' number of \code{sigma} values
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
56
57
58
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}}
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
59
#' @details
Armin Rauschenberger's avatar
Armin Rauschenberger committed
60
61
62
63
#' - INCLUDE note on deviance (not comparable between lin and log models)
#' - alpha: elastic net parameter\strong{:}
#' numeric between \eqn{0} (ridge) and \eqn{1} (lasso)
#' - do not use "family"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
64
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
65
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
66
67
68
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
69
#' net <- cornet(y=y,cutoff=0,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
70
#' ### Add ... to all glmnet::glmnet calls !!! ###
Armin Rauschenberger's avatar
Armin Rauschenberger committed
71
cornet <- function(y,cutoff,X,alpha=1,npi=101,pi=NULL,nsigma=99,sigma=NULL,nfolds=10,foldid=NULL,type.measure="deviance",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
72
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
73
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
74
  # cutoff <- 0; npi <- 101; pi <- NULL; nsigma <- 99; sigma <- NULL; nfolds <- 10;  foldid <- NULL; type.measure <- "deviance"; logistic <- TRUE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
75
  test <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
76
  test$sigma <- test$pi <- FALSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
77
  test$grid <- TRUE
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
  cornet:::.check(x=y,type="vector")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
81
  if(all(y %in% c(0,1))){warning("Binary response.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
82
83
  cornet:::.check(x=cutoff,type="scalar",min=min(y),max=max(y))
  cornet:::.check(x=X,type="matrix")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
84
  cornet:::.check(x=alpha,type="scalar",min=0,max=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
85
  if(length(y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
86
87
88
89
90
91
92
  cornet:::.check(x=npi,type="scalar",min=1)
  cornet:::.check(x=pi,type="vector",min=0,max=1,null=TRUE)
  cornet:::.check(x=nsigma,type="scalar",min=1)
  cornet:::.check(x=sigma,type="vector",min=.Machine$double.eps,null=TRUE)
  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)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
93
94
  if(!is.null(list(...)$family)){stop("Reserved argument \"family\".",call.=FALSE)}
  n <- length(y)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
95
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
96
  # binarisation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
97
98
  z <- 1*(y > cutoff)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
99
  # fold identifiers
Armin Rauschenberger's avatar
Armin Rauschenberger committed
100
101
102
103
  if(is.null(foldid)){
    foldid <- palasso:::.folds(y=z,nfolds=nfolds)
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
104
  #--- model fitting ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
105
  fit <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
106
107
  fit$gaussian <- glmnet::glmnet(y=y,x=X,family="gaussian",alpha=alpha)
  fit$binomial <- glmnet::glmnet(y=z,x=X,family="binomial",alpha=alpha)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
108
109
   
  #--- sigma sequence ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
110
111
  if(is.null(sigma)){
    fit$sigma <- exp(seq(from=log(0.05*stats::sd(y)),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
112
                  to=log(10*stats::sd(y)),length.out=nsigma))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
113
114
  } else {
    fit$sigma <- sigma
Armin Rauschenberger's avatar
Armin Rauschenberger committed
115
    nsigma <- length(sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
116
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
117
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
118
119
  #--- pi sequence ---
  if(is.null(pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
120
    fit$pi <- seq(from=0,to=1,length.out=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
121
122
123
  } else {
    fit$pi <- pi
    npi <- length(pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
124
125
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
126
127
128
129
130
  #--- tuning parameters ---
  lab.pi <- paste0("pi",seq_len(npi))
  lab.sigma <- paste0("si",seq_len(nsigma))
  names(fit$sigma) <- lab.sigma
  names(fit$pi) <- lab.pi
Armin Rauschenberger's avatar
Armin Rauschenberger committed
131
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
132
  #--- cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
133
  pred <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
134
135
136
  pred$y <- matrix(data=NA,nrow=n,ncol=length(fit$gaussian$lambda))
  pred$z <- matrix(data=NA,nrow=n,ncol=length(fit$binomial$lambda))

Armin Rauschenberger's avatar
Armin Rauschenberger committed
137
138
139
140
  if(test$sigma){
    pred$sigma <- matrix(data=NA,nrow=n,ncol=nsigma)
  }
  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
141
    pred$pi <- matrix(data=NA,nrow=n,ncol=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
142
143
  }
  if(test$grid){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
144
    dimnames <- list(NULL,lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
145
    pred$grid <- array(data=NA,dim=c(n,nsigma,npi),dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
146
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
147

Armin Rauschenberger's avatar
Armin Rauschenberger committed
148
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
149

Armin Rauschenberger's avatar
Armin Rauschenberger committed
150
151
152
153
154
155
156
    y0 <- y[foldid!=k]
    y1 <- y[foldid==k]
    z0 <- z[foldid!=k]
    z1 <- z[foldid==k]
    X0 <- X[foldid!=k,,drop=FALSE]
    X1 <- X[foldid==k,,drop=FALSE]
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
157
    # linear regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
158
    net <- glmnet::glmnet(y=y0,x=X0,family="gaussian",alpha=alpha)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
159
160
    temp_y <- stats::predict(object=net,newx=X1,type="response",s=fit$gaussian$lambda)
    pred$y[foldid==k,seq_len(ncol(temp_y))] <- temp_y
Armin Rauschenberger's avatar
Armin Rauschenberger committed
161
    cvm <- cornet:::.loss(y=y1,fit=temp_y,family="gaussian",type.measure="deviance")[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
162
    y_hat <- temp_y[,which.min(cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
163
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
164
    # logistic regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
165
    net <- glmnet::glmnet(y=z0,x=X0,family="binomial",alpha=alpha)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
166
167
    temp_z <- stats::predict(object=net,newx=X1,type="response",s=fit$binomial$lambda)
    pred$z[foldid==k,seq_len(ncol(temp_z))] <- temp_z
Armin Rauschenberger's avatar
Armin Rauschenberger committed
168
    cvm <- cornet:::.loss(y=z1,fit=temp_z,family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
169
170
    z_hat <- temp_z[,which.min(cvm)]
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
171
    # fusion (sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
172
173
    if(test$sigma){
      for(i in seq_along(fit$sigma)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
174
        #pred$sigma[foldid==k,i] <- stats::pnorm(q=y_hat,mean=cutoff,sd=fit$sigma[i])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
175
      }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
176
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
177
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
178
    # fusion (pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
179
180
    if(test$pi){
      for(i in seq_along(fit$pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
181
182
        #cont <- stats::pnorm(q=y_hat,mean=cutoff,sd=stats::sd(y))
        #pred$pi[foldid==k,i] <- fit$pi[i]*cont + (1-fit$pi[i])*z_hat
Armin Rauschenberger's avatar
Armin Rauschenberger committed
183
184
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
185
186
187

    # fusion (grid)
    if(test$grid){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
188
189
190
      for(i in seq_along(fit$sigma)){
        for(j in seq_along(fit$pi)){
          cont <- stats::pnorm(q=y_hat,mean=cutoff,sd=fit$sigma[i])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
191
          pred$grid[foldid==k,i,j] <- fit$pi[j]*cont + (1-fit$pi[j])*z_hat
Armin Rauschenberger's avatar
Armin Rauschenberger committed
192
193
194
        }
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
195
196
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
197
198
  #--- evaluation ---
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
199
  # deviance (not comparable between Gaussian and binomial families)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
200
  fit$gaussian$cvm <- cornet:::.loss(y=y,fit=pred$y,family="gaussian",type.measure="deviance")[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
201
  fit$gaussian$lambda.min <- fit$gaussian$lambda[which.min(fit$gaussian$cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
202
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
203
  fit$binomial$cvm <- cornet:::.loss(y=z,fit=pred$z,family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
204
  fit$binomial$lambda.min <- fit$binomial$lambda[which.min(fit$binomial$cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
205

Armin Rauschenberger's avatar
Armin Rauschenberger committed
206
  if(test$sigma){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
207
    #fit$sigma.cvm <- cornet:::.loss(y=z,fit=pred$sigma,family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
208
    #fit$sigma.min1 <- fit$sigma[which.min(fit$sigma.cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
209
210
211
  }

  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
212
    #fit$pi.cvm <- cornet:::.loss(y=z,fit=pred$pi,family="binomial",type.measure=type.measure)[[1]] # trial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
213
    #fit$pi.min1 <- fit$pi[which.min(fit$pi.cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
214
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
215

Armin Rauschenberger's avatar
Armin Rauschenberger committed
216
  if(test$grid){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
217
    dimnames <- list(lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
218
    fit$cvm <- matrix(data=NA,nrow=nsigma,ncol=npi,dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
219
    for(i in seq_len(nsigma)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
220
      for(j in seq_len(npi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
221
        fit$cvm[i,j] <- cornet:::.loss(y=z,fit=pred$grid[,i,j],family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
222
223
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
224
    temp <- which(fit$cvm==min(fit$cvm),arr.ind=TRUE,useNames=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
225
    if(nrow(temp)>1){warning("MULTIPLE!",call.=FALSE);temp <- temp[1,,drop=FALSE]}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
226
227
228
    fit$sigma.min <- fit$sigma[temp[1]]
    fit$pi.min <- fit$pi[temp[2]]
    if(fit$cvm[names(fit$sigma.min),names(fit$pi.min)]!=min(fit$cvm)){stop("Internal error.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
229
230
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
231
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
232
  fit$cutoff <- cutoff
Armin Rauschenberger's avatar
Armin Rauschenberger committed
233
234
  fit$info <- list(type.measure=type.measure,
                   sd.y=stats::sd(y),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
235
                   table=table(z),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
236
                   test=as.data.frame(test))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
237

Armin Rauschenberger's avatar
Armin Rauschenberger committed
238
  class(fit) <- "cornet"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
239
240
241
  return(fit)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
242
#--- Methods -------------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
243

Armin Rauschenberger's avatar
Armin Rauschenberger committed
244
245
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
246
#' Extract estimated coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
247
248
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
249
250
#' Extracts estimated coefficients from linear and logistic regression,
#' under the penalty parameter that minimises the cross-validated loss.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
251
252
#'
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
253
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
254
255
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
256
257
258
259
260
261
262
#' further arguments (not applicable)
#' 
#' @return
#' This function returns a matrix with \eqn{n} rows and two columns,
#' where \eqn{n} is the sample size. It includes the estimated coefficients
#' from linear (first column, \code{"beta"})
#' and logistic (second column, \code{"gamma"}) regression.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
263
264
265
266
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
267
coef.cornet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
268
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
269
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
270
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
271
272
273
274
275
  s <- object$gaussian$lambda.min
  beta <- glmnet::coef.glmnet(object=object$gaussian,s=s)
  
  s <- object$binomial$lambda.min
  gamma <- glmnet::coef.glmnet(object=object$binomial,s=s)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
276
277
278
279
280
281
  
  coef <- cbind(beta,gamma)
  colnames(coef) <- c("beta","gamma")
  return(coef)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
282
283
284

#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
285
#' Plot loss matrix
Armin Rauschenberger's avatar
Armin Rauschenberger committed
286
287
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
288
289
#' Plots the loss for difference combinations of
#' the weight (pi) and scale (sigma) paramters.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
290
291
#'
#' @param x
Armin Rauschenberger's avatar
Armin Rauschenberger committed
292
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
293
294
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
295
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
296
297
298
299
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
300
plot.cornet <- function(x,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
301
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
302
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
303

Armin Rauschenberger's avatar
Armin Rauschenberger committed
304
  k <- 100
Armin Rauschenberger's avatar
Armin Rauschenberger committed
305
  levels <- stats::quantile(x$cvm,probs=seq(from=0,to=1,length.out=k+1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
306
307
308
309
310
311
312
313
314
  
  ## RColorBrewer
  if("RColorBrewer" %in% .packages(all.available=TRUE)){
    pal <- rev(c("white",RColorBrewer::brewer.pal(n=9,name="Blues")))
    col <- grDevices::colorRampPalette(colors=pal)(k)
  } else {
    col <- grDevices::heat.colors(n=k)
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
315
  nsigma <- length(x$sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
316
  npi <- length(x$pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
317
318
319
  
  graphics::plot.new()
  graphics::par(xaxs="i",yaxs="i")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
320
  graphics::plot.window(xlim=c(1-0.5,nsigma+0.5),ylim=c(1-0.5,npi+0.5))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
321
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
322
  graphics::title(xlab=expression(sigma),ylab=expression(pi),cex.lab=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
323
  #graphics::.filled.contour(x=seq_along(x$sigma),y=seq_along(x$pi),z=x$cvm,levels=levels,col=col)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
324
  graphics::image(x=seq_along(x$sigma),y=seq_along(x$pi),z=x$cvm,breaks=levels,col=col,add=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
325
  graphics::box()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
326
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
327
328
329
330
  ssigma <- which(x$sigma %in% x$sigma.min)
  spi <- which(x$pi %in% x$pi.min)
  
  if(length(ssigma)==1 & length(spi)==1){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
331
    # axes with labels for tuned parameters
Armin Rauschenberger's avatar
Armin Rauschenberger committed
332
333
    graphics::axis(side=1,at=c(1,ssigma,nsigma),labels=signif(x$sigma[c(1,ssigma,nsigma)],digits=2))
    graphics::axis(side=2,at=c(1,spi,npi),labels=signif(x$pi[c(1,spi,npi)],digits=2))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
334
    # point for tuned parameters
Armin Rauschenberger's avatar
Armin Rauschenberger committed
335
    graphics::points(x=ssigma,y=spi,pch=4,col="white",cex=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
336
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
337
    # axes with standard labels
Armin Rauschenberger's avatar
Armin Rauschenberger committed
338
339
340
341
    at <- seq(from=1,to=nsigma,length.out=5)
    graphics::axis(side=1,at=at,labels=signif(x$sigma,digits=2)[at])
    at <- seq(from=1,to=nsigma,length.out=5)
    graphics::axis(side=2,at=at,labels=signif(x$pi,digits=2)[at])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
342
343
344
    # points for selected parameters
    isigma <- sapply(x$sigma.min,function(y) which(x$sigma==y))
    ipi <- sapply(x$pi.min,function(y) which(x$pi==y))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
345
    graphics::points(x=isigma,y=ipi,pch=4,col="white",cex=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
346
347
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
348
349
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
350
351
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
352
#' Predict binary outcome
Armin Rauschenberger's avatar
Armin Rauschenberger committed
353
354
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
355
356
357
358
359
360
361
362
#' Predicts the binary outcome with linear, logistic, and combined regression.
#' 
#' @details
#' For linear regression, this function tentatively transforms
#' the predicted values to predicted probabilities,
#' using a Gaussian distribution with a fixed mean (threshold)
#' and a fixed variance (estimated variance of the numeric outcome).
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
363
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
364
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
365
366
367
368
369
370
371
372
373
374
#' 
#' @param newx
#' covariates\strong{:}
#' numeric matrix with \eqn{n} rows (samples)
#' and \eqn{p} columns (variables)
#' 
#' @param type
#' \code{"probability"}, \code{"odds"}, \code{"log-odds"}
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
375
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
376
377
378
379
#' 
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
380
predict.cornet <- function(object,newx,type="probability",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
381
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
382
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
383
384
  
  x <- object; rm(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
385
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
386
387
  test <- x$info$test
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
388
389
  .check(x=newx,type="matrix")
  .check(x=type,type="string",values=c("probability","odds","log-odds"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
390
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391
392
393
394
395
396
397
  # linear, logistic and mixed
  prob <- list()
  link <- as.numeric(stats::predict(object=x$gaussian,
                  newx=newx,s=x$gaussian$lambda.min,type="response"))
  prob$gaussian <- stats::pnorm(q=link,mean=x$cutoff,sd=x$info$sd.y)
  prob$binomial <- as.numeric(stats::predict(object=x$binomial,
                  newx=newx,s=x$binomial$lambda.min,type="response"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
398
399
  
  if(test$sigma){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
400
    #prob$sigma <- stats::pnorm(q=link,mean=x$cutoff,sd=x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
401
402
403
  }
  
  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
404
     #prob$pi <- x$pi.min*prob$gaussian + (1-x$pi.min)*prob$binomial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
405
406
407
  }
 
  if(test$grid){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
408
    cont <- stats::pnorm(q=link,mean=x$cutoff,sd=x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
409
    prob$grid <- x$pi.min*cont + (1-x$pi.min)*prob$binomial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
410
411
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
412
413
  # consistency tests
  lapply(X=prob,FUN=function(p) .check(x=p,type="vector",min=0,max=1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
414
  .equal(link>x$cutoff,prob$gaussian>0.5)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
415
416
417
418
419
420
421
422
423
424

  # transformation
  if(type=="probability"){
    frame <- prob
  } else if(type=="odds"){
    frame <- lapply(X=prob,FUN=function(x) x/(1-x))
  } else if(type=="log-odds"){
    frame <- lapply(X=prob,FUN=function(x) log(x/(1-x)))
  } else {
    stop("Invalid type.",call.=FALSE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
425
426
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
427
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
428
429
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
430
#--- Application ---------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
431

Armin Rauschenberger's avatar
Armin Rauschenberger committed
432
433
434
435
436
437
438
#' @export
#' @title
#' Comparison
#'
#' @description
#' Compares models for a continuous response with a cutoff value
#'
Armin Rauschenberger's avatar
Armin Rauschenberger committed
439
#' @inheritParams  cornet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
440
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
441
442
443
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
444
.compare <- function(y,cutoff,X,alpha=1,nfolds=5,foldid=NULL,type.measure="deviance"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
445
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
446
  z <- 1*(y > cutoff)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
447
448
449
450
451
  if(is.null(foldid)){
    fold <- palasso:::.folds(y=z,nfolds=nfolds)
  } else {
    fold <- foldid
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
452
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
453
  cols <- c("gaussian","binomial","grid")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
454
455
  pred <- matrix(data=NA,nrow=length(y),ncol=length(cols),
                 dimnames=list(NULL,cols))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
456
457
  
  select <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
458
  for(i in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
459
    fit <- cornet::cornet(y=y[fold!=i],cutoff=cutoff,X=X[fold!=i,],alpha=alpha,type.measure=type.measure)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
460
461
    tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL)
    temp <- cornet:::predict.cornet(fit,newx=X[fold==i,])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
462
    if(any(temp<0|temp>1)){stop("Outside unit interval.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
463
464
465
466
    model <- colnames(pred)
    for(j in seq_along(model)){
      pred[fold==i,model[j]] <- temp[[model[j]]]
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
467
468
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
469
  type <- c("deviance","class","mse","mae","auc")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
470
  loss <- lapply(X=type,FUN=function(x) cornet:::.loss(y=z,fit=pred,family="binomial",type.measure=x,foldid=fold)[[1]])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
471
  names(loss) <- type
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
497
498
499
500
501
502
503
  ###################
  ### start trial ###
  
  # squared deviance residuals
  limit <- 1e-05
  pred[pred < limit] <- limit
  pred[pred > 1 - limit] <- 1 - limit
  res <- -2 * (z * log(pred) + (1 - z) * log(1 - pred))
  rxs <- res[,"binomial"]
  rys <- res[,"grid"]
  
  ## examine differences per fold
  loss$cv.diff <- loss$cv.size <- loss$cv.pval <- numeric()
  for(i in seq_len(nfolds)){
    cond <- fold==i
    loss$cv.size[i] <- mean((rxs>rys)[cond])
    loss$cv.diff[i] <- stats::median(((rys-rxs)/rxs)[cond])
    loss$cv.pval[i] <- stats::wilcox.test(x=rxs[cond],y=rys[cond],paired=TRUE,
                           alternative="greater")$p.value
  }
  
  # examine all differences
  loss$all.size <- mean(rxs>rys)
  loss$all.diff <- stats::median((rys-rxs)/rxs)
  loss$all.pval <- stats::wilcox.test(x=rxs,y=rys,paired=TRUE,
                                 alternative="greater")$p.value
  # The overall p-value is anti-conservative!
  
  ### end trial ###
  #################
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
504
505
506
507
508
509
510
511
512
513
514
  #if(trial){
  #  list <- list(diff=(pred-z)^2,fold=fold,loss=loss)
  #  return(list)
  #} else {
  #  return(loss)
  #}
  
  return(loss)

}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
515

Armin Rauschenberger's avatar
Armin Rauschenberger committed
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
#' @export
#' @title
#' Testing
#'
#' @description
#' Compares models for a continuous response with a cutoff value (testing)
#'
#' @inheritParams cornet
#' 
#' @examples
#' NA
#' 
.test <- function(y,cutoff,X,alpha=1,type.measure="deviance"){
  
  z <- 1*(y > cutoff)
  fold <- palasso:::.folds(y=z,nfolds=5)
  fold <- ifelse(fold==1,1,0)
  
  fit <- cornet::cornet(y=y[fold==0],cutoff=cutoff,X=X[fold==0,],alpha=alpha)
  tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL)
  pred <- cornet:::predict.cornet(fit,newx=X[fold==1,])
  if(any(pred<0|pred>1)){stop("Outside unit interval.",call.=FALSE)}
  
  #res <- (pred-z[fold==1])^2 # MSE
  #pvalue <- wilcox.test(x=res[,"binomial"],y=res[,"grid"],paired=TRUE,alternative="greater")$p.value
  #colMeans(abs(pred-0.5)) # distance from 0.5
Armin Rauschenberger's avatar
Armin Rauschenberger committed
542
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
543
544
545
  limit <- 1e-05
  pred[pred < limit] <- limit
  pred[pred > 1 - limit] <- 1 - limit
Armin Rauschenberger's avatar
Armin Rauschenberger committed
546
547
  res <- -2 * (z[fold==1] * log(pred) + (1 - z[fold==1]) * log(1 - pred))
  # Changed y to z (2019-02-08)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
548
549
550
551
552
553
554
555
556
557
558
559
560
561
  pvalue <- stats::wilcox.test(x=res[,"binomial"],y=res[,"grid"],paired=TRUE,alternative="greater")$p.value
  
  return(pvalue)
}

# Simulates y and X.
.simulate <- function(n,p,prob=0.2,fac=1){
  beta <- stats::rnorm(n=p)
  cond <- stats::rbinom(n=p,size=1,prob=prob)
  beta[cond==0] <- 0
  X <- matrix(stats::rnorm(n=n*p),nrow=n,ncol=p)
  mean <- X %*% beta
  y <- stats::rnorm(n=n,mean=mean,sd=fac*stats::sd(mean))
  return(list(y=y,X=X))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
562
563
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
564
565
566
567
568
569
570
571
572
573
574
575
576
#--- start trial ---
if(FALSE){
  n <- 1000
  y_hat <- runif(n)
  y <- y_hat > 0.9
  y <- rbinom(n=n,size=1,prob=0.5)
  foldid <- rep(1:10,length.out=n)
  .loss(y=y,fit=y_hat,family="binomial",type.measure="auc",foldid=foldid)
}
#--- end trial ---

#--- Internal functions --------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
#' @title
#' Arguments
#'
#' @description
#' Verifies whether two or more arguments are identical.
#'
#' @param ...
#' scalars, vectors, or matrices of equal dimensions
#' 
#' @param na.rm
#' remove missing values\strong{:}
#' logical
#' 
#' @examples 
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
.equal <- function(...,na.rm=FALSE){
  list <- list(...)
  cond <- vapply(X=list,
                 FUN=function(x) all(x==list[[1]],na.rm=na.rm),
                 FUN.VALUE=logical(1))
  if(any(!cond)){
    stop("Unequal elements.",call.=FALSE)
  }
  return(invisible(NULL))
}

#' @title
#' Arguments
#'
#' @description
#' Verifies whether an argument matches formal requirements.
#'
#' @param x
#' argument
#' 
#' @param type
#' character \code{"string"}, \code{"scalar"}, \code{"vector"}, \code{"matrix"}
#' 
#' @param miss
#' accept missing values\strong{:}
#' logical
#' 
#' @param min
#' lower limit\strong{:}
#' numeric
#' 
#' @param max
#' upper limit\strong{:}
#' numeric
#' 
#' @param values
#' only accept specific values\strong{:}
#' vector 
#' 
#' @param inf
#' accept infinite (\code{Inf} or \code{-Inf}) values\strong{:}
#' logical
#' 
#' @param null
#' accept \code{NULL}\strong{:}
#' logical
#'
#' @examples
#' NA
#' 
.check <- function(x,type,miss=FALSE,min=NULL,max=NULL,values=NULL,inf=FALSE,null=FALSE){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
644
  name <- deparse(substitute(x))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
645
646
647
  if(null && is.null(x)){
    return(invisible(NULL)) 
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
  if(type=="string"){
    cond <- is.vector(x) & is.character(x) & length(x)==1
  } else if(type=="scalar"){
    cond <- is.vector(x) & is.numeric(x) & length(x)==1
  } else if(type=="vector"){
    cond <- is.vector(x) & is.numeric(x)
  } else if(type=="matrix"){
    cond <- is.matrix(x) & is.numeric(x)
  } else {
    warning("Unknown type.")
  }
  if(!cond){
    stop(paste0("Argument \"",name,"\" does not match formal requirements."),call.=FALSE)
  }
  if(!miss && any(is.na(x))){
    stop(paste0("Argument \"",name,"\" contains missing values."),call.=FALSE)
  }
  if(!is.null(min) && any(x<min)){
    stop(paste0("expecting ",name," >= ",min),call.=FALSE)
  }
  if(!is.null(max) && any(x>max)){
    stop(paste0("expecting ",name," <= ",max),call.=FALSE)
  }
  if(!is.null(values) && !(x %in% values)){
    stop(paste0("Argument \"",name,"\" contains invalid values."),call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
674
675
676
  if(!inf && any(is.infinite(values))){
    stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
677
678
679
  return(invisible(NULL))
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
680
681
682



Armin Rauschenberger's avatar
Armin Rauschenberger committed
683
684
# Correct this function in the palasso package (search twice for "# typo").
.loss <- function (y,fit,family,type.measure,foldid=NULL){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
  if (!is.list(fit)) {
    fit <- list(fit)
  }
  loss <- list()
  for (i in seq_along(fit)) {
    if (is.vector(fit[[i]])) {
      fit[[i]] <- as.matrix(fit[[i]])
    }
    if (is.null(foldid) & (family == "cox" | type.measure == 
                           "auc")) {
      stop("Missing foldid.", call. = FALSE)
    }
    if (family == "gaussian") {
      if (type.measure %in% c("deviance", "mse")) {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) mean((x - y)^2))
      }
      else if (type.measure == "mae") {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) mean(abs(x - y)))
      }
      else {
        stop("Invalid type measure.", call. = FALSE)
      }
    }
    else if (family == "binomial") {
      if (type.measure == "deviance") {
        limit <- 1e-05
        fit[[i]][fit[[i]] < limit] <- limit
        fit[[i]][fit[[i]] > 1 - limit] <- 1 - limit
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) mean(-2 * (y * log(x) + (1 - 
                                                                        y) * log(1 - x))))
      }
      else if (type.measure == "mse") {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) 2 * mean((x - y)^2))
      }
      else if (type.measure == "mae") {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) 2 * mean(abs(x - y)))
      }
      else if (type.measure == "class") {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) mean(abs(round(x) - y)))
      }
      else if (type.measure == "auc") {
        weights <- table(foldid)
        cvraw <- matrix(data = NA, nrow = length(weights), 
                        ncol = ncol(fit[[i]])) # typo in palasso package !
        for (k in seq_along(weights)) {
          cvraw[k, ] <- apply(X = fit[[i]], MARGIN = 2, 
                              FUN = function(x) glmnet::auc(y = y[foldid == 
                                                                    k], prob = x[foldid == k]))
        }
        loss[[i]] <- apply(X = cvraw, MARGIN = 2, FUN = function(x) stats::weighted.mean(x = x, 
                                                                                         w = weights, na.rm = TRUE))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
742
        names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
      }
      else {
        stop("Invalid type measure.", call. = FALSE)
      }
    }
    else if (family == "poisson") {
      if (type.measure == "deviance") {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) mean(2 * (ifelse(y == 0, 
                                                              0, y * log(y/x)) - y + x), na.rm = TRUE))
      }
      else if (type.measure == "mse") {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) mean((x - y)^2))
      }
      else if (type.measure == "mae") {
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) mean(abs(x - y)))
      }
      else {
        stop("Invalid type measure.", call. = FALSE)
      }
    }
    else if (family == "cox") {
      if (type.measure == "deviance") {
        weights <- tapply(X = y[, "status"], INDEX = foldid, 
                          FUN = sum)
        loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, 
                           FUN = function(x) stats::weighted.mean(x = x/weights, 
                                                                  w = weights, na.rm = TRUE))
      }
      else {
        stop("Invalid type measure.", call. = FALSE)
      }
    }
    else {
      stop("Invalid family.", call. = FALSE)
    }
    if (sum(diff(is.na(loss[[i]]))) == 1) {
      loss[[i]] <- loss[[i]][!is.na(loss[[i]])]
    }
  }
  return(loss)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
788
789
#--- Lost and found ------------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
790
# calibrate (for cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
791
792
793
#if(test$calibrate){
#  fit$calibrate <- CalibratR::calibrate(actual=z,predicted=pred$y[,which.min(fit$gaussian$cvm)],nCores=1,model_idx=5)$calibration_models
#}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
794
# calibrate (for predict.cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
#if(test$calibrate){
#  prob$calibrate <- CalibratR::predict_calibratR(calibration_models=x$calibrate,new=link,nCores=1)$GUESS_2
#}

.args <- function(...){
  args <- list(...)
  names <- names(formals(glmnet::glmnet))
  if(!is.null(args$family)){
    warning("Unexpected argument \"family\".",call.=FALSE) 
  }
  if(any(!names(args) %in% names)){
    stop("Unexpected argument.",call.=FALSE)
  }
  if(is.null(args$alpha)) {
    args$alpha <- 1
  }
  if(is.null(args$nlambda)){
    args$nlambda <- 100
  }
  if(is.null(args$lambda)){
    if(is.null(args$nlambda)){
      args$nlambda <- 100
    }
  } else {
    args$nlambda <- length(args$lambda)
  }
  return(args)
}




Armin Rauschenberger's avatar
Armin Rauschenberger committed
827