##
## plots.q
##
## Plots for Effect Measures Lecture
##
## 1-CI Plot
t	<- c(0:10)
IR1	<- 0.2
IR2	<- 0.1
Risk1	<- exp(-IR1*t)
Risk2	<- exp(-IR2*t)
time	<- c(0:10)
##
#
postscript("IncidenceRates.ps")
par(cex=1.2)
plot(time, Risk1, type="l", ylim=c(0,1), xlim=c(0,10), 
     ylab="Proportion Disease-Free (1-CI)",
     xlab="Time (years)", lty=1)
lines(time, Risk2, lty=2)
text(2.2,0.4,"IR=0.2 / person-year")
text(6.5,0.65,"IR=0.1 / person-year")
dev.off()

##Incidence / Prevalence Plots
#
postscript("incidence.ps")
par( cex=1.25)
time	<-  c(0:6)
Subject	<-  c(0:6)
plot(time,Subject, pch=".", axes=F)
abline(h=c(6,5,3))

rect(4,5.9,6.5,6.1, col=heat.colors(11))
rect(-1,3.9,6.5,4.1, col=heat.colors(11))
rect(2,2.9,4,3.1, col=heat.colors(11))
rect(2,1.9,3,2.1, col=heat.colors(11))

segments(-1,2, 2, 2)
segments(1,1, 2.5, 1)
segments(4,1, 5, 1)

rect(2.5,0.9,4,1.1, col=heat.colors(11))

text(1.5, .67, "With disease", adj=0)
rect(.5,.62,1.3,.72, col=heat.colors(11))

text(1.5, .33, "At risk of becoming disease", adj=0)
segments(.5,.32,1.3,.32)

axis(1, labels=c("Initiation", "Termination"), at=c(0,6),cex=3)
axis(2, labels=c(1,2,3,4,5,6), at=c(1,2,3,4,5,6))

box()
dev.off()

%%%%%%%%% Additive versus Multiplicative Effect Plots%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Multiplicative Effects Example
n	 <- 100
SDLogE	 <- 0.35
SDLogX	 <- 0.5
Delta	 <- 1
MeanLogX <- 2

LogE	 <- rnorm(n, 0, SDLogE) 
LogX	 <- rnorm(n, MeanLogX, SDLogX)
LogY	 <- LogX + Delta + LogE

postscript("multiplicative.ps")
par(mfrow=c(1,2), cex=1.5)
plot( (exp(LogX)+exp(LogY))/2, exp(LogY)-exp(LogX), xlab=quote((Y[2]+Y[1])/2),
      ylab=quote(Y[2]-Y[1]), xlim=c(10,30), ylim=c(10,30), pch="+")
plot(  (LogX + LogY)/2, LogY-LogX, xlab=quote((log(Y[2])+log(Y[1]))/2),
      ylab=quote(log(Y[2])-log(Y[1])), pch="+")
dev.off()     
  
## Additive Effects Example
#
n	<- 100
SDE	<- 0.4
SDX	<- 1
Delta	<- 1
MeanX	<- 2

E	<- rnorm(n, 0, SDE) 
X	<- rnorm(n, MeanX, SDX)
Y	<- X + Delta + E

postscript("additive.ps")
par(mfrow=c(1,1), cex=1.5)
plot(  (X + Y)/2, Y-X, xlab=quote((Y[2]+Y[1])/2),
      ylab=quote(Y[2]-Y[1]), pch="+")
dev.off()    
 
%%%%%%%%%%%%Correlation Plots%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


GenRho<-function(MeanX, SDX, SDE, a, b, c, Xpoint=NULL){
	E <- rnorm(200,0,SDE)
	X <- rnorm(200,MeanX, SDX)

	Y<-a+b*X + c*X^2+ E
	X<-c(X, Xpoint)
	Y<-c(Y, Xpoint)
	rho<-cor(X,Y)
	print(rho)
	plot(X,Y, #main=quote(paste(rho,"=", sep="")), 
     		xlim=range(X), ylim=range(X), pch="+")
}

postscript("Rhos.ps")
par(mfrow=c(2,3))

GenRho(0,1,.5,0,1,0)
GenRho(0,1,.5,0,0,0)
GenRho(0,1,.5,0,-1,0)
GenRho(0,1,.1,0,.25,0)
GenRho(0,1,.05,2,0,-2)
GenRho(0,1,1,0,0,0, Xpoint=20)
dev.off()