#Plot Residuals
library(ggplot2)

#Plot residuals and observed against predicted (obsdat, preddat, varname)
#Previously called the Show function.
resprplt <- function(obs,pred,varname) {
	res <- obs - pred
	Type1 <- rep("Residuals", length(res))
	Type2 <- rep("Observed", length(res))
	Data <- c(res, obs)
	Type <- c(Type1,Type2)
	Predicted <- rep(pred,2)
	graph <- data.frame(Data, Type, Predicted)
	unity <- function(x) {x}
	res_unity <- function(x) 0
	s1 <- "Predicted"
	s2 <- "Observed or Residual"
	s3 <- varname
	xlabel <- paste(s1, s3, sep = " ")
	ylabel <- paste(s2, s3, sep = " ")
	xlab <- xlabel
	ylab <- ylabel
	m1 <- ggplot(graph, aes(y=Data, x=Predicted, color=Type, shape=Type)) +
		xlab(xlab)+ylab(ylab)+
		scale_color_manual(values=c("black", "red"))+
		geom_point()+
		geom_smooth(method=lm, se=FALSE, linetype="dotted")+
		stat_function(fun=unity, color="black", size=1)+
		stat_function(fun=res_unity, color="red", size=1)+
		theme(axis.title=element_text(family= "serif", colour="black", size=16),
			axis.text=element_text(family= "serif", colour="black", size=16),
			legend.title=element_blank(),
			legend.text=element_text(family="serif"), 
			legend.position="top")
	dev.new()
	return(m1)
}


#Plot residuals and observed against predicted with lines by group (ObsDat, PredDat, GroupDat, VarName)
#ObsDat, PredDat, and GroupDat are vectors of observed, predicted, and grouping data. VarName is the variable name to be used on the axes.
resprgrpplt <- function(obs,pred,grp,varname) {
	res <- obs - pred
	Type1 <- rep("Residuals", length(res))
	Type2 <- rep("Observed", length(res))
	Data <- c(res, obs)
	Type <- c(Type1,Type2)
	Predicted <- rep(pred,2)
	Group <- rep(grp, 2)
	graph0 <- data.frame(Data, Type, Predicted, Group)
	graph <- subset(graph0, !is.na(graph0$Data))
	unity <- function(x) {x}
	res_unity <- function(x) 0
	s1 <- "Predicted"
	s2 <- "Observed or Residual"
	s3 <- varname
	xlabel <- paste(s1, s3, sep = " ")
	ylabel <- paste(s2, s3, sep = " ")
	xlab <- xlabel
	ylab <- ylabel
	m1 <- ggplot(graph, aes(y=Data, x=Predicted, group=interaction(Type,Group), color=Type, shape=Type)) +
		xlab(xlab)+ylab(ylab)+
		scale_color_manual(values=c("black", "red"))+
		geom_point()+
		#geom_line() +
		geom_smooth(method=lm, se=FALSE, linetype="solid", size=0.75)+
#		the following will add a quadratic prediction by PubID.  need 3 points/PubID.
#		geom_smooth(method=lm, formula = y ~ poly(x, 2), se=FALSE, linetype="solid", size=0.75)+
		stat_function(fun=unity, color="black", size=1.5)+
		stat_function(fun=res_unity, color="red", size=1.5)+
		theme(axis.title=element_text(family= "serif", colour="black", size=16),
			axis.text=element_text(family= "serif", colour="black", size=16),
			legend.title=element_blank(),
			legend.text=element_text(family="serif"), 
			legend.position="top")
	dev.new()
	return(m1)
}


