functions.R 20.4 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
246
247
248
249
250
251
#' @export
#' @title
#' to do
#'
#' @description
#' to do
#'
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
252
#' cornet object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
253
254
255
256
257
258
259
#' 
#' @param ...
#' to do
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
260
coef.cornet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
261
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
262
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
263
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
264
265
266
267
268
  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
269
270
271
272
273
274
  
  coef <- cbind(beta,gamma)
  colnames(coef) <- c("beta","gamma")
  return(coef)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291

#' @export
#' @title
#' to do
#'
#' @description
#' to do
#'
#' @param x
#' to do
#' 
#' @param ...
#' further arguments
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
292
plot.cornet <- function(x,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
293
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
294
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
295

Armin Rauschenberger's avatar
Armin Rauschenberger committed
296
  k <- 100
Armin Rauschenberger's avatar
Armin Rauschenberger committed
297
  levels <- stats::quantile(x$cvm,probs=seq(from=0,to=1,length.out=k+1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
298
299
  col <- colorspace::diverge_hsv(n=k)
  nsigma <- length(x$sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
300
  npi <- length(x$pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
301
302
303
  
  graphics::plot.new()
  graphics::par(xaxs="i",yaxs="i")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
304
  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
305
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
306
307
  ssigma <- which(x$sigma==x$sigma.min)
  graphics::axis(side=1,at=c(1,ssigma,nsigma),labels=signif(x$sigma[c(1,ssigma,nsigma)],digits=2))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
308
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
309
310
  spi <- which(x$pi==x$pi.min)
  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
311
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
312
313
  graphics::title(xlab=expression(sigma),ylab=expression(pi))
  #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
314
  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
315
  graphics::box()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
316
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
317
318
319
320
321
  #graphics::abline(v=ssigma,lty=2,col="grey")
  #graphics::abline(h=spi,lty=2,col="grey")
  
  graphics::points(x=ssigma,y=spi,pch=4,col="black",cex=1)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
322
323
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
324
325
326
327
328
329
330
331
#' @export
#' @title
#' to do
#'
#' @description
#' to do
#'
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
332
#' cornet object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
#' 
#' @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 ...
#' to do
#' 
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
348
predict.cornet <- function(object,newx,type="probability",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
349
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
350
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
351
352
  
  x <- object; rm(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
353
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
354
355
  test <- x$info$test
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
356
357
  .check(x=newx,type="matrix")
  .check(x=type,type="string",values=c("probability","odds","log-odds"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
358
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
359
360
361
362
363
364
365
  # 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
366
367
  
  if(test$sigma){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
368
    #prob$sigma <- stats::pnorm(q=link,mean=x$cutoff,sd=x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
369
370
371
  }
  
  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
372
     #prob$pi <- x$pi.min*prob$gaussian + (1-x$pi.min)*prob$binomial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
373
374
375
  }
 
  if(test$grid){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
376
    cont <- stats::pnorm(q=link,mean=x$cutoff,sd=x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
377
    prob$grid <- x$pi.min*cont + (1-x$pi.min)*prob$binomial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
378
379
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
380
381
  # consistency tests
  lapply(X=prob,FUN=function(p) .check(x=p,type="vector",min=0,max=1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
382
  .equal(link>x$cutoff,prob$gaussian>0.5)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
383
384
385
386
387
388
389
390
391
392

  # 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
393
394
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
395
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
396
397
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
398
399
#--- Internal functions --------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
400
401
402
403
404
405
406
#' @export
#' @title
#' Comparison
#'
#' @description
#' Compares models for a continuous response with a cutoff value
#'
Armin Rauschenberger's avatar
Armin Rauschenberger committed
407
#' @inheritParams  cornet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
408
409
410
411
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
412
.compare <- function(y,cutoff,X,alpha=1,nfolds=5,foldid=NULL,type.measure="deviance"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
413
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
414
  z <- 1*(y > cutoff)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
415
416
417
418
419
  if(is.null(foldid)){
    fold <- palasso:::.folds(y=z,nfolds=nfolds)
  } else {
    fold <- foldid
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
420
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
421
422
  #cols <- c("gaussian","binomial","sigma","pi","trial","grid","unit")
  cols <- c("gaussian","binomial","grid")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
423
424
  pred <- matrix(data=NA,nrow=length(y),ncol=length(cols),
                 dimnames=list(NULL,cols))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
425
426
  
  select <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
427
  for(i in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
428
    fit <- cornet::cornet(y=y[fold!=i],cutoff=cutoff,X=X[fold!=i,],alpha=alpha,logistic=TRUE,type.measure=type.measure)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
429
430
431
    tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL)
    #cornet:::plot.cornet(fit)
    temp <- cornet:::predict.cornet(fit,newx=X[fold==i,])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
432
    if(any(temp<0|temp>1)){stop("Outside unit interval.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
433
434
435
436
    model <- colnames(pred)
    for(j in seq_along(model)){
      pred[fold==i,model[j]] <- temp[[model[j]]]
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
437
438
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
439
  type <- c("deviance","class","mse","mae","auc")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
440
  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
441
  names(loss) <- type
Armin Rauschenberger's avatar
Armin Rauschenberger committed
442

Armin Rauschenberger's avatar
Armin Rauschenberger committed
443
444
445
  return(loss)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
#' @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
462
463
464
465
466
467
468
469
470
471
472
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
504
505
506
507
508
509
510
511
512
.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
513
  name <- deparse(substitute(x))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
514
515
516
  if(null && is.null(x)){
    return(invisible(NULL)) 
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
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
542
  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
543
544
545
  if(!inf && any(is.infinite(values))){
    stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
546
547
548
  return(invisible(NULL))
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
549
550
551
552
553
554
555
556
557
558
# 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
559

Armin Rauschenberger's avatar
Armin Rauschenberger committed
560
561
562
563
564
565
566
567
568
569
570
571
572
#--- 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 ---



Armin Rauschenberger's avatar
Armin Rauschenberger committed
573
574
# 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
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
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
  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
632
        names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
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
674
675
676
677
      }
      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
678
679
#--- Lost and found ------------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
680
# calibrate (for cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
681
682
683
#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
684
# calibrate (for predict.cornet)
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
#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
717