functions.R 23.5 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
#' @export
Armin Rauschenberger's avatar
Armin Rauschenberger committed
5
#' @aliases cornet-package
Armin Rauschenberger's avatar
Armin Rauschenberger committed
6
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
7
#' Logistic regression with a continuous response
Armin Rauschenberger's avatar
Armin Rauschenberger committed
8 9
#' 
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
10
#' Implements logistic regression with a continuous response.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
11 12 13 14 15 16
#'  
#' @param y
#' continuous response\strong{:}
#' vector of length \eqn{n}
#' 
#' @param cutoff
Armin Rauschenberger's avatar
Armin Rauschenberger committed
17
#' cutoff point for dichotomising response into classes\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
18 19 20 21
#' value between \code{min(y)} and \code{max(y)}
#' 
#' @param X
#' covariates\strong{:}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
22
#' numeric matrix with \eqn{n} rows (samples)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
23
#' and \eqn{p} columns (variables)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
24
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
25 26 27 28
#' @param alpha
#' elastic net mixing parameter\strong{:}
#' numeric between \eqn{0} (ridge) and \eqn{1} (lasso)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
29 30 31 32
#' @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
33 34 35 36
#' 
#' @param nfolds
#' number of folds
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
37
#' @param type.measure
Armin Rauschenberger's avatar
Armin Rauschenberger committed
38 39 40
#' loss function for binary classification
#' (linear regression uses the deviance)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
41 42 43 44 45 46 47 48
#' @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
49 50 51 52 53 54 55 56
#' @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
57 58 59
#' @param ...
#' further arguments passed to \code{\link[glmnet]{glmnet}}
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
60
#' @details
Armin Rauschenberger's avatar
Armin Rauschenberger committed
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
#' This function fits a \code{"gaussian"} model for the numeric response,
#' and a \code{"binomial"} model for the binary response,
#' meaning that the \code{glmnet} argument \code{family} is unavailable.
#' Also if \code{type.measure} equals \code{"deviance"},
#' the loss is uncomparable between linear and logistic regression.
#' 
#' @return
#' Returns an object of class \code{cornet}, a list with multiple slots:
#' \itemize{
#'    \item \code{"gaussian"}: fitted linear model, class \code{glmnet}
#'    \item \code{"binomial"}: fitted logistic model, class \code{glmnet}
#'    \item \code{"sigma"}: scaling parameters \code{sigma},
#'           vector of length \code{nsigma}
#'    \item \code{"pi"}: weighting parameters \code{pi},
#'           vector of length \code{npi}
#'    \item \code{cvm}: evaluation loss,
#'           matrix with \code{nsigma} rows and \code{npi} columns
#'    \item \code{sigma.min}: optimal scaling parameter,
#'           positive scalar
#'    \item \code{pi.min}: optimal weighting parameter,
#'           scalar in unit interval
#'    \item \code{cutoff}: threshold for dichotomisation
#' }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
84
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
85
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
86 87 88
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
89
#' net <- cornet(y=y,cutoff=0,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
90
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
91
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
92
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
93
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
94
  # 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
95
  test <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
96
  test$combined <- TRUE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
97
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
98
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
99
  cornet:::.check(x=y,type="vector")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
