Over the past few weeks you have explored linear models and how to interpret model summary output. Again we have stepped up the complexity, now venturing into the world of non-linear models.
This practical focuses on fitting non-linear models to data with an emphasis on 3 important classes of functions that all budding biologists and environmental scientists should know
polynomials,
exponential models, and
logistic models.
A question before we begin:
What are some advantages and disadvantages of non-linear models as compared to polynomials?
Polynomials
Quadratic
\(y = \beta_0 + \beta_1 x + \beta_2 x^2\)
where the parameters are the y-intercept (b0), the linear component (b1) and the quadratic component (b2).
If b2 is negative then the shape of the function is convex upwards, i.e. y increases with x until reaches a peak and then y decreases.
It is easy to understand so it has been commonly used for modelling the response of yield to inputs such as fertiliser, seeding rates. This is despite much criticism for being unrealistic.
Limitations:
rate of increase to peak is same as rate of decrease past peak
does not level off as x becomes small or very large, y just keeps increasing or decreasing.
Compared to the quadratic model which has 1 turning point, a cubic model has 2 turning points.
Exercise 1: Interpreting polynomials
A study was performed to examine the soil properties that control the within-field variation in crop yield. The focus of this question is on soil pH which (among other things) controls the availability of nutrients to plants.
This exercise does not require you to read in any data, but rather focus on interpreting the model output and comparing the models.
The figure below shows the raw observations of yield plotted against pH with three models fitted; a linear model, quadratic polynomial and a cubic polynomial.
which line corresponds to which model?
based on the output from the 3 models below, which model fits the data best? Note: no hypothesis testing yet, just how well the model fits the data (R2).
Linear model:
Quadratic model:
Cubic model:
Use the R output to perform hypothesis testing to find the best model. Write out the hypotheses you are testing.
Exercise 2: Fitting polynomials in R
This exercise will use real data from a yield-fertiliser trial in Bedfordshire, United Kingdom.
First thing we can do is fit a linear model to the fertiliser data:
# create fertiliser and yield objectsfert <-c(0, 100, 170, 225)yield <-c(3.32, 5.23, 5.41, 5.02)# Fits a linear model and saves it to an object called lin.modlin.mod <-lm(yield ~ fert)# Summarises key features of modelsummary(lin.mod)
Fit and plot a quadratic polynomial in R. In R a quadratic polynomial can be fitted using the following lines of code:
# create a new variable which is the square of the fertilizer ratesfert2 <- fert^2# fit the quadratic model incorporating fert2quad.mod <-lm(yield ~ fert + fert2)summary(quad.mod)
What is the fit like for our quadratic model? is it better than our linear model?
In Excel it is easy to fit a line, by creating a scatterplot, then add Trendline… and selecting Polynomial (2nd order).
To fit our polynomial line in R, we need to obtain model predictions first.
To plot model predictions you first need to predict at fine intervals of the predictor to make a continuous plot that is not jagged or stepped. To create a new prediction dataset you can use the following code:
# creates a sequence of numbers from 0 to 225 going up in increments of 1new.fert <-seq(0, 225, 1)
We can use our model to predict at the values in the new prediction dataset, in this case new.fert.
The general form of the predict function is predict(model object, list object).
The list object tells R what object contains the data we will use to predict. For example in our case the model was built on fert and fert2 so we have to tell the predict function what object contains the new values for each of these, in our case new.fert.
Now we plot the raw observation as points and add an overlay of the model fit as lines:
plot(fert, yield, xlab ="Fertilizer", ylab ="Yield")lines(new.fert, new.pred) # Adds lines to original plot
Exponential function
\(y=y_0e^{kx}\)
where the parameters are y0 which is the multiplier which expresses the starting or final value, and k which is negative for exponential decay and positive for the exponential growth.
The half life (for decay) or doubling time (for growth) can be calculated as
\(\frac{log_e 2}{k}\)
Limitations:
harder to fit than polynomials
exponential growth has no horizontal asymptote; keeps going up.
Exercise 3: Initial estimates for exponential growth function
In this exercise we will find initial estimates of the parameters of an exponential growth model by visual assessment of plots of the data and/or rough calculations. The initial estimates of the parameters are needed as starting points for the iterative fitting methods we will use in the practicals, e.g. SOLVER in Excel and the nls() function in R.
The plot and table below presents the population of the world from 1650-1965.
We wish to model the data with an exponential growth function of the form;
\(y=y_0e^{kx}\)
where
y is the population in the year x,
y0 is the population in 1650 and
k is the rate constant.
Provide an initial estimate of y0.
The parameter k can be estimated from a linear model fitted to loge population against year.
Rather than formally fitting a linear model you could estimate the slope approximately by using the smallest and largest value to estimate the slope and therefore k.
For an exponential growth model the doubling time of a population can be estimated by loge2 /k.
Examine the graph and/or table to estimate the doubling time and use this to estimate k. You will have to make k the subject in the equation for estimating the doubling time.
How similar were the estimates of k?
Exercise 4 : Exponential growth models
This data is from Jenkins & Adams (2010) who studied soil respiration rates against temperature for different vegetation communities in the Snowy Mountains. They fitted an exponential growth model to the data.
The purpose of this exercise is to illustrate the dangers of using Excel’s in-built functions for statistics more complex than calculating means and fitting simple models.
Plot the data in Excel and using the Add Trendline… option. Make sure tick the option for displaying the equation in the graph.
The researchers performed the experiment up to a temperature of 40 degrees C, would you expect exponential growth in the respiration rate to continue if high temperatures were considered? Is there a better model?
Now fit the same model in R using the nls function. Code to get you started is:
Now you can fit a line to the plot. Does this look similar to your trendline in Excel?
#Plots raw dataplot(temp,respiration,xlab='Temperature',ylab='Respiration')#Creates new dataset for predictions ( 5 to 40 at an interval of 1)temp.new<-seq(5,40,1)#Makes predictions onto temp.newpred.exp<-predict(res.exp,list(temp=temp.new))#Adds model fit to existing plotslines (temp.new,pred.exp)
Compare the parameters values between Excel and R. You can extract the RSS value from an nls object by using the code below:
deviance(res.exp)
Calculate the RSS value for the Excel exponential model. Based on this, which is the better model?
When faced with the need to fit an exponential function, one approach that was used before computing power became readily accessible was to log the y values which linearises the relationship with x, enabling the modeller to use a simple linear model.
If we linearise, model would be \(log_e(y) = b_0 + b_1x\), where \(e^{b_0}\) is the y0 parameter in an exponential model, and b1 is the k parameter in the exponential model. This is similar to what was demonstrated in the lecture this week.
In Excel, log the soil respiration data and fit a linear model. You will see that the fitted model gives the same values as the exponential model fitted to the untransformed data.
If you compare the R2 values for both you will see they are the same. This means that Excel reports the R2 of the linear model fitted to logged respiration as the R2 of the exponential model fitted to the raw data. This is naughty of Excel.
For the dataset used here the exponential model fits it so well the Excel approach is only slightly different to the correct approach used in R.
In cases where the model does not fit the data so well the differences would be larger. Logarithm makes smaller values larger and larger values smaller.
WHY IS THIS SUB-OPTIMAL?
Regression modelling assumes that the residuals are normally distributed so logging normally distributed data will change the distribution to a non-normal one – it is best to analyse the data without transformation.
Modelling data on the logged scale reduces the impact that larger values have on minimising the RSS but when you plot the fitted model with the original data you may observe large discrepancies for larger values. In other words, using a linear model on the log(data) can result in a higher discrepancy for larger values when plotting the fitted model on the original data. Therefore, the model fitted to the logged data is not necessarily the best on the original data.
The reporting of the R-squared on the logged data on a model purported to be fitted to untransformed data is just wrong as the R-squared on the log scale will be better as the variation in the data has been reduced but we really want to know how well an exponential model fits the raw data.
Logistic function
There are several versions of the logistic function. We will use the following:
\[ y = \frac{Asym}{1+e^{\frac{xmid-x}{scal}}} \]
where
\(Asym\) is the maximum value of \(y\) (upper limit, horizontal asymptote).
\(xmid\) is the value of \(x\) when \(y\) is halfway between the lower and upper limits (inflection point, \(y = 0.5 \times Asym\)).
\(scal\) is the rate of change: the rate at which \(y\) approaches the upper limit.
Commonly used to model growth that has a sigmoid shape, i.e. where growth is initially slow, then picks up to a maximum, then slows down again as the system reaches a maximum.
Limitations:
Harder to fit than polynomials.
Exercise 5: Logistic models
In this exercise you will model the yield of pasture over time (since sowing).
Note, we will assume the yield at sowing = 0 is 0 t/ha, which allows us to use the equation above (that is pre-defined in SSlogis()).
If this were not the case, we would need to use a slightly different equation; \(y = y_0 + \frac{Asym}{1+e^{-scal(x-xmid)}}\), where \(y_0\) is the yield at sowing.
Fit the model in R using the nls function. Estimate the starting
parameters manually (using the previous exercise and the equation above as a guide). Then use the SSlogis() function to automatically estimate the parameters (see L12, Tut12 or the handbook).
Plot the fitted model with the observations.
Compare the final model parameters when we provide starting estimates, and when we use SSlogis(). Are they the same?
STARTING VALUES
When fitting non-linear functions (i.e. logistic or exponential) using iterative procedures such as nls or SOLVER the starting estimates of the parameters need to be approximated. If the values are too far, the model will not run and return an error.
The best way to ensure that you have suitable starting values is to plot the data with the predictions overlaid for your starting parameters. You can then see how close your initial model is to the data.
The reason we go through the process of estimating parameters is because nls and SOLVER can fit any nonlinear equation. They are more versatile
Realistically, most nonlinear relationships you would fit are covered by a self-starting function; i.e. SSlogis(), SSasymp(), or nlraa::SSexpf(). These estimate the parameters for you and are more efficient to run. We recommend using these when suitable.
That’s it for Module 3! Great work exploring non-linear models!
Thank you all (students and demonstrators!) for your hard work and enthusiasm throughout this Module. Good luck with Project 3 and the final exam!