# 
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

save_plots = F

#
# EPage 140
#

DF = read.csv( "../../Data/Puromycin.csv", header=TRUE )

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter6/soar_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot( DF$Initial_Concentration, DF$Rate, col=DF$Puromycin )
if( save_plots ){ dev.off() }

m1 = lm( Rate ~ Initial_Concentration, data=DF )
m2 = lm( log(Rate) ~ Initial_Concentration, data=DF )
m3 = lm( sqrt(Rate) ~ Initial_Concentration, data=DF )

# Plot each model:
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter6/soar_model_plots.eps",onefile=FALSE, horizontal=FALSE) }
plot( DF$Initial_Concentration, DF$Rate, col=DF$Puromycin )
abline( m1 )
lines( DF$Initial_Concentration, exp(predict(m2, newdata=DF)), type='l', col='green' )
lines( DF$Initial_Concentration, predict(m3, newdata=DF)^2, type='l', col='blue' ) 
legend( x=0.8, y=100, legend=c("linear", "logarithmic", "sqrt"), col=c("black","green","blue"), lty=c(1,1,1) )
if( save_plots ){ dev.off() }

# A square root transformation of the "y=concentration" looks better than the above two transforamtions:
#
m4 = lm( Rate ~ sqrt(Initial_Concentration), data=DF )
lines( DF$Initial_Concentration, predict(m4, newdata=DF), type='l', col='purple' )

# How does puromycin affect the reaction:
#
m5 = lm( Rate ~ Initial_Concentration + Puromycin, data=DF )
summary(m5)

m6 = lm( Rate ~ sqrt(Initial_Concentration) + Puromycin, data=DF )
summary(m6)

# Generate some data to predict rate from:
#
ics = seq( min(DF$Initial_Concentration), max(DF$Initial_Concentration), length.out=100 )
ndf_T = data.frame( Initial_Concentration=ics, Puromycin=rep( 'True' ) )
ndf_F = data.frame( Initial_Concentration=ics, Puromycin=rep( 'False' ) )

# Make some predictions of the models built so far:
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter6/soar_puromycin_models.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
plot( DF$Initial_Concentration, DF$Rate, col=DF$Puromycin )
abline( m1, col='gray' ) # the no puromycin model
lines( ics, predict(m5, newdata=ndf_T), type='l', col='red' ) # linear model with puromycin
lines( ics, predict(m5, newdata=ndf_F), type='l', col='black' )
legend( x=0.2, y=75, legend=c("global linear", "linear(Puro=T)", "linear (Puro=F)"), col=c("gray","red","black"), lty=c(1,1,1) )
plot( DF$Initial_Concentration, DF$Rate, col=DF$Puro )
abline( m1, col='gray' ) # the no puromycin model
lines( ics, predict(m6, newdata=ndf_T), type='l', col='red' ) # with puromycin and sqrt(Initial_Concentration)
lines( ics, predict(m6, newdata=ndf_F), type='l', col='black' )
legend( x=0.2, y=75, legend=c("global linear", "sqrt (Puro=T)","sqrt (Puro=F)"), col=c("gray","red","black"), lty=c(1,1,1) )
if( save_plots ){ dev.off() }