R/plot_pool.R

Defines functions interflex.plot.pool

Documented in interflex.plot.pool

interflex.plot.pool <- function( # only for discrete treatments 
  out,
  diff.values = NULL,
  order = NULL,
  subtitles = NULL,
  show.subtitles = NULL,
  legend.title = NULL,
  CI = TRUE,
  Xdistr = "histogram",
  main = NULL,
  Ylabel = NULL,
  Dlabel = NULL,
  Xlabel = NULL,
  xlab = NULL,
  ylab = NULL,
  xlim = NULL,
  ylim = NULL,
  theme.bw = FALSE,
  show.grid = TRUE,     
  cex.main = NULL,
  cex.lab = NULL,
  cex.axis = NULL,
  cex.sub = NULL,
  bin.labs = TRUE, # bin labels    
  interval = NULL, # interval in replicated papers
  color = NULL,
  file = NULL,
  scale = 1.1,
  height = 7,
  width = 10
){
	X <- NULL
	TE <- NULL
	Treatment <- NULL
	CI_lower <- NULL
	CI_upper <- NULL
	ME <- NULL
	x0 <- NULL
	CI.lower <- NULL
	CI.upper <- NULL
	x <- NULL
	y <- NULL
	end_level <- NULL
	xmin <- NULL
	xmax <- NULL
	count1 <- NULL
	ymax <- NULL
	D <- NULL
	name <- NULL
	r <- NULL
	

	if (!class(out) %in% c("interflex")) {
		stop("Not an \"interflex\" object.")
	} 
						
	treat.info <- out$treat.info
	treat.type <- treat.info[["treat.type"]]
	if(treat.type=='discrete'){
		other.treat <- names(treat.info[["other.treat"]])
		all.treat <- names(treat.info[["all.treat"]])
		base <- names(treat.info[["base"]])
	}
	if(treat.type=='continuous'){
		D.sample <- treat.info[["D.sample"]]
		label.name <- names(D.sample)
	}
	
	de <- out$de
	de.tr <- out$de.tr
	hist.out <- out$hist.out
	count.tr <- out$count.tr
	estimator <- out$estimator
	
	if(is.null(show.subtitles)==FALSE){
		if (is.logical(show.subtitles) == FALSE & is.numeric(show.subtitles)==FALSE) {
			stop("\"show.subtitles\" is not a logical flag.")
		}
	}else{
		show.subtitles <- TRUE
	}
	
	# CI
	if(is.null(CI)==FALSE){
		if(is.logical(CI) == FALSE & is.numeric(CI)==FALSE) {
			stop("\"CI\" is not a logical flag.")
		}
		if(estimator=='kernel'){
			if(CI==TRUE & out$CI==FALSE){
				stop("Confidence intervals are not estimated, please set CI to FALSE.")
			}
		}		
	}
	
	if(estimator=='kernel'){
		if(is.null(CI)==TRUE){
			CI <- out$CI
		}
	}
	
	if(estimator=='binning'|estimator=='linear'){
		if(is.null(CI)==TRUE){
			CI <- TRUE
		}
	}

	# Xdistr
	if (!Xdistr %in% c("hist","histogram","density","none")){
		stop("\"Xdistr\" must be \"histogram\", \"density\", or \"none\".")
	}

	#main
	if (is.null(main)==FALSE) {
		main <- as.character(main)[1]
	}
  
	#Ylabel
	if (is.null(Ylabel)==TRUE) {
		Ylabel <- out$Ylabel
	} else {
		if (is.character(Ylabel) == FALSE) {
			stop("\"Ylabel\" is not a string.")
		} else {
			Ylabel <- Ylabel[1]
		}   
	}
  
	#Dlabel  
	if (is.null(Dlabel)==TRUE) {
		Dlabel <- out$Dlabel   
	} else {
		if (is.character(Dlabel) == FALSE) {
			stop("\"Dlabel\" is not a string.")
		} else {
			Dlabel <- Dlabel[1]
		}   
	}

	#Xlabel
	if (is.null(Xlabel)==TRUE) {
		Xlabel <- out$Xlabel  
	} else {
		if (is.character(Xlabel) == FALSE) {
			stop("\"Xlabel\" is not a string.")
		} else {
			Xlabel <- Xlabel[1]
		}   
	}
  
	## axis labels
	if(is.null(xlab)==FALSE){
		if (is.character(xlab) == FALSE) {
			stop("\"xlab\" is not a string.")
		}        
	}
	if(is.null(ylab)==FALSE){
		if (is.character(ylab) == FALSE) {
			stop("\"ylab\" is not a string.")
		}        
	}
  
	if(is.null(xlab)==TRUE){
		xlab<-c(paste("Moderator: ", Xlabel, sep=""))
	} else {
		if (is.character(xlab) == FALSE) {
			stop("\"xlab\" is not a string.")
		}        
	}
	if(is.null(ylab)==TRUE){
		ylab<-c(paste("Marginal Effect of ",Dlabel," on ",Ylabel,sep=""))
	} else {
		if (is.character(ylab) == FALSE) {
			stop("\"ylab\" is not a string.")
		}        
	}
  
	## xlim ylim
	if (is.null(xlim)==FALSE) {
		if (is.numeric(xlim)==FALSE) {
			stop("Some element in \"xlim\" is not numeric.")
		} else {
			if (length(xlim)!=2) {
				stop("\"xlim\" must be of length 2.")
			}
		}
	}
  
	if (is.null(ylim)==FALSE) {
		if (is.numeric(ylim)==FALSE) {
			stop("Some element in \"ylim\" is not numeric.")
		} else {
			if (length(ylim)!=2) {
				stop("\"ylim\" must be of length 2.")
			}
		}
	}
  
	## theme.bw
	if (is.logical(theme.bw) == FALSE & is.numeric(theme.bw)==FALSE) {
		stop("\"theme.bw\" is not a logical flag.")
	}
  
	## show.grid
	if (is.logical(show.grid) == FALSE & is.numeric(show.grid)==FALSE) {
		stop("\"show.grid\" is not a logical flag.")
	}
  
	## font size
	if (is.null(cex.main)==FALSE) {
		if (is.numeric(cex.main)==FALSE) {
			stop("\"cex.main\" is not numeric.")
		}
	}
	if (is.null(cex.sub)==FALSE) {
		if (is.numeric(cex.sub)==FALSE) {
			stop("\"cex.sub\" is not numeric.")
		}
	}
	if (is.null(cex.lab)==FALSE) {
		if (is.numeric(cex.lab)==FALSE) {
			stop("\"cex.lab\" is not numeric.")
		}
	}
	if (is.null(cex.axis)==FALSE) {
		if (is.numeric(cex.axis)==FALSE) {
			stop("\"cex.axis\" is not numeric.")
		}    
	}

	## bin.labs
	if (is.logical(bin.labs) == FALSE & is.numeric(bin.labs)==FALSE) {
		stop("\"bin.labs\" is not a logical flag.")
	}
  
	## interval
	if (is.null(interval)==FALSE) {
		if (is.numeric(interval)==FALSE) {
			stop("Some element in \"interval\" is not numeric.")
		} 
	}

	## file
	if (is.null(file)==FALSE) {
		if (is.character(file)==FALSE) {
			stop("Wrong file name.")
		} 
	}
	
	## color
	if(is.null(color)==FALSE){
		color <- as.character(color)
		color.in <- c()
		for(char in color){
			res <- try(col2rgb(char),silent=TRUE)
			if(!"try-error"%in%class(res)){
				color.in <- c(color.in,char)
			}else{stop(paste0(char," is not one name for a color.\n"))}
		}
		color <- color.in
	}
	
	if (is.null(legend.title)==FALSE){
		legend.title <- as.character(legend.title)[1]
	}
	
	if(treat.type=='discrete' & (estimator=='linear'|estimator=='binning')){
		tempxx <- out$est.lin[[other.treat[1]]][,'X']
	}
	if(treat.type=='discrete' & estimator=='kernel'){
		tempxx <- out$est.kernel[[other.treat[1]]][,'X']
	}
	if(treat.type=='continuous' & (estimator=='linear'|estimator=='binning')){
		tempxx <- out$est.lin[[label.name[1]]][,'X']
	}
	if(treat.type=='continuous' & estimator=='kernel'){
		tempxx <- out$est.kernel[[label.name[1]]][,'X']
	}
	min.XX <- min(tempxx)
	max.XX <- max(tempxx)
	

	

	## order/subtitles
	if(treat.type=='discrete') {
		other.treat <- sort(all.treat[which(all.treat!=base)])
		if(is.null(order)==FALSE){
			order <- as.character(order)
			if(length(order)!=length(unique(order))){
				stop("\"order\" should not contain repeated values.")
			}
	
			if(length(order)!=length(other.treat)){
				stop("\"order\" should include all kinds of treatment arms except for the baseline group.")
			}

			if(sum(!is.element(order,other.treat))!=0 | sum(!is.element(other.treat,order))!=0){
				stop("\"order\" should include all kinds of treatment arms except for the baseline group.")
			}
			other.treat <- order
		}
	
		if(is.null(show.subtitles)==TRUE){
			show.subtitles <- TRUE
		}
    
		if(is.null(subtitles)==FALSE){
			if(length(subtitles)!=length(all.treat)){
				stop("The number of elements in \"subtitles\" should be m(m is the number of different treatment arms including the baseline group).")
			}
		}
		
		if(is.null(subtitles)==TRUE){
			base.name <- paste0("Base Group (",base,")")
			subtitles <- c(base.name,other.treat)
		}
		else{base.name <- subtitles[1]}
		
		subtitles.all <- as.character(subtitles)
		subtitles <- subtitles.all[2:length(subtitles.all)]
		
	}
		
	if(treat.type=='continuous'){
		if(is.null(order)==FALSE){
			if(is.numeric(order)==FALSE){
				stop("\"order\" should be numeric.")
			}
			if(length(order)!=length(unique(order))){
				stop("\"order\" should not contain repeated values.")
			}
			if(length(order)!=length(D.sample)){
				stop("\"order\" should contain all reference values of D.")
			}
			if(sum(!is.element(order,D.sample))!=0 | sum(!is.element(D.sample,order))!=0){
				stop("\"order\" should contain all reference values of D.")
			}
			label.name.order <- c()
			for(a in order){
				label.name.order <- c(label.name.order,names(D.sample[which(D.sample==a)]))
			}
			label.name <- label.name.order
		}
		
		if(is.null(show.subtitles)==TRUE){
			show.subtitles <- TRUE
		}
    
		if(is.null(subtitles)==FALSE){
			if(length(subtitles)!=length(label.name)){
				stop("The number of elements in \"subtitles\" should equal to the number of values in D.ref.")
			}
		}
		
		if(is.null(subtitles)==TRUE){
			subtitles <- label.name
		}
		
		subtitles.all <- subtitles
		
	}
	
	if(is.null(diff.values)==FALSE){
		if(estimator=='binning'){
			stop("\"diff.values\" can only work after linear or kernel model is applied.")
		}
		if(is.numeric(diff.values)==FALSE){
			stop("\"diff.values\" is not numeric.")
		}
		if(length(diff.values)<2){
			stop("\"diff.values\" must be of length 2 or more.")
		}
		if(treat.type=='discrete' & estimator=='linear'){
			tempxx <- out$est.lin[[other.treat[1]]][,'X']
		}
		if(treat.type=='discrete' & estimator=='kernel'){
			tempxx <- out$est.kernel[[other.treat[1]]][,'X']
		}
		if(treat.type=='continuous' & estimator=='linear'){
			tempxx <- out$est.lin[[label.name[1]]][,'X']
		}
		if(treat.type=='continuous' & estimator=='kernel'){
			tempxx <- out$est.kernel[[label.name[1]]][,'X']
		}
		min.XX <- min(tempxx)
		max.XX <- max(tempxx)
		for(a in diff.values){
			if(a<min.XX|a>max.XX){
				stop("Elements in \"diff.values\" should be within the range of the moderator.")
			}
		}
	}else{
		if(estimator=='binning'){
			diff.values <- NULL
		}
		else{
			#diff.values <- out$diff.info[["diff.values.plot"]]
			diff.values <- NULL
		}
	}
  
	#yrange
	if(estimator=='binning'){
		nbins <- out$nbins
		if(treat.type=='discrete'){
			est.lin <- out$est.lin
			est.bin <- out$est.bin
			est.bin2 <- list() ## non missing part
			est.bin3 <- list() ## missing part
			yrange <- c(0)
			for(char in other.treat) {
				est.bin2[[char]] <- as.matrix(est.bin[[char]][which(is.na(est.bin[[char]][,2])==FALSE),]) 
				est.bin3[[char]] <- as.matrix(est.bin[[char]][which(is.na(est.bin[[char]][,2])==TRUE),])
				if(dim(est.bin2[[char]])[2]==1){
					est.bin2[[char]] <- t(est.bin2[[char]])
				}
				if(dim(est.bin3[[char]])[2]==1){
					est.bin3[[char]] <- t(est.bin3[[char]])
				}
				if(CI==TRUE){
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[char]][,c(4,5)],est.bin[[char]][,c(4,5)]))))
				}else{
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[char]][,2],est.bin[[char]][,2]))))
				}
			}
      
			if (is.null(ylim)==FALSE) {yrange<-c(ylim[2],ylim[1]+(ylim[2]-ylim[1])*1/8)}
			X.lvls <- est.lin[[other.treat[1]]][,1]
			errorbar.width<-(max(X.lvls)-min(X.lvls))/20
			maxdiff<-(max(yrange)-min(yrange))
			pos<-max(yrange)-maxdiff/20
		}
		
		if (treat.type=='continuous'){
			est.lin <- out$est.lin
			est.bin<-out$est.bin
			est.bin2 <- list() ## non missing part
			est.bin3 <- list() ## missing part
			yrange <- c(0)
			for(label in label.name){
				est.bin2[[label]] <- as.matrix(est.bin[[label]][which(is.na(est.bin[[label]][,2])==FALSE),]) 
				est.bin3[[label]] <- as.matrix(est.bin[[label]][which(is.na(est.bin[[label]][,2])==TRUE),])
				if(dim(est.bin2[[label]])[2]==1){
					est.bin2[[label]] <- t(est.bin2[[label]])
				}
				if(dim(est.bin3[[label]])[2]==1){
					est.bin3[[label]] <- t(est.bin3[[label]])
				}
				if(CI==TRUE){
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[label]][,c(4,5)],est.bin[[label]][,c(4,5)]))))
				}else{
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[label]][,2],est.bin[[label]][,2]))))
				}
			}
			
			X.lvls <- est.lin[[label.name[1]]][,1]
			errorbar.width<-(max(X.lvls)-min(X.lvls))/20
			if (is.null(ylim)==FALSE) {yrange<-c(ylim[2],ylim[1]+(ylim[2]-ylim[1])*1/8)}
			maxdiff<-(max(yrange)-min(yrange))
			pos<-max(yrange)-maxdiff/20
		}
		ymin <- min(yrange)-maxdiff/5
	}
	
	

	if(estimator=='linear'){
		if(treat.type=='discrete'){
			est.lin <- out$est.lin
			yrange <- c(0)
			for(char in other.treat) {
				if(CI==TRUE){
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[char]][,c(4,5)]))))
				}else{
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[char]][,2]))))
				}
			}
			X.lvls <- est.lin[[other.treat[1]]][,1]
		}
		
		if(treat.type=='continuous'){
			est.lin <- out$est.lin
			yrange <- c(0)
			for(label in label.name){
				if(CI==TRUE){
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[label]][,c(4,5)]))))
				}else{
					yrange <- c(yrange,na.omit(unlist(c(est.lin[[label]][,2]))))
				}
			}
			X.lvls <- est.lin[[label.name[1]]][,1]
		}
		
		errorbar.width <- (max(X.lvls)-min(X.lvls))/20
		if (is.null(ylim)==FALSE) {yrange<-c(ylim[2],ylim[1]+(ylim[2]-ylim[1])*1/8)}
		maxdiff<-(max(yrange)-min(yrange))
		pos<-max(yrange)-maxdiff/20
		ymin <- min(yrange)-maxdiff/5
	}
	
	if(estimator=='kernel'){
		est.kernel <- out$est.kernel
		yrange <- c(0)
		if(CI==FALSE){
			if(treat.type=='discrete'){
				for(char in other.treat) {
					yrange <- c(yrange,na.omit(unlist(c(est.kernel[[char]][,2]))))
				}
				X.lvls <- est.kernel[[other.treat[1]]][,1]
			}
			
			if (treat.type=='continuous'){
				for(label in label.name){
					yrange <- c(yrange,na.omit(unlist(c(est.kernel[[label]][,2]))))
				}
				X.lvls <- est.kernel[[label.name[1]]][,1]
			}	
		}
		
		if(CI==TRUE){
			if(treat.type=='discrete'){
				for(char in other.treat) {
					yrange <- c(yrange,na.omit(unlist(c(est.kernel[[char]][,c(4,5)]))))
				}
				X.lvls <- est.kernel[[other.treat[1]]][,1]
			}
			
			if (treat.type=='continuous'){
				for(label in label.name){
					yrange <- c(yrange,na.omit(unlist(c(est.kernel[[label]][,c(4,5)]))))
				}
				X.lvls <- est.kernel[[label.name[1]]][,1]
			}
		}

		if (is.null(ylim)==FALSE) {yrange<-c(ylim[2],ylim[1]+(ylim[2]-ylim[1])*1/8)}
		errorbar.width<-(max(X.lvls)-min(X.lvls))/20
		maxdiff<-(max(yrange)-min(yrange))
		pos<-max(yrange)-maxdiff/20	
		ymin <- min(yrange)-maxdiff/5
	}
	
	#color
	#get base.color & platte(for discrete data)
	requireNamespace("RColorBrewer")
	platte <- brewer.pal(n=8, "Set2")
	
	if(is.null(color)==TRUE){
		base.color <- 'gray50'
	}
	
	if(is.null(color)==FALSE){
		if(treat.type=='discrete'){
			base.color <- color[1]
			if(length(color)==1){
				platte <- platte
			}else{platte <- c(color[2:length(color)],platte)}
		}
		
		if(treat.type=='continuous'){
			platte <- c(color,platte)
		}
	}
	
	if(treat.type=='discrete'){
		num.treat <- length(other.treat)
	}
	
	if(treat.type=='continuous'){
		num.treat <- length(label.name)
	}
	
	platte <- platte[1:num.treat]
	# initialize
	p1 <- ggplot()
  
	## black white theme and mark zero
	if (theme.bw == FALSE) {
		p1 <- p1 + geom_hline(yintercept=0,colour="white",size=2)
	} else {
		p1 <- p1 + theme_bw() + geom_hline(yintercept=0,colour="#AAAAAA50",size=2)
	}
  
	if (show.grid == FALSE) {
		p1 <- p1 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
	}
	
	if(estimator=='kernel'|estimator=='linear'){
		if(estimator=='kernel'){
			est <- est.kernel
		}else{est <- est.lin}
		
		if(treat.type=='discrete'){	  
			for(char in other.treat){
				est.touse <- est[[char]]
				if(CI==TRUE){ 
					colnames(est.touse) <- c("X","TE","sd","CI_lower","CI_upper")
				}
				if(CI==FALSE){
					est.touse <- est.touse[,c(1,2)]
					colnames(est.touse) <- c("X","TE")
				}
				est.touse <- as.data.frame(est.touse)
				est.touse[['Treatment']] <- rep(char,dim(est.touse)[1])
				
				if(char==other.treat[1]){
					tempest <- est.touse
				}else{
					tempest <- rbind(tempest,est.touse)
				}
			}
			tempest <- as.data.frame(tempest)
			tempest$Treatment <- factor(tempest$Treatment, levels = other.treat)
			p1 <- p1 + geom_line(data=tempest,aes(x = X,y = TE,color = Treatment),show.legend = FALSE)
			p1 <- p1 + scale_color_manual(values = platte,labels = subtitles)
			if (CI == TRUE) {
				p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper,fill = Treatment),alpha=0.2,show.legend = F)
				p1 <- p1 + scale_fill_manual(values = platte,labels = subtitles)
			}
	  
			if(is.null(diff.values)==FALSE){
				k <- 1
				for(char in other.treat){
					tempest <- est[[char]]
					for(target.value in diff.values){
						Xnew<-abs(tempest[,'X']-target.value)
						d1<-min(Xnew)     
						label1<-which.min(Xnew)
						Xnew[label1]<-Inf
						d2<-min(Xnew)     
						label2<-which.min(Xnew)
						if(d1==0){
							est.mark <- tempest[label1,2]
							if(CI==TRUE){
								lb.mark <- tempest[label1,4]
								ub.mark <- tempest[label1,5]
							}
						}  
						else if(d2==0){
							est.mark <- tempest[label2,2]
							if(CI==TRUE){
								lb.mark <- tempest[label2,4]
								ub.mark <- tempest[label2,5]
							}
						} 
						else{ ## weighted average
							est.mark1 <- tempest[label1,2]
							est.mark2 <- tempest[label2,2]
							est.mark <- ((est.mark1 * d2 + est.mark2 * d1)/(d1 + d2))
							if(CI==TRUE){
								lb.mark1 <- tempest[label1,4]
								ub.mark1 <- tempest[label1,5]
								lb.mark2 <- tempest[label2,4]
								ub.mark2 <- tempest[label2,5]
								lb.mark <- ((lb.mark1 * d2 + lb.mark2 * d1)/(d1 + d2))
								ub.mark <- ((ub.mark1 * d2 + ub.mark2 * d1)/(d1 + d2))
							}
						}
						
						p1 <- p1 + annotate("point",x=target.value,y=est.mark,size=1,colour=platte[k])
						if(CI==TRUE){
							p1 <- p1+ annotate("errorbar",x=target.value,ymin=lb.mark,ymax=ub.mark,colour=platte[k],size=0.5,width= (max(tempxx)-min(tempxx))/20)
						}
					}
					k <- k + 1
				}
			}
		}
		
		if(treat.type=='continuous'){
			for(label in label.name){
				est.touse <- est[[label]]
				if(CI==TRUE){ 
					colnames(est.touse) <- c("X","ME","sd","CI_lower","CI_upper")
				}
				if(CI==FALSE){
					est.touse <- est.touse[,c(1,2)]
					colnames(est.touse) <- c("X","ME")
				}
				
				est.touse <- as.data.frame(est.touse)
				est.touse[['Treatment']] <- rep(label,dim(est.touse)[1])
				
				if(label==label.name[1]){
					tempest <- est.touse
				}else{
					tempest <- rbind(tempest,est.touse)
				}
			}
			tempest <- as.data.frame(tempest)
			tempest$Treatment <- factor(tempest$Treatment, levels = label.name)
			p1 <- p1 + geom_line(data=tempest,aes(x = X,y = ME,color = Treatment),show.legend = FALSE)
			p1 <- p1 + scale_color_manual(values = platte,labels = subtitles)
			if (CI == TRUE) {
				p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper,fill = Treatment),alpha=0.2,show.legend = FALSE)
				p1 <- p1 + scale_fill_manual(values = platte,labels = subtitles)
			}
			
			if(is.null(diff.values)==FALSE){
				k <- 1
				for(label in label.name){
					tempest <- est[[label]]
					for(target.value in diff.values){
						Xnew<-abs(tempest[,'X']-target.value)
						d1<-min(Xnew)     
						label1<-which.min(Xnew)
						Xnew[label1]<-Inf
						d2<-min(Xnew)     
						label2<-which.min(Xnew)
						if(d1==0){
							est.mark <- tempest[label1,2]
							if(CI==TRUE){
								lb.mark <- tempest[label1,4]
								ub.mark <- tempest[label1,5]
							}
						}  
						else if(d2==0){
							est.mark <- tempest[label2,2]
							if(CI==TRUE){
								lb.mark <- tempest[label2,4]
								ub.mark <- tempest[label2,5]
							}
						} 
						else{ ## weighted average
							est.mark1 <- tempest[label1,2]
							est.mark2 <- tempest[label2,2]
							est.mark <- ((est.mark1 * d2 + est.mark2 * d1)/(d1 + d2))
							if(CI==TRUE){
								lb.mark1 <- tempest[label1,4]
								ub.mark1 <- tempest[label1,5]
								lb.mark2 <- tempest[label2,4]
								ub.mark2 <- tempest[label2,5]
								lb.mark <- ((lb.mark1 * d2 + lb.mark2 * d1)/(d1 + d2))
								ub.mark <- ((ub.mark1 * d2 + ub.mark2 * d1)/(d1 + d2))
							}
						}
						
						p1 <- p1 + annotate("point",x=target.value,y=est.mark,size=1,colour=platte[k])
						if(CI==TRUE){
							p1 <- p1+ annotate("errorbar",x=target.value,ymin=lb.mark,ymax=ub.mark,colour=platte[k],size=0.5,width= (max(tempxx)-min(tempxx))/20)
						}
					}
					k <- k + 1
				}
			}
		}
	
	}
	
	if(estimator=='binning'){
		#est <- est.lin
		if(treat.type=='discrete'){
			for(char in other.treat){
				est.touse <- est.lin[[char]]
				if(CI==TRUE){ 
					colnames(est.touse) <- c("X","TE","sd","CI_lower","CI_upper")
				}
				if(CI==FALSE){
					est.touse <- est.touse[,c(1,2)]
					colnames(est.touse) <- c("X","TE")
				}
				est.touse <- as.data.frame(est.touse)
				est.touse[['Treatment']] <- rep(char,dim(est.touse)[1])
				
				if(char==other.treat[1]){
					tempest <- est.touse
				}else{
					tempest <- rbind(tempest,est.touse)
				}
			}
			tempest$Treatment <- factor(tempest$Treatment, levels = other.treat)
			p1 <- p1 + geom_line(data=tempest,aes(x = X,y = TE,color = Treatment),show.legend = FALSE)
			p1 <- p1 + scale_color_manual(values = platte, labels = subtitles)
			if (CI == TRUE){
				p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper,fill = Treatment),alpha=0.2,show.legend = FALSE)
				p1 <- p1 + scale_fill_manual(values = platte,labels = subtitles)
			}
			
			k <- 1
			for(char in other.treat){
				tempest2 <- as.data.frame(est.bin2[[char]])
				tempest3 <- as.data.frame(est.bin3[[char]])
				p1 <- p1+ geom_errorbar(data=tempest2, aes(x=x0, ymin=CI.lower, ymax=CI.upper),color = platte[k],
										width= errorbar.width/3)+
						  geom_point(data=tempest2,aes(x=x0,y=coef),size=3,shape=21,fill = platte[k],color = platte[k])
      
				if(dim(tempest3)[1]!=0){
					p1 <- p1+geom_text(data=tempest3,aes(x=x0,y=0),label="NaN",color = platte[k])
				}
				k <- k+1
			}	
		}
		
		if(treat.type=='continuous'){
			for(label in label.name){
				est.touse <- est.lin[[label]]
				if(CI==TRUE){ 
					colnames(est.touse) <- c("X","ME","sd","CI_lower","CI_upper")
				}
				if(CI==FALSE){
					est.touse <- est.touse[,c(1,2)]
					colnames(est.touse) <- c("X","ME")
				}
				est.touse <- as.data.frame(est.touse)
				est.touse[['Treatment']] <- rep(label,dim(est.touse)[1])
				
				if(label==label.name[1]){
					tempest <- est.touse
				}else{
					tempest <- rbind(tempest,est.touse)
				}
			}
			tempest$Treatment <- factor(tempest$Treatment, levels = label.name)
			p1 <- p1 + geom_line(data=tempest,aes(x = X,y = ME,color = Treatment),show.legend = FALSE)
			p1 <- p1 + scale_color_manual(values = platte, labels = subtitles)
			if (CI == TRUE){
				p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper,fill = Treatment),alpha=0.2,show.legend = FALSE)
				p1 <- p1 + scale_fill_manual(values = platte,labels = subtitles)
			}
			
			k <- 1
			for(label in label.name){
				tempest2 <- as.data.frame(est.bin2[[label]])
				tempest3 <- as.data.frame(est.bin3[[label]])
     
				p1 <- p1+ geom_errorbar(data=tempest2, aes(x=x0, ymin=CI.lower, ymax=CI.upper),color = platte[k],
										width= errorbar.width/3)+
						  geom_point(data=tempest2,aes(x=x0,y=coef),size=3,shape=21,fill = platte[k],color = platte[k])
      
				if(dim(tempest3)[1]!=0){
					p1 <- p1+geom_text(data=tempest3,aes(x=x0,y=0),label="NaN",color = platte[k])
				}
				k <- k+1
			}
		}
		
		if (bin.labs == TRUE){
			if(treat.type=='discrete'){
				char0 <- other.treat[1]
			}
			if(treat.type=='continuous'){
				char0 <- label.name[1]
			}
			if (nbins==3){
				p1 <- p1 + annotate(geom="text", x=est.bin[[char0]][1,1], y=pos,
									label="L",colour="gray50",size=10) +
						   annotate(geom="text", x=est.bin[[char0]][2,1], y=pos,
									label="M",colour="gray50",size=10) +
						   annotate(geom="text", x=est.bin[[char0]][3,1], y=pos,
									label="H",colour="gray50",size=10)
			} 
			else if (nbins==4){
				p1 <- p1 + annotate(geom="text", x=est.bin[[char0]][1,1], y=pos,
									label="L",colour="gray50",size=10) +
						   annotate(geom="text", x=est.bin[[char0]][2,1], y=pos,
								    label="M1",colour="gray50",size=10) +
						   annotate(geom="text", x=est.bin[[char0]][3,1], y=pos,
								    label="M2",colour="gray50",size=10) +
						   annotate(geom="text", x=est.bin[[char0]][4,1], y=pos,
								    label="H",colour="gray50",size=10)
			} 
			else if (nbins==2){
					p1 <- p1 + annotate(geom="text", x=est.bin[[char0]][1,1], y=pos,
										label="L",colour="gray50",size=10) +
							   annotate(geom="text", x=est.bin[[char0]][2,1], y=pos,
									    label="H",colour="gray50",size=10)
			} 
		}
	}
	
	if(Xdistr == "density") { # density plot
		if(treat.type=='discrete'){
			## put in data frames
			dist<-hist.out$mids[2]-hist.out$mids[1]
			deX.ymin <- min(yrange)-maxdiff/5
			deX.co <- data.frame(x = de.tr[[base]]$x,
								y = de.tr[[base]]$y/max(de.tr[[base]]$y) * maxdiff/10 + min(yrange) - maxdiff/5)
			## color
			p1 <- p1 + geom_ribbon(data = deX.co, aes(x = x, ymax = y, ymin = deX.ymin),color=base.color,
								fill = base.color, alpha = 0.0, size=0.3)
			k <- 1
			char0 <- other.treat[1]
			start_level <- rep(deX.ymin,length(de.tr[[char0]]$x))
			for(char in other.treat){
				dex.tr.plot <- data.frame(x = de.tr[[char]]$x,
										start_level = start_level,
										end_level = de.tr[[char]]$y/max(de.tr[[char0]]$y)*maxdiff/10+start_level)
      
				p1 <- p1+geom_ribbon(data = dex.tr.plot, aes(x = x, ymax = end_level, ymin = start_level), color=platte[k],
									alpha = 0.0,fill = platte[k],size=0.3)
      
				k <- k+1
			}
			p1 <- p1+geom_line(data = dex.tr.plot, aes(x = x, y = ymin), color='gray50',size=0.3)
		}
		
		if(treat.type=='continuous'){
			deX.ymin <- min(yrange)-maxdiff/5
			deX <- data.frame(x = de$x,
							  y = de$y/max(de$y) * maxdiff/5 + min(yrange) - maxdiff/5) 
      
			## color
			feed.col<-col2rgb("gray50")
			col<-rgb(feed.col[1]/1000, feed.col[2]/1000,feed.col[3]/1000)
			p1 <- p1 + geom_ribbon(data = deX, aes(x = x, ymax = y, ymin = deX.ymin),
								   fill = col, alpha = 0.2)

		}
	}
	
	if (Xdistr %in% c("histogram","hist")){ # histogram plot
		if(treat.type=='discrete'){
			n.hist<-length(hist.out$mids)
			dist<-hist.out$mids[2]-hist.out$mids[1]
			hist.max<-max(hist.out$counts)
			hist.col<-data.frame(ymin=rep(min(yrange)-maxdiff/5,n.hist),
								ymax=hist.out$counts/hist.max*maxdiff/5+min(yrange)-maxdiff/5,
								xmin=hist.out$mids-dist/2,
								xmax=hist.out$mids+dist/2,
								count1=count.tr[[base]]/hist.max*maxdiff/5+min(yrange)-maxdiff/5) 
			
			p1 <- p1 + geom_rect(data=hist.col,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=count1),fill=base.color,color='gray50',
								 alpha=0.3,size=0.5) 
      
			k <- 1
			start_level <- count.tr[[base]]/hist.max*maxdiff/5+min(yrange)-maxdiff/5
			for (char in other.treat){
				hist.treat<-data.frame(ymin=start_level,
									   ymax=count.tr[[char]]/hist.max*maxdiff/5+start_level,
									   xmin=hist.out$mids-dist/2,
									   xmax=hist.out$mids+dist/2)
        
				start_level <- count.tr[[char]]/hist.max*maxdiff/5+start_level
        
				p1 <- p1 + geom_rect(data=hist.treat,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),fill=platte[k],color='gray50',
									 alpha=0.5,size=0.5)
				k <- k + 1
			}
		}
		
		if(treat.type=='continuous'){
			n.hist <- length(hist.out$mids)
			dist <- hist.out$mids[2]-hist.out$mids[1]
			hist.max <- max(hist.out$counts)            
			histX <- data.frame(ymin=rep(min(yrange)-maxdiff/5,n.hist),
							  ymax=hist.out$counts/hist.max*maxdiff/5+min(yrange)-maxdiff/5,
							  xmin=hist.out$mids-dist/2,
							  xmax=hist.out$mids+dist/2)
			p1 <- p1 + geom_rect(data=histX,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
									 colour="gray50",alpha=0.3,size=0.5)
		}	
    }
 
	
  ## other properties
  if(is.null(legend.title)==FALSE){
    p1 <- p1 + labs(fill = legend.title,color = legend.title)
  }
 
  ## mark the original interval (in replicated papers)
  if (is.null(interval)==FALSE) {
    p1<- p1 + geom_vline(xintercept=interval,colour="steelblue", linetype=2,size=1.5)
  }
  
  ## Other universal options
  ## axis labels
  if (is.null(cex.lab)==TRUE) {
    cex.lab <- 15
  } else {
    cex.lab <- 15 * cex.lab
  }
  if (is.null(cex.axis)==TRUE) {
    cex.axis <- 15
  } else {
    cex.axis <- 15 * cex.axis
  }
  p1 <- p1 + xlab(xlab) + ylab(ylab) + 
    theme(axis.text = element_text(size=cex.axis), axis.title = element_text(size=cex.lab))
  
  ## title
  if (is.null(cex.main)==TRUE) {
    cex.main <- 18
  } else {
    cex.main <- 18 * cex.main
  }
  
  if (is.null(cex.sub)==TRUE) {
    cex.sub <- 10
  } else {
    cex.sub <- 10 * cex.sub
  }
  
  if (is.null(main)==FALSE) {
    p1<-p1 + ggtitle(main) +
      theme(plot.title = element_text(hjust = 0.5, size=cex.main,
                                      lineheight=.8, face="bold"))
  } 
  
  ## xlim and ylim
  if (is.null(ylim)==FALSE) {
    ylim2 = c(ylim[1]-(ylim[2]-ylim[1])*0.25/6, ylim[2]+(ylim[2]-ylim[1])*0.4/6)
  }
  if (is.null(xlim)==FALSE & is.null(ylim)==FALSE) {
    p1<-p1+coord_cartesian(xlim = xlim, ylim = ylim2)
  }
  if (is.null(xlim)==TRUE & is.null(ylim)==FALSE) {
    p1<-p1+coord_cartesian(ylim = ylim2)
  }
  if (is.null(xlim)==FALSE & is.null(ylim)==TRUE) {
    p1<-p1+coord_cartesian(xlim = xlim)
  } 
	
  #legend
  if(show.subtitles==TRUE){
	if(treat.type=='discrete'){
		p1_table <- ggplot_gtable(ggplot_build(p1))
		data.touse3 <- data.frame(X=rep(1,length(all.treat)),ymin=-1,ymax=1,D=all.treat)
		data.touse3$D <- factor(data.touse3$D,levels = all.treat)
		p0 <- ggplot() + geom_ribbon(data=data.touse3, aes(x=X,ymin=ymin,ymax=ymax,fill=D),
									alpha=0.3)
		p0 <-  p0 + scale_fill_manual(values = c(base.color,platte[1:length(other.treat)]),
									  labels = as.character(subtitles.all))
		if(is.null(legend.title)==FALSE){
			p0 <- p0 + labs(fill = legend.title,color = legend.title)
		}else{p0 <- p0 + labs(fill = "Treatment",color = "Treatment")}
		
		p0 <- p0 + theme(legend.title = element_text(colour="black", size=cex.sub),
						 legend.text = element_text(color = "black", size = cex.sub*0.95))
	
		p0 <- p0 + xlab(xlab) + ylab(ylab) + theme(axis.text = element_text(size=cex.axis), axis.title = element_text(size=cex.lab))
	
		if(is.null(main)==FALSE){
			p0 <- p0 + ggtitle(main) + theme(plot.title = element_text(hjust = 0.5, size=cex.main,
                                       lineheight=.8, face="bold"))
		}
	
		if (is.null(ylim)==FALSE) {
			ylim2 = c(ylim[1]-(ylim[2]-ylim[1])*0.25/6, ylim[2]+(ylim[2]-ylim[1])*0.4/6)
		}
		if (is.null(xlim)==FALSE & is.null(ylim)==FALSE) {
			p0<-p0+coord_cartesian(xlim = xlim, ylim = ylim2)
		}
		if (is.null(xlim)==TRUE & is.null(ylim)==FALSE) {
			p0<-p0+coord_cartesian(ylim = ylim2)
		}
		if (is.null(xlim)==FALSE & is.null(ylim)==TRUE) {
			p0<-p0+coord_cartesian(xlim = xlim)
		}
		y.limits <- layer_scales(p1)$y$range$range
		x.limits <- layer_scales(p1)$x$range$range
		ymaxmax <- y.limits[2]
		yminmin <- y.limits[1]
		xmaxmax <- x.limits[2]
		xminmin <- x.limits[1]

		suppressWarnings(
			p0 <- p0+ylim(c(yminmin,ymaxmax))+xlim(c(xminmin,xmaxmax))
		)
		suppressWarnings(
			p0 <- ggplot_gtable(ggplot_build(p0))
		)
		suppressWarnings(
			pp <-c(subset(p0$layout, name == "panel", se=t:r))
		)

		gt <- gtable_add_grob(p0,
							p1_table$grobs[[which(p1_table$layout$name == "panel")]],
							pp$t,pp$l,pp$b,pp$l)
	
		gt <- as.ggplot(gt)
		p1 <- gt
		
		
	}
	
	if(treat.type=='continuous'){
		p1_table <- ggplot_gtable(ggplot_build(p1))
		data.touse3 <- data.frame(X=rep(1,length(label.name)),ymin=-1,ymax=1,D=label.name)
		data.touse3$D <- factor(data.touse3$D,levels = label.name)
		p0 <- ggplot() + geom_ribbon(data=data.touse3, aes(x=X,ymin=ymin,ymax=ymax,fill=D),
									alpha=0.3)
		p0 <-  p0 + scale_fill_manual(values = platte[1:length(label.name)],
									  labels = as.character(subtitles.all))
		if(is.null(legend.title)==FALSE){
			p0 <- p0 + labs(fill = legend.title,color = legend.title)
		}else{p0 <- p0 + labs(fill = "Treatment",color = "Treatment")}
		
		p0 <- p0 + theme(legend.title = element_text(colour="black", size=cex.sub),
						 legend.text = element_text(color = "black", size = cex.sub*0.95))
	
		p0 <- p0 + xlab(xlab) + ylab(ylab) + theme(axis.text = element_text(size=cex.axis), axis.title = element_text(size=cex.lab))
	
		if(is.null(main)==FALSE){
			p0 <- p0 + ggtitle(main) + theme(plot.title = element_text(hjust = 0.5, size=cex.main,
                                       lineheight=.8, face="bold"))
		}
	
		## xlim and ylim
		if (is.null(ylim)==FALSE) {
			ylim2 = c(ylim[1]-(ylim[2]-ylim[1])*0.25/6, ylim[2]+(ylim[2]-ylim[1])*0.4/6)
		}
		if (is.null(xlim)==FALSE & is.null(ylim)==FALSE) {
			p0<-p0+coord_cartesian(xlim = xlim, ylim = ylim2)
		}
		if (is.null(xlim)==TRUE & is.null(ylim)==FALSE) {
			p0<-p0+coord_cartesian(ylim = ylim2)
		}
		if (is.null(xlim)==FALSE & is.null(ylim)==TRUE) {
			p0<-p0+coord_cartesian(xlim = xlim)
		}
		y.limits <- layer_scales(p1)$y$range$range
		x.limits <- layer_scales(p1)$x$range$range
		ymaxmax <- y.limits[2]
		yminmin <- y.limits[1]
		xmaxmax <- x.limits[2]
		xminmin <- x.limits[1]
		suppressWarnings(
			p0 <- p0+ylim(c(yminmin,ymaxmax))+xlim(c(xminmin,xmaxmax))
		)
		suppressWarnings(
			p0 <- ggplot_gtable(ggplot_build(p0))
		)
		suppressWarnings(
			pp <-c(subset(p0$layout, name == "panel", se=t:r))
		)

		gt <- gtable_add_grob(p0,
							p1_table$grobs[[which(p1_table$layout$name == "panel")]],
							pp$t,pp$l,pp$b,pp$l)
	
		gt <- as.ggplot(gt)
		p1 <- gt
	}
  }  

  ## save to file
	if (is.null(file)==FALSE) {
		ggsave(file, p1,scale = scale,width=width,height = height)
    }

  return(p1)


}
xuyiqing/interflex documentation built on Feb. 12, 2022, 10:10 p.m.