functions.R 25.6 KB
Newer Older
Armin Rauschenberger's avatar
Armin Rauschenberger committed
1

Armin Rauschenberger's avatar
Armin Rauschenberger committed
2
#--- Workhorse function --------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
3

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
167
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
168

Armin Rauschenberger's avatar
Armin Rauschenberger committed
169 170 171 172 173 174 175
    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
176
    # linear regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
177
    net <- glmnet::glmnet(y=y0,x=X0,family="gaussian",alpha=alpha,...)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
178 179
    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
180
    cvm <- palasso:::.loss(y=y1,fit=temp_y,family="gaussian",type.measure="deviance")[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
181
    y_hat <- temp_y[,which.min(cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
182
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
183
    # logistic regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
184
    net <- glmnet::glmnet(y=z0,x=X0,family="binomial",alpha=alpha,...)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
185 186
    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
187
    cvm <- palasso:::.loss(y=z1,fit=temp_z,family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
188 189
    z_hat <- temp_z[,which.min(cvm)]
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
190
    # combined regression
Armin Rauschenberger's avatar
Armin Rauschenberger committed
191
    if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
192 193 194
      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
195
          pred$combined[foldid==k,i,j] <- fit$pi[j]*cont + (1-fit$pi[j])*z_hat
Armin Rauschenberger's avatar
Armin Rauschenberger committed
196 197 198
        }
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
199
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
200 201
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
202 203
  #--- evaluation ---
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
204
  # linear loss
Armin Rauschenberger's avatar
Armin Rauschenberger committed
205
  fit$gaussian$cvm <- palasso:::.loss(y=y,fit=pred$y,family="gaussian",type.measure="deviance")[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
206
  fit$gaussian$lambda.min <- fit$gaussian$lambda[which.min(fit$gaussian$cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
207
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
208
  # logistic loss
Armin Rauschenberger's avatar
Armin Rauschenberger committed
209
  fit$binomial$cvm <- palasso:::.loss(y=z,fit=pred$z,family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
210
  fit$binomial$lambda.min <- fit$binomial$lambda[which.min(fit$binomial$cvm)]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
211

Armin Rauschenberger's avatar
Armin Rauschenberger committed
212
  # combined loss
Armin Rauschenberger's avatar
Armin Rauschenberger committed
213
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
214
    dimnames <- list(lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
215
    fit$cvm <- matrix(data=NA,nrow=nsigma,ncol=npi,dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
216
    for(i in seq_len(nsigma)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
217
      for(j in seq_len(npi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
218
        fit$cvm[i,j] <- palasso:::.loss(y=z,fit=pred$combined[,i,j],family="binomial",type.measure=type.measure)[[1]]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
219 220
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
221
    temp <- which(fit$cvm==min(fit$cvm),arr.ind=TRUE,useNames=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
222
    if(nrow(temp)>1){warning("MULTIPLE!",call.=FALSE);temp <- temp[1,,drop=FALSE]}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
223 224 225
    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
226 227
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
228
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
229
  fit$cutoff <- cutoff
Armin Rauschenberger's avatar
Armin Rauschenberger committed
230 231
  fit$info <- list(type.measure=type.measure,
                   sd.y=stats::sd(y),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
232 233
                   "+"=sum(z==1),
                   "-"=sum(z==0),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
234
                   table=table(z),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
235
                   n=n,p=ncol(X),
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 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
#' @export
#' @title
#' Combined regression
#'
#' @description
#' Prints summary of cornet object.
#'
#' @param x
#' \link[cornet]{cornet} object
#' 
#' @param ...
#' further arguments (not applicable)
#' 
#' @return
#' Returns sample size \eqn{n},
#' number of covariates \eqn{p},
#' information on dichotomisation,
#' tuned scaling parameter (sigma),
#' tuned weighting parameter (pi),
#' and corresponding loss.
#'
#' @examples
#' NA
#' 
print.cornet <- function(x,...){
  cat("cornet object:\n")
  cat(paste0("n = ",x$info$n,","," p = ",x$info$p),"\n")
  cat(paste0("z = I(y > ",signif(x$cutoff,digits=2),"): "))
  cat(paste0(x$info$"+","+"," vs ",x$info$"-","-"),"\n")
  cat(paste0("sigma.min = ",signif(x$sigma.min,digits=1)),"\n")
  cat(paste0("pi.min = ",round(x$pi.min,digits=2)),"\n")
  type <- x$info$type.measure
  value <- signif(x$cvm[names(x$sigma.min),names(x$pi.min)],digits=2)
  cat(paste0(type," = ",value))
  return(invisible(NULL))
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
281 282
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
283
#' Extract estimated coefficients
Armin Rauschenberger's avatar
Armin Rauschenberger committed
284 285
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
286 287
#' 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
288 289
#'
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
290
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
291 292
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
293 294 295 296 297 298 299
#' 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
300 301 302 303
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
304
coef.cornet <- function(object,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
305
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
306
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
307
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
308 309 310 311 312
  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
313 314 315 316 317 318
  
  coef <- cbind(beta,gamma)
  colnames(coef) <- c("beta","gamma")
  return(coef)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
319 320
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
321
#' Plot loss matrix
Armin Rauschenberger's avatar
Armin Rauschenberger committed
322 323
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
324 325
#' Plots the loss for different combinations of
#' scaling (sigma) and weighting (pi) parameters.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
326 327
#'
#' @param x
Armin Rauschenberger's avatar
Armin Rauschenberger committed
328
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
329 330
#' 
#' @param ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
331
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
332 333 334 335
#'
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
336
plot.cornet <- function(x,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
337
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
338
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
339

Armin Rauschenberger's avatar
Armin Rauschenberger committed
340
  k <- 100
Armin Rauschenberger's avatar
Armin Rauschenberger committed
341
  levels <- stats::quantile(x$cvm,probs=seq(from=0,to=1,length.out=k+1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
342
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
343
  # colours
Armin Rauschenberger's avatar
Armin Rauschenberger committed
344 345 346 347 348 349 350
  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
351
  nsigma <- length(x$sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
352
  npi <- length(x$pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
353 354 355
  
  graphics::plot.new()
  graphics::par(xaxs="i",yaxs="i")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
356
  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
357
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
358
  graphics::title(xlab=expression(sigma),ylab=expression(pi),cex.lab=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
359
  #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
360
  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
361
  graphics::box()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
362
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
363 364 365 366
  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
367
    # axes with labels for tuned parameters
Armin Rauschenberger's avatar
Armin Rauschenberger committed
368 369
    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
370
    # point for tuned parameters
Armin Rauschenberger's avatar
Armin Rauschenberger committed
371
    graphics::points(x=ssigma,y=spi,pch=4,col="white",cex=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
372
  } else {
Armin Rauschenberger's avatar
Armin Rauschenberger committed
373
    # axes with standard labels
Armin Rauschenberger's avatar
Armin Rauschenberger committed
374 375 376 377
    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
378 379 380
    # 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
381
    graphics::points(x=isigma,y=ipi,pch=4,col="white",cex=1)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
382 383
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
384 385
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
386 387
#' @export
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
388
#' Predict binary outcome
Armin Rauschenberger's avatar
Armin Rauschenberger committed
389 390
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
391 392 393 394 395 396 397 398
#' 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
399
#' @param object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
400
#' \link[cornet]{cornet} object
Armin Rauschenberger's avatar
Armin Rauschenberger committed
401 402 403 404 405 406 407 408 409 410
#' 
#' @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
411
#' further arguments (not applicable)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
412 413 414 415
#' 
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
416
predict.cornet <- function(object,newx,type="probability",...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
417
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
418
  if(length(list(...))!=0){warning("Ignoring arguments.")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
419 420
  
  x <- object; rm(object)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
421
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
422 423
  test <- x$info$test
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
424 425
  .check(x=newx,type="matrix")
  .check(x=type,type="string",values=c("probability","odds","log-odds"))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
426
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
427
  # linear and logistic
Armin Rauschenberger's avatar
Armin Rauschenberger committed
428 429 430 431 432 433
  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
434
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
435
  # combined
Armin Rauschenberger's avatar
Armin Rauschenberger committed
436
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
437
    cont <- stats::pnorm(q=link,mean=x$cutoff,sd=x$sigma.min)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
438
    prob$combined <- x$pi.min*cont + (1-x$pi.min)*prob$binomial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
439 440
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
441 442
  # consistency tests
  lapply(X=prob,FUN=function(p) .check(x=p,type="vector",min=0,max=1))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
443
  .equal(link>x$cutoff,prob$gaussian>0.5)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
444 445 446 447 448 449 450 451 452 453

  # 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
454 455
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
456
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
457 458
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
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 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
#--- Internal functions --------------------------------------------------------

#' @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 
#' cornet:::.equal(1,1,1)
#' 
.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
#' cornet:::.check(0.5,type="scalar",min=0,max=1)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
527 528
#' #.check(x=c(1,2,3,4,45),type="vector",values=1:10)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556
.check <- function(x,type,miss=FALSE,min=NULL,max=NULL,values=NULL,inf=FALSE,null=FALSE){
  name <- deparse(substitute(x))
  if(null && is.null(x)){
    return(invisible(NULL)) 
  }
  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)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
557
  if(!is.null(values) && any(!x %in% values)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
558 559 560 561 562 563 564 565
    stop(paste0("Argument \"",name,"\" contains invalid values."),call.=FALSE)
  }
  if(!inf && any(is.infinite(values))){
    stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE)
  }
  return(invisible(NULL))
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
566
#--- Application ---------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
567

Armin Rauschenberger's avatar
Armin Rauschenberger committed
568
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
569
#' Performance measurement by cross-validation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
570 571
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
572 573 574 575 576 577 578 579
#' 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
580
#' @inheritParams  cornet
Armin Rauschenberger's avatar
Armin Rauschenberger committed
581
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
582 583 584
#' @examples
#' NA
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
585
.compare <- function(y,cutoff,X,alpha=1,nfolds=5,foldid=NULL,type.measure="deviance"){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
586
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
587
  z <- 1*(y > cutoff)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
588
  if(is.null(foldid)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
589
    fold <- palasso:::.folds(y=z,nfolds=nfolds)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
590 591 592
  } else {
    fold <- foldid
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
593
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
594 595 596
  #--- cross-validated loss ---
  
  cols <- c("gaussian","binomial","combined")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
597 598
  pred <- matrix(data=NA,nrow=length(y),ncol=length(cols),
                 dimnames=list(NULL,cols))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
599
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
600
  for(i in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
601
    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
602 603
    tryCatch(expr=plot.cornet(fit),error=function(x) NULL)
    temp <- predict.cornet(fit,newx=X[fold==i,])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
604
    if(any(temp<0|temp>1)){stop("Outside unit interval.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
605 606 607 608
    model <- colnames(pred)
    for(j in seq_along(model)){
      pred[fold==i,model[j]] <- temp[[model[j]]]
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
609 610
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
611
  type <- c("deviance","class","mse","mae","auc")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
612
  loss <- lapply(X=type,FUN=function(x) palasso:::.loss(y=z,fit=pred,family="binomial",type.measure=x,foldid=fold)[[1]])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
613
  names(loss) <- type
Armin Rauschenberger's avatar
Armin Rauschenberger committed
614
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
615
  #--- deviance residuals ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
616 617 618 619 620 621 622
  
  # 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
623 624 625 626
  rys <- res[,"combined"]
  
  # residual increase/decrease
  loss$resid.factor <- stats::median((rys-rxs)/rxs)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
627
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
628 629
  # paired test for each fold
  loss$resid.pvalue <- numeric()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
630 631
  for(i in seq_len(nfolds)){
    cond <- fold==i
Armin Rauschenberger's avatar
Armin Rauschenberger committed
632
    loss$resid.pvalue[i] <- stats::wilcox.test(x=rxs[cond],y=rys[cond],
Armin Rauschenberger's avatar
Armin Rauschenberger committed
633
                                               paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
634 635
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
636
  return(loss)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
637
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
638 639 640
}

#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
641
#' Single-split test
Armin Rauschenberger's avatar
Armin Rauschenberger committed
642 643
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
644 645 646 647 648 649 650 651
#' 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
652 653 654 655 656 657 658 659
#' @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
660
  fold <- palasso:::.folds(y=z,nfolds=5)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
661 662 663
  fold <- ifelse(fold==1,1,0)
  
  fit <- cornet::cornet(y=y[fold==0],cutoff=cutoff,X=X[fold==0,],alpha=alpha)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
664 665
  tryCatch(expr=plot.cornet(fit),error=function(x) NULL)
  pred <- predict.cornet(fit,newx=X[fold==1,])
Armin Rauschenberger's avatar
Armin Rauschenberger committed
666 667 668
  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
669
  #pvalue <- wilcox.test(x=res[,"binomial"],y=res[,"combined"],paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
670
  #colMeans(abs(pred-0.5)) # distance from 0.5
Armin Rauschenberger's avatar
Armin Rauschenberger committed
671
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
672 673 674
  limit <- 1e-05
  pred[pred < limit] <- limit
  pred[pred > 1 - limit] <- 1 - limit
Armin Rauschenberger's avatar
Armin Rauschenberger committed
675
  res <- -2 * (z[fold==1] * log(pred) + (1 - z[fold==1]) * log(1 - pred))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
676
  pvalue <- stats::wilcox.test(x=res[,"binomial"],y=res[,"combined"],paired=TRUE,alternative="greater")$p.value
Armin Rauschenberger's avatar
Armin Rauschenberger committed
677 678 679 680
  
  return(pvalue)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
681
#' @title
Armin Rauschenberger's avatar
Armin Rauschenberger committed
682
#' Data simulation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
683 684
#'
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
685
#' Simulates data for unit tests
Armin Rauschenberger's avatar
Armin Rauschenberger committed
686
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
687 688 689
#' @param n
#' sample size\strong{:}
#' positive integer
Armin Rauschenberger's avatar
Armin Rauschenberger committed
690
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
691 692 693
#' @param p
#' covariate space\strong{:}
#' positive integer
Armin Rauschenberger's avatar
Armin Rauschenberger committed
694
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
695 696 697
#' @param prob
#' (approximate) proportion of causal covariates\strong{:}
#' numeric between \eqn{0} and \eqn{1}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
698
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
699 700 701
#' @param fac
#' noise factor\strong{:}
#' positive real number
Armin Rauschenberger's avatar
Armin Rauschenberger committed
702
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
703 704
#' @return
#' Returns invisible list with elements \code{y} and \code{X}.
Armin Rauschenberger's avatar
Armin Rauschenberger committed
705 706
#' 
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
707 708
#' data <- cornet:::.simulate(n=10,p=20,prob=0.2,fac=2)
#' names(data)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
709
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
710 711 712 713 714 715 716 717
.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(invisible(list(y=y,X=X)))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
718 719
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
720
#--- Legacy --------------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
721

Armin Rauschenberger's avatar
Armin Rauschenberger committed
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 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845
# # Import this function from the palasso package.
# .loss <- function (y,fit,family,type.measure,foldid=NULL){
#   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))
#         names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package!
#       }
#       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)
# }
# 
# # 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)
# }