if( !require('drc') ){
    install.packages('drc', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(drc)

# The example from the book:
#
RScompetition.book = nls(biomass ~ a/(1 + b * (x + c * z)), data=RScompetition, start = list(a=20, b=1, c=1))

# Try to fit the suggested models:
#
#RScompetition.m3 = nls(biomass ~ a/(1 + b * (x + c * z)), data=RScompetition, start = list(a=20, b=0, c=1))
#RScompetition.m3 = nls(biomass ~ a/(1 + b * (x + c * z)), data=RScompetition, start = list(a=200, b=-1, c=100))
RScompetition.m3 = nls(biomass ~ a/(1 + b * (x + c * z)), data=RScompetition, start = list(a=20, b=-1, c=1))


# Compare the two models:
#
virDensity.book = with(RScompetition, x + coef(RScompetition.book)[3] * z)
virDensitygrid.book = seq(0, max(virDensity.book), length.out=100)
biomassVal.book = predict(RScompetition.book, newdata=data.frame(x=virDensitygrid.book, z=0))

virDensity.m3 = with(RScompetition, x + coef(RScompetition.m3)[3]*2)
virDensitygrid.m3 = seq(0, max(virDensity.m3), length.out=100)
biomassVal.m3 = predict(RScompetition.m3, newdata=data.frame(x=virDensitygrid.m3, z=0))

#postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_4.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1, 2))
plot(biomassVal.book ~ virDensitygrid.book, type='l', xlab='virtual density', ylab='biomass', main='book')
with(RScompetition, points(biomass ~ virDensity.book))
grid()

plot(biomassVal.m3 ~ virDensitygrid.m3, type='l', xlab='virtual density', ylab='biomass', main='m3')
with(RScompetition, points(biomass ~ virDensity.m3))
grid()
par(mfrow=c(1, 1))
#dev.off()