functions.R 19.8 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 foldid
#' fold identifiers\strong{:}
#' vector with entries between \eqn{1} and \code{nfolds};
#' or \code{NULL} (balance)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
28
29
30
31
#' 
#' @param nfolds
#' number of folds
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
32
#' @param type.measure
Armin Rauschenberger's avatar
Armin Rauschenberger committed
33
34
35
#' loss function for binary classification
#' (linear regression uses the deviance)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
36
37
38
39
40
41
42
43
#' @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
44
45
46
47
48
49
50
51
#' @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
52
53
54
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}}
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
55
#' @details
Armin Rauschenberger's avatar
Armin Rauschenberger committed
56
57
58
59
#' - 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
60
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
61
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
62
63
64
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
65
#' net <- cornet(y=y,cutoff=0,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
66
#' ### Add ... to all glmnet::glmnet calls !!! ###
Armin Rauschenberger's avatar
Armin Rauschenberger committed
67
cornet <- function(y,cutoff,X,npi=101,pi=NULL,nsigma=99,sigma=NULL,nfolds=10,foldid=NULL,type.measure="deviance",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
68
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
69
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
70
  # 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
71
  test <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
72
  test$sigma <- test$pi <- FALSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
73
  test$grid <- TRUE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
74
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
75
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
76
  cornet:::.check(x=y,type="vector")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
77
  if(all(y %in% c(0,1))){warning("Binary response.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
78
79
  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
80
  if(length(y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
81
82
83
84
85
86
87
  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
88
89
  if(!is.null(list(...)$family)){stop("Reserved argument \"family\".",call.=FALSE)}
  n <- length(y)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
90
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
91
  # binarisation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
92
93
  z <- 1*(y > cutoff)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
94
  # fold identifiers
Armin Rauschenberger's avatar
Armin Rauschenberger committed
95
96
97
98
  if(is.null(foldid)){
    foldid <- palasso:::.folds(y=z,nfolds=nfolds)
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
99
  #--- model fitting ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
100
  fit <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
101
  fit$gaussian <- glmnet::glmnet(y=y,x=X,family="gaussian")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
102
103
104
  fit$binomial <- glmnet::glmnet(y=z,x=X,family="binomial")
   
  #--- sigma sequence ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
105
106
  if(is.null(sigma)){
    fit$sigma <- exp(seq(from=log(0.05*stats::sd(y)),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
107
                  to=log(10*stats::sd(y)),length.out=nsigma))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
108
109
  } else {
    fit$sigma <- sigma
Armin Rauschenberger's avatar
Armin Rauschenberger committed
110
    nsigma <- length(sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
111
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
112
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
113
114
  #--- pi sequence ---
  if(is.null(pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
115
    fit$pi <- seq(from=0,to=1,length.out=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
116
117
118
  } else {
    fit$pi <- pi
    npi <- length(pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
119
120
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
121
122
123
124
125
  #--- 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
126
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
127
  #--- cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
128
  pred <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
129
130
131
  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
132
133
134
135
  if(test$sigma){
    pred$sigma <- matrix(data=NA,nrow=n,ncol=nsigma)
  }
  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
136
    pred$pi <- matrix(data=NA,nrow=n,ncol=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
137
138
  }
  if(test$grid){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
139
    dimnames <- list(NULL,lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
140
    pred$grid <- array(data=NA,dim=c(n,nsigma,npi),dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
141
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
142

Armin Rauschenberger's avatar
Armin Rauschenberger committed
143
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
144

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

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

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

  if(test$pi){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
207
    #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
208
    #fit$pi.min1 <- fit$pi[which.min(fit$pi.cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
209
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
210

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
233
  class(fit) <- "cornet"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
234
235
236
  return(fit)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
237
#--- Methods -------------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
238

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
291
  k <- 100
Armin Rauschenberger's avatar
Armin Rauschenberger committed
292
  levels <- stats::quantile(x$cvm,probs=seq(from=0,to=1,length.out=k+1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
293
294
  col <- colorspace::diverge_hsv(n=k)
  nsigma <- length(x$sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
295
  npi <- length(x$pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
296
297
298
  
  graphics::plot.new()
  graphics::par(xaxs="i",yaxs="i")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
299
  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
300
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
301
  sel <- which(x$sigma==x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
302
303
  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
304
  sel <- which(x$pi==x$pi.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
305
  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
306
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
307
308
  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
309
  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
310
  graphics::box()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
311
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
312
313
}

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

  # 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
383
384
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
385
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
386
387
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
388
389
#--- Internal functions --------------------------------------------------------

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
433
434
435
  return(loss)
}

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
539
540
541
542
543
544
545
546
547
548
# 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
549
550
551

# 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
552
553
554
555
556
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
  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
609
        names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
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
652
653
654
      }
      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
655
656
#--- Lost and found ------------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
657
# calibrate (for cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
658
659
660
#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
661
# calibrate (for predict.cornet)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
662
663
664
665
666
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
#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
694