R/plot.R

Defines functions plot.interflex

Documented in plot.interflex

plot.interflex <- function(x,
						  order = NULL,
						  subtitles = NULL,
						  show.subtitles = NULL,
						  CI = NULL,
					      diff.values = NULL,
						  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.sub = NULL,
						  cex.lab = NULL,
						  cex.axis = NULL,
						  bin.labs = TRUE, # bin labels    
						  interval = NULL, # interval in replicated papers
						  file = NULL,
						  ncols = NULL,
						  #pool plot
						  pool = FALSE,
						  legend.title = NULL,
						  color = NULL,
						  show.all = FALSE,
						  scale = 1.1,
  						  height = 7,
  						  width = 10,
						  ...
){	
	y <- NULL
	xmin <- NULL
	xmax <- NULL
	ymin <- NULL
	count1 <- NULL
	ymax <- NULL
	X <- NULL
	TE <- NULL
	CI_lower <- NULL
	CI_upper <- NULL
	ME <- NULL
	x0 <- NULL



	if(pool==TRUE){
		p <- interflex.plot.pool(out=x,
								diff.values = diff.values,
								order = order,
								subtitles = subtitles,
								show.subtitles = show.subtitles,
								legend.title = legend.title,
								CI = CI,
								Xdistr = Xdistr,
								main = main,
								Ylabel = Ylabel,
								Dlabel = Dlabel,
								Xlabel = Xlabel,
								xlab = xlab,
								ylab = ylab,
								xlim = xlim,
								ylim = ylim,
								theme.bw = theme.bw,
								show.grid = show.grid,     
								cex.main = cex.main,
								cex.lab = cex.lab,
								cex.axis = cex.axis,
								cex.sub = cex.sub,
								bin.labs = bin.labs, # bin labels    
								interval = interval, # interval in replicated papers
								color = color,
								file = file)
		return(p)
	}

	
	out <- x
	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.")
		}
	}
	
	# 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("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.")
		} 
	}
	
	## 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){
			if(length(other.treat)==1){
				show.subtitles <- FALSE
			}
			else{
				show.subtitles <- TRUE
			}
      
			if(is.null(subtitles)==FALSE){
				show.subtitles <- TRUE
			}
		}
    
		if(is.null(subtitles)==FALSE){
			if(length(subtitles)!=length(other.treat)){
				stop("The number of elements in \"subtitles\" should be m-1(m is the number of different treatment arms).")
			}
		}
	}
		
	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){
			if(length(label.name)==1){
				show.subtitles <- FALSE
			}
			else{
				show.subtitles <- TRUE
			}
      
			if(is.null(subtitles)==FALSE){
				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.")
			}
		}
	}
	
	## ncols
	if (is.null(ncols) == FALSE) {
		if (ncols%%1 != 0) {
			stop("\"ncols\" is not a positive integer.")
		} else {
			ncols <- ncols[1]
		}
		if (ncols < 1) {
			stop("\"ncols\" is not a positive integer.")
		}
	} else{
		if(treat.type=="discrete"){
			ncols <- length(other.treat)
		}
		if(treat.type=="continuous"){
			ncols <- length(label.name)
		}
	}
	treat_sc <- max(1,ncols-1)
	
	## diff.values
	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)
	
	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
		}
	}
	
	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
	}
	
	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	
	}
	
	
	#plot initialization
	p.group <- list()
	if(treat.type=='discrete'){
		for(char in other.treat){
			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())
			}
			p.group[[char]] <- p1
		}
	}
	
	if(treat.type=='continuous'){
		for(label in label.name){
			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())
			}
			p.group[[label]] <- p1
		}
	}
						
	#density
	if (Xdistr == "density") { 
		if (treat.type == 'discrete') { 
			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/5 + min(yrange) - maxdiff/5)
      
			## color
			feed.col <- col2rgb("gray50")
			col.co <- rgb(feed.col[1]/1000, feed.col[2]/1000,feed.col[3]/1000)
			col.tr <- rgb(red=1, blue=0, green=0)
    
			for(char in other.treat){
				deX.tr <- data.frame(x = de.tr[[char]]$x,
									 y = de.tr[[char]]$y/max(de.tr[[char]]$y) * maxdiff/5 + min(yrange) - maxdiff/5)
        
				p1 <- p.group[[char]] + geom_ribbon(data = deX.co, aes(x = x, ymax = y, ymin = deX.ymin),
													fill = col.co, alpha = 0.2) +
										geom_ribbon(data = deX.tr, aes(x = x, ymax = y, ymin = deX.ymin),
													fill = col.tr, alpha = 0.2)
				p.group[[char]] <- p1
			}
		} 
    
		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)
      
			for(label in label.name){
				## plotting
				p1 <- p.group[[label]]
				p1 <- p1 + geom_ribbon(data = deX, aes(x = x, ymax = y, ymin = deX.ymin),
									   fill = col, alpha = 0.2)
				p.group[[label]] <- p1
			}
		}
	}
  
	#hist
	if(Xdistr %in% c("histogram","hist")){
		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) 
      
			for (char in other.treat){
				hist.treat<-data.frame(ymin=hist.col[,'count1'],
									   #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[[char]]/hist.max*maxdiff/5+hist.col[,'count1']) 
        
				p1 <- p.group[[char]] + geom_rect(data=hist.col,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=count1),fill="gray50",
											colour="gray50",alpha=0.3,size=0.3) + # control
										geom_rect(data=hist.treat,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=count1),
											fill="red",colour="gray50",alpha=0.3,size=0.3)
				p.group[[char]] <- p1
			}
		}
		
		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)
			for(label in label.name){
				p1 <- p.group[[label]]
				p1 <- p1 + geom_rect(data=histX,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
									 colour="gray50",alpha=0.3,size=0.5)
				p.group[[label]] <- p1
			}
		}
	
	}
	
	#ME/TE in kernel/linear
	if(estimator=='kernel'|estimator=='linear'){
		if(estimator=='kernel'){
			est <- est.kernel
		}else{est <- est.lin}
		
		if(treat.type=='discrete'){
			for(char in other.treat){
				p1 <- p.group[[char]]
				tempest <- est[[char]]
				if(CI==TRUE){ 
					colnames(tempest) <- c("X","TE","sd","CI_lower","CI_upper")
					tempest <- as.data.frame(tempest)
				}
				if(CI==FALSE){
					tempest <- tempest[,c(1,2)]
					colnames(tempest) <- c("X","TE")
					tempest <- as.data.frame(tempest)
				}
				p1 <-  p1 + geom_line(data=tempest,aes(X,TE),size=1.2)
				if (CI == TRUE) {
					p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper),alpha=0.2)
				}
				#ymin=min(yrange)-maxdiff/5
		
				if(is.null(diff.values)==FALSE){
					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,"TE"]
							if(CI==TRUE){
								lb.mark <- tempest[label1,"CI_lower"]
								ub.mark <- tempest[label1,"CI_upper"]
							}
						}  
						else if(d2==0){
							est.mark <- tempest[label2,"TE"]
							if(CI==TRUE){
								lb.mark <- tempest[label2,"CI_lower"]
								ub.mark <- tempest[label2,"CI_upper"]
							}
						} 
						else{ ## weighted average
							est.mark1 <- tempest[label1,"TE"]
							est.mark2 <- tempest[label2,"TE"]
							est.mark <- ((est.mark1 * d2 + est.mark2 * d1)/(d1 + d2))
							if(CI==TRUE){
								lb.mark1 <- tempest[label1,"CI_lower"]
								ub.mark1 <- tempest[label1,"CI_upper"]
								lb.mark2 <- tempest[label2,"CI_lower"]
								ub.mark2 <- tempest[label2,"CI_upper"]
								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='red')
						if(CI==TRUE){
							p1 <- p1+ annotate("errorbar",x=target.value,ymin=lb.mark,ymax=ub.mark,colour='red',size=0.5,width= (max(tempxx)-min(tempxx))/30)
						}
					}
				}
				p.group[[char]] <- p1
			}
		}
		
		if(treat.type=='continuous'){
			for(label in label.name){
				p1 <- p.group[[label]]
				tempest <- est[[label]]
				if(CI==TRUE){
					colnames(tempest) <- c("X","ME","sd","CI_lower","CI_upper")
					tempest <- as.data.frame(tempest)
				}
				if(CI==FALSE){
					tempest <- tempest[,c(1,2)]
					colnames(tempest) <- c("X","ME")
					tempest <- as.data.frame(tempest)
				}
				p1 <-  p1 + geom_line(data=tempest,aes(X,ME),size=1.2)
				if (CI == TRUE) {
					p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper),alpha=0.2)
				}
				#ymin=min(yrange)-maxdiff/5
		
				if(is.null(diff.values)==FALSE){
					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,"ME"]
							if(CI==TRUE){
								lb.mark <- tempest[label1,"CI_lower"]
								ub.mark <- tempest[label1,"CI_upper"]
							}
						}  
						else if(d2==0){
							est.mark <- tempest[label2,"ME"]
							if(CI==TRUE){
								lb.mark <- tempest[label2,"CI_lower"]
								ub.mark <- tempest[label2,"CI_upper"]
							}
						} 
						else{ ## weighted average
							est.mark1 <- tempest[label1,"ME"]
							est.mark2 <- tempest[label2,"ME"]
							est.mark <- ((est.mark1 * d2 + est.mark2 * d1)/(d1 + d2))
							if(CI==TRUE){
								lb.mark1 <- tempest[label1,"CI_lower"]
								ub.mark1 <- tempest[label1,"CI_upper"]
								lb.mark2 <- tempest[label2,"CI_lower"]
								ub.mark2 <- tempest[label2,"CI_upper"]
								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='red')
						if(CI==TRUE){
							p1 <- p1+ annotate("errorbar",x=target.value,ymin=lb.mark,ymax=ub.mark,colour='red',size=0.5,width= (max(tempxx)-min(tempxx))/30)
						}
					}
				}
				
				p.group[[label]] <- p1
			}

		}
	
	}
	
	if(estimator=='binning'){
		if(treat.type=='discrete'){ 
			for(char in other.treat){	
				p1 <- p.group[[char]]
				tempest <- est.lin[[char]]
				tempest.bin2 <- est.bin2[[char]]
				if(CI==TRUE){
					colnames(tempest) <- c("X","TE","sd","CI_lower","CI_upper")
					colnames(tempest.bin2) <- c("x0","TE","sd","CI_lower","CI_upper")
					tempest <- as.data.frame(tempest)
					tempest.bin2 <- as.data.frame(tempest.bin2)
				}
				if(CI==FALSE){
					tempest <- tempest[,c(1,2)]
					colnames(tempest) <- c("X","TE")
					tempest.bin2 <- tempest.bin2[,c(1,2)]
					colnames(tempest.bin2) <- c("x0","TE")
					tempest <- as.data.frame(tempest)
					tempest.bin2 <- as.data.frame(tempest.bin2)
				}	
					
				p1 <- p1 + geom_line(data=tempest,aes(X,TE),size=1.2)
				if(CI==TRUE){
					p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper),alpha=0.2)
				}
				## bin estimates
				p1 <- p1 + geom_point(data=tempest.bin2,aes(x0,TE),size=4/treat_sc,shape=21,fill="white",colour="red")
				if(CI==TRUE){
					p1<-p1+ geom_errorbar(data=tempest.bin2, aes(x=x0, ymin=CI_lower, ymax=CI_upper),colour="red",size = 1,
										  width= errorbar.width)
				}
									   
				## in case there's non-overlap
				
				
				p1<-p1+annotate(geom="text", x=est.bin3[[char]][,1], y=rep(0,dim(est.bin3[[char]])[1]),
								label="NaN",colour="red")  
				
				## labels: L, M, H and so on 
				if (bin.labs == TRUE) {
					if (nbins==3) {
						p1<-p1 + annotate(geom="text", x=est.bin[[char]][1,1], y=pos,
										  label="L",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[char]][2,1], y=pos,
										  label="M",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[char]][3,1], y=pos,
										  label="H",colour="gray50",size=10/treat_sc)
					} 
					else if (nbins==4) {
						p1<-p1 + annotate(geom="text", x=est.bin[[char]][1,1], y=pos,
										  label="L",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[char]][2,1], y=pos,
										  label="M1",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[char]][3,1], y=pos,
										  label="M2",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[char]][4,1], y=pos,
										  label="H",colour="gray50",size=10/treat_sc)
					}
					else if (nbins==2) {
						p1<-p1 + annotate(geom="text", x=est.bin[[char]][1,1], y=pos,
										  label="L",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[char]][2,1], y=pos,
										  label="H",colour="gray50",size=10/treat_sc)
					}
				}
				p.group[[char]] <- p1
			}
		}
		
		if(treat.type=='continuous'){ 
			for(label in label.name){	
				p1 <- p.group[[label]]
				tempest <- est.lin[[label]]
				tempest.bin2 <- est.bin2[[label]]
				if(CI==TRUE){
					colnames(tempest) <- c("X","ME","sd","CI_lower","CI_upper")
					colnames(tempest.bin2) <- c("x0","ME","sd","CI_lower","CI_upper")
					tempest <- as.data.frame(tempest)
					tempest.bin2 <- as.data.frame(tempest.bin2)
				}
				if(CI==FALSE){
					tempest <- tempest[,c(1,2)]
					colnames(tempest) <- c("X","ME")
					tempest.bin2 <- tempest.bin2[,c(1,2)]
					colnames(tempest.bin2) <- c("x0","ME")
					tempest <- as.data.frame(tempest)
					tempest.bin2 <- as.data.frame(tempest.bin2)	
				}	
					
				p1 <- p1 + geom_line(data=tempest,aes(X,ME),size=1.2)
				if(CI==TRUE){
					p1 <- p1 + geom_ribbon(data=tempest, aes(x=X,ymin=CI_lower,ymax=CI_upper),alpha=0.2)
				}
				## bin estimates
				p1 <- p1 + geom_point(data=tempest.bin2,aes(x0,ME),size=4/treat_sc,shape=21,fill="white",colour="red")
				if(CI==TRUE){
					p1<-p1+ geom_errorbar(data=tempest.bin2, aes(x=x0, ymin=CI_lower, ymax=CI_upper),colour="red",size = 1,
										  width= errorbar.width)
				}	   
				## in case there's non-overlap
				p1<-p1+annotate(geom="text", x=est.bin3[[label]][,1], y=rep(0,dim(est.bin3[[label]])[1]),
								label="NaN",colour="red")  
        
				## labels: L, M, H and so on 
				if (bin.labs == TRUE) {
					if (nbins==3) {
						p1<-p1 + annotate(geom="text", x=est.bin[[label]][1,1], y=pos,
										  label="L",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[label]][2,1], y=pos,
										  label="M",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[label]][3,1], y=pos,
										  label="H",colour="gray50",size=10/treat_sc)
					} 
					else if (nbins==4) {
						p1<-p1 + annotate(geom="text", x=est.bin[[label]][1,1], y=pos,
										  label="L",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[label]][2,1], y=pos,
										  label="M1",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[label]][3,1], y=pos,
										  label="M2",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[label]][4,1], y=pos,
										  label="H",colour="gray50",size=10/treat_sc)
					}
					else if (nbins==2) {
						p1<-p1 + annotate(geom="text", x=est.bin[[label]][1,1], y=pos,
										  label="L",colour="gray50",size=10/treat_sc) +
								 annotate(geom="text", x=est.bin[[label]][2,1], y=pos,
										  label="H",colour="gray50",size=10/treat_sc)
					}
				}
				p.group[[label]] <- p1
			}
		}
	
	}
	
	#cex/title...
    if (is.null(cex.lab)==TRUE) {
      cex.lab <- 15
    } else {
      cex.lab <- 15 * cex.lab
    }
    if (is.null(cex.axis)==TRUE) {
      cex.axis <- 15/treat_sc
    } else {
      cex.axis <- 15 * cex.axis/treat_sc
    }
    ## title
    if (is.null(cex.main)==TRUE) {
      cex.main <- 18
    } else {
      cex.main <- 18 * cex.main
    }
    
    if (is.null(cex.sub)==TRUE) {
      cex.sub <- 12
    } else {
      cex.sub <- 12 * cex.sub
    }
    
    ## 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(treat.type=='discrete'){
		k <- 1
		for(char in other.treat) {
			p1 <- p.group[[char]]
			## 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
			p1 <- p1 + xlab(NULL) + ylab(NULL) + 
				  theme(axis.text = element_text(size=cex.axis))

			if(show.subtitles==TRUE){
				if(is.null(subtitles)==TRUE){
					subtitle.temp <- paste0("Treated = ",char,", Baseline = ",base)
					p1 <- p1 + labs(subtitle = subtitle.temp) + theme(plot.subtitle = element_text(hjust = 0.5, size=cex.sub,lineheight=.8))
				}
			
				if(is.null(subtitles)==FALSE){
					subtitle.temp<- subtitles[k]
					p1 <- p1 + labs(subtitle = subtitle.temp) + theme(plot.subtitle = element_text(hjust = 0.5, size=cex.sub, lineheight=.8))
				}
			}

			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)
			} 

			p.group[[char]] <- p1
			k <- k+1
		}
    
		## ylim
		for(char in other.treat){
			y.limits <- layer_scales(p.group[[char]])$y$range$range
			if(char == other.treat[1]){
				ymaxmax <- y.limits[2]
				yminmin <- y.limits[1]
			}
			else{
				ymaxmax <- max(ymaxmax,y.limits[2])
				yminmin <- min(yminmin,y.limits[1])
			}
		}
		for(char in other.treat){
			p.group[[char]] <- p.group[[char]]+ylim(c(yminmin,ymaxmax))
		}
    
		requireNamespace("gridExtra")
		requireNamespace("grid")
		requireNamespace("ggplotify")
		suppressMessages(
			graph <- arrangeGrob(grobs=p.group,ncol=ncols,align = 'h',
								 top = textGrob(main,gp=gpar(fontsize=cex.main,face="bold")),
								 left = textGrob(ylab, rot = 90, vjust=1,gp=gpar(fontsize=cex.lab)),
								 bottom = textGrob(xlab, gp=gpar(fontsize=cex.lab))
								)
		)
	}	
    
	if(treat.type=='continuous'){
		k <- 1
		for(label in label.name) {
			p1 <- p.group[[label]]
			## 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
			p1 <- p1 + xlab(NULL) + ylab(NULL) + 
				  theme(axis.text = element_text(size=cex.axis))

			if(show.subtitles==TRUE){
				if(is.null(subtitles)==TRUE){
					subtitle.temp <- label
					p1 <- p1 + labs(subtitle = subtitle.temp) + theme(plot.subtitle = element_text(hjust = 0.5, size=cex.sub,lineheight=.8))
				}
			
				if(is.null(subtitles)==FALSE){
					subtitle.temp<- subtitles[k]
					p1 <- p1 + labs(subtitle = subtitle.temp) + theme(plot.subtitle = element_text(hjust = 0.5, size=cex.sub, lineheight=.8))
				}
			}

			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)
			} 

			p.group[[label]] <- p1
			k <- k+1
		}
    
		## ylim
		for(label in label.name){
			y.limits <- layer_scales(p.group[[label]])$y$range$range
			if(label == label.name[1]){
				ymaxmax <- y.limits[2]
				yminmin <- y.limits[1]
			}
			else{
				ymaxmax <- max(ymaxmax,y.limits[2])
				yminmin <- min(yminmin,y.limits[1])
			}
		}
		for(label in label.name){
			p.group[[label]] <- p.group[[label]]+ylim(c(yminmin,ymaxmax))
		}
    
		requireNamespace("gridExtra")
		requireNamespace("grid")
		requireNamespace("ggplotify")
		suppressMessages(
			graph <- arrangeGrob(grobs=p.group,ncol=ncols,align = 'h',
								 top = textGrob(main,gp=gpar(fontsize=cex.main,face="bold")),
								 left = textGrob(ylab, rot = 90, vjust=1,gp=gpar(fontsize=cex.lab)),
								 bottom = textGrob(xlab, gp=gpar(fontsize=cex.lab))
								)
		)
	}	

	
	
	## save to file
	if (is.null(file)==FALSE) {
		graph <- as.ggplot(graph)
		ggsave(file, graph,scale = scale,width=width,height = height)
    }

	if(show.all==FALSE){
		graph <- as.ggplot(graph)
		return(graph)
	}
	else{
		cex.axis <- cex.axis
		if(treat.type=='discrete'){
			for(char in other.treat){
			## Other universal options
				p1 <- p.group[[char]] + xlab(xlab) + ylab(ylab) + 
					  theme(axis.text = element_text(size=cex.axis),
							axis.title.x = element_text(size=cex.lab),
							axis.title.y = element_text(size=cex.lab))
				p1 <- p1 + ggtitle(main) +
					  theme(plot.title = element_text(hjust = 0.5, size=cex.main,lineheight=.8, face="bold"))
				p.group[[char]] <- p1
			}
		}
		
		if(treat.type=='continuous'){
			for(label in label.name){
			## Other universal options
				p1 <- p.group[[label]] + xlab(xlab) + ylab(ylab) + 
					  theme(axis.text = element_text(size=cex.axis),
							axis.title.x = element_text(size=cex.lab),
							axis.title.y = element_text(size=cex.lab))	  
				p1 <- p1 + ggtitle(main) +
					  theme(plot.title = element_text(hjust = 0.5, size=cex.main,lineheight=.8, face="bold"))
				p.group[[label]] <- p1
			}
		}
	return(p.group)
	}
  }

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