load(file='../../BookCode/nlrwr/data/M.merluccius.rda') # Plot the data to get an idea of what it looks like: # plot(num.fish ~ spawn.biomass, data=M.merluccius, xlab='spawn', ylab='number of fish') grid() # Define the fitted function we will use: # BHfct = function(S, alpha, k){ (alpha * S)/(1 + S/k) } # Guess at some initial values for the parameters (alpha, k): # alpha = 40/20 # close to the slope at S=0 k = with(M.merluccius, max(num.fish))/alpha # limit_{S -> inf} f(S;(alpha, k)) = alpha k fit = nls(num.fish ~ BHfct(spawn.biomass, alpha, k), data=M.merluccius, start=c(alpha=alpha, k=k)) print(summary(fit)) if( !require('investr') ){ install.packages('investr', dependencies=TRUE, repos='http://cran.rstudio.com/') } library(investr) plotFit(fit, interval='confidence', xlab='spawn', ylab='number fish') grid() # Lets compare the two suggested models on this data: RSS_m1 = sum( residuals(fit)^2 ) RSS_m2 = sum( (M.merluccius$num.fish - BHfct(M.merluccius$spawn.biomass, 4.91, 45.39))^2 ) print(sprintf('RSS_m1= %f; RSS_m2= %f', RSS_m1, RSS_m2)) # Plot the two models: # plot(num.fish ~ spawn.biomass, data=M.merluccius, xlab='spawn', ylab='number of fish') bm_grid = with(M.merluccius, seq(min(spawn.biomass), max(spawn.biomass), length.out=100)) nfish_m1 = predict(fit, data.frame(spawn.biomass=bm_grid)) lines(bm_grid, nfish_m1, col='blue') nfish_m2 = BHfct(bm_grid, 4.91, 45.39) lines(bm_grid, nfish_m2, col='red') legend_str = c( sprintf('nls fit: RSS=%.1f', RSS_m1), sprintf('Cadima fit: RSS=%.1f', RSS_m2) ) legend('topleft', legend_str, lty=c(1, 1), col=c('blue','red')) grid()