Compare manually estimated Cv, Cb ~ gnls estimations

Mid-point of H2O2-dependent viability can be consistently estimated manually and by gnls. 

These are the fitting results of H2O2-LOH data.

Comparison of manually fitted Cb ~ gnls results. Red line is y=x, blue line is the model of gnls_Cb = a*manualCb + b. 

Gnls seem to give larger Cb than manual estimations. So, I should use manual estimations as initial values for 'gnls' and see whether it helps. Function 'gnls' treats all points equally correct, but I know data points in higher H2O2 concentration are more noisy due to fewer colonies. So, manual fit has some advantage, though 'gnls' can be viewed as kind of 'objective' method (well, if initial values are not considered as human biases).



# I am concerned how initial values can influence Cb estimation.
# The black colonies data are very noisy, so initial values can trapped fitting to different local optima.

rm(list=ls())
setwd("~/github/LOH_H2O2_2012-master/analysis")
files = list.files(path="../", pattern="_*csv")

tb29 = read.csv("../_ctl.tb_out.20130529.csv")
tb30 = read.csv("../_ctl.tb_out.20130530a.csv")
tb30$files = as.character(tb30$files)

head(tb29)
head(tb30)
tb30$change = tb30$Cb / tb29$Cb
summary(tb30)

tb30$color = ifelse(tb30$change>1.5 | tb30$change<0.5, "red", "black"   )
tb30$pch = ifelse(tb30$change>1.5 | tb30$change<0.5, 19, 1   )
plot( tb29$Cb ~ tb30$Cb, xlim=c(0,0.4),pch=tb30$pch, col=tb30$color)

tb30$label = ifelse(tb30$change>2 | tb30$change<0.5, as.character(tb30$files), NA   )
text( tb30$Cb, tb29$Cb, tb30$label)

red= tb30[tb30$color=='red', ]
head(tb30)
tb30$Cv.vs.Cb = tb30$Cv / tb30$Cb

######load the manual fitting results
tbManual = read.csv("H2O2_Log_Plot_Summarized_data,2013May30.csv")
head(tbManual)
names(tbManual)=c("files", "Date", "strains", "Cv", "Cb", "OD600", "note")
tbManual$files = as.character(tbManual$files)

tbManual$Cv.vs.Cb = tbManual$Cv / tbManual$Cb

summary(tbManual$Cv.vs.Cb)
summary(tb30$Cv.vs.Cb)
tbManual[tbManual$Cv.vs.Cb>5, ]
str(tbManual)

### merge
intersect( tbManual$files, tb30$files )
head(tbManual)
head(tb30)

tmp = tbManual[match(tb30$files, tbManual$files), ]
cbind( tmp[,1], tb30[,1]) #visual check, passed. 

tb30M = data.frame( cbind( tb30, tmp[, c(2:9)]) )
head(tb30M)
plot( Cv.vs.Cb ~ Cv.vs.Cb.1, data=tb30M)

plot( Cv ~ Cv.1, data=tb30M, xlab="manual fit Cv", ylab="gnls Cv")

plot( Cb ~ Cb.1, data=tb30M, xlab= "manual Cb", ylab="gnls Cb")
x = seq(0,1, by=0.1); y = x; 
lines( y ~ x, col='red')
m = lm( Cb ~ Cb.1, data=tb30M)
summary(m)
abline(m, col='blue')