#Plot residuals against another variable with lines by group (obs,pred,oth,grp, obsname,othname)
#obst, pred, oth, and grp are vectors of observed, predicted, other, and grouping data. obsname is the obs vare name, 
#	othname is the oth var name for axes labels.
resothgrpplt <- function(obs,pred,oth,grp,obsname,othname) {
	s1 <- "Predicted"
	s2 <- "Residual"
	s3 <- obsname
	s4 <- othname
	xlabel <- s4
	ylabel <- paste(s2, s3, sep = " ")
	xlab <- xlabel
	ylab <- ylabel
	res_unity <- function(x) 0
	graphdata <- data.frame(obs, pred, oth, as.factor(grp))
	graphdata$res <- graphdata$obs - graphdata$pred

	m1 <- ggplot(graphdata, aes(y=res, x=oth, group=grp, color="black")) +
		xlab(xlab)+ylab(ylab) + guides(fill=FALSE) +
		scale_color_manual(values="black")+	geom_point()+
		geom_smooth(method=lm, se=FALSE, linetype="solid", size=0.75)+
		stat_function(fun=res_unity, color="black", size=1)+
		theme(axis.title=element_text(family= "serif", colour="black", size=16),
			axis.text=element_text(family= "serif", colour="black", size=16),
			legend.title=element_text(family="serif"),
			legend.text=element_text(family="serif"), 
			legend.position="none")
	dev.new()
	return(m1)
}


	
#Plot observed and predicted against another variable (obs,pred,othdat,obsname,othname)
obsprothplt <- function(obs,pred,oth,obsname,othname) {
	Type1 <- rep("Predicted", length(pred))
	Type2 <- rep("Observed", length(oth))
	Data <- c(pred, obs)
	Type <- c(Type1,Type2)
	Oth <- rep(oth,2)
	graphdata <- data.frame(Oth,Data, Type)
	s1 <- "Predicted and Observed"
	s2 <- "Predicted"
	s3 <- obsname
	s4 <- othname
	xlabel <- paste(s2, s4, sep = " ")
	ylabel <- paste(s1, s3, sep = " ")
	xlab <- xlabel
	ylab <- ylabel
	m1 <- ggplot(graphdata, aes(y=Data, x=Oth, color=Type, shape=Type)) +
		xlab(xlab)+ylab(ylab)+
		scale_color_manual(values=c("black", "red"))+
		geom_point()+
		geom_smooth(method=lm, se=FALSE, linetype="dotted")+
		theme(axis.title=element_text(family= "serif", colour="black", size=16),
			axis.text=element_text(family= "serif", colour="black", size=16),
			legend.title=element_blank(),
			legend.text=element_text(family="serif"), 
			legend.position="top")
	dev.new()
	return(m1)
}


#Plot residuals against another variable (obs,pred,oth,obsname,othname)
resothplt <- function(obs,pred,oth,obsname,othname) {
	res <- obs - pred
	res_unity <- function(x) 0
	s1 <- "Predicted"
	s2 <- "Residual"
	s3 <- obsname
	s4 <- othname
	xlabel <- s4
	ylabel <- paste(s2, s3, sep = " ")
	xlab <- xlabel
	ylab <- ylabel
	graphdata <- data.frame(oth, res)
	m1 <- ggplot(graphdata, aes(y=res, x=oth, color="black")) +
		xlab(xlab)+ylab(ylab) + guides(fill=FALSE) +
		scale_color_manual(values="black")+	geom_point()+
		geom_smooth(method=lm, se=FALSE, linetype="dotted")+
		stat_function(fun=res_unity, color="black", size=1)+
		theme(axis.title=element_text(family= "serif", colour="black", size=16),
			axis.text=element_text(family= "serif", colour="black", size=16),
			legend.title=element_text(family="serif"),
			legend.text=element_text(family="serif"), 
			legend.position="none")
	dev.new()
	return(m1)
}

showPubID <- function(obs,pred,name) {
	res <- obs-pred
	PMeans <- t.test(res)$p.value
	SlopeBias <- coef(lm(res~pred))[2]
	Type1 <- rep("Residuals", length(res))
	Type2 <- rep("Observed", length(res))
	Data <- c(res, obs)
	Type <- c(Type1,Type2)
	Predicted <- rep(pred,2)
	Trt <- rep(rownames(pred),2)
	PubIDRef <- p[c("TID", "PubID")]
	q <- subset(PubIDRef, PubIDRef$TID %in% Trt)
	q$Trt <- q$TID
	unity <- function(x) {x}
	zero <- function(x) {0}
	s1 <- "Predicted"
	s2 <- "Observed or Residual"
	s3 <- name
	xlabel <- paste(s1, s3, sep = " ")
	ylabel <- paste(s2, s3, sep = " ")
	xlab <- xlabel
	ylab <- ylabel

	graph <- data.frame(Data, Type, Predicted, Trt)
	graph2 <- merge(graph, q, by="Trt", all.x=T)
	graph2$x <- mean(Predicted, na.rm=TRUE)
	graph2$y <- min(res, na.rm=TRUE) * 1.2
	graph2$x2 <- mean(Predicted, na.rm=TRUE)
	graph2$y2 <- min(res, na.rm=TRUE) * 1.4
	graph2$x3 <- mean(Predicted, na.rm=TRUE)
	graph2$y3 <- min(res, na.rm=TRUE) * 1.6
	graph2$x4 <- mean(Predicted, na.rm=TRUE)
	graph2$y4 <- min(res, na.rm=TRUE) * 1.8

	for(i in unique(graph2$PubID)) {
		pvalue <- subset(graph2, graph2$PubID == i & graph2$Type=="Residuals")
		pvalue$Check <- ifelse(is.na(pvalue$Predicted)==TRUE | is.na(pvalue$Data)==TRUE, 1, 0)
		if(sum(pvalue$Check>0)){
			slope <- NA} else{
			slope <- coefficients(lm(pvalue$Data~pvalue$Predicted))[[2]]
		}
		if(i == unique(graph2$PubID)[1]) {
			out <- data.frame(slope) 
		} else {
			out <- rbind(out, slope)
		}
	}
	out <- out[is.na(out$slope)==FALSE,]
	pval <- t.test(out)$p.value


	if (pval <=0.05) {
		label3 <- paste("Within study p-slope, P =", signif(pval,3)) 		#Within-study
	} else {
		label3 <- paste("Within study p-slope, P =", signif(pval,3))
	}

	if (PMeans <=0.05) {
		label2 <- paste("Across slope p-slope, P =", signif(PMeans,3)) 		#Global (t test of the residuals - same as RMSE bias)
	} else {
		label2 <- paste("Across slope p-slope, P =", signif(PMeans,3))		#Global
	}

	Meandat <- subset(graph2, graph2$Type=="Residuals")
	Means <- tapply(Meandat$Data, Meandat$PubID, mean)
	Pvalmean <- t.test(Means)$p.value

	if(Pvalmean <= 0.05) {
		label1 <- paste("Mean bias p-slope, P=", signif(Pvalmean, 3))		#Mean Bias
	} else {
		label1 <- paste("Mean bias p-slope, P=", signif(Pvalmean, 3))
	}


	SlpDiff <- out - SlopeBias
	Pvaldiff <- t.test(SlpDiff)$p.value

	if(Pvaldiff <= 0.05) {
		label4 <- paste("Difference p-slope, P=", signif(Pvaldiff, 3))		#Global - within
	} else {
		label4 <- paste("Difference p-slope, P=", signif(Pvaldiff, 3))
	}


	colorvec <- rep("black", length(unique(graph2$TID)))

	m1 <- ggplot(graph2, aes(y=Data, x=Predicted, color=PubID, shape=Type))+
		xlab(xlab)+ylab(ylab)+geom_point()+
		geom_text(aes(x=graph2$x, y=graph2$y, label=label1)) +
		geom_text(aes(x=graph2$x2, y=graph2$y2, label=label2)) + 
		geom_text(aes(x=graph2$x3, y=graph2$y3, label=label3)) +
		geom_text(aes(x=graph2$x4, y=graph2$y4, label=label4)) +
		#scale_color_manual(values=colorvec)+
		geom_smooth(method=lm, se=FALSE) + 
		stat_function(fun=unity, linetype="solid", color="darkgrey", size=2)+
		stat_function(fun=zero, linetype="solid", color="darkgrey", size=2) +
		theme_bw()+
		theme(axis.title=element_text(family= "serif", colour="black", size=12),
			axis.text=element_text(family= "serif", colour="black", size=12),legend.position="none")
	dev.new()
	return(m1)
}




showMod <- function(m,name) {
pred <- fitted(m)
res <- residuals(m)
o <- res + pred
Type1 <- rep("Residuals", length(res))
Type2 <- rep("Observed", length(res))
Data <- c(res, o)
Type <- c(Type1,Type2)
Predicted <- rep(pred,2)
Trt <- rep(names(residuals(m)),2)
PubIDRef <- p[c("TID", "PubID")]
q <- subset(PubIDRef, rownames(PubIDRef) %in% Trt)
q$Trt <- rownames(q)
unity <- function(x) {x}
s1 <- "Predicted"
s2 <- "Observed or Residual"
s3 <- name
xlabel <- paste(s1, s3, sep = " ")
ylabel <- paste(s2, s3, sep = " ")
xlab <- xlabel
#xlab <- "Predicted Values"
ylab <- ylabel
#ylab <- "Observed Values or Residuals"

graph <- data.frame(Data, Type, Predicted, Trt)
graph2 <- merge(graph, q, by="Trt")
m1 <- ggplot(graph2, aes(y=Data, x=Predicted, color=PubID, shape=Type))+xlab(xlab)+ylab(ylab)+geom_point()+geom_smooth(method=lm, se=FALSE)+stat_function(fun=unity, linetype="dotted", color="black", size=1)+theme(axis.title=element_text(family= "serif", colour="black", size=12),axis.text=element_text(family= "serif", colour="black", size=12),legend.title=element_text(family="serif"),legend.text=element_text(family="serif"), legend.position="top")
dev.new()
return(m1)
}

resc <- function(o,pred,f) {
if(is.factor(f)==FALSE) f <- as.factor(f)
res <- o-pred
graph <- data.frame(f,res)
m1 <- ggplot(graph, aes(y=res, x=pred, color=f, shape=f))+geom_point()+geom_smooth(method=lm, se=FALSE, aes(color=f))+theme(axis.title=element_text(family="serif", colour="black", size=12))
return(m1)
}