100
  if(all(y %in% c(0,1))){warning("Binary response.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
101 102
  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
103
  cornet:::.check(x=alpha,type="scalar",min=0,max=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
104
  if(length(y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
105 106 107 108 109 110 111
  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
112 113
  if(!is.null(list(...)$family)){stop("Reserved argument \"family\".",call.=FALSE)}
  n <- length(y)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
114
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
115
  # binarisation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
116 117
  z <- 1*(y > cutoff)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
118
  # fold identifiers
Armin Rauschenberger's avatar
Armin Rauschenberger committed
119
  if(is.null(foldid)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
120
    foldid <- cornet:::.folds(y=z,nfolds=nfolds)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
121 122
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
123
  #--- model fitting ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
124
  fit <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
125 126
  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
127 128
   
  #--- sigma sequence ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
129 130
  if(is.null(sigma)){
    fit$sigma <- exp(seq(from=log(0.05*stats::sd(y)),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
131
                  to=log(10*stats::sd(y)),length.out=nsigma))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
132 133
  } else {
    fit$sigma <- sigma
Armin Rauschenberger's avatar
Armin Rauschenberger committed
134
    nsigma <- length(sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
135
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
136
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
137 138
  #--- pi sequence ---
  if(is.null(pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
139
    fit$pi <- seq(from=0,to=1,length.out=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
140 141 142
  } else {
    fit$pi <- pi
    npi <- length(pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
143 144
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
145 146 147 148 149
  #--- 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
150
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
151
  #--- cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
152
  pred <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
153 154 155
  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
156
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
157
    dimnames <- list(NULL,lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
158
    pred$combined <- array(data=NA,dim=c(n,nsigma,npi),dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
159
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
160

Armin Rauschenberger's avatar
Armin Rauschenberger committed
161
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
162

Armin Rauschenberger's avatar
Armin Rauschenberger committed
163 164 165 166 167 168 169
    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
170
    # linear regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
171
    net <- glmnet::glmnet(y=y0,x=X0,family="gaussian",alpha=alpha,...)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
172 173
    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
174
    cvm <- cornet:::.loss(y=y1,fit=temp_y,family="gaussian",type.measure="deviance")[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
175
    y_hat <- temp_y[,which.min(cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
176
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
177
    # logistic regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
178
    net <- glmnet::glmnet(y=z0,x=X0,family="binomial",alpha=alpha,...)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
179 180
    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
181
    cvm <- cornet:::.loss(y=z1,fit=temp_z,family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
182 183
    z_hat <- temp_z[,which.min(cvm)]
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
184
    # combined regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
185
    if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
186 187 188
      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
189
          pred$combined[foldid==k,i,j] <- fit$pi[j]*cont + (1-fit$pi[j])*z_hat
Armin Rauschenberger's avatar
Armin Rauschenberger committed
190 191 192
        }
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
193
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
194 195
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
196 197
  #--- evaluation ---
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
198
  # linear loss
Armin Rauschenberger's avatar
Armin Rauschenberger committed
199
  fit$gaussian$cvm <- cornet:::.loss(y=y,fit=pred$y,family="gaussian",type.measure="deviance")[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
200
  fit$gaussian$lambda.min <- fit$gaussian$lambda[which.min(fit$gaussian$cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
201
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
202
  # logistic loss
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
  # combined loss
Armin Rauschenberger's avatar
Armin Rauschenberger committed
207
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
208
    dimnames <- list(lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
209
    fit$cvm <- matrix(data=NA,nrow=nsigma,ncol=npi,dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
210
    for(i in seq_len(nsigma)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
211
      for(j in seq_len(npi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
212
        fit$cvm[i,j] <- cornet:::.loss(y=z,fit=pred$combined[,i,j],family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
213 214
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
215
    temp <- which(fit$cvm==min(fit$cvm),arr.ind=TRUE,useNames=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
216
    if(nrow(temp)>1){warning("MULTIPLE!",call.=FALSE);temp <- temp[1,,drop=FALSE]}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
217 218 219
    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
220 221
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
222
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
223
  fit$cutoff <- cutoff
Armin Rauschenberger's avatar
Armin Rauschenberger committed
224 225
  fit$info <- list(type.measure=type.measure,
                   sd.y=stats::sd(y),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
226
                   table=table(z),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
227
                   test=as.data.frame(test))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
228

Armin Rauschenberger's avatar
Armin Rauschenberger committed
229
  class(fit) <- "cornet"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
230 231 232
  return(fit)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
233
#--- Methods -------------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
234

Armin Rauschenberger's avatar
Armin Rauschenberger committed
235 236
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
237
#' Extract estimated coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
238 239
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
240 241
#' Extracts estimated coefficients from linear and logistic regression,
#' under the penalty parameter that minimises the cross-validated loss.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
242 243
#'
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
244
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
245 246
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
247 248 249 250 251 252 253
#' further arguments (not applicable)
#' 
#' @return
#' This function returns a matrix with \eqn{n} rows and two columns,
#' where \eqn{n} is the sample size. It includes the estimated coefficients
#' from linear (first column, \code{"beta"})
#' and logistic (second column, \code{"gamma"}) regression.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
254 255 256 257
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
258
coef.cornet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
259
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
260
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
261
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
262 263 264 265 266
  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
267 268 269 270 271 272
  
  coef <- cbind(beta,gamma)
  colnames(coef) <- c("beta","gamma")
  return(coef)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
273 274 275

#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
276
#' Plot loss matrix
Armin Rauschenberger's avatar
Armin Rauschenberger committed
277 278
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
279 280
#' Plots the loss for difference combinations of
#' the weight (pi) and scale (sigma) paramters.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
281 282
#'
#' @param x
Armin Rauschenberger's avatar
Armin Rauschenberger committed
283
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
284 285
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
286
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
287 288 289 290
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
291
plot.cornet <- function(x,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
292
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
293
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
294

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

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

  # 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
409 410
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
411
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
412 413
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
414
#--- Application ---------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
415

Armin Rauschenberger's avatar
Armin Rauschenberger committed
416 417
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
418
#' Performance measurement by cross-validation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
419 420
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
421 422 423 424 425 426 427 428
#' Compares models for a continuous response with a cutoff value.
#' 
#' @details
#' Uses k-fold cross-validation,
#' fits linear, logistic, and combined regression,
#' calculates different loss functions,
#' and examines squared deviance residuals.
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
429
#' @inheritParams  cornet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
430
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
431 432 433
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
434
.compare <- function(y,cutoff,X,alpha=1,nfolds=5,foldid=NULL,type.measure="deviance"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
435
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
436
  z <- 1*(y > cutoff)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
437
  if(is.null(foldid)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
438
    fold <- cornet:::.folds(y=z,nfolds=nfolds)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
439 440 441
  } else {
    fold <- foldid
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
442
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
443 444 445
  #--- cross-validated loss ---
  
  cols <- c("gaussian","binomial","combined")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
446 447
  pred <- matrix(data=NA,nrow=length(y),ncol=length(cols),
                 dimnames=list(NULL,cols))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
448

Armin Rauschenberger's avatar
Armin Rauschenberger committed
449
  for(i in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
450
    fit <- cornet::cornet(y=y[fold!=i],cutoff=cutoff,X=X[fold!=i,],alpha=alpha,type.measure=type.measure)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
451 452
    tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL)
    temp <- cornet:::predict.cornet(fit,newx=X[fold==i,])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
453
    if(any(temp<0|temp>1)){stop("Outside unit interval.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
454 455 456 457
    model <- colnames(pred)
    for(j in seq_along(model)){
      pred[fold==i,model[j]] <- temp[[model[j]]]
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
458 459
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
460
  type <- c("deviance","class","mse","mae","auc")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
461
  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
462
  names(loss) <- type
Armin Rauschenberger's avatar
Armin Rauschenberger committed
463
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
464
  #--- deviance residuals ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
465 466 467 468 469 470 471
  
  # squared deviance residuals
  limit <- 1e-05
  pred[pred < limit] <- limit
  pred[pred > 1 - limit] <- 1 - limit
  res <- -2 * (z * log(pred) + (1 - z) * log(1 - pred))
  rxs <- res[,"binomial"]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
472 473 474 475
  rys <- res[,"combined"]
  
  # residual increase/decrease
  loss$resid.factor <- stats::median((rys-rxs)/rxs)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
476
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
477 478
  # paired test for each fold
  loss$resid.pvalue <- numeric()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
479 480
  for(i in seq_len(nfolds)){
    cond <- fold==i
Armin Rauschenberger's avatar
Armin Rauschenberger committed
481 482
    loss$resid.pvalue[i] <- stats::wilcox.test(x=rxs[cond],y=rys[cond],
          paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
483 484
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
485 486 487 488 489 490
  return(loss)

}

#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
491
#' Single-split test
Armin Rauschenberger's avatar
Armin Rauschenberger committed
492 493
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
494 495 496 497 498 499 500 501
#' Compares models for a continuous response with a cutoff value.
#' 
#' @details
#' Splits samples into 80% for training and 20% for testing,
#' calculates squared deviance residuals of logistic and combined regression,
#' conducts the paired one-sided Wilcoxon signed rank test,
#' and returns the p-value.
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
502 503 504 505 506 507 508 509
#' @inheritParams cornet
#' 
#' @examples
#' NA
#' 
.test <- function(y,cutoff,X,alpha=1,type.measure="deviance"){
  
  z <- 1*(y > cutoff)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
510
  fold <- cornet:::.folds(y=z,nfolds=5)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
511 512 513 514 515 516 517 518
  fold <- ifelse(fold==1,1,0)
  
  fit <- cornet::cornet(y=y[fold==0],cutoff=cutoff,X=X[fold==0,],alpha=alpha)
  tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL)
  pred <- cornet:::predict.cornet(fit,newx=X[fold==1,])
  if(any(pred<0|pred>1)){stop("Outside unit interval.",call.=FALSE)}
  
  #res <- (pred-z[fold==1])^2 # MSE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
519
  #pvalue <- wilcox.test(x=res[,"binomial"],y=res[,"combined"],paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
520
  #colMeans(abs(pred-0.5)) # distance from 0.5
Armin Rauschenberger's avatar
Armin Rauschenberger committed
521
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
522 523 524
  limit <- 1e-05
  pred[pred < limit] <- limit
  pred[pred > 1 - limit] <- 1 - limit
Armin Rauschenberger's avatar
Armin Rauschenberger committed
525
  res <- -2 * (z[fold==1] * log(pred) + (1 - z[fold==1]) * log(1 - pred))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
526
  pvalue <- stats::wilcox.test(x=res[,"binomial"],y=res[,"combined"],paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
527 528 529 530 531 532 533 534 535 536 537 538 539
  
  return(pvalue)
}

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
542 543
#--- Internal functions --------------------------------------------------------

Armin Rauschenberger's avatar
Armin Rauschenberger committed
544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559
#' @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
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
.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
611
  name <- deparse(substitute(x))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
612 613 614
  if(null && is.null(x)){
    return(invisible(NULL)) 
  }
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
  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
641 642 643
  if(!inf && any(is.infinite(values))){
    stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
644 645 646
  return(invisible(NULL))
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
647
# Import this function from the palasso package.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
648
.loss <- function (y,fit,family,type.measure,foldid=NULL){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
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 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705
  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
706
        names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package!
Armin Rauschenberger's avatar
Armin Rauschenberger committed
707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751
      }
      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
752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770
# Import this function from the palasso package.
.folds <- function (y, nfolds, foldid = NULL){
  if(!is.null(foldid)){
    return(foldid)
  }
  #if (survival::is.Surv(y)){ # active in palasso
  #  y <- y[, "status"] # active in palasso
  #} # active in palasso
  if(all(y %in% c(0, 1))){
    foldid <- rep(x = NA, times = length(y))
    foldid[y == 0] <- sample(x = rep(x = seq_len(nfolds), 
                                     length.out = sum(y == 0)))
    foldid[y == 1] <- sample(x = rep(x = seq_len(nfolds), 
                                     length.out = sum(y == 1)))
  } else {
    foldid <- sample(x = rep(x = seq_len(nfolds), length.out = length(y)))
  }
  return(foldid)
}