#=================================================================
# Figure 18.3.R.
# R script to calculate nonlinear least squares estimates for
# parameters b1 (maximum feeding rate) and b2 (half saturation
# constant) in the rectangular hyperbolic equation for feeding
# rate (Holling type 2 functional response, Monod nutrient uptake
# rate, Michaelis-Menten enzyme equation).  The equation is defined
# by
#           b1*prey
#  rate = ----------   .
#          b2 + prey
#
# Here 0<b1, 0<b2, and "prey" is the density, concentration, or
# abundance of prey, substrate, or nutrient.
#
# Data in example are moose density (# per 1000 km^2) and number
# killed per wolf in 100 d from:
# Messier, F.  1994.  Ungulate Population Models with Predation: A
# Case Study with the North American Moose.  Ecology 75:478-488.
#=================================================================

#-----------------------------------------------------------------
# 1.  Input the data.  Text file of data has two columns, one labeled
# "prey" and one labeled "rate".
#-----------------------------------------------------------------
Wolf.Moose=read.table("wolf_moose.txt",header=TRUE) # Change file name
                                                    #  if necessary.
attach(Wolf.Moose)

#-----------------------------------------------------------------
#  2.  Calculate initial values using a linearization transform.  The
#  transform allows initial values to be calculated with a multiple
#  regression without intercept: 
#                rate*prey = b1*prey - b2*rate.
#-----------------------------------------------------------------
yy=rate*prey;
xx=cbind(prey,rate);
bb=solve(t(xx)%*%xx,t(xx)%*%t(t(yy)));
b1.0=bb[1];
b2.0=-bb[2];

#-----------------------------------------------------------------
# 3.  Use nls( ) function to minimize sum of squares with an
# iterative numerial method.  The object "curve.fit" will be a list
# of results from the calculations.
#-----------------------------------------------------------------
curve.fit=nls(rate~b1*prey/(b2+prey),data=Wolf.Moose,
  start=list(b1=b1.0,b2=b2.0))

#-----------------------------------------------------------------
# 4.  Print the results of the calculations.  Assign values to b1
# and b2 for plotting the fitted model.
#----------------------------------------------------------------- 
summary(curve.fit)
b=coef(curve.fit)
b1=b[1]
b2=b[2]

#-----------------------------------------------------------------
# 5.  Calculate the fitted rate curve for a range of values of prey,
# and store the values in "prey.vals".
# Range is from 0 to slightly beyond max(prey).  Change range if 
# desired.  Values of fitted rate curve are in "fitted.curve".
#-----------------------------------------------------------------
prey.vals=(0:100)*1.01*max(prey)/100;
fitted.curve=prey.vals*b1/(b2+prey.vals);

#-----------------------------------------------------------------
# 6.  Plot the data in a scatterplot.  Overlay the fitted rate
# equation.
#-----------------------------------------------------------------
plot(prey,rate,ylim=c(0,4),xlim=c(0,2.5));
points(prey.vals,fitted.curve,type="l");

#-----------------------------------------------------------------
# 7.  Tidy up a bit.
#-----------------------------------------------------------------
detach(Wolf.Moose)