functions.R 23.3 KB
Newer Older
 Armin Rauschenberger committed Dec 12, 2018 1 `````` `````` Armin Rauschenberger committed Jan 08, 2019 2 ``````#--- Workhorse function -------------------------------------------------------- `````` Armin Rauschenberger committed Dec 21, 2018 3 `````` `````` Armin Rauschenberger committed Dec 12, 2018 4 5 ``````#' @export #' @title `````` Armin Rauschenberger committed Dec 13, 2018 6 ``````#' Logistic regression with a continuous response `````` Armin Rauschenberger committed Dec 12, 2018 7 8 ``````#' #' @description `````` Armin Rauschenberger committed Dec 19, 2018 9 ``````#' Implements logistic regression with a continuous response. `````` Armin Rauschenberger committed Dec 12, 2018 10 11 12 13 14 15 ``````#' #' @param y #' continuous response\strong{:} #' vector of length \eqn{n} #' #' @param cutoff `````` Armin Rauschenberger committed Dec 19, 2018 16 ``````#' cutoff point for dichotomising response into classes\strong{:} `````` Armin Rauschenberger committed Dec 12, 2018 17 18 19 20 ``````#' value between \code{min(y)} and \code{max(y)} #' #' @param X #' covariates\strong{:} `````` Armin Rauschenberger committed Dec 19, 2018 21 ``````#' numeric matrix with \eqn{n} rows (samples) `````` Armin Rauschenberger committed Dec 13, 2018 22 ``````#' and \eqn{p} columns (variables) `````` Armin Rauschenberger committed Dec 12, 2018 23 ``````#' `````` Armin Rauschenberger committed Jan 17, 2019 24 25 26 27 ``````#' @param alpha #' elastic net mixing parameter\strong{:} #' numeric between \eqn{0} (ridge) and \eqn{1} (lasso) #' `````` Armin Rauschenberger committed Dec 19, 2018 28 29 30 31 ``````#' @param foldid #' fold identifiers\strong{:} #' vector with entries between \eqn{1} and \code{nfolds}; #' or \code{NULL} (balance) `````` Armin Rauschenberger committed Dec 12, 2018 32 33 34 35 ``````#' #' @param nfolds #' number of folds #' `````` Armin Rauschenberger committed Dec 13, 2018 36 ``````#' @param type.measure `````` Armin Rauschenberger committed Dec 19, 2018 37 38 39 ``````#' loss function for binary classification #' (linear regression uses the deviance) #' `````` Armin Rauschenberger committed Jan 08, 2019 40 41 42 43 44 45 46 47 ``````#' @param pi #' pi sequence\strong{:} #' vector of values in the unit interval; #' or \code{NULL} (default sequence) #' #' @param npi #' number of \code{pi} values #' `````` Armin Rauschenberger committed Dec 19, 2018 48 49 50 51 52 53 54 55 ``````#' @param sigma #' sigma sequence\strong{:} #' vector of increasing positive values; #' or \code{NULL} (default sequence) #' #' @param nsigma #' number of \code{sigma} values #' `````` Armin Rauschenberger committed Dec 20, 2018 56 57 58 ``````#' @param ... #' further arguments passed to \code{\link[glmnet]{glmnet}} #' `````` Armin Rauschenberger committed Dec 19, 2018 59 ``````#' @details `````` Armin Rauschenberger committed Dec 20, 2018 60 61 62 63 ``````#' - INCLUDE note on deviance (not comparable between lin and log models) #' - alpha: elastic net parameter\strong{:} #' numeric between \eqn{0} (ridge) and \eqn{1} (lasso) #' - do not use "family" `````` Armin Rauschenberger committed Dec 14, 2018 64 ``````#' `````` Armin Rauschenberger committed Dec 12, 2018 65 ``````#' @examples `````` Armin Rauschenberger committed Dec 19, 2018 66 67 68 ``````#' n <- 100; p <- 200 #' y <- rnorm(n) #' X <- matrix(rnorm(n*p),nrow=n,ncol=p) `````` Armin Rauschenberger committed Jan 17, 2019 69 ``````#' net <- cornet(y=y,cutoff=0,X=X) `````` Armin Rauschenberger committed Dec 21, 2018 70 ``````#' ### Add ... to all glmnet::glmnet calls !!! ### `````` Armin Rauschenberger committed Jan 17, 2019 71 ``````cornet <- function(y,cutoff,X,alpha=1,npi=101,pi=NULL,nsigma=99,sigma=NULL,nfolds=10,foldid=NULL,type.measure="deviance",...){ `````` Armin Rauschenberger committed Dec 19, 2018 72 `````` `````` Armin Rauschenberger committed Dec 20, 2018 73 `````` #--- temporary --- `````` Armin Rauschenberger committed Jan 08, 2019 74 `````` # cutoff <- 0; npi <- 101; pi <- NULL; nsigma <- 99; sigma <- NULL; nfolds <- 10; foldid <- NULL; type.measure <- "deviance"; logistic <- TRUE `````` Armin Rauschenberger committed Dec 20, 2018 75 `````` test <- list() `````` Armin Rauschenberger committed Jan 08, 2019 76 `````` test\$sigma <- test\$pi <- FALSE `````` Armin Rauschenberger committed Dec 20, 2018 77 `````` test\$grid <- TRUE `````` Armin Rauschenberger committed Jan 08, 2019 78 `````` `````` Armin Rauschenberger committed Dec 20, 2018 79 `````` #--- checks --- `````` Armin Rauschenberger committed Jan 17, 2019 80 `````` cornet:::.check(x=y,type="vector") `````` Armin Rauschenberger committed Dec 21, 2018 81 `````` if(all(y %in% c(0,1))){warning("Binary response.",call.=FALSE)} `````` Armin Rauschenberger committed Jan 17, 2019 82 83 `````` cornet:::.check(x=cutoff,type="scalar",min=min(y),max=max(y)) cornet:::.check(x=X,type="matrix") `````` Armin Rauschenberger committed Jan 17, 2019 84 `````` cornet:::.check(x=alpha,type="scalar",min=0,max=1) `````` Armin Rauschenberger committed Dec 19, 2018 85 `````` if(length(y)!=nrow(X)){stop("Contradictory sample size.",call.=FALSE)} `````` Armin Rauschenberger committed Jan 17, 2019 86 87 88 89 90 91 92 `````` cornet:::.check(x=npi,type="scalar",min=1) cornet:::.check(x=pi,type="vector",min=0,max=1,null=TRUE) cornet:::.check(x=nsigma,type="scalar",min=1) cornet:::.check(x=sigma,type="vector",min=.Machine\$double.eps,null=TRUE) cornet:::.check(x=nfolds,type="scalar",min=3) cornet:::.check(x=foldid,type="vector",values=seq_len(nfolds),null=TRUE) cornet:::.check(x=type.measure,type="string",values=c("deviance","class","mse","mae")) # not auc (min/max confusion) `````` Armin Rauschenberger committed Dec 20, 2018 93 94 `````` if(!is.null(list(...)\$family)){stop("Reserved argument \"family\".",call.=FALSE)} n <- length(y) `````` Armin Rauschenberger committed Dec 12, 2018 95 `````` `````` Armin Rauschenberger committed Dec 13, 2018 96 `````` # binarisation `````` Armin Rauschenberger committed Dec 12, 2018 97 98 `````` z <- 1*(y > cutoff) `````` Armin Rauschenberger committed Dec 13, 2018 99 `````` # fold identifiers `````` Armin Rauschenberger committed Dec 19, 2018 100 101 102 103 `````` if(is.null(foldid)){ foldid <- palasso:::.folds(y=z,nfolds=nfolds) } `````` Armin Rauschenberger committed Dec 20, 2018 104 `````` #--- model fitting --- `````` Armin Rauschenberger committed Dec 12, 2018 105 `````` fit <- list() `````` Armin Rauschenberger committed Jan 17, 2019 106 107 `````` fit\$gaussian <- glmnet::glmnet(y=y,x=X,family="gaussian",alpha=alpha) fit\$binomial <- glmnet::glmnet(y=z,x=X,family="binomial",alpha=alpha) `````` Armin Rauschenberger committed Jan 08, 2019 108 109 `````` #--- sigma sequence --- `````` Armin Rauschenberger committed Dec 19, 2018 110 111 `````` if(is.null(sigma)){ fit\$sigma <- exp(seq(from=log(0.05*stats::sd(y)), `````` Armin Rauschenberger committed Dec 20, 2018 112 `````` to=log(10*stats::sd(y)),length.out=nsigma)) `````` Armin Rauschenberger committed Dec 19, 2018 113 114 `````` } else { fit\$sigma <- sigma `````` Armin Rauschenberger committed Dec 21, 2018 115 `````` nsigma <- length(sigma) `````` Armin Rauschenberger committed Dec 19, 2018 116 `````` } `````` Armin Rauschenberger committed Dec 20, 2018 117 `````` `````` Armin Rauschenberger committed Jan 08, 2019 118 119 `````` #--- pi sequence --- if(is.null(pi)){ `````` Armin Rauschenberger committed Jan 07, 2019 120 `````` fit\$pi <- seq(from=0,to=1,length.out=npi) `````` Armin Rauschenberger committed Jan 08, 2019 121 122 123 `````` } else { fit\$pi <- pi npi <- length(pi) `````` Armin Rauschenberger committed Dec 20, 2018 124 125 `````` } `````` Armin Rauschenberger committed Jan 08, 2019 126 127 128 129 130 `````` #--- tuning parameters --- lab.pi <- paste0("pi",seq_len(npi)) lab.sigma <- paste0("si",seq_len(nsigma)) names(fit\$sigma) <- lab.sigma names(fit\$pi) <- lab.pi `````` Armin Rauschenberger committed Dec 21, 2018 131 `````` `````` Armin Rauschenberger committed Dec 20, 2018 132 `````` #--- cross-validation --- `````` Armin Rauschenberger committed Dec 13, 2018 133 `````` pred <- list() `````` Armin Rauschenberger committed Jan 08, 2019 134 135 136 `````` pred\$y <- matrix(data=NA,nrow=n,ncol=length(fit\$gaussian\$lambda)) pred\$z <- matrix(data=NA,nrow=n,ncol=length(fit\$binomial\$lambda)) `````` Armin Rauschenberger committed Dec 20, 2018 137 138 139 140 `````` if(test\$sigma){ pred\$sigma <- matrix(data=NA,nrow=n,ncol=nsigma) } if(test\$pi){ `````` Armin Rauschenberger committed Jan 08, 2019 141 `````` pred\$pi <- matrix(data=NA,nrow=n,ncol=npi) `````` Armin Rauschenberger committed Dec 20, 2018 142 143 `````` } if(test\$grid){ `````` Armin Rauschenberger committed Dec 21, 2018 144 `````` dimnames <- list(NULL,lab.sigma,lab.pi) `````` Armin Rauschenberger committed Jan 08, 2019 145 `````` pred\$grid <- array(data=NA,dim=c(n,nsigma,npi),dimnames=dimnames) `````` Armin Rauschenberger committed Dec 21, 2018 146 `````` } `````` Armin Rauschenberger committed Jan 08, 2019 147 `````` `````` Armin Rauschenberger committed Dec 20, 2018 148 `````` for(k in seq_len(nfolds)){ `````` Armin Rauschenberger committed Dec 13, 2018 149 `````` `````` Armin Rauschenberger committed Dec 12, 2018 150 151 152 153 154 155 156 `````` y0 <- y[foldid!=k] y1 <- y[foldid==k] z0 <- z[foldid!=k] z1 <- z[foldid==k] X0 <- X[foldid!=k,,drop=FALSE] X1 <- X[foldid==k,,drop=FALSE] `````` Armin Rauschenberger committed Dec 13, 2018 157 `````` # linear regression `````` Armin Rauschenberger committed Jan 17, 2019 158 `````` net <- glmnet::glmnet(y=y0,x=X0,family="gaussian",alpha=alpha) `````` Armin Rauschenberger committed Dec 20, 2018 159 160 `````` temp_y <- stats::predict(object=net,newx=X1,type="response",s=fit\$gaussian\$lambda) pred\$y[foldid==k,seq_len(ncol(temp_y))] <- temp_y `````` Armin Rauschenberger committed Jan 17, 2019 161 `````` cvm <- cornet:::.loss(y=y1,fit=temp_y,family="gaussian",type.measure="deviance")[[1]] `````` Armin Rauschenberger committed Jan 08, 2019 162 `````` y_hat <- temp_y[,which.min(cvm)] `````` Armin Rauschenberger committed Dec 12, 2018 163 `````` `````` Armin Rauschenberger committed Dec 13, 2018 164 `````` # logistic regression `````` Armin Rauschenberger committed Jan 17, 2019 165 `````` net <- glmnet::glmnet(y=z0,x=X0,family="binomial",alpha=alpha) `````` Armin Rauschenberger committed Jan 08, 2019 166 167 `````` temp_z <- stats::predict(object=net,newx=X1,type="response",s=fit\$binomial\$lambda) pred\$z[foldid==k,seq_len(ncol(temp_z))] <- temp_z `````` Armin Rauschenberger committed Jan 17, 2019 168 `````` cvm <- cornet:::.loss(y=z1,fit=temp_z,family="binomial",type.measure=type.measure)[[1]] `````` Armin Rauschenberger committed Jan 08, 2019 169 170 `````` z_hat <- temp_z[,which.min(cvm)] `````` Armin Rauschenberger committed Dec 19, 2018 171 `````` # fusion (sigma) `````` Armin Rauschenberger committed Dec 20, 2018 172 173 `````` if(test\$sigma){ for(i in seq_along(fit\$sigma)){ `````` Armin Rauschenberger committed Jan 08, 2019 174 `````` #pred\$sigma[foldid==k,i] <- stats::pnorm(q=y_hat,mean=cutoff,sd=fit\$sigma[i]) `````` Armin Rauschenberger committed Dec 20, 2018 175 `````` } `````` Armin Rauschenberger committed Dec 19, 2018 176 `````` } `````` Armin Rauschenberger committed Dec 12, 2018 177 `````` `````` Armin Rauschenberger committed Dec 14, 2018 178 `````` # fusion (pi) `````` Armin Rauschenberger committed Dec 20, 2018 179 180 `````` if(test\$pi){ for(i in seq_along(fit\$pi)){ `````` Armin Rauschenberger committed Jan 08, 2019 181 182 `````` #cont <- stats::pnorm(q=y_hat,mean=cutoff,sd=stats::sd(y)) #pred\$pi[foldid==k,i] <- fit\$pi[i]*cont + (1-fit\$pi[i])*z_hat `````` Armin Rauschenberger committed Dec 20, 2018 183 184 `````` } } `````` Armin Rauschenberger committed Jan 08, 2019 185 186 187 `````` # fusion (grid) if(test\$grid){ `````` Armin Rauschenberger committed Dec 21, 2018 188 189 190 `````` for(i in seq_along(fit\$sigma)){ for(j in seq_along(fit\$pi)){ cont <- stats::pnorm(q=y_hat,mean=cutoff,sd=fit\$sigma[i]) `````` Armin Rauschenberger committed Jan 08, 2019 191 `````` pred\$grid[foldid==k,i,j] <- fit\$pi[j]*cont + (1-fit\$pi[j])*z_hat `````` Armin Rauschenberger committed Dec 21, 2018 192 193 194 `````` } } } `````` Armin Rauschenberger committed Dec 12, 2018 195 196 `````` } `````` Armin Rauschenberger committed Dec 20, 2018 197 198 `````` #--- evaluation --- `````` Armin Rauschenberger committed Dec 13, 2018 199 `````` # deviance (not comparable between Gaussian and binomial families) `````` Armin Rauschenberger committed Jan 17, 2019 200 `````` fit\$gaussian\$cvm <- cornet:::.loss(y=y,fit=pred\$y,family="gaussian",type.measure="deviance")[[1]] `````` Armin Rauschenberger committed Dec 12, 2018 201 `````` fit\$gaussian\$lambda.min <- fit\$gaussian\$lambda[which.min(fit\$gaussian\$cvm)] `````` Armin Rauschenberger committed Dec 19, 2018 202 `````` `````` Armin Rauschenberger committed Jan 17, 2019 203 `````` fit\$binomial\$cvm <- cornet:::.loss(y=z,fit=pred\$z,family="binomial",type.measure=type.measure)[[1]] `````` Armin Rauschenberger committed Jan 08, 2019 204 `````` fit\$binomial\$lambda.min <- fit\$binomial\$lambda[which.min(fit\$binomial\$cvm)] `````` Armin Rauschenberger committed Dec 19, 2018 205 `````` `````` Armin Rauschenberger committed Dec 20, 2018 206 `````` if(test\$sigma){ `````` Armin Rauschenberger committed Jan 17, 2019 207 `````` #fit\$sigma.cvm <- cornet:::.loss(y=z,fit=pred\$sigma,family="binomial",type.measure=type.measure)[[1]] `````` Armin Rauschenberger committed Jan 08, 2019 208 `````` #fit\$sigma.min1 <- fit\$sigma[which.min(fit\$sigma.cvm)] `````` Armin Rauschenberger committed Dec 20, 2018 209 210 211 `````` } if(test\$pi){ `````` Armin Rauschenberger committed Jan 17, 2019 212 `````` #fit\$pi.cvm <- cornet:::.loss(y=z,fit=pred\$pi,family="binomial",type.measure=type.measure)[[1]] # trial `````` Armin Rauschenberger committed Jan 08, 2019 213 `````` #fit\$pi.min1 <- fit\$pi[which.min(fit\$pi.cvm)] `````` Armin Rauschenberger committed Dec 20, 2018 214 `````` } `````` Armin Rauschenberger committed Dec 19, 2018 215 `````` `````` Armin Rauschenberger committed Dec 20, 2018 216 `````` if(test\$grid){ `````` Armin Rauschenberger committed Dec 21, 2018 217 `````` dimnames <- list(lab.sigma,lab.pi) `````` Armin Rauschenberger committed Jan 08, 2019 218 `````` fit\$cvm <- matrix(data=NA,nrow=nsigma,ncol=npi,dimnames=dimnames) `````` Armin Rauschenberger committed Dec 21, 2018 219 `````` for(i in seq_len(nsigma)){ `````` Armin Rauschenberger committed Jan 08, 2019 220 `````` for(j in seq_len(npi)){ `````` Armin Rauschenberger committed Jan 17, 2019 221 `````` fit\$cvm[i,j] <- cornet:::.loss(y=z,fit=pred\$grid[,i,j],family="binomial",type.measure=type.measure)[[1]] `````` Armin Rauschenberger committed Dec 21, 2018 222 223 `````` } } `````` Armin Rauschenberger committed Jan 08, 2019 224 `````` temp <- which(fit\$cvm==min(fit\$cvm),arr.ind=TRUE,useNames=TRUE) `````` Armin Rauschenberger committed Dec 21, 2018 225 `````` if(nrow(temp)>1){warning("MULTIPLE!",call.=FALSE);temp <- temp[1,,drop=FALSE]} `````` Armin Rauschenberger committed Jan 08, 2019 226 227 228 `````` fit\$sigma.min <- fit\$sigma[temp[1]] fit\$pi.min <- fit\$pi[temp[2]] if(fit\$cvm[names(fit\$sigma.min),names(fit\$pi.min)]!=min(fit\$cvm)){stop("Internal error.")} `````` Armin Rauschenberger committed Jan 08, 2019 229 230 `````` } `````` Armin Rauschenberger committed Dec 20, 2018 231 `````` #--- return --- `````` Armin Rauschenberger committed Dec 13, 2018 232 `````` fit\$cutoff <- cutoff `````` Armin Rauschenberger committed Dec 19, 2018 233 234 `````` fit\$info <- list(type.measure=type.measure, sd.y=stats::sd(y), `````` Armin Rauschenberger committed Dec 20, 2018 235 `````` table=table(z), `````` Armin Rauschenberger committed Jan 08, 2019 236 `````` test=as.data.frame(test)) `````` Armin Rauschenberger committed Dec 13, 2018 237 `````` `````` Armin Rauschenberger committed Jan 17, 2019 238 `````` class(fit) <- "cornet" `````` Armin Rauschenberger committed Dec 12, 2018 239 240 241 `````` return(fit) } `````` Armin Rauschenberger committed Jan 08, 2019 242 ``````#--- Methods ------------------------------------------------------------------- `````` Armin Rauschenberger committed Dec 20, 2018 243 `````` `````` Armin Rauschenberger committed Jan 08, 2019 244 245 ``````#' @export #' @title `````` Armin Rauschenberger committed Feb 05, 2019 246 ``````#' Extract estimated coefficients `````` Armin Rauschenberger committed Jan 08, 2019 247 248 ``````#' #' @description `````` Armin Rauschenberger committed Feb 05, 2019 249 250 ``````#' Extracts estimated coefficients from linear and logistic regression, #' under the penalty parameter that minimises the cross-validated loss. `````` Armin Rauschenberger committed Jan 08, 2019 251 252 ``````#' #' @param object `````` Armin Rauschenberger committed Feb 05, 2019 253 ``````#' \link[cornet]{cornet} object `````` Armin Rauschenberger committed Jan 08, 2019 254 255 ``````#' #' @param ... `````` Armin Rauschenberger committed Feb 05, 2019 256 257 258 259 260 261 262 ``````#' 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 committed Jan 08, 2019 263 264 265 266 ``````#' #' @examples #' NA #' `````` Armin Rauschenberger committed Jan 17, 2019 267 ``````coef.cornet <- function(object,...){ `````` Armin Rauschenberger committed Jan 08, 2019 268 `````` `````` Armin Rauschenberger committed Jan 08, 2019 269 `````` if(length(list(...))!=0){warning("Ignoring arguments.")} `````` Armin Rauschenberger committed Dec 12, 2018 270 `````` `````` Armin Rauschenberger committed Jan 08, 2019 271 272 273 274 275 `````` 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 committed Dec 13, 2018 276 277 278 279 280 281 `````` coef <- cbind(beta,gamma) colnames(coef) <- c("beta","gamma") return(coef) } `````` Armin Rauschenberger committed Jan 08, 2019 282 283 284 `````` #' @export #' @title `````` Armin Rauschenberger committed Feb 05, 2019 285 ``````#' Plot loss matrix `````` Armin Rauschenberger committed Jan 08, 2019 286 287 ``````#' #' @description `````` Armin Rauschenberger committed Feb 05, 2019 288 289 ``````#' Plots the loss for difference combinations of #' the weight (pi) and scale (sigma) paramters. `````` Armin Rauschenberger committed Jan 08, 2019 290 291 ``````#' #' @param x `````` Armin Rauschenberger committed Feb 05, 2019 292 ``````#' \link[cornet]{cornet} object `````` Armin Rauschenberger committed Jan 08, 2019 293 294 ``````#' #' @param ... `````` Armin Rauschenberger committed Feb 05, 2019 295 ``````#' further arguments (not applicable) `````` Armin Rauschenberger committed Jan 08, 2019 296 297 298 299 ``````#' #' @examples #' NA #' `````` Armin Rauschenberger committed Jan 17, 2019 300 ``````plot.cornet <- function(x,...){ `````` Armin Rauschenberger committed Jan 08, 2019 301 `````` `````` Armin Rauschenberger committed Jan 08, 2019 302 `````` if(length(list(...))!=0){warning("Ignoring arguments.")} `````` Armin Rauschenberger committed Dec 21, 2018 303 `````` `````` Armin Rauschenberger committed Dec 20, 2018 304 `````` k <- 100 `````` Armin Rauschenberger committed Jan 08, 2019 305 `````` levels <- stats::quantile(x\$cvm,probs=seq(from=0,to=1,length.out=k+1)) `````` Armin Rauschenberger committed Feb 05, 2019 306 307 308 309 310 311 312 313 314 `````` ## RColorBrewer 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 committed Dec 20, 2018 315 `````` nsigma <- length(x\$sigma) `````` Armin Rauschenberger committed Dec 21, 2018 316 `````` npi <- length(x\$pi) `````` Armin Rauschenberger committed Dec 20, 2018 317 318 319 `````` graphics::plot.new() graphics::par(xaxs="i",yaxs="i") `````` Armin Rauschenberger committed Dec 21, 2018 320 `````` graphics::plot.window(xlim=c(1-0.5,nsigma+0.5),ylim=c(1-0.5,npi+0.5)) `````` Armin Rauschenberger committed Dec 20, 2018 321 `````` `````` Armin Rauschenberger committed Feb 05, 2019 322 `````` graphics::title(xlab=expression(sigma),ylab=expression(pi),cex.lab=1) `````` Armin Rauschenberger committed Dec 21, 2018 323 `````` #graphics::.filled.contour(x=seq_along(x\$sigma),y=seq_along(x\$pi),z=x\$cvm,levels=levels,col=col) `````` Armin Rauschenberger committed Jan 08, 2019 324 `````` graphics::image(x=seq_along(x\$sigma),y=seq_along(x\$pi),z=x\$cvm,breaks=levels,col=col,add=TRUE) `````` Armin Rauschenberger committed Dec 20, 2018 325 `````` graphics::box() `````` Armin Rauschenberger committed Dec 21, 2018 326 `````` `````` Armin Rauschenberger committed Feb 04, 2019 327 328 329 330 `````` 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 committed Feb 05, 2019 331 `````` # axes with labels for tuned parameters `````` Armin Rauschenberger committed Feb 04, 2019 332 333 `````` 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 committed Feb 05, 2019 334 `````` # point for tuned parameters `````` Armin Rauschenberger committed Feb 05, 2019 335 `````` graphics::points(x=ssigma,y=spi,pch=4,col="white",cex=1) `````` Armin Rauschenberger committed Feb 04, 2019 336 `````` } else { `````` Armin Rauschenberger committed Feb 05, 2019 337 `````` # axes with standard labels `````` Armin Rauschenberger committed Feb 04, 2019 338 339 340 341 `````` 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 committed Feb 05, 2019 342 343 344 `````` # 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 committed Feb 05, 2019 345 `````` graphics::points(x=isigma,y=ipi,pch=4,col="white",cex=1) `````` Armin Rauschenberger committed Feb 04, 2019 346 347 `````` } `````` Armin Rauschenberger committed Dec 19, 2018 348 349 ``````} `````` Armin Rauschenberger committed Jan 08, 2019 350 351 ``````#' @export #' @title `````` Armin Rauschenberger committed Feb 05, 2019 352 ``````#' Predict binary outcome `````` Armin Rauschenberger committed Jan 08, 2019 353 354 ``````#' #' @description `````` Armin Rauschenberger committed Feb 05, 2019 355 356 357 358 359 360 361 362 ``````#' 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 committed Jan 08, 2019 363 ``````#' @param object `````` Armin Rauschenberger committed Feb 05, 2019 364 ``````#' \link[cornet]{cornet} object `````` Armin Rauschenberger committed Jan 08, 2019 365 366 367 368 369 370 371 372 373 374 ``````#' #' @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 committed Feb 05, 2019 375 ``````#' further arguments (not applicable) `````` Armin Rauschenberger committed Jan 08, 2019 376 377 378 379 ``````#' #' @examples #' NA #' `````` Armin Rauschenberger committed Jan 17, 2019 380 ``````predict.cornet <- function(object,newx,type="probability",...){ `````` Armin Rauschenberger committed Jan 08, 2019 381 `````` `````` Armin Rauschenberger committed Jan 08, 2019 382 `````` if(length(list(...))!=0){warning("Ignoring arguments.")} `````` Armin Rauschenberger committed Jan 08, 2019 383 384 `````` x <- object; rm(object) `````` Armin Rauschenberger committed Dec 14, 2018 385 `````` `````` Armin Rauschenberger committed Dec 20, 2018 386 387 `````` test <- x\$info\$test `````` Armin Rauschenberger committed Dec 19, 2018 388 389 `````` .check(x=newx,type="matrix") .check(x=type,type="string",values=c("probability","odds","log-odds")) `````` Armin Rauschenberger committed Dec 13, 2018 390 `````` `````` Armin Rauschenberger committed Dec 19, 2018 391 392 393 394 395 396 397 `````` # linear, logistic and mixed prob <- list() link <- as.numeric(stats::predict(object=x\$gaussian, newx=newx,s=x\$gaussian\$lambda.min,type="response")) prob\$gaussian <- stats::pnorm(q=link,mean=x\$cutoff,sd=x\$info\$sd.y) prob\$binomial <- as.numeric(stats::predict(object=x\$binomial, newx=newx,s=x\$binomial\$lambda.min,type="response")) `````` Armin Rauschenberger committed Dec 20, 2018 398 399 `````` if(test\$sigma){ `````` Armin Rauschenberger committed Jan 08, 2019 400 `````` #prob\$sigma <- stats::pnorm(q=link,mean=x\$cutoff,sd=x\$sigma.min) `````` Armin Rauschenberger committed Dec 20, 2018 401 402 403 `````` } if(test\$pi){ `````` Armin Rauschenberger committed Jan 08, 2019 404 `````` #prob\$pi <- x\$pi.min*prob\$gaussian + (1-x\$pi.min)*prob\$binomial `````` Armin Rauschenberger committed Dec 20, 2018 405 406 407 `````` } if(test\$grid){ `````` Armin Rauschenberger committed Jan 08, 2019 408 `````` cont <- stats::pnorm(q=link,mean=x\$cutoff,sd=x\$sigma.min) `````` Armin Rauschenberger committed Jan 08, 2019 409 `````` prob\$grid <- x\$pi.min*cont + (1-x\$pi.min)*prob\$binomial `````` Armin Rauschenberger committed Dec 21, 2018 410 411 `````` } `````` Armin Rauschenberger committed Dec 19, 2018 412 413 `````` # consistency tests lapply(X=prob,FUN=function(p) .check(x=p,type="vector",min=0,max=1)) `````` Armin Rauschenberger committed Dec 20, 2018 414 `````` .equal(link>x\$cutoff,prob\$gaussian>0.5) `````` Armin Rauschenberger committed Dec 19, 2018 415 416 417 418 419 420 421 422 423 424 `````` # 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 committed Dec 13, 2018 425 426 `````` } `````` Armin Rauschenberger committed Dec 19, 2018 427 `````` return(as.data.frame(frame)) `````` Armin Rauschenberger committed Dec 13, 2018 428 429 ``````} `````` Armin Rauschenberger committed Feb 04, 2019 430 ``````#--- Application --------------------------------------------------------------- `````` Armin Rauschenberger committed Jan 08, 2019 431 `````` `````` Armin Rauschenberger committed Dec 19, 2018 432 433 434 435 436 437 438 ``````#' @export #' @title #' Comparison #' #' @description #' Compares models for a continuous response with a cutoff value #' `````` Armin Rauschenberger committed Jan 17, 2019 439 ``````#' @inheritParams cornet `````` Armin Rauschenberger committed Jan 21, 2019 440 ``````#' `````` Armin Rauschenberger committed Dec 19, 2018 441 442 443 ``````#' @examples #' NA #' `````` Armin Rauschenberger committed Feb 04, 2019 444 ``````.compare <- function(y,cutoff,X,alpha=1,nfolds=5,foldid=NULL,type.measure="deviance"){ `````` Armin Rauschenberger committed Dec 12, 2018 445 `````` `````` Armin Rauschenberger committed Dec 13, 2018 446 `````` z <- 1*(y > cutoff) `````` Armin Rauschenberger committed Jan 08, 2019 447 448 449 450 451 `````` if(is.null(foldid)){ fold <- palasso:::.folds(y=z,nfolds=nfolds) } else { fold <- foldid } `````` Armin Rauschenberger committed Dec 13, 2018 452 `````` `````` Armin Rauschenberger committed Jan 08, 2019 453 `````` cols <- c("gaussian","binomial","grid") `````` Armin Rauschenberger committed Dec 14, 2018 454 455 `````` pred <- matrix(data=NA,nrow=length(y),ncol=length(cols), dimnames=list(NULL,cols)) `````` Armin Rauschenberger committed Dec 12, 2018 456 457 `````` select <- list() `````` Armin Rauschenberger committed Dec 20, 2018 458 `````` for(i in seq_len(nfolds)){ `````` Armin Rauschenberger committed Feb 04, 2019 459 `````` fit <- cornet::cornet(y=y[fold!=i],cutoff=cutoff,X=X[fold!=i,],alpha=alpha,type.measure=type.measure) `````` Armin Rauschenberger committed Jan 17, 2019 460 461 `````` tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL) temp <- cornet:::predict.cornet(fit,newx=X[fold==i,]) `````` Armin Rauschenberger committed Dec 21, 2018 462 `````` if(any(temp<0|temp>1)){stop("Outside unit interval.",call.=FALSE)} `````` Armin Rauschenberger committed Dec 14, 2018 463 464 465 466 `````` model <- colnames(pred) for(j in seq_along(model)){ pred[fold==i,model[j]] <- temp[[model[j]]] } `````` Armin Rauschenberger committed Dec 12, 2018 467 468 `````` } `````` Armin Rauschenberger committed Dec 14, 2018 469 `````` type <- c("deviance","class","mse","mae","auc") `````` Armin Rauschenberger committed Jan 17, 2019 470 `````` loss <- lapply(X=type,FUN=function(x) cornet:::.loss(y=z,fit=pred,family="binomial",type.measure=x,foldid=fold)[[1]]) `````` Armin Rauschenberger committed Dec 14, 2018 471 `````` names(loss) <- type `````` Armin Rauschenberger committed Jan 21, 2019 472 `````` `````` Armin Rauschenberger committed Feb 04, 2019 473 474 475 476 477 478 479 480 481 482 483 `````` #if(trial){ # list <- list(diff=(pred-z)^2,fold=fold,loss=loss) # return(list) #} else { # return(loss) #} return(loss) } `````` Armin Rauschenberger committed Dec 19, 2018 484 `````` `````` Armin Rauschenberger committed Feb 04, 2019 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 ``````#' @export #' @title #' Testing #' #' @description #' Compares models for a continuous response with a cutoff value (testing) #' #' @inheritParams cornet #' #' @examples #' NA #' .test <- function(y,cutoff,X,alpha=1,type.measure="deviance"){ z <- 1*(y > cutoff) fold <- palasso:::.folds(y=z,nfolds=5) fold <- ifelse(fold==1,1,0) fit <- cornet::cornet(y=y[fold==0],cutoff=cutoff,X=X[fold==0,],alpha=alpha) tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL) pred <- cornet:::predict.cornet(fit,newx=X[fold==1,]) if(any(pred<0|pred>1)){stop("Outside unit interval.",call.=FALSE)} #res <- (pred-z[fold==1])^2 # MSE #pvalue <- wilcox.test(x=res[,"binomial"],y=res[,"grid"],paired=TRUE,alternative="greater")\$p.value #colMeans(abs(pred-0.5)) # distance from 0.5 `````` Armin Rauschenberger committed Jan 21, 2019 511 `````` `````` Armin Rauschenberger committed Feb 04, 2019 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 `````` limit <- 1e-05 pred[pred < limit] <- limit pred[pred > 1 - limit] <- 1 - limit res <- -2 * (y[fold==1] * log(pred) + (1 - y[fold==1]) * log(1 - pred)) pvalue <- stats::wilcox.test(x=res[,"binomial"],y=res[,"grid"],paired=TRUE,alternative="greater")\$p.value return(pvalue) } # Simulates y and X. .simulate <- function(n,p,prob=0.2,fac=1){ beta <- stats::rnorm(n=p) cond <- stats::rbinom(n=p,size=1,prob=prob) beta[cond==0] <- 0 X <- matrix(stats::rnorm(n=n*p),nrow=n,ncol=p) mean <- X %*% beta y <- stats::rnorm(n=n,mean=mean,sd=fac*stats::sd(mean)) return(list(y=y,X=X)) `````` Armin Rauschenberger committed Dec 12, 2018 530 531 ``````} `````` Armin Rauschenberger committed Feb 04, 2019 532 533 534 535 536 537 538 539 540 541 542 543 544 ``````#--- start trial --- if(FALSE){ n <- 1000 y_hat <- runif(n) y <- y_hat > 0.9 y <- rbinom(n=n,size=1,prob=0.5) foldid <- rep(1:10,length.out=n) .loss(y=y,fit=y_hat,family="binomial",type.measure="auc",foldid=foldid) } #--- end trial --- #--- Internal functions -------------------------------------------------------- `````` Armin Rauschenberger committed Jan 08, 2019 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 ``````#' @title #' Arguments #' #' @description #' Verifies whether two or more arguments are identical. #' #' @param ... #' scalars, vectors, or matrices of equal dimensions #' #' @param na.rm #' remove missing values\strong{:} #' logical #' #' @examples #' NA #' `````` Armin Rauschenberger committed Dec 19, 2018 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 ``````.equal <- function(...,na.rm=FALSE){ list <- list(...) cond <- vapply(X=list, FUN=function(x) all(x==list[[1]],na.rm=na.rm), FUN.VALUE=logical(1)) if(any(!cond)){ stop("Unequal elements.",call.=FALSE) } return(invisible(NULL)) } #' @title #' Arguments #' #' @description #' Verifies whether an argument matches formal requirements. #' #' @param x #' argument #' #' @param type #' character \code{"string"}, \code{"scalar"}, \code{"vector"}, \code{"matrix"} #' #' @param miss #' accept missing values\strong{:} #' logical #' #' @param min #' lower limit\strong{:} #' numeric #' #' @param max #' upper limit\strong{:} #' numeric #' #' @param values #' only accept specific values\strong{:} #' vector #' #' @param inf #' accept infinite (\code{Inf} or \code{-Inf}) values\strong{:} #' logical #' #' @param null #' accept \code{NULL}\strong{:} #' logical #' #' @examples #' NA #' .check <- function(x,type,miss=FALSE,min=NULL,max=NULL,values=NULL,inf=FALSE,null=FALSE){ `````` Armin Rauschenberger committed Dec 13, 2018 612 `````` name <- deparse(substitute(x)) `````` Armin Rauschenberger committed Dec 19, 2018 613 614 615 `````` if(null && is.null(x)){ return(invisible(NULL)) } `````` Armin Rauschenberger committed Dec 13, 2018 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 `````` 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),call.=FALSE) } if(!is.null(max) && any(x>max)){ stop(paste0("expecting ",name," <= ",max),call.=FALSE) } if(!is.null(values) && !(x %in% values)){ stop(paste0("Argument \"",name,"\" contains invalid values."),call.=FALSE) } `````` Armin Rauschenberger committed Dec 19, 2018 642 643 644 `````` if(!inf && any(is.infinite(values))){ stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE) } `````` Armin Rauschenberger committed Dec 13, 2018 645 646 647 `````` return(invisible(NULL)) } `````` Armin Rauschenberger committed Jan 18, 2019 648 649 650 `````` `````` Armin Rauschenberger committed Dec 13, 2018 651 652 ``````# Correct this function in the palasso package (search twice for "# typo"). .loss <- function (y,fit,family,type.measure,foldid=NULL){ `````` Armin Rauschenberger committed Dec 12, 2018 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 `````` if (!is.list(fit)) { fit <- list(fit) } loss <- list() for (i in seq_along(fit)) { if (is.vector(fit[[i]])) { fit[[i]] <- as.matrix(fit[[i]]) } if (is.null(foldid) & (family == "cox" | type.measure == "auc")) { stop("Missing foldid.", call. = FALSE) } if (family == "gaussian") { if (type.measure %in% c("deviance", "mse")) { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) mean((x - y)^2)) } else if (type.measure == "mae") { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) mean(abs(x - y))) } else { stop("Invalid type measure.", call. = FALSE) } } else if (family == "binomial") { if (type.measure == "deviance") { limit <- 1e-05 fit[[i]][fit[[i]] < limit] <- limit fit[[i]][fit[[i]] > 1 - limit] <- 1 - limit loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) mean(-2 * (y * log(x) + (1 - y) * log(1 - x)))) } else if (type.measure == "mse") { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) 2 * mean((x - y)^2)) } else if (type.measure == "mae") { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) 2 * mean(abs(x - y))) } else if (type.measure == "class") { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) mean(abs(round(x) - y))) } else if (type.measure == "auc") { weights <- table(foldid) cvraw <- matrix(data = NA, nrow = length(weights), ncol = ncol(fit[[i]])) # typo in palasso package ! for (k in seq_along(weights)) { cvraw[k, ] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) glmnet::auc(y = y[foldid == k], prob = x[foldid == k])) } loss[[i]] <- apply(X = cvraw, MARGIN = 2, FUN = function(x) stats::weighted.mean(x = x, w = weights, na.rm = TRUE)) `````` Armin Rauschenberger committed Dec 13, 2018 710 `````` names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package! `````` Armin Rauschenberger committed Dec 12, 2018 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 `````` } else { stop("Invalid type measure.", call. = FALSE) } } else if (family == "poisson") { if (type.measure == "deviance") { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) mean(2 * (ifelse(y == 0, 0, y * log(y/x)) - y + x), na.rm = TRUE)) } else if (type.measure == "mse") { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) mean((x - y)^2)) } else if (type.measure == "mae") { loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) mean(abs(x - y))) } else { stop("Invalid type measure.", call. = FALSE) } } else if (family == "cox") { if (type.measure == "deviance") { weights <- tapply(X = y[, "status"], INDEX = foldid, FUN = sum) loss[[i]] <- apply(X = fit[[i]], MARGIN = 2, FUN = function(x) stats::weighted.mean(x = x/weights, w = weights, na.rm = TRUE)) } else { stop("Invalid type measure.", call. = FALSE) } } else { stop("Invalid family.", call. = FALSE) } if (sum(diff(is.na(loss[[i]]))) == 1) { loss[[i]] <- loss[[i]][!is.na(loss[[i]])] } } return(loss) } `````` Armin Rauschenberger committed Jan 08, 2019 756 757 ``````#--- Lost and found ------------------------------------------------------------ `````` Armin Rauschenberger committed Jan 17, 2019 758 ``````# calibrate (for cornet) `````` Armin Rauschenberger committed Jan 08, 2019 759 760 761 ``````#if(test\$calibrate){ # fit\$calibrate <- CalibratR::calibrate(actual=z,predicted=pred\$y[,which.min(fit\$gaussian\$cvm)],nCores=1,model_idx=5)\$calibration_models #} `````` Armin Rauschenberger committed Jan 17, 2019 762 ``````# calibrate (for predict.cornet) `````` Armin Rauschenberger committed Jan 08, 2019 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 ``````#if(test\$calibrate){ # prob\$calibrate <- CalibratR::predict_calibratR(calibration_models=x\$calibrate,new=link,nCores=1)\$GUESS_2 #} .args <- function(...){ args <- list(...) names <- names(formals(glmnet::glmnet)) if(!is.null(args\$family)){ warning("Unexpected argument \"family\".",call.=FALSE) } if(any(!names(args) %in% names)){ stop("Unexpected argument.",call.=FALSE) } if(is.null(args\$alpha)) { args\$alpha <- 1 } if(is.null(args\$nlambda)){ args\$nlambda <- 100 } if(is.null(args\$lambda)){ if(is.null(args\$nlambda)){ args\$nlambda <- 100 } } else { args\$nlambda <- length(args\$lambda) } return(args) } `````` Armin Rauschenberger committed Dec 13, 2018 795