Nothing
inst/doc/script.R
In joinet: Multivariate Elastic Net Regression
## ----start,include=FALSE------------------------------------------------------
knitr::opts_chunk$set(echo=TRUE,eval=FALSE)
#setwd("~/Desktop/joinet")
#source("~/Desktop/joinet/package/R/functions.R")
#pkgs <- c("devtools","missRanger","mice","mvtnorm","glmnet","earth","spls","MTPS","RMTL","MultivariateRandomForest","mcen","MRCE","remMap","SiER","GPM","MLPUGS","benchmarkme")
#install.packages(pkgs)
#devtools::install_github("rauschenberger/joinet")
#devtools::install_github("cran/QUIC")
#devtools::install_github("cran/MRCE")
#devtools::install_github("cran/remMap")
#install.packages("joinet_0.0.X.tar.gz",repos=NULL,type="source")
library(joinet)
## ----figure_ABS---------------------------------------------------------------
# box <- function(x,y,w=0.2,h=0.2,labels="",col="black",turn=FALSE,...){
# xs <- x + 0.5*c(-1,-1,1,1)*w
# ys <- y + 0.5*c(-1,1,1,-1)*h
# graphics::polygon(x=xs,y=ys,border=col,col=grey(ifelse(col=="grey",0.99,0.99)),lwd=2,...)
# graphics::text(x=x,y=y,labels=labels,col=col,srt=ifelse(turn,90,0))
# }
#
# #grDevices::pdf(file="manuscript/figure_ABS.pdf",width=4,height=2.5)
# grDevices::setEPS()
# grDevices::postscript(file="manuscript/figure_ABS.eps",width=4,height=2.5)
# graphics::par(mar=c(0,0,0,0))
# graphics::plot.new()
# graphics::plot.window(xlim=c(0,2),ylim=c(-0.1,1.1))
#
# n <- 0.45; p <- 0.4; q <- 0.4
# h <- w <- 0.4
#
# y0 <- 0.10; y1 <- 0.50; y2 <- 0.90
# x0 <- 0.20; x1 <- 1.00; x2 <- 1.80
# y <- seq(from=1,to=0,length.out=5)
#
# w <- 0.15
# d <- seq(from=0.06,to=-0.06,length.out=4)[c(1,2,3,Inf,4)]
# for(i in 1:5){
# if(i==4){graphics::text(x=c(x1,x2),y=y[4],labels="...",srt=90,font=2,cex=1.2);next}
# graphics::arrows(x0=x0+p/2,x1=x1-q/2-0.02,y0=y1,y1=y[i],lwd=2,length=0.1,col="grey")
# labels <- paste("target",ifelse(i==5,"q",i))
# for(j in 1:5){
# if(j==4){next}
# graphics::arrows(x0=x1+p/2,x1=x2-q/2-0.02,y0=y[i],y1=y[j]+d[i],lwd=2,length=0.1,col=ifelse(i==j,"grey","grey"))
# }
# box(x=x1,y=y[i],h=0.5*h,w=q,labels=labels)
# box(x=x2,y=y[i],h=0.5*h,w=q,labels=labels)
# }
# box(x=x0,y=y1,h=0.5*h,w=p,labels="features")
#
# grDevices::dev.off()
## ----figure_OUT---------------------------------------------------------------
# #grDevices::pdf(file="manuscript/figure_OUT.pdf",width=5,height=3)
# grDevices::postscript(file="manuscript/figure_OUT.eps",width=5,height=3)
#
# ellipse <- function(x,y,text=NULL,a=0.5,b=0.5,border="black",col=NULL,txt.col="black",...){
# n <- max(c(length(x),length(y)))
# if(is.null(col)){col <- rep(grey(0.9),times=n)}
# if(length(col)==1){col <- rep(col,times=n)}
# if(length(x)==1){x <- rep(x,times=n)}
# if(length(y)==1){y <- rep(y,times=n)}
# if(length(text)==1){text <- rep(text,times=n)}
# if(length(border)==1){border <- rep(border,times=n)}
# for(i in seq_len(n)){
# angle <- seq(from=0,to=2*pi,length=100)
# xs <- x[i] + a * cos(angle)
# ys <- y[i] + b * sin(angle)
# graphics::polygon(x=xs,y=ys,col=col[i],border=border[i])
# graphics::text(labels=text[i],x=x[i],y=y[i],col=txt.col,...)
# }
# }
#
# txt <- list()
# txt$x <- expression(x[j])
# txt$y <- eval(parse(text=paste0("expression(",paste0("y[",c(1:3,"k","q"),"]",collapse=","),")")))
# txt$beta <- eval(parse(text=paste0("expression(",paste0("beta[j",c(1:3,"k","q"),"]",collapse=","),")")))
# txt$omega <- eval(parse(text=paste0("expression(",paste0("omega[\"",c(1:3,"k","q"),"k\"]",collapse=","),")")))
# txt$dots <- expression(cdots)
#
# pos <- list()
# pos$x <- 4
# pos$y <- c(1,2,3,5,7)
# pos$beta <- median(pos$x)+(pos$y-median(pos$x))/2
# pos$omega <- 5+(pos$y-5)/2
#
# a <- b <- 0.3
# graphics::plot.new()
# graphics::par(mar=c(0,0,0,0))
# graphics::plot.window(xlim=c(0.5,7.5),ylim=c(0.5,5.5))
#
# # beta
# segments(x0=4,y0=5-a,x1=pos$y,y1=3+a,lwd=2,col="blue")
# ellipse(x=pos$beta,y=4,text=txt$beta,a=0.21,b=0.21,cex=1.2,col="white",border="white",txt.col="blue")
#
# # omega
# segments(x0=rep(pos$y,each=4),y0=3-a,x1=rep(pos$y,times=4),y1=1+a,lwd=2,col="grey")
# segments(x0=pos$y,y0=3-a,y1=1+a,x1=5,lwd=2,col="red")
# ellipse(x=pos$omega,y=2,text=txt$omega,a=0.25,b=0.18,cex=1.2,col="white",border="white",txt.col="red")
#
# # x and y
# ellipse(x=pos$x,y=5,text=txt$x,a=a,b=b,cex=1.2)
# text(x=c(3,5),y=5,labels=txt$dots,cex=1.2)
# ellipse(x=pos$y,y=3,text=txt$y,a=a,b=b,cex=1.2)
# text(x=c(4,6),y=3,labels=txt$dots,cex=1.2)
# ellipse(x=pos$y,y=1,text=txt$y,a=a,b=b,cex=1.2)
# text(x=c(4,6),y=1,labels=txt$dots,cex=1.2)
#
# grDevices::dev.off()
## ----simulation,eval=FALSE----------------------------------------------------
# #<<start>>
# library(joinet)
#
# grid <- list()
# grid$rho_x <- c(0.00,0.10,0.30)
# grid$rho_b <- c(0.00,0.50,0.90)
# delta <- 0.8
# grid <- expand.grid(grid)
# grid <- rbind(grid,grid,grid)
# grid$p <- rep(c(10,500,500),each=nrow(grid)/3)
# grid$nzero <- rep(c(5,5,100),each=nrow(grid)/3)
# grid <- grid[rep(1:nrow(grid),times=10),]
#
# n0 <- 100; n1 <- 10000
# n <- n0 + n1
# q <- 3
# foldid.ext <- rep(c(0,1),times=c(n0,n1))
#
# loss <- list(); cor <- numeric(); seed <- list()
# set.seed(1) # new
# for(i in seq_len(nrow(grid))){
# p <- grid$p[i]
# #set.seed(i) # old
# seed[[i]] <- .Random.seed # new
# cat("i =",i,"\n")
#
# #--- features ---
# mean <- rep(0,times=p)
# sigma <- matrix(grid$rho_x[i],nrow=p,ncol=p)
# diag(sigma) <- 1
# X <- mvtnorm::rmvnorm(n=n,mean=mean,sigma=sigma)
#
# #--- effects --- (multivariate Gaussian)
# mean <- rep(0,times=q)
# sigma <- matrix(data=grid$rho_b[i],nrow=q,ncol=q)
# diag(sigma) <- 1
# beta <- mvtnorm::rmvnorm(n=p,mean=mean,sigma=sigma)
# #beta <- 1*apply(beta,2,function(x) x>(sort(x,decreasing=TRUE)[grid$nzero[i]])) # old (either zero or one)
# beta <- 1*apply(beta,2,function(x) ifelse(x>sort(x,decreasing=TRUE)[grid$nzero[i]],x,0)) # new (either zero or non-zero)
#
# #-- effects --- (multivariate binomial)
# #sigma <- matrix(grid$rho_b[i],nrow=q,ncol=q); diag(sigma) <- 1
# #beta <- bindata::rmvbin(n=p,margprob=rep(grid$prop[i],times=q),bincorr=sigma)
#
# #--- outcomes ---
# signal <- scale(X%*%beta)
# signal[is.na(signal)] <- 0
# noise <- matrix(rnorm(n*q),nrow=n,ncol=q)
# Y <- sqrt(delta)*signal + sqrt(1-delta)*noise
# # binomial: Y <- round(exp(Y)/(1+exp(Y)))
# cors <- stats::cor(Y)
# diag(cors) <- NA
# cor[i] <- mean(cors,na.rm=TRUE)
#
# #--- holdout ---
# alpha.base <- 1*(grid$nzero[i] <= 10) # sparse vs dense
# compare <- TRUE
# loss[[i]] <- tryCatch(expr=cv.joinet(Y=Y,X=X,family="gaussian",compare=compare,foldid.ext=foldid.ext,alpha.base=alpha.base,alpha.meta=1,times=TRUE,sign=1),error=function(x) NA)
# }
# save(grid,loss,cor,seed,file="results/simulation.RData")
# writeLines(text=capture.output(utils::sessionInfo(),cat("\n"),sessioninfo::session_info()),con="results/info_sim.txt")
## ----figure_SIM,results="hide"------------------------------------------------
# #<<start>>
# load("results/simulation.RData")
#
# cond <- lapply(loss,length)==2
# if(any(!cond)){
# warning("At least one error.",call.=FALSE)
# grid <- grid[cond,]; loss <- loss[cond]
# }
#
# #--- computation time ---
# time <- sapply(loss,function(x) unlist(x$time))
# #round(sort(colMeans(apply(time,1,function(x) x/time["meta",]))),digits=1)
# sort(round(rowMeans(time),digits=1))
#
# #--- average ---
# loss <- lapply(loss,function(x) x$loss)
#
# #prop <- sapply(loss[cond],function(x) rowMeans(100*x/matrix(x["none",],nrow=nrow(x),ncol=ncol(x),byrow=TRUE))[-nrow(x)]) # old (first re-scale, then average)
# mean <- sapply(loss,function(x) rowMeans(x)) # new (first average)
# prop <- 100*mean[rownames(mean)!="none",]/matrix(mean["none",],nrow=nrow(mean)-1,ncol=ncol(mean),byrow=TRUE) # new (then re-scale)
#
# mode <- ifelse(grid$p==10,"ld",ifelse(grid$nzero==5,"hd-s","hd-d"))
# set <- as.numeric(sapply(rownames(grid),function(x) strsplit(x,split="\\.")[[1]][[1]]))
#
# #--- mean rank ---
# mult <- rownames(prop)!="base"
# sort(round(rowMeans(apply(prop[mult,mode=="ld"],2,rank)),digits=1))
# sort(round(rowMeans(apply(prop[mult,mode=="hd-s"],2,rank)),digits=1))
# sort(round(rowMeans(apply(prop[mult,mode=="hd-d"],2,rank)),digits=1))
#
# #--- testing ---
# apply(prop[mult,],1,function(x) sum(tapply(X=x-prop["base",],INDEX=set,FUN=function(x) wilcox.test(x,alternative="less",exact=FALSE)$p.value<0.05),na.rm=TRUE))
#
# colSums(tapply(X=prop["meta",]-prop["base",],INDEX=list(set=set,mode=mode),FUN=function(x) wilcox.test(x,alternative="less",exact=FALSE)$p.value<0.05),na.rm=TRUE)
# colSums(tapply(X=prop["spls",]-prop["base",],INDEX=list(set=set,mode=mode),FUN=function(x) wilcox.test(x,alternative="less",exact=FALSE)$p.value<0.05),na.rm=TRUE)
## ----table_SIM,results="asis"-------------------------------------------------
# beta <- sapply(unique(set),function(i) rowMeans(prop[,set==i]))
# cor <- sapply(unique(set),function(i) mean(cor[set==i]))
#
# rownames(beta)[rownames(beta)=="mnorm"] <- "mvn"
# sign <- apply(beta,2,function(x) sign(x["base"]-x))
# #min <- apply(beta,2,function(x) which.min(x)) # incorrect (old)
# min <- apply(beta,2,function(x) which(x==min(x))) # correct (new)
# beta <- format(round(beta,digits=1),trim=TRUE)
# beta[sign<=0] <- paste0("\\textcolor{gray}{",beta[sign<=0],"}")
# index <- cbind(min,1:ncol(beta))
# beta[index] <- paste0("\\underline{",beta[index],"}")
# unique <- unique(grid)
# info <- format(round(cbind("$\\rho_x$"=unique$rho_x,"$\\rho_b$"=unique$rho_b,"$\\rho_y$"=cor),digits=1))
# temp <- paste0("\\cran{",sapply(rownames(beta),function(x) switch(x,base="glmnet",meta="joinet",mvn="glmnet",mars="earth",spls="spls",mrce="MRCE",map="remMap",mrf="MultivariateRandomForest",sier="SiER",mcen="mcen",gpm="GPM",rmtl="RMTL",mtps="MTPS",NULL)),"}")
# temp[1] <- paste0(temp[1],"$^1$")
# temp[3] <- paste0(temp[3],"$^2$")
# temp[8] <- "\\href{https://CRAN.R-project.org/package=MultivariateRandomForest}{\\texttt{MRF}}$^3$"
# rownames(beta) <- paste0("\\begin{sideways}",temp,"\\end{sideways}")
# xtable <- xtable::xtable(cbind(info,t(beta)),align=paste0("rccc",paste0(rep("c",times=nrow(beta)),collapse=""),collapse=""),caption="")
# xtable::print.xtable(xtable,comment=FALSE,floating=TRUE,sanitize.text.function=identity,hline.after=c(0,9,18,ncol(beta)),include.rownames=FALSE,size="\\footnotesize",caption.placement="top")
## ----data,eval=FALSE----------------------------------------------------------
# # clinical features
# X <- read.csv("data/PPMI_Baseline_Data_02Jul2018.csv",row.names="PATNO",na.strings=c(".",""))
# colnames(X) <- tolower(colnames(X))
# X <- X[X$apprdx==1,] # Parkinson's disease
# X[c("site","apprdx","event_id","symptom5_comment")] <- NULL
# for(i in seq_len(ncol(X))){
# if(is.factor(X[[i]])){levels(X[[i]]) <- paste0("-",levels(X[[i]]))}
# }
# 100*mean(is.na(X)) # proportion missingness
# x <- lapply(seq_len(1),function(x) missRanger::missRanger(data=X,pmm.k=3,
# num.trees=100,verbose=0,seed=1))
# x <- lapply(x,function(x) model.matrix(~.-1,data=x))
#
# # genomic features
# load("data/ppmi_rnaseq_bl_pd_hc-2019-01-11.Rdata",verbose=TRUE)
# counts <- t(ppmi_rnaseq_bl_pdhc)
# mean(grepl(pattern="ENSG0000|ENSGR0000",x=colnames(counts)))
# cond <- apply(counts,2,function(x) sd(x)>0) & !grepl(pattern="ENSG0000|ENSGR0000",x=colnames(counts))
# Z <- palasso:::.prepare(counts[,cond],filter=10,cutoff="knee")$X
#
# # outcome
# Y <- read.csv("data/PPMI_Year_1-3_Data_02Jul2018.csv",na.strings=".")
# Y <- Y[Y$APPRDX==1 & Y$EVENT_ID %in% c("V04","V06","V08"),]
# colnames(Y)[colnames(Y)=="updrs_totscore"] <- "updrs"
# vars <- c("moca","quip","updrs","gds","scopa","ess","bjlot","rem")
# # too few levels: "NP1HALL","NP1DPRS"
# Y <- Y[,c("EVENT_ID","PATNO",vars)]
# Y <- reshape(Y,idvar="PATNO",timevar="EVENT_ID",direction="wide")
# rownames(Y) <- Y$PATNO; Y$PATNO <- NULL
#
# # overlap
# names <- Reduce(intersect,list(rownames(X),rownames(Y),rownames(Z)))
# Z <- Z[names,]
# Y <- Y[names,]
# Y <- sapply(vars,function(x) Y[,grepl(pattern=x,x=colnames(Y))],simplify=FALSE)
# for(i in seq_along(Y)){
# colnames(Y[[i]]) <- c("V04","V06","V08")
# }
# x <- lapply(x,function(x) x[names,]); rm(names)
# X <- x[[1]]; rm(x) # impute multiple times!
#
# # inversion for positive correlation
# Y$moca <- -Y$moca # "wrong" sign
# Y$bjlot <- -Y$bjlot # "wrong" sign
# sapply(Y,function(x) range(unlist(x),na.rm=TRUE))
#
# save(Y,X,Z,file="results/data.RData")
# writeLines(text=capture.output(utils::sessionInfo(),cat("\n"),sessioninfo::session_info()),con="results/info_dat.txt")
## ----figure_COR---------------------------------------------------------------
# load("results/data.RData",verbose=TRUE)
#
# #grDevices::pdf(file="manuscript/figure_COR.pdf",width=6,height=3)
# grDevices::postscript(file="manuscript/figure_COR.eps",width=6,height=3)
# graphics::par(mar=c(0.5,3,2,0.5))
# graphics::layout(mat=matrix(c(1,2),nrow=1,ncol=2),width=c(0.2,0.8))
#
# # correlation between years
# cor <- cbind(sapply(Y,function(x) cor(x[,1],x[,2],use="complete.obs",method="spearman")),
# sapply(Y,function(x) cor(x[,2],x[,3],use="complete.obs",method="spearman")),
# sapply(Y,function(x) cor(x[,1],x[,3],use="complete.obs",method="spearman")))
# colnames(cor) <- c("1-2","2-3","1-3")
# cor <- rowMeans(cor)
# joinet:::plot.matrix(cor,range=c(-3,3),margin=1,cex=0.7)
#
# # correlation between variables
# cor <- 1/3*cor(sapply(Y,function(x) x[,1]),use="complete.obs",method="spearman")+
# 1/3*cor(sapply(Y,function(x) x[,2]),use="complete.obs",method="spearman")+
# 1/3*cor(sapply(Y,function(x) x[,3]),use="complete.obs",method="spearman")
# joinet:::plot.matrix(cor,range=c(-3,3),margin=c(1,2),cex=0.7)
# grDevices::dev.off()
#
# # other information
# sapply(Y,colMeans,na.rm=TRUE) # increasing values
# sapply(Y,function(x) apply(x,2,sd,na.rm=TRUE)) # increasing variance
# sapply(Y,function(x) colSums(is.na(x))) # increasing numbers of NAs
## ----application,eval=FALSE---------------------------------------------------
# #<<start>>
# library(joinet)
#
# set.seed(1)
# load("results/data.RData",verbose=TRUE)
#
# set.seed(1)
# foldid.ext <- rep(1:5,length.out=nrow(Y$moca))
# foldid.int <- rep(rep(1:10,each=5),length.out=nrow(Y$moca))
# table(foldid.ext,foldid.int)
#
# #- - - - - - - - - - - - -
# #- - internal coaching - -
# #- - - - - - - - - - - - -
#
# table <- list()
# table$alpha <- c("lasso","ridge")
# table$data <- c("clinic","omics","both")
# table$var <- names(Y)
# table <- rev(expand.grid(table,stringsAsFactors=FALSE))
#
# loss <- fit <- list()
# for(i in seq_len(nrow(table))){
# cat(rep("*",times=5),"setting",i,rep("*",times=5),"\n")
# y <- Y[[table$var[i]]]
# x <- list(clinic=X,omics=Z,both=cbind(X,Z))[[table$data[i]]]
# alpha <- 1*(table$alpha[i]=="lasso")
# loss[[i]] <- cv.joinet(Y=y,X=x,alpha.base=alpha,foldid.ext=foldid.ext,
# foldid.int=foldid.int,sign=1) # add joinet::
# #fit[[i]] <- joinet(Y=y,X=x,alpha.base=alpha,foldid=foldid.int,sign=1)
# }
#
# save(table,loss,file="results/internal.RData")
#
# #- - - - - - - - - - - - -
# #- - external coaching - -
# #- - - - - - - - - - - - -
#
# table <- list()
# temp <- utils::combn(x=names(Y),m=2)
# table$comb <- paste0(temp[1,],"-",temp[2,])
# table$step <- c("V04","V06","V08")
# table$alpha <- c("lasso","ridge")
# table$data <- c("clinic","omics","both")
# table <- rev(expand.grid(table,stringsAsFactors=FALSE))
# temp <- strsplit(table$comb,split="-"); table$comb <- NULL
# table$var1 <- sapply(temp,function(x) x[[1]])
# table$var2 <- sapply(temp,function(x) x[[2]])
#
# loss <- fit <- list()
# for(i in seq_len(nrow(table))){
# cat(rep("*",times=5),"setting",i,rep("*",times=5),"\n")
# y <- cbind(Y[[table$var1[i]]][,table$step[i]],
# Y[[table$var2[i]]][,table$step[i]])
# x <- list(clinic=X,omics=Z,both=cbind(X,Z))[[table$data[i]]]
# alpha <- 1*(table$alpha[i]=="lasso")
# loss[[i]] <- cv.joinet(Y=y,X=x,alpha.base=alpha,
# foldid.ext=foldid.ext,foldid.int=foldid.int,sign=1) # add joinet::
# #fit[[i]] <- joinet(Y=y,X=x,alpha.base=alpha,foldid=foldid.int,sign=1)
# }
#
# save(table,loss,file="results/external.RData")
# writeLines(text=capture.output(utils::sessionInfo(),cat("\n"),sessioninfo::session_info()),con="results/info_app.txt")
## ----figure_INT---------------------------------------------------------------
# load("results/internal.RData")
#
# # standardised loss
# vars <- unique(table$var)
# former <- t(sapply(loss,function(x) x["base",]))
# min <- sapply(vars,function(x) min(former[table$var==x,]))
# max <- sapply(vars,function(x) max(former[table$var==x,]))
# index <- match(x=table$var,table=vars)
# former <- (former-min[index])/(max[index]-min[index])
# dimnames(former) <- list(table$var,seq_len(3))
#
# # percentage change
# change <- t(sapply(loss,function(x) 100*(x["meta",]-x["base",])/x["base",]))
# dimnames(change) <- list(table$var,c("1st","2nd","3rd"))
#
# # overview
# #grDevices::pdf(file="manuscript/figure_INT.pdf",width=6,height=3,pointsize=12)
# grDevices::postscript(file="manuscript/figure_INT.eps",width=6,height=3,pointsize=12)
# graphics::par(mfrow=c(2,3),mar=c(0.1,3,2,0.1),oma=c(0,1.1,1,0))
# for(alpha in c("lasso","ridge")){
# for(data in c("clinic","omics","both")){
# cond <- table$alpha==alpha & table$data==data
# joinet:::plot.matrix(X=change[cond,],range=c(-50,50),cex=0.7)
# #graphics::title(main=paste0(alpha,"-",data),col.main="red",line=0) # check
# if(alpha=="lasso"){graphics::mtext(text=data,side=3,line=1.5,cex=0.8)}
# if(data=="clinic"){graphics::mtext(text=alpha,side=2,line=3,cex=0.8)}
# }
# }
# grDevices::dev.off()
#
# TEMP <- tapply(X=rowMeans(change),INDEX=table$var,FUN=mean)[vars]
# mean(change<0)
# round(tapply(X=rowMeans(change),INDEX=table$data,FUN=mean),digits=2)
# round(tapply(X=rowMeans(change),INDEX=table$alpha,FUN=mean),digits=2)
# round(colMeans(change),digits=2)
## ----figure_EXT---------------------------------------------------------------
# load("results/external.RData")
#
# data <- c("clinic","omics","both")
# alpha <- c("lasso","ridge")
# step <- c("V04","V06","V08")
#
# # percentage change
# change <- t(sapply(loss,function(x) 100*(x["meta",]-x["base",])/x["base",]))
#
# # overview
# vars <- unique(c(table$var1,table$var2))
# temp <- matrix(NA,nrow=length(vars),ncol=length(vars),dimnames=list(vars,vars))
# array <- array(data=list(temp),dim=c(3,2,3),dimnames=list(data,alpha,step))
# #grDevices::pdf(file="manuscript/figure_EXT.pdf",width=7.5,height=10,pointsize=14)
# grDevices::postscript(file="manuscript/figure_EXT.eps",width=7.5,height=10,pointsize=14)
# graphics::par(mfrow=c(6,3),mar=c(0.1,2.5,2.5,0.1),oma=c(0,1,2,0))
# for(i in data){
# for(j in alpha){
# for(k in step){
# cond <- table$data==i & table$alpha==j & table$step==k
# array[i,j,k][[1]][cbind(table$var1,table$var2)[cond,]] <- change[cond,1]
# array[i,j,k][[1]][cbind(table$var2,table$var1)[cond,]] <- change[cond,2]
# joinet:::plot.matrix(array[i,j,k][[1]],margin=0,las=2,range=c(-20,20),cex=0.6)
# #graphics::title(main=paste0(i,"-",j,"-",k),col.main="red",line=0) # check
# if(i=="clinic" & j=="lasso"){graphics::mtext(text=ifelse(k=="V04","1st",ifelse(k=="V06","2nd","3rd")),side=3,line=2.5,cex=0.8)}
# if(k=="V04"){graphics::mtext(text=paste0(i,"-",j),side=2,line=2.5,cex=0.8)}
# }
# }
# }
# grDevices::dev.off()
#
# # check
# i <- sample(seq_len(nrow(table)),size=1)
# table[i,]
# x <- loss[[i]]
# 100*(x["meta",]-x["base",])/x["base",]
## ----figure_ALL---------------------------------------------------------------
# #grDevices::pdf(file="manuscript/figure_ALL.pdf",height=3,width=6)
# grDevices::postscript(file="manuscript/figure_ALL.eps",height=3,width=6)
#
# graphics::par(mar=c(0.5,3,2,0.5))
# graphics::layout(mat=matrix(c(1,2),nrow=1,ncol=2),width=c(0.2,0.8))
# joinet:::plot.matrix(as.matrix(TEMP),margin=1,las=1,range=c(-20,20),cex=0.7)
#
# sum(unlist(array)<0,na.rm=TRUE)/sum(!is.na(unlist(array)))
# means <- apply(array,c(1,2,3),function(x) mean(x[[1]],na.rm=TRUE))
# lapply(seq_len(3),function(x) apply(means,x,mean))
# mean <- 1/length(array)*Reduce(f="+",x=array)
# joinet:::plot.matrix(mean,margin=1,las=1,range=c(-20,20),cex=0.7)
# # rows: target variable, columns: coaching variable
#
# grDevices::dev.off()
## ----figure_DIF---------------------------------------------------------------
# #grDevices::pdf(file="manuscript/figure_DIF.pdf",height=1.2,width=5)
# grDevices::postscript(file="manuscript/figure_DIF.eps",height=1.2,width=5)
#
# load("results/internal.RData")
# vars <- unique(table$var)
# base <- t(sapply(loss,function(x) 100*(x["base",]-x["none",])/x["none",]))
# meta <- t(sapply(loss,function(x) 100*(x["meta",]-x["none",])/x["none",]))
# dimnames(meta) <- dimnames(base) <- list(table$var,c("1st","2nd","3rd"))
# standard <- tapply(X=rowMeans(base),INDEX=table$var,FUN=mean)[vars]
# internal <- tapply(X=rowMeans(meta),INDEX=table$var,FUN=mean)[vars]
#
# load("results/external.RData")
# vars <- unique(c(table$var1,table$var2))
# base <- meta <- list()
# for(i in seq_len(2)){
# base[[i]] <- sapply(loss,function(x) 100*(x["base",i]-x["none",i])/x["none",i])
# meta[[i]] <- sapply(loss,function(x) 100*(x["meta",i]-x["none",i])/x["none",i])
# }
# index <- c(table$var1,table$var2)
# base <- unlist(base); meta <- unlist(meta)
# #standard <- tapply(X=base,INDEX=index,FUN=mean)[vars]
# external <- tapply(X=meta,INDEX=index,FUN=mean)[vars]
#
# matrix <- round(rbind(standard,internal,external),digits=2)
# rownames(matrix) <- c("","","")
# graphics::par(mfrow=c(1,1),mar=c(0.5,3,1.5,1))
# joinet:::plot.matrix(matrix,margin=c(1,2),las=1,range=c(-100,0),cex=0.7,digits=3)
#
# grDevices::dev.off()
Try the joinet package in your browser
Any scripts or data that you put into this service are public.
joinet documentation built on Aug. 9, 2021, 9:13 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----start,include=FALSE------------------------------------------------------
knitr::opts_chunk$set(echo=TRUE,eval=FALSE)
#setwd("~/Desktop/joinet")
#source("~/Desktop/joinet/package/R/functions.R")
#pkgs <- c("devtools","missRanger","mice","mvtnorm","glmnet","earth","spls","MTPS","RMTL","MultivariateRandomForest","mcen","MRCE","remMap","SiER","GPM","MLPUGS","benchmarkme")
#install.packages(pkgs)
#devtools::install_github("rauschenberger/joinet")
#devtools::install_github("cran/QUIC")
#devtools::install_github("cran/MRCE")
#devtools::install_github("cran/remMap")
#install.packages("joinet_0.0.X.tar.gz",repos=NULL,type="source")
library(joinet)
## ----figure_ABS---------------------------------------------------------------
# box <- function(x,y,w=0.2,h=0.2,labels="",col="black",turn=FALSE,...){
# xs <- x + 0.5*c(-1,-1,1,1)*w
# ys <- y + 0.5*c(-1,1,1,-1)*h
# graphics::polygon(x=xs,y=ys,border=col,col=grey(ifelse(col=="grey",0.99,0.99)),lwd=2,...)
# graphics::text(x=x,y=y,labels=labels,col=col,srt=ifelse(turn,90,0))
# }
#
# #grDevices::pdf(file="manuscript/figure_ABS.pdf",width=4,height=2.5)
# grDevices::setEPS()
# grDevices::postscript(file="manuscript/figure_ABS.eps",width=4,height=2.5)
# graphics::par(mar=c(0,0,0,0))
# graphics::plot.new()
# graphics::plot.window(xlim=c(0,2),ylim=c(-0.1,1.1))
#
# n <- 0.45; p <- 0.4; q <- 0.4
# h <- w <- 0.4
#
# y0 <- 0.10; y1 <- 0.50; y2 <- 0.90
# x0 <- 0.20; x1 <- 1.00; x2 <- 1.80
# y <- seq(from=1,to=0,length.out=5)
#
# w <- 0.15
# d <- seq(from=0.06,to=-0.06,length.out=4)[c(1,2,3,Inf,4)]
# for(i in 1:5){
# if(i==4){graphics::text(x=c(x1,x2),y=y[4],labels="...",srt=90,font=2,cex=1.2);next}
# graphics::arrows(x0=x0+p/2,x1=x1-q/2-0.02,y0=y1,y1=y[i],lwd=2,length=0.1,col="grey")
# labels <- paste("target",ifelse(i==5,"q",i))
# for(j in 1:5){
# if(j==4){next}
# graphics::arrows(x0=x1+p/2,x1=x2-q/2-0.02,y0=y[i],y1=y[j]+d[i],lwd=2,length=0.1,col=ifelse(i==j,"grey","grey"))
# }
# box(x=x1,y=y[i],h=0.5*h,w=q,labels=labels)
# box(x=x2,y=y[i],h=0.5*h,w=q,labels=labels)
# }
# box(x=x0,y=y1,h=0.5*h,w=p,labels="features")
#
# grDevices::dev.off()
## ----figure_OUT---------------------------------------------------------------
# #grDevices::pdf(file="manuscript/figure_OUT.pdf",width=5,height=3)
# grDevices::postscript(file="manuscript/figure_OUT.eps",width=5,height=3)
#
# ellipse <- function(x,y,text=NULL,a=0.5,b=0.5,border="black",col=NULL,txt.col="black",...){
# n <- max(c(length(x),length(y)))
# if(is.null(col)){col <- rep(grey(0.9),times=n)}
# if(length(col)==1){col <- rep(col,times=n)}
# if(length(x)==1){x <- rep(x,times=n)}
# if(length(y)==1){y <- rep(y,times=n)}
# if(length(text)==1){text <- rep(text,times=n)}
# if(length(border)==1){border <- rep(border,times=n)}
# for(i in seq_len(n)){
# angle <- seq(from=0,to=2*pi,length=100)
# xs <- x[i] + a * cos(angle)
# ys <- y[i] + b * sin(angle)
# graphics::polygon(x=xs,y=ys,col=col[i],border=border[i])
# graphics::text(labels=text[i],x=x[i],y=y[i],col=txt.col,...)
# }
# }
#
# txt <- list()
# txt$x <- expression(x[j])
# txt$y <- eval(parse(text=paste0("expression(",paste0("y[",c(1:3,"k","q"),"]",collapse=","),")")))
# txt$beta <- eval(parse(text=paste0("expression(",paste0("beta[j",c(1:3,"k","q"),"]",collapse=","),")")))
# txt$omega <- eval(parse(text=paste0("expression(",paste0("omega[\"",c(1:3,"k","q"),"k\"]",collapse=","),")")))
# txt$dots <- expression(cdots)
#
# pos <- list()
# pos$x <- 4
# pos$y <- c(1,2,3,5,7)
# pos$beta <- median(pos$x)+(pos$y-median(pos$x))/2
# pos$omega <- 5+(pos$y-5)/2
#
# a <- b <- 0.3
# graphics::plot.new()
# graphics::par(mar=c(0,0,0,0))
# graphics::plot.window(xlim=c(0.5,7.5),ylim=c(0.5,5.5))
#
# # beta
# segments(x0=4,y0=5-a,x1=pos$y,y1=3+a,lwd=2,col="blue")
# ellipse(x=pos$beta,y=4,text=txt$beta,a=0.21,b=0.21,cex=1.2,col="white",border="white",txt.col="blue")
#
# # omega
# segments(x0=rep(pos$y,each=4),y0=3-a,x1=rep(pos$y,times=4),y1=1+a,lwd=2,col="grey")
# segments(x0=pos$y,y0=3-a,y1=1+a,x1=5,lwd=2,col="red")
# ellipse(x=pos$omega,y=2,text=txt$omega,a=0.25,b=0.18,cex=1.2,col="white",border="white",txt.col="red")
#
# # x and y
# ellipse(x=pos$x,y=5,text=txt$x,a=a,b=b,cex=1.2)
# text(x=c(3,5),y=5,labels=txt$dots,cex=1.2)
# ellipse(x=pos$y,y=3,text=txt$y,a=a,b=b,cex=1.2)
# text(x=c(4,6),y=3,labels=txt$dots,cex=1.2)
# ellipse(x=pos$y,y=1,text=txt$y,a=a,b=b,cex=1.2)
# text(x=c(4,6),y=1,labels=txt$dots,cex=1.2)
#
# grDevices::dev.off()
## ----simulation,eval=FALSE----------------------------------------------------
# #<<start>>
# library(joinet)
#
# grid <- list()
# grid$rho_x <- c(0.00,0.10,0.30)
# grid$rho_b <- c(0.00,0.50,0.90)
# delta <- 0.8
# grid <- expand.grid(grid)
# grid <- rbind(grid,grid,grid)
# grid$p <- rep(c(10,500,500),each=nrow(grid)/3)
# grid$nzero <- rep(c(5,5,100),each=nrow(grid)/3)
# grid <- grid[rep(1:nrow(grid),times=10),]
#
# n0 <- 100; n1 <- 10000
# n <- n0 + n1
# q <- 3
# foldid.ext <- rep(c(0,1),times=c(n0,n1))
#
# loss <- list(); cor <- numeric(); seed <- list()
# set.seed(1) # new
# for(i in seq_len(nrow(grid))){
# p <- grid$p[i]
# #set.seed(i) # old
# seed[[i]] <- .Random.seed # new
# cat("i =",i,"\n")
#
# #--- features ---
# mean <- rep(0,times=p)
# sigma <- matrix(grid$rho_x[i],nrow=p,ncol=p)
# diag(sigma) <- 1
# X <- mvtnorm::rmvnorm(n=n,mean=mean,sigma=sigma)
#
# #--- effects --- (multivariate Gaussian)
# mean <- rep(0,times=q)
# sigma <- matrix(data=grid$rho_b[i],nrow=q,ncol=q)
# diag(sigma) <- 1
# beta <- mvtnorm::rmvnorm(n=p,mean=mean,sigma=sigma)
# #beta <- 1*apply(beta,2,function(x) x>(sort(x,decreasing=TRUE)[grid$nzero[i]])) # old (either zero or one)
# beta <- 1*apply(beta,2,function(x) ifelse(x>sort(x,decreasing=TRUE)[grid$nzero[i]],x,0)) # new (either zero or non-zero)
#
# #-- effects --- (multivariate binomial)
# #sigma <- matrix(grid$rho_b[i],nrow=q,ncol=q); diag(sigma) <- 1
# #beta <- bindata::rmvbin(n=p,margprob=rep(grid$prop[i],times=q),bincorr=sigma)
#
# #--- outcomes ---
# signal <- scale(X%*%beta)
# signal[is.na(signal)] <- 0
# noise <- matrix(rnorm(n*q),nrow=n,ncol=q)
# Y <- sqrt(delta)*signal + sqrt(1-delta)*noise
# # binomial: Y <- round(exp(Y)/(1+exp(Y)))
# cors <- stats::cor(Y)
# diag(cors) <- NA
# cor[i] <- mean(cors,na.rm=TRUE)
#
# #--- holdout ---
# alpha.base <- 1*(grid$nzero[i] <= 10) # sparse vs dense
# compare <- TRUE
# loss[[i]] <- tryCatch(expr=cv.joinet(Y=Y,X=X,family="gaussian",compare=compare,foldid.ext=foldid.ext,alpha.base=alpha.base,alpha.meta=1,times=TRUE,sign=1),error=function(x) NA)
# }
# save(grid,loss,cor,seed,file="results/simulation.RData")
# writeLines(text=capture.output(utils::sessionInfo(),cat("\n"),sessioninfo::session_info()),con="results/info_sim.txt")
## ----figure_SIM,results="hide"------------------------------------------------
# #<<start>>
# load("results/simulation.RData")
#
# cond <- lapply(loss,length)==2
# if(any(!cond)){
# warning("At least one error.",call.=FALSE)
# grid <- grid[cond,]; loss <- loss[cond]
# }
#
# #--- computation time ---
# time <- sapply(loss,function(x) unlist(x$time))
# #round(sort(colMeans(apply(time,1,function(x) x/time["meta",]))),digits=1)
# sort(round(rowMeans(time),digits=1))
#
# #--- average ---
# loss <- lapply(loss,function(x) x$loss)
#
# #prop <- sapply(loss[cond],function(x) rowMeans(100*x/matrix(x["none",],nrow=nrow(x),ncol=ncol(x),byrow=TRUE))[-nrow(x)]) # old (first re-scale, then average)
# mean <- sapply(loss,function(x) rowMeans(x)) # new (first average)
# prop <- 100*mean[rownames(mean)!="none",]/matrix(mean["none",],nrow=nrow(mean)-1,ncol=ncol(mean),byrow=TRUE) # new (then re-scale)
#
# mode <- ifelse(grid$p==10,"ld",ifelse(grid$nzero==5,"hd-s","hd-d"))
# set <- as.numeric(sapply(rownames(grid),function(x) strsplit(x,split="\\.")[[1]][[1]]))
#
# #--- mean rank ---
# mult <- rownames(prop)!="base"
# sort(round(rowMeans(apply(prop[mult,mode=="ld"],2,rank)),digits=1))
# sort(round(rowMeans(apply(prop[mult,mode=="hd-s"],2,rank)),digits=1))
# sort(round(rowMeans(apply(prop[mult,mode=="hd-d"],2,rank)),digits=1))
#
# #--- testing ---
# apply(prop[mult,],1,function(x) sum(tapply(X=x-prop["base",],INDEX=set,FUN=function(x) wilcox.test(x,alternative="less",exact=FALSE)$p.value<0.05),na.rm=TRUE))
#
# colSums(tapply(X=prop["meta",]-prop["base",],INDEX=list(set=set,mode=mode),FUN=function(x) wilcox.test(x,alternative="less",exact=FALSE)$p.value<0.05),na.rm=TRUE)
# colSums(tapply(X=prop["spls",]-prop["base",],INDEX=list(set=set,mode=mode),FUN=function(x) wilcox.test(x,alternative="less",exact=FALSE)$p.value<0.05),na.rm=TRUE)
## ----table_SIM,results="asis"-------------------------------------------------
# beta <- sapply(unique(set),function(i) rowMeans(prop[,set==i]))
# cor <- sapply(unique(set),function(i) mean(cor[set==i]))
#
# rownames(beta)[rownames(beta)=="mnorm"] <- "mvn"
# sign <- apply(beta,2,function(x) sign(x["base"]-x))
# #min <- apply(beta,2,function(x) which.min(x)) # incorrect (old)
# min <- apply(beta,2,function(x) which(x==min(x))) # correct (new)
# beta <- format(round(beta,digits=1),trim=TRUE)
# beta[sign<=0] <- paste0("\\textcolor{gray}{",beta[sign<=0],"}")
# index <- cbind(min,1:ncol(beta))
# beta[index] <- paste0("\\underline{",beta[index],"}")
# unique <- unique(grid)
# info <- format(round(cbind("$\\rho_x$"=unique$rho_x,"$\\rho_b$"=unique$rho_b,"$\\rho_y$"=cor),digits=1))
# temp <- paste0("\\cran{",sapply(rownames(beta),function(x) switch(x,base="glmnet",meta="joinet",mvn="glmnet",mars="earth",spls="spls",mrce="MRCE",map="remMap",mrf="MultivariateRandomForest",sier="SiER",mcen="mcen",gpm="GPM",rmtl="RMTL",mtps="MTPS",NULL)),"}")
# temp[1] <- paste0(temp[1],"$^1$")
# temp[3] <- paste0(temp[3],"$^2$")
# temp[8] <- "\\href{https://CRAN.R-project.org/package=MultivariateRandomForest}{\\texttt{MRF}}$^3$"
# rownames(beta) <- paste0("\\begin{sideways}",temp,"\\end{sideways}")
# xtable <- xtable::xtable(cbind(info,t(beta)),align=paste0("rccc",paste0(rep("c",times=nrow(beta)),collapse=""),collapse=""),caption="")
# xtable::print.xtable(xtable,comment=FALSE,floating=TRUE,sanitize.text.function=identity,hline.after=c(0,9,18,ncol(beta)),include.rownames=FALSE,size="\\footnotesize",caption.placement="top")
## ----data,eval=FALSE----------------------------------------------------------
# # clinical features
# X <- read.csv("data/PPMI_Baseline_Data_02Jul2018.csv",row.names="PATNO",na.strings=c(".",""))
# colnames(X) <- tolower(colnames(X))
# X <- X[X$apprdx==1,] # Parkinson's disease
# X[c("site","apprdx","event_id","symptom5_comment")] <- NULL
# for(i in seq_len(ncol(X))){
# if(is.factor(X[[i]])){levels(X[[i]]) <- paste0("-",levels(X[[i]]))}
# }
# 100*mean(is.na(X)) # proportion missingness
# x <- lapply(seq_len(1),function(x) missRanger::missRanger(data=X,pmm.k=3,
# num.trees=100,verbose=0,seed=1))
# x <- lapply(x,function(x) model.matrix(~.-1,data=x))
#
# # genomic features
# load("data/ppmi_rnaseq_bl_pd_hc-2019-01-11.Rdata",verbose=TRUE)
# counts <- t(ppmi_rnaseq_bl_pdhc)
# mean(grepl(pattern="ENSG0000|ENSGR0000",x=colnames(counts)))
# cond <- apply(counts,2,function(x) sd(x)>0) & !grepl(pattern="ENSG0000|ENSGR0000",x=colnames(counts))
# Z <- palasso:::.prepare(counts[,cond],filter=10,cutoff="knee")$X
#
# # outcome
# Y <- read.csv("data/PPMI_Year_1-3_Data_02Jul2018.csv",na.strings=".")
# Y <- Y[Y$APPRDX==1 & Y$EVENT_ID %in% c("V04","V06","V08"),]
# colnames(Y)[colnames(Y)=="updrs_totscore"] <- "updrs"
# vars <- c("moca","quip","updrs","gds","scopa","ess","bjlot","rem")
# # too few levels: "NP1HALL","NP1DPRS"
# Y <- Y[,c("EVENT_ID","PATNO",vars)]
# Y <- reshape(Y,idvar="PATNO",timevar="EVENT_ID",direction="wide")
# rownames(Y) <- Y$PATNO; Y$PATNO <- NULL
#
# # overlap
# names <- Reduce(intersect,list(rownames(X),rownames(Y),rownames(Z)))
# Z <- Z[names,]
# Y <- Y[names,]
# Y <- sapply(vars,function(x) Y[,grepl(pattern=x,x=colnames(Y))],simplify=FALSE)
# for(i in seq_along(Y)){
# colnames(Y[[i]]) <- c("V04","V06","V08")
# }
# x <- lapply(x,function(x) x[names,]); rm(names)
# X <- x[[1]]; rm(x) # impute multiple times!
#
# # inversion for positive correlation
# Y$moca <- -Y$moca # "wrong" sign
# Y$bjlot <- -Y$bjlot # "wrong" sign
# sapply(Y,function(x) range(unlist(x),na.rm=TRUE))
#
# save(Y,X,Z,file="results/data.RData")
# writeLines(text=capture.output(utils::sessionInfo(),cat("\n"),sessioninfo::session_info()),con="results/info_dat.txt")
## ----figure_COR---------------------------------------------------------------
# load("results/data.RData",verbose=TRUE)
#
# #grDevices::pdf(file="manuscript/figure_COR.pdf",width=6,height=3)
# grDevices::postscript(file="manuscript/figure_COR.eps",width=6,height=3)
# graphics::par(mar=c(0.5,3,2,0.5))
# graphics::layout(mat=matrix(c(1,2),nrow=1,ncol=2),width=c(0.2,0.8))
#
# # correlation between years
# cor <- cbind(sapply(Y,function(x) cor(x[,1],x[,2],use="complete.obs",method="spearman")),
# sapply(Y,function(x) cor(x[,2],x[,3],use="complete.obs",method="spearman")),
# sapply(Y,function(x) cor(x[,1],x[,3],use="complete.obs",method="spearman")))
# colnames(cor) <- c("1-2","2-3","1-3")
# cor <- rowMeans(cor)
# joinet:::plot.matrix(cor,range=c(-3,3),margin=1,cex=0.7)
#
# # correlation between variables
# cor <- 1/3*cor(sapply(Y,function(x) x[,1]),use="complete.obs",method="spearman")+
# 1/3*cor(sapply(Y,function(x) x[,2]),use="complete.obs",method="spearman")+
# 1/3*cor(sapply(Y,function(x) x[,3]),use="complete.obs",method="spearman")
# joinet:::plot.matrix(cor,range=c(-3,3),margin=c(1,2),cex=0.7)
# grDevices::dev.off()
#
# # other information
# sapply(Y,colMeans,na.rm=TRUE) # increasing values
# sapply(Y,function(x) apply(x,2,sd,na.rm=TRUE)) # increasing variance
# sapply(Y,function(x) colSums(is.na(x))) # increasing numbers of NAs
## ----application,eval=FALSE---------------------------------------------------
# #<<start>>
# library(joinet)
#
# set.seed(1)
# load("results/data.RData",verbose=TRUE)
#
# set.seed(1)
# foldid.ext <- rep(1:5,length.out=nrow(Y$moca))
# foldid.int <- rep(rep(1:10,each=5),length.out=nrow(Y$moca))
# table(foldid.ext,foldid.int)
#
# #- - - - - - - - - - - - -
# #- - internal coaching - -
# #- - - - - - - - - - - - -
#
# table <- list()
# table$alpha <- c("lasso","ridge")
# table$data <- c("clinic","omics","both")
# table$var <- names(Y)
# table <- rev(expand.grid(table,stringsAsFactors=FALSE))
#
# loss <- fit <- list()
# for(i in seq_len(nrow(table))){
# cat(rep("*",times=5),"setting",i,rep("*",times=5),"\n")
# y <- Y[[table$var[i]]]
# x <- list(clinic=X,omics=Z,both=cbind(X,Z))[[table$data[i]]]
# alpha <- 1*(table$alpha[i]=="lasso")
# loss[[i]] <- cv.joinet(Y=y,X=x,alpha.base=alpha,foldid.ext=foldid.ext,
# foldid.int=foldid.int,sign=1) # add joinet::
# #fit[[i]] <- joinet(Y=y,X=x,alpha.base=alpha,foldid=foldid.int,sign=1)
# }
#
# save(table,loss,file="results/internal.RData")
#
# #- - - - - - - - - - - - -
# #- - external coaching - -
# #- - - - - - - - - - - - -
#
# table <- list()
# temp <- utils::combn(x=names(Y),m=2)
# table$comb <- paste0(temp[1,],"-",temp[2,])
# table$step <- c("V04","V06","V08")
# table$alpha <- c("lasso","ridge")
# table$data <- c("clinic","omics","both")
# table <- rev(expand.grid(table,stringsAsFactors=FALSE))
# temp <- strsplit(table$comb,split="-"); table$comb <- NULL
# table$var1 <- sapply(temp,function(x) x[[1]])
# table$var2 <- sapply(temp,function(x) x[[2]])
#
# loss <- fit <- list()
# for(i in seq_len(nrow(table))){
# cat(rep("*",times=5),"setting",i,rep("*",times=5),"\n")
# y <- cbind(Y[[table$var1[i]]][,table$step[i]],
# Y[[table$var2[i]]][,table$step[i]])
# x <- list(clinic=X,omics=Z,both=cbind(X,Z))[[table$data[i]]]
# alpha <- 1*(table$alpha[i]=="lasso")
# loss[[i]] <- cv.joinet(Y=y,X=x,alpha.base=alpha,
# foldid.ext=foldid.ext,foldid.int=foldid.int,sign=1) # add joinet::
# #fit[[i]] <- joinet(Y=y,X=x,alpha.base=alpha,foldid=foldid.int,sign=1)
# }
#
# save(table,loss,file="results/external.RData")
# writeLines(text=capture.output(utils::sessionInfo(),cat("\n"),sessioninfo::session_info()),con="results/info_app.txt")
## ----figure_INT---------------------------------------------------------------
# load("results/internal.RData")
#
# # standardised loss
# vars <- unique(table$var)
# former <- t(sapply(loss,function(x) x["base",]))
# min <- sapply(vars,function(x) min(former[table$var==x,]))
# max <- sapply(vars,function(x) max(former[table$var==x,]))
# index <- match(x=table$var,table=vars)
# former <- (former-min[index])/(max[index]-min[index])
# dimnames(former) <- list(table$var,seq_len(3))
#
# # percentage change
# change <- t(sapply(loss,function(x) 100*(x["meta",]-x["base",])/x["base",]))
# dimnames(change) <- list(table$var,c("1st","2nd","3rd"))
#
# # overview
# #grDevices::pdf(file="manuscript/figure_INT.pdf",width=6,height=3,pointsize=12)
# grDevices::postscript(file="manuscript/figure_INT.eps",width=6,height=3,pointsize=12)
# graphics::par(mfrow=c(2,3),mar=c(0.1,3,2,0.1),oma=c(0,1.1,1,0))
# for(alpha in c("lasso","ridge")){
# for(data in c("clinic","omics","both")){
# cond <- table$alpha==alpha & table$data==data
# joinet:::plot.matrix(X=change[cond,],range=c(-50,50),cex=0.7)
# #graphics::title(main=paste0(alpha,"-",data),col.main="red",line=0) # check
# if(alpha=="lasso"){graphics::mtext(text=data,side=3,line=1.5,cex=0.8)}
# if(data=="clinic"){graphics::mtext(text=alpha,side=2,line=3,cex=0.8)}
# }
# }
# grDevices::dev.off()
#
# TEMP <- tapply(X=rowMeans(change),INDEX=table$var,FUN=mean)[vars]
# mean(change<0)
# round(tapply(X=rowMeans(change),INDEX=table$data,FUN=mean),digits=2)
# round(tapply(X=rowMeans(change),INDEX=table$alpha,FUN=mean),digits=2)
# round(colMeans(change),digits=2)
## ----figure_EXT---------------------------------------------------------------
# load("results/external.RData")
#
# data <- c("clinic","omics","both")
# alpha <- c("lasso","ridge")
# step <- c("V04","V06","V08")
#
# # percentage change
# change <- t(sapply(loss,function(x) 100*(x["meta",]-x["base",])/x["base",]))
#
# # overview
# vars <- unique(c(table$var1,table$var2))
# temp <- matrix(NA,nrow=length(vars),ncol=length(vars),dimnames=list(vars,vars))
# array <- array(data=list(temp),dim=c(3,2,3),dimnames=list(data,alpha,step))
# #grDevices::pdf(file="manuscript/figure_EXT.pdf",width=7.5,height=10,pointsize=14)
# grDevices::postscript(file="manuscript/figure_EXT.eps",width=7.5,height=10,pointsize=14)
# graphics::par(mfrow=c(6,3),mar=c(0.1,2.5,2.5,0.1),oma=c(0,1,2,0))
# for(i in data){
# for(j in alpha){
# for(k in step){
# cond <- table$data==i & table$alpha==j & table$step==k
# array[i,j,k][[1]][cbind(table$var1,table$var2)[cond,]] <- change[cond,1]
# array[i,j,k][[1]][cbind(table$var2,table$var1)[cond,]] <- change[cond,2]
# joinet:::plot.matrix(array[i,j,k][[1]],margin=0,las=2,range=c(-20,20),cex=0.6)
# #graphics::title(main=paste0(i,"-",j,"-",k),col.main="red",line=0) # check
# if(i=="clinic" & j=="lasso"){graphics::mtext(text=ifelse(k=="V04","1st",ifelse(k=="V06","2nd","3rd")),side=3,line=2.5,cex=0.8)}
# if(k=="V04"){graphics::mtext(text=paste0(i,"-",j),side=2,line=2.5,cex=0.8)}
# }
# }
# }
# grDevices::dev.off()
#
# # check
# i <- sample(seq_len(nrow(table)),size=1)
# table[i,]
# x <- loss[[i]]
# 100*(x["meta",]-x["base",])/x["base",]
## ----figure_ALL---------------------------------------------------------------
# #grDevices::pdf(file="manuscript/figure_ALL.pdf",height=3,width=6)
# grDevices::postscript(file="manuscript/figure_ALL.eps",height=3,width=6)
#
# graphics::par(mar=c(0.5,3,2,0.5))
# graphics::layout(mat=matrix(c(1,2),nrow=1,ncol=2),width=c(0.2,0.8))
# joinet:::plot.matrix(as.matrix(TEMP),margin=1,las=1,range=c(-20,20),cex=0.7)
#
# sum(unlist(array)<0,na.rm=TRUE)/sum(!is.na(unlist(array)))
# means <- apply(array,c(1,2,3),function(x) mean(x[[1]],na.rm=TRUE))
# lapply(seq_len(3),function(x) apply(means,x,mean))
# mean <- 1/length(array)*Reduce(f="+",x=array)
# joinet:::plot.matrix(mean,margin=1,las=1,range=c(-20,20),cex=0.7)
# # rows: target variable, columns: coaching variable
#
# grDevices::dev.off()
## ----figure_DIF---------------------------------------------------------------
# #grDevices::pdf(file="manuscript/figure_DIF.pdf",height=1.2,width=5)
# grDevices::postscript(file="manuscript/figure_DIF.eps",height=1.2,width=5)
#
# load("results/internal.RData")
# vars <- unique(table$var)
# base <- t(sapply(loss,function(x) 100*(x["base",]-x["none",])/x["none",]))
# meta <- t(sapply(loss,function(x) 100*(x["meta",]-x["none",])/x["none",]))
# dimnames(meta) <- dimnames(base) <- list(table$var,c("1st","2nd","3rd"))
# standard <- tapply(X=rowMeans(base),INDEX=table$var,FUN=mean)[vars]
# internal <- tapply(X=rowMeans(meta),INDEX=table$var,FUN=mean)[vars]
#
# load("results/external.RData")
# vars <- unique(c(table$var1,table$var2))
# base <- meta <- list()
# for(i in seq_len(2)){
# base[[i]] <- sapply(loss,function(x) 100*(x["base",i]-x["none",i])/x["none",i])
# meta[[i]] <- sapply(loss,function(x) 100*(x["meta",i]-x["none",i])/x["none",i])
# }
# index <- c(table$var1,table$var2)
# base <- unlist(base); meta <- unlist(meta)
# #standard <- tapply(X=base,INDEX=index,FUN=mean)[vars]
# external <- tapply(X=meta,INDEX=index,FUN=mean)[vars]
#
# matrix <- round(rbind(standard,internal,external),digits=2)
# rownames(matrix) <- c("","","")
# graphics::par(mfrow=c(1,1),mar=c(0.5,3,1.5,1))
# joinet:::plot.matrix(matrix,margin=c(1,2),las=1,range=c(-100,0),cex=0.7,digits=3)
#
# grDevices::dev.off()
Try the joinet package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.