# Goal: To do nonlinear regression, in three ways # By just supplying the function to be fit, # By also supplying the analytical derivatives, and # By having him analytically differentiate the function to be fit. # # John Fox has a book "An R and S+ companion to applied regression" # (abbreviated CAR). # An appendix associated with this book, titled # "Nonlinear regression and NLS" # is up on the web, and I strongly recommend that you go read it. # # This file is essentially from there (I have made slight changes). # First take some data - from the CAR book -- library(car) data(US.pop) attach(US.pop) plot(year, population, type="l", col="blue") # So you see, we have a time-series of the US population. We want to # fit a nonlinear model to it. library(stats) # Contains nonlinear regression time <- 0:20 pop.mod <- nls(population ~ beta1/(1 + exp(beta2 + beta3*time)), start=list(beta1=350, beta2=4.5, beta3=-0.3), trace=TRUE) # You just write in the formula that you want to fit, and supply # starting values. "trace=TRUE" makes him show iterations go by. summary(pop.mod) # Add in predicted values into the plot lines(year, fitted.values(pop.mod), lwd=3, col="red") # Look at residuals plot(year, residuals(pop.mod), type="b") abline(h=0, lty=2) # Using analytical derivatives: model <- function(beta1, beta2, beta3, time) { m <- beta1/(1+exp(beta2+beta3*time)) term <- exp(beta2 + beta3*time) gradient <- cbind((1+term)^-1, -beta1*(1+term)^-2 * term, -beta1*(1+term)^-2 * term * time) attr(m, 'gradient') <- gradient return(m) } summary(nls(population ~ model(beta1, beta2, beta3, time), start=list(beta1=350, beta2=4.5, beta3=-0.3))) # Using analytical derivatives, using automatic differentiation (!!!): model <- deriv(~ beta1/(1 + exp(beta2+beta3*time)), # rhs of model c('beta1', 'beta2', 'beta3'), # parameter names function(beta1, beta2, beta3, time){} # arguments for result ) summary(nls(population ~ model(beta1, beta2, beta3, time), start=list(beta1=350, beta2=4.5, beta3=-0.3)))