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

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
4
#' @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
#' cut-off 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
92
93
94
95
#' @references 
#' A Rauschenberger, E Glaab (2019).
#' "Lasso and ridge regression with dichotomised outcomes".
#' \emph{Manuscript in preparation}.
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
96
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
97
98
99
#' n <- 100; p <- 200
#' y <- rnorm(n)
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
100
#' net <- cornet(y=y,cutoff=0,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
101
#' net
Armin Rauschenberger's avatar
Armin Rauschenberger committed
102
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
103
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
104
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
105
  #--- temporary ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
106
  # 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
107
  test <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
108
  test$combined <- TRUE
Armin Rauschenberger's avatar
Armin Rauschenberger committed
109
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
110
  #--- checks ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
111
  .check(x=y,type="vector")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
112
  if(all(y %in% c(0,1))){warning("Binary response.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
113
114
115
  .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
116
  if(length(y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
117
118
119
120
121
122
123
  .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
124
125
  if(!is.null(list(...)$family)){stop("Reserved argument \"family\".",call.=FALSE)}
  n <- length(y)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
126
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
127
  # binarisation
Armin Rauschenberger's avatar
Armin Rauschenberger committed
128
129
  z <- 1*(y > cutoff)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
130
  # fold identifiers
Armin Rauschenberger's avatar
Armin Rauschenberger committed
131
  if(is.null(foldid)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
132
    foldid <- palasso:::.folds(y=z,nfolds=nfolds)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
133
134
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
135
  #--- model fitting ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
136
  fit <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
137
138
  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
139
140
   
  #--- sigma sequence ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
141
142
  if(is.null(sigma)){
    fit$sigma <- exp(seq(from=log(0.05*stats::sd(y)),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
143
                  to=log(10*stats::sd(y)),length.out=nsigma))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
144
145
  } else {
    fit$sigma <- sigma
Armin Rauschenberger's avatar
Armin Rauschenberger committed
146
    nsigma <- length(sigma)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
147
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
148
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
149
150
  #--- pi sequence ---
  if(is.null(pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
151
    fit$pi <- seq(from=0,to=1,length.out=npi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
152
153
154
  } else {
    fit$pi <- pi
    npi <- length(pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
155
156
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
157
158
159
160
161
  #--- 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
162
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
163
  #--- cross-validation ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
164
  pred <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
165
166
167
  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
168
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
169
    dimnames <- list(NULL,lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
170
    pred$combined <- array(data=NA,dim=c(n,nsigma,npi),dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
171
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
172

Armin Rauschenberger's avatar
Armin Rauschenberger committed
173
  for(k in seq_len(nfolds)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
174

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
218
  # combined loss
Armin Rauschenberger's avatar
Armin Rauschenberger committed
219
  if(test$combined){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
220
    dimnames <- list(lab.sigma,lab.pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
221
    fit$cvm <- matrix(data=NA,nrow=nsigma,ncol=npi,dimnames=dimnames)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
222
    for(i in seq_len(nsigma)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
223
      for(j in seq_len(npi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
224
        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
225
226
      }
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
227
    temp <- which(fit$cvm==min(fit$cvm),arr.ind=TRUE,useNames=TRUE)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
228
    if(nrow(temp)>1){warning("Multiple!",call.=FALSE);temp <- temp[1,,drop=FALSE]}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
229
230
231
    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
232
233
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
234
  #--- return ---
Armin Rauschenberger's avatar
Armin Rauschenberger committed
235
  fit$cutoff <- cutoff
Armin Rauschenberger's avatar
Armin Rauschenberger committed
236
237
  fit$info <- list(type.measure=type.measure,
                   sd.y=stats::sd(y),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
238
239
                   "+"=sum(z==1),
                   "-"=sum(z==0),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
240
                   table=table(z),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
241
                   n=n,p=ncol(X),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
242
                   test=as.data.frame(test))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
243

Armin Rauschenberger's avatar
Armin Rauschenberger committed
244
  class(fit) <- "cornet"
Armin Rauschenberger's avatar
Armin Rauschenberger committed
245
246
247
  return(fit)
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
248
#--- Methods -------------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
249

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

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

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

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

  # 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
460
461
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
462
  return(as.data.frame(frame))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
463
464
}

Armin Rauschenberger's avatar
Armin Rauschenberger committed
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
527
528
529
530
531
532
#--- 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
533
534
#' #.check(x=c(1,2,3,4,45),type="vector",values=1:10)
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
.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
563
  if(!is.null(values) && any(!x %in% values)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
564
565
566
567
568
569
570
571
    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
572
#--- Application ---------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
573

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

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

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

Armin Rauschenberger's avatar
Armin Rauschenberger committed
726
#--- Legacy --------------------------------------------------------------------
Armin Rauschenberger's avatar
Armin Rauschenberger committed
727

Armin Rauschenberger's avatar
Armin Rauschenberger committed
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
846
847
848
849
850
851
# # 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)
# }