functions.R 20.4 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 246 247 248 249 250 251 ``````#' @export #' @title #' to do #' #' @description #' to do #' #' @param object `````` Armin Rauschenberger committed Jan 17, 2019 252 ``````#' cornet object `````` Armin Rauschenberger committed Jan 08, 2019 253 254 255 256 257 258 259 ``````#' #' @param ... #' to do #' #' @examples #' NA #' `````` Armin Rauschenberger committed Jan 17, 2019 260 ``````coef.cornet <- function(object,...){ `````` Armin Rauschenberger committed Jan 08, 2019 261 `````` `````` Armin Rauschenberger committed Jan 08, 2019 262 `````` if(length(list(...))!=0){warning("Ignoring arguments.")} `````` Armin Rauschenberger committed Dec 12, 2018 263 `````` `````` Armin Rauschenberger committed Jan 08, 2019 264 265 266 267 268 `````` 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 269 270 271 272 273 274 `````` coef <- cbind(beta,gamma) colnames(coef) <- c("beta","gamma") return(coef) } `````` Armin Rauschenberger committed Jan 08, 2019 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 `````` #' @export #' @title #' to do #' #' @description #' to do #' #' @param x #' to do #' #' @param ... #' further arguments #' #' @examples #' NA #' `````` Armin Rauschenberger committed Jan 17, 2019 292 ``````plot.cornet <- function(x,...){ `````` Armin Rauschenberger committed Jan 08, 2019 293 `````` `````` Armin Rauschenberger committed Jan 08, 2019 294 `````` if(length(list(...))!=0){warning("Ignoring arguments.")} `````` Armin Rauschenberger committed Dec 21, 2018 295 `````` `````` Armin Rauschenberger committed Dec 20, 2018 296 `````` k <- 100 `````` Armin Rauschenberger committed Jan 08, 2019 297 `````` levels <- stats::quantile(x\$cvm,probs=seq(from=0,to=1,length.out=k+1)) `````` Armin Rauschenberger committed Dec 20, 2018 298 299 `````` col <- colorspace::diverge_hsv(n=k) nsigma <- length(x\$sigma) `````` Armin Rauschenberger committed Dec 21, 2018 300 `````` npi <- length(x\$pi) `````` Armin Rauschenberger committed Dec 20, 2018 301 302 303 `````` graphics::plot.new() graphics::par(xaxs="i",yaxs="i") `````` Armin Rauschenberger committed Dec 21, 2018 304 `````` 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 305 `````` `````` Armin Rauschenberger committed Jan 18, 2019 306 307 `````` ssigma <- which(x\$sigma==x\$sigma.min) graphics::axis(side=1,at=c(1,ssigma,nsigma),labels=signif(x\$sigma[c(1,ssigma,nsigma)],digits=2)) `````` Armin Rauschenberger committed Dec 20, 2018 308 `````` `````` Armin Rauschenberger committed Jan 18, 2019 309 310 `````` spi <- which(x\$pi==x\$pi.min) graphics::axis(side=2,at=c(1,spi,npi),labels=signif(x\$pi[c(1,spi,npi)],digits=2)) `````` Armin Rauschenberger committed Dec 20, 2018 311 `````` `````` Armin Rauschenberger committed Dec 21, 2018 312 313 `````` graphics::title(xlab=expression(sigma),ylab=expression(pi)) #graphics::.filled.contour(x=seq_along(x\$sigma),y=seq_along(x\$pi),z=x\$cvm,levels=levels,col=col) `````` Armin Rauschenberger committed Jan 08, 2019 314 `````` 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 315 `````` graphics::box() `````` Armin Rauschenberger committed Dec 21, 2018 316 `````` `````` Armin Rauschenberger committed Jan 18, 2019 317 318 319 320 321 `````` #graphics::abline(v=ssigma,lty=2,col="grey") #graphics::abline(h=spi,lty=2,col="grey") graphics::points(x=ssigma,y=spi,pch=4,col="black",cex=1) `````` Armin Rauschenberger committed Dec 19, 2018 322 323 ``````} `````` Armin Rauschenberger committed Jan 08, 2019 324 325 326 327 328 329 330 331 ``````#' @export #' @title #' to do #' #' @description #' to do #' #' @param object `````` Armin Rauschenberger committed Jan 17, 2019 332 ``````#' cornet object `````` Armin Rauschenberger committed Jan 08, 2019 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 ``````#' #' @param newx #' covariates\strong{:} #' numeric matrix with \eqn{n} rows (samples) #' and \eqn{p} columns (variables) #' #' @param type #' \code{"probability"}, \code{"odds"}, \code{"log-odds"} #' #' @param ... #' to do #' #' @examples #' NA #' `````` Armin Rauschenberger committed Jan 17, 2019 348 ``````predict.cornet <- function(object,newx,type="probability",...){ `````` Armin Rauschenberger committed Jan 08, 2019 349 `````` `````` Armin Rauschenberger committed Jan 08, 2019 350 `````` if(length(list(...))!=0){warning("Ignoring arguments.")} `````` Armin Rauschenberger committed Jan 08, 2019 351 352 `````` x <- object; rm(object) `````` Armin Rauschenberger committed Dec 14, 2018 353 `````` `````` Armin Rauschenberger committed Dec 20, 2018 354 355 `````` test <- x\$info\$test `````` Armin Rauschenberger committed Dec 19, 2018 356 357 `````` .check(x=newx,type="matrix") .check(x=type,type="string",values=c("probability","odds","log-odds")) `````` Armin Rauschenberger committed Dec 13, 2018 358 `````` `````` Armin Rauschenberger committed Dec 19, 2018 359 360 361 362 363 364 365 `````` # 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 366 367 `````` if(test\$sigma){ `````` Armin Rauschenberger committed Jan 08, 2019 368 `````` #prob\$sigma <- stats::pnorm(q=link,mean=x\$cutoff,sd=x\$sigma.min) `````` Armin Rauschenberger committed Dec 20, 2018 369 370 371 `````` } if(test\$pi){ `````` Armin Rauschenberger committed Jan 08, 2019 372 `````` #prob\$pi <- x\$pi.min*prob\$gaussian + (1-x\$pi.min)*prob\$binomial `````` Armin Rauschenberger committed Dec 20, 2018 373 374 375 `````` } if(test\$grid){ `````` Armin Rauschenberger committed Jan 08, 2019 376 `````` cont <- stats::pnorm(q=link,mean=x\$cutoff,sd=x\$sigma.min) `````` Armin Rauschenberger committed Jan 08, 2019 377 `````` prob\$grid <- x\$pi.min*cont + (1-x\$pi.min)*prob\$binomial `````` Armin Rauschenberger committed Dec 21, 2018 378 379 `````` } `````` Armin Rauschenberger committed Dec 19, 2018 380 381 `````` # consistency tests lapply(X=prob,FUN=function(p) .check(x=p,type="vector",min=0,max=1)) `````` Armin Rauschenberger committed Dec 20, 2018 382 `````` .equal(link>x\$cutoff,prob\$gaussian>0.5) `````` Armin Rauschenberger committed Dec 19, 2018 383 384 385 386 387 388 389 390 391 392 `````` # 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 393 394 `````` } `````` Armin Rauschenberger committed Dec 19, 2018 395 `````` return(as.data.frame(frame)) `````` Armin Rauschenberger committed Dec 13, 2018 396 397 ``````} `````` Armin Rauschenberger committed Jan 08, 2019 398 399 ``````#--- Internal functions -------------------------------------------------------- `````` Armin Rauschenberger committed Dec 19, 2018 400 401 402 403 404 405 406 ``````#' @export #' @title #' Comparison #' #' @description #' Compares models for a continuous response with a cutoff value #' `````` Armin Rauschenberger committed Jan 17, 2019 407 ``````#' @inheritParams cornet `````` Armin Rauschenberger committed Dec 19, 2018 408 409 410 411 ``````#' #' @examples #' NA #' `````` Armin Rauschenberger committed Jan 17, 2019 412 ``````.compare <- function(y,cutoff,X,alpha=1,nfolds=5,foldid=NULL,type.measure="deviance"){ `````` Armin Rauschenberger committed Dec 12, 2018 413 `````` `````` Armin Rauschenberger committed Dec 13, 2018 414 `````` z <- 1*(y > cutoff) `````` Armin Rauschenberger committed Jan 08, 2019 415 416 417 418 419 `````` if(is.null(foldid)){ fold <- palasso:::.folds(y=z,nfolds=nfolds) } else { fold <- foldid } `````` Armin Rauschenberger committed Dec 13, 2018 420 `````` `````` Armin Rauschenberger committed Jan 08, 2019 421 422 `````` #cols <- c("gaussian","binomial","sigma","pi","trial","grid","unit") cols <- c("gaussian","binomial","grid") `````` Armin Rauschenberger committed Dec 14, 2018 423 424 `````` pred <- matrix(data=NA,nrow=length(y),ncol=length(cols), dimnames=list(NULL,cols)) `````` Armin Rauschenberger committed Dec 12, 2018 425 426 `````` select <- list() `````` Armin Rauschenberger committed Dec 20, 2018 427 `````` for(i in seq_len(nfolds)){ `````` Armin Rauschenberger committed Jan 17, 2019 428 `````` fit <- cornet::cornet(y=y[fold!=i],cutoff=cutoff,X=X[fold!=i,],alpha=alpha,logistic=TRUE,type.measure=type.measure) `````` Armin Rauschenberger committed Jan 17, 2019 429 430 431 `````` tryCatch(expr=cornet:::plot.cornet(fit),error=function(x) NULL) #cornet:::plot.cornet(fit) temp <- cornet:::predict.cornet(fit,newx=X[fold==i,]) `````` Armin Rauschenberger committed Dec 21, 2018 432 `````` if(any(temp<0|temp>1)){stop("Outside unit interval.",call.=FALSE)} `````` Armin Rauschenberger committed Dec 14, 2018 433 434 435 436 `````` model <- colnames(pred) for(j in seq_along(model)){ pred[fold==i,model[j]] <- temp[[model[j]]] } `````` Armin Rauschenberger committed Dec 12, 2018 437 438 `````` } `````` Armin Rauschenberger committed Dec 14, 2018 439 `````` type <- c("deviance","class","mse","mae","auc") `````` Armin Rauschenberger committed Jan 17, 2019 440 `````` 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 441 `````` names(loss) <- type `````` Armin Rauschenberger committed Dec 19, 2018 442 `````` `````` Armin Rauschenberger committed Dec 12, 2018 443 444 445 `````` return(loss) } `````` Armin Rauschenberger committed Jan 08, 2019 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 ``````#' @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 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 ``````.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 513 `````` name <- deparse(substitute(x)) `````` Armin Rauschenberger committed Dec 19, 2018 514 515 516 `````` if(null && is.null(x)){ return(invisible(NULL)) } `````` Armin Rauschenberger committed Dec 13, 2018 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 `````` 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 543 544 545 `````` if(!inf && any(is.infinite(values))){ stop(paste0("Argument \"",name,"\ contains infinite values."),call.=FALSE) } `````` Armin Rauschenberger committed Dec 13, 2018 546 547 548 `````` return(invisible(NULL)) } `````` Armin Rauschenberger committed Jan 08, 2019 549 550 551 552 553 554 555 556 557 558 ``````# 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 13, 2018 559 `````` `````` Armin Rauschenberger committed Jan 18, 2019 560 561 562 563 564 565 566 567 568 569 570 571 572 ``````#--- 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 --- `````` Armin Rauschenberger committed Dec 13, 2018 573 574 ``````# 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 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 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 `````` 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 632 `````` names(loss[[i]]) <- colnames(fit[[i]]) # typo in palasso package! `````` Armin Rauschenberger committed Dec 12, 2018 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 `````` } 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 678 679 ``````#--- Lost and found ------------------------------------------------------------ `````` Armin Rauschenberger committed Jan 17, 2019 680 ``````# calibrate (for cornet) `````` Armin Rauschenberger committed Jan 08, 2019 681 682 683 ``````#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 684 ``````# calibrate (for predict.cornet) `````` Armin Rauschenberger committed Jan 08, 2019 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 710 711 712 713 714 715 716 ``````#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 717