functions.R 20 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
  sel <- which(x$sigma==x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
307
308
  graphics::axis(side=1,at=c(1,sel,nsigma),labels=signif(x$sigma[c(1,sel,nsigma)],digits=2))
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
309
  sel <- which(x$pi==x$pi.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
310
  graphics::axis(side=2,at=c(1,sel,npi),labels=signif(x$pi[c(1,sel,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
}

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

  # 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
388
389
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
390
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391
392
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
393
394
#--- Internal functions --------------------------------------------------------

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
438
439
440
  return(loss)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
#' @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
457
458
459
460
461
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
.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
508
  name <- deparse(substitute(x))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
509
510
511
  if(null && is.null(x)){
    return(invisible(NULL)) 
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
  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
538
539
540
  if(!inf && any(is.infinite(values))){
    stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
541
542
543
  return(invisible(NULL))
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
544
545
546
547
548
549
550
551
552
553
# 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
554
555
556

# 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
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
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
  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
614
        names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
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
652
653
654
655
656
657
658
659
      }
      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
660
661
#--- Lost and found ------------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
662
# calibrate (for cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
663
664
665
#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
666
# calibrate (for predict.cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
#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
699