functions.R 24.6 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
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
71
#' ### Add ... to all glmnet::glmnet calls !!! ###
Armin Rauschenberger's avatar
Armin Rauschenberger committed
72
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
73
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
74
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
75
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
76
  # 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
77
  test <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
78
  test$sigma <- test$pi <- FALSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
79
  test$combined <- TRUE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
80
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
81
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
82
  cornet:::.check(x=y,type="vector")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
83
  if(all(y %in% c(0,1))){warning("Binary response.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
84
85
  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
86
  cornet:::.check(x=alpha,type="scalar",min=0,max=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
87
  if(length(y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
88
89
90
91
92
93
94
  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
95
96
  if(!is.null(list(...)$family)){stop("Reserved argument \"family\".",call.=FALSE)}
  n <- length(y)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
97
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
98
  # binarisation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
99
100
  z <- 1*(y > cutoff)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
101
  # fold identifiers
Armin Rauschenberger's avatar
Armin Rauschenberger committed
102
103
104
105
  if(is.null(foldid)){
    foldid <- palasso:::.folds(y=z,nfolds=nfolds)
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
106
  #--- model fitting ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
107
  fit <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
108
109
  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
110
111
   
  #--- sigma sequence ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
112
113
  if(is.null(sigma)){
    fit$sigma <- exp(seq(from=log(0.05*stats::sd(y)),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
114
                  to=log(10*stats::sd(y)),length.out=nsigma))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
115
116
  } else {
    fit$sigma <- sigma
Armin Rauschenberger's avatar
Armin Rauschenberger committed
117
    nsigma <- length(sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
118
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
119
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
120
121
  #--- pi sequence ---
  if(is.null(pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
122
    fit$pi <- seq(from=0,to=1,length.out=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
123
124
125
  } else {
    fit$pi <- pi
    npi <- length(pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
126
127
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
128
129
130
131
132
  #--- 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
133
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
134
  #--- cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
135
  pred <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
136
137
138
  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
139
140
141
142
  if(test$sigma){
    pred$sigma <- matrix(data=NA,nrow=n,ncol=nsigma)
  }
  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
143
    pred$pi <- matrix(data=NA,nrow=n,ncol=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
144
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
145
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
146
    dimnames <- list(NULL,lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
147
    pred$combined <- array(data=NA,dim=c(n,nsigma,npi),dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
148
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
149

Armin Rauschenberger's avatar
Armin Rauschenberger committed
150
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
151

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
188
189
    # fusion (combined)
    if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
190
191
192
      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
193
          pred$combined[foldid==k,i,j] <- fit$pi[j]*cont + (1-fit$pi[j])*z_hat
Armin Rauschenberger's avatar
Armin Rauschenberger committed
194
195
196
        }
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
197
198
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
199
200
  #--- evaluation ---
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
201
  # deviance (not comparable between Gaussian and binomial families)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
202
  fit$gaussian$cvm <- cornet:::.loss(y=y,fit=pred$y,family="gaussian",type.measure="deviance")[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
203
  fit$gaussian$lambda.min <- fit$gaussian$lambda[which.min(fit$gaussian$cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
204
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
205
  fit$binomial$cvm <- cornet:::.loss(y=z,fit=pred$z,family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
206
  fit$binomial$lambda.min <- fit$binomial$lambda[which.min(fit$binomial$cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
207

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

  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
214
    #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
215
    #fit$pi.min1 <- fit$pi[which.min(fit$pi.cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
216
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
217

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

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
244
#--- Methods -------------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
245

Armin Rauschenberger's avatar
Armin Rauschenberger committed
246
247
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
248
#' Extract estimated coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
249
250
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
251
252
#' 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
253
254
#'
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
255
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
256
257
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
258
259
260
261
262
263
264
#' 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
265
266
267
268
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
269
coef.cornet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
270
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
271
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
272
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
273
274
275
276
277
  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
278
279
280
281
282
283
  
  coef <- cbind(beta,gamma)
  colnames(coef) <- c("beta","gamma")
  return(coef)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
284
285
286

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
306
  k <- 100
Armin Rauschenberger's avatar
Armin Rauschenberger committed
307
  levels <- stats::quantile(x$cvm,probs=seq(from=0,to=1,length.out=k+1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
308
309
310
311
312
313
314
315
316
  
  ## 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
317
  nsigma <- length(x$sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
318
  npi <- length(x$pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
319
320
321
  
  graphics::plot.new()
  graphics::par(xaxs="i",yaxs="i")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
322
  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
323
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
324
  graphics::title(xlab=expression(sigma),ylab=expression(pi),cex.lab=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
325
  #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
326
  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
327
  graphics::box()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
328
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
329
330
331
332
  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
333
    # axes with labels for tuned parameters
Armin Rauschenberger's avatar
Armin Rauschenberger committed
334
335
    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
336
    # point for tuned parameters
Armin Rauschenberger's avatar
Armin Rauschenberger committed
337
    graphics::points(x=ssigma,y=spi,pch=4,col="white",cex=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
338
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
339
    # axes with standard labels
Armin Rauschenberger's avatar
Armin Rauschenberger committed
340
341
342
343
    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
344
345
346
    # 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
347
    graphics::points(x=isigma,y=ipi,pch=4,col="white",cex=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
348
349
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
350
351
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
352
353
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
354
#' Predict binary outcome
Armin Rauschenberger's avatar
Armin Rauschenberger committed
355
356
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
357
358
359
360
361
362
363
364
#' 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
365
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
366
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
367
368
369
370
371
372
373
374
375
376
#' 
#' @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
377
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
378
379
380
381
#' 
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
382
predict.cornet <- function(object,newx,type="probability",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
383
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
384
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
385
386
  
  x <- object; rm(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
387
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
388
389
  test <- x$info$test
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
390
391
  .check(x=newx,type="matrix")
  .check(x=type,type="string",values=c("probability","odds","log-odds"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
392
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
393
394
395
396
397
398
399
  # 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
400
401
  
  if(test$sigma){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
402
    #prob$sigma <- stats::pnorm(q=link,mean=x$cutoff,sd=x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
403
404
405
  }
  
  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
406
     #prob$pi <- x$pi.min*prob$gaussian + (1-x$pi.min)*prob$binomial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
407
408
  }
 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
409
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
410
    cont <- stats::pnorm(q=link,mean=x$cutoff,sd=x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
411
    prob$combined <- x$pi.min*cont + (1-x$pi.min)*prob$binomial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
412
413
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
414
415
  # consistency tests
  lapply(X=prob,FUN=function(p) .check(x=p,type="vector",min=0,max=1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
416
  .equal(link>x$cutoff,prob$gaussian>0.5)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
417
418
419
420
421
422
423
424
425
426

  # 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
427
428
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
429
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
430
431
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
432
#--- Application ---------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
433

Armin Rauschenberger's avatar
Armin Rauschenberger committed
434
435
436
437
438
#' @export
#' @title
#' Comparison
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
439
440
441
442
443
444
445
446
#' Compares models for a continuous response with a cutoff value.
#' 
#' @details
#' Uses k-fold cross-validation,
#' fits linear, logistic, and combined regression,
#' calculates different loss functions,
#' and examines squared deviance residuals.
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
447
#' @inheritParams  cornet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
448
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
449
450
451
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
452
.compare <- function(y,cutoff,X,alpha=1,nfolds=5,foldid=NULL,type.measure="deviance"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
453
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
454
  z <- 1*(y > cutoff)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
455
456
457
458
459
  if(is.null(foldid)){
    fold <- palasso:::.folds(y=z,nfolds=nfolds)
  } else {
    fold <- foldid
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
460
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
461
462
463
  #--- cross-validated loss ---
  
  cols <- c("gaussian","binomial","combined")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
464
465
  pred <- matrix(data=NA,nrow=length(y),ncol=length(cols),
                 dimnames=list(NULL,cols))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
466

Armin Rauschenberger's avatar
Armin Rauschenberger committed
467
  for(i in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
468
    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
469
470
    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
471
    if(any(temp<0|temp>1)){stop("Outside unit interval.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
472
473
474
475
    model <- colnames(pred)
    for(j in seq_along(model)){
      pred[fold==i,model[j]] <- temp[[model[j]]]
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
476
477
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
478
  type <- c("deviance","class","mse","mae","auc")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
479
  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
480
  names(loss) <- type
Armin Rauschenberger's avatar
Armin Rauschenberger committed
481
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
482
  #--- deviance residuals ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
483
484
485
486
487
488
489
  
  # 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"]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
490
491
492
493
  rys <- res[,"combined"]
  
  # residual increase/decrease
  loss$resid.factor <- stats::median((rys-rxs)/rxs)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
494
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
495
496
497
498
499
500
501
502
503
504
505
506
  if(FALSE){# tests
    # equality deviance
    loss$deviance["binomial"]==mean(res[,"binomial"])
    loss$deviance["combined"]==mean(res[,"combined"])
    # percentage decrease
    #range((rys-rxs)/rxs)
    stats::median((rys-rxs)/rxs)
    mean((rys-rxs)/rxs)
    (sum(rys)-sum(rxs))/sum(rxs)
    (loss$deviance["combined"]-loss$deviance["binomial"])/loss$deviance["binomial"]
  } 
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
507
508
  # paired test for each fold
  loss$resid.pvalue <- numeric()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
509
510
  for(i in seq_len(nfolds)){
    cond <- fold==i
Armin Rauschenberger's avatar
Armin Rauschenberger committed
511
512
    loss$resid.pvalue[i] <- stats::wilcox.test(x=rxs[cond],y=rys[cond],
          paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
513
514
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
515
516
517
518
519
520
  return(loss)

}

#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
521
#' Single-split test
Armin Rauschenberger's avatar
Armin Rauschenberger committed
522
523
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
524
525
526
527
528
529
530
531
#' Compares models for a continuous response with a cutoff value.
#' 
#' @details
#' Splits samples into 80% for training and 20% for testing,
#' calculates squared deviance residuals of logistic and combined regression,
#' conducts the paired one-sided Wilcoxon signed rank test,
#' and returns the p-value.
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
#' @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
Armin Rauschenberger's avatar
Armin Rauschenberger committed
549
  #pvalue <- wilcox.test(x=res[,"binomial"],y=res[,"combined"],paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
550
  #colMeans(abs(pred-0.5)) # distance from 0.5
Armin Rauschenberger's avatar
Armin Rauschenberger committed
551
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
552
553
554
  limit <- 1e-05
  pred[pred < limit] <- limit
  pred[pred > 1 - limit] <- 1 - limit
Armin Rauschenberger's avatar
Armin Rauschenberger committed
555
  res <- -2 * (z[fold==1] * log(pred) + (1 - z[fold==1]) * log(1 - pred))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
556
  pvalue <- stats::wilcox.test(x=res[,"binomial"],y=res[,"combined"],paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
557
558
559
560
561
562
563
564
565
566
567
568
569
  
  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
570
571
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
572
573
574
575
576
577
578
579
580
581
582
583
584
#--- 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
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
#' @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
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
644
645
646
647
648
649
650
651
.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
652
  name <- deparse(substitute(x))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
653
654
655
  if(null && is.null(x)){
    return(invisible(NULL)) 
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
  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
682
683
684
  if(!inf && any(is.infinite(values))){
    stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
685
686
687
  return(invisible(NULL))
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
688
689
690



Armin Rauschenberger's avatar
Armin Rauschenberger committed
691
692
# 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
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
742
743
744
745
746
747
748
749
  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
750
        names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
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
788
789
790
791
792
793
794
795
      }
      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
796
797
#--- Lost and found ------------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
798
# calibrate (for cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
799
800
801
#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
802
# calibrate (for predict.cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
#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
835