4
$\begingroup$

I have measured growth rates over a range of temperatures (temperature response curve) and would like to fit an already established equation/model to it. I'm very new to R and have trouble coding it right.

This is my data:

    x <- c(37.8,
   34.8,
   32.25,
   29.65,
   26.75,
   24.4,
   21.45,
   19.15,
   16.2,
   13.7,
   10.95,
   8.2,
   5.4)
    y <- c(0.0381,
   0.018133333,
   1.057033333,
   2.3395,
   2.5188,
   2.6132,
   2.398466667,
   2.082966667,
   1.687366667,
   1.389233333,
   0.9835,
   0.639833333,
   0.3757)

The equation I'm trying to fit is as follows:

$y=(3.5-x/10)(x/25)^{(5/2)}$

There is also another data set:

x2 <- c(34.8, 32.25, 29.65, 26.75, 24.4, 21.45, 19.15, 16.2, 5.4)

and

y2 <- c(0.033966667, 1.5153, 2.226266667, 2.3462, 2.096766667, 2.0986, 1.941866667, 1.5183, 0.2717)

I tried modifing code that I found here to my equation but have been unsuccessful so far. Can someone help with this with a step by step guide for dummies?

Updated model for my dataset:

$$y = Y_{max}\frac{T_{max}-x}{T_{max}-T_{opt}}\left(\frac{x}{T_{opt}}\right)^{(T_{opt}/(T_{max}-T_{opt}))}$$

$\endgroup$
  • 6
    $\begingroup$ Your function doesn't appear to contain any unknowns. In what sense is there an "fitting" going on? Also ... when you write the "(x/25)exp(5/2)" part do you actually mean (x/25)^(5/2)? That looks a lot more like the shape of your data (but the entire function is only about 37% of what it should be). Even with that change the function isn't quite right, but it's sort of closeish. $\endgroup$ – Glen_b Sep 4 '15 at 18:28
  • $\begingroup$ Hi Glen, fitting is the wrong word I guess then. I basically want to see how well my data fits to the equation. There are some unknowns in the original equation that can be filled with measurements. The following equation can be used too: y=((35−x/(35-z))(x/z)^(z/(10)) $\endgroup$ – Kay Sep 4 '15 at 20:23
  • $\begingroup$ With the previous version of your "updated model," y can never reach a value greater than 1. The peak for y occurs at $x=T_{opt}$, and at that value of x, y=1. I edited your question to add an additional parameter $Y_{max}$ to cover this. I also put the additional data set (which I now use for examining differences between parameter values) into your question directly, as comments (unlike questions) can be lost. Make sure that you are checking plots by eye and that you have reasonable fits before you start worrying about differences between fitted parameters from different experiments. $\endgroup$ – EdM Sep 27 '15 at 17:51
19
$\begingroup$

Fitting data to an equation without free parameters

First, let's clarify what it means "to see how well my data fits to the equation" when you have a fixed equation as in your original question. You have data on $x$ and $y$, and your equation: $$y=(3.5-(x/10))(x/25)^{5/2}$$ has no free parameters. Here's a plot; data points are circles, the equation a solid curve:

original curve and data

So to answer your original question, your data don't fit that equation very well.

Fitting parameters of an equation to match data

Second, "fitting an equation" generally means having a general form of a relation between $y$ and $x$ in mind, but the details of the form depend on some additional unknown parameters of the equation. "Fitting" then means using the data to estimate the values of the parameters that best match the data, say by minimizing the sum of the squares of the discrepancies between the curve and the data. That is what's implicit in the form of the equation in your comment. (I corrected an apparently missing parenthesis so that it matches your original parameter-free equation at $z=25$): $$y = \frac{35-x}{35-z}\left(\frac{x}{z}\right)^{(z/10)}$$ Now there is something to "fit": find the value of $z$ that comes closest to your data. The simplest R function to do this is nls. (The R code I used is at the end of the answer.) The best fit to your $x$ and $y$ data is provided by a $z$ value of 33.1, not the value of 25 that would match your first parameter-free equation. It comes closer, but doesn't do very well, either, substantially below most of the true values:

fit z only

It looks like your equation might have an even more general form, like: $$y = \frac{a-x}{a-z}\left(\frac{x}{z}\right)^{(z/b)}$$ where $z$, $a$ and $b$ are parameters whose values might be estimated from the data. With that more general form, you can come fairly close to your data, at $z=34.0$, $a=36.7$, and $b=17.8$: enter image description here

Dangers in fitting equations

Third, you must be very careful if you head down this fitting path. As John von Neumann allegedly said, "With four parameters I can fit an elephant, and with five I can make him wiggle his trunk." Overfitting is a serious temptation when you start playing with data.

In a case like yours you have to think about what the parameters might mean in the real world. It seems that there are competing processes that lead to a peak in growth rate ($y$) at a certain temperature ($x$), and if the parameters $z$, $a$ and $b$ have definite meanings in terms of those competing processes then you might be justified in doing this type of fit. You might, say, learn something about how those processes differ from one strain of bacteria to another. But if such parameters (or the fixed values in your original equation) were simply chosen by someone to make some earlier set of data fit a particular curve then you are headed for big trouble.

Comparing fits to different data sets

If you want to compare the parameter values fitted to two different data sets, it's best to start by examining whether the overall fits differ between the data sets. There is a "useful and rather easily overlooked" way to do this with nls as Peter Dalgaard put it in a recent answer on R-help. (It can be found at the very bottom of the nls help page.)

Instead of working with separate vectors of $x$ and $y$ values from different experiments, construct a pooled data frame containing all the data. Each row should contain an $x$ value, the associated observed $y$ value, and a factor variable indicating the specific experiment. (Make sure the experiments are coded as factors, not as numbers.)

nls then allows you to fit different sets of parameter values to different experiments. Instead of specifying general parameters like Topt in your call, you can subscript the parameters by experiment, e.g. Topt[expt], where "expt" is the name of the factor distinguishing experiments.

Thus you can use anova to compare a fit to all the data, ignoring individual experiments, against a fit that allows different parameter values for different experiments. If these fits aren't statistically distinguishable, then you shouldn't be looking at differences of individual fitted parameter values among experiments. There's an example below where I present R code.

If the fits do differ, then you can consider comparing fitted parameter values. In the limit of very large samples the fitted values should have normal distributions, but it's not clear what "large" means for nonlinear fits and the meaning might differ depending on the equation. Thus there will be questions about the validity of p-values.

If you are willing to assume normality, nls reports standard errors for the parameter values, which you can use to obtain p-values under that assumption. A z-test can provide a rough overview, but it would probably be better to use the approach of Welch's t-test, which allows for different standard errors and degrees of freedom between two parameter estimates. To apply that approach, note that terms like $\frac{s_i^2}{N_i}$ in the formula for Welch's test are equivalent to $SE_i^2$ in the nls output.

R code

In terms of R code, it was simplest to define a general function for your temperature-response curve:

trcFunc <- function(x,z,a,b){((a-x)/(a-z))*((x/z)^(z/b))}

then give specific values for the parameters that matched your original equation, to compare your data to the curve:

plot(x,y)
curve(trcFunc(x,25,35,10),add=TRUE)

For fitting, I used the nls function, for example for $z$ alone

Zfit <- nls(y~trcFunc(x,z,35,10),start=list(z=30))

or for all 3 parameter values

ZABfit <- nls(y~trcFunc(x,z,a,b),start=list(z=30,a=35,b=10))

Get the fitted parameter values with

ZABfit
Nonlinear regression model
  model: y ~ trcFunc(x, z, a, b)
   data: parent.frame()
    z     a     b 
33.98 36.74 17.84 
 residual sum-of-squares: 1.228

Number of iterations to convergence: 5 
Achieved convergence tolerance: 1.35e-06

and compare data and curve:

plot(x,y,main="fit z, a, b")
curve(trcFunc(x,33.98,36.74,17.84),add=TRUE)

New form for relation; comparing fits With the new general form for the relation of $y$ to $x$ in an edit to the OP, a new general function for fitting is defined as:

newFunc <- function(x,Topt,Tmax,Ymax){Ymax*((Tmax-x)/(Tmax-Topt))*(x/Topt)^(Topt/(Tmax-Topt))}

For comparing two experiments, take expt1 as the data at the beginning of the question and expt2 as the second data set (x2,y2) toward the end, and construct a pooled data frame as suggested above. Then the fit ignoring the differences between experiments is:

fitPool <- nls(y~newFunc(x,Topt,Tmax,Ymax),data=newDF,start=list(Topt=25,Tmax=35,Ymax=2.5))

and that allowing separate parameter values for the two experiments is:

fitSeparate <- nls(y~newFunc(x,Topt[expt],Tmax[expt],Ymax[expt]),data=newDF,start=list(Topt=c(25,25),Tmax=c(35,35),Ymax=c(2.5,2.5)))

The pooled and separate fits are not significantly different:

> anova(fitPool,fitSeparate)
Analysis of Variance Table
Model 1: y ~ newFunc(x, Topt, Tmax, Ymax)
Model 2: y ~ newFunc(x, Topt[expt], Tmax[expt], Ymax[expt])
  Res.Df Res.Sum Sq Df  Sum Sq F value Pr(>F)
1     19     1.9215                          
2     16     1.4410  3 0.48046  1.7782 0.1918

Warning on choice of starting values

Note that, as can happen in non-linear curve fitting, what you get as the "best" fit can depend on what you provide as an initial guess for the parameter values. With your data, if I provided an initial guess of 26 for $z$ and fixed values of 35 and 10 for $a$ and $b$:

Zfit26start <- nls(y~trcFunc(x,z,35,10),start=list(z=26))

the value returned for $z$ is 13.1,which doesn't fit very well at all:

fit with different initial parameter choice

This happened because the program wasn't set up to evaluate an entire range of $z$ values and got stuck in what's called a "local optimum" when it started searching from $z=26$. If I tried starting values of $z$ from 20 to 25 then the nls program returned an error:

Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

while starting values of 19 or less gave the poorly fitting estimate of 13.1 for $z$. With a non-linear curve fit you can't count on getting the "global" optimum without doing a lot of exploration.

$\endgroup$
  • $\begingroup$ Hi Ed, thanks a lot for your detailed answer and comments! This helps me a lot! $\endgroup$ – Kay Sep 5 '15 at 19:24
  • $\begingroup$ @Kay I added an additional caution about choice of initial parameter values in non-linear curve fitting that was illustrated nicely by your data, and corrected a couple of errors in my earlier post. $\endgroup$ – EdM Sep 6 '15 at 15:14
  • $\begingroup$ Hi Ed, thanks a lot for following up on that. I've noticed the same trends when I change the starting values. I'm comparing 3 different temperature settings and I can make an educated guess for the starting values to not run into the mentioned problems. There is a third parameter that can be added to the equation which gives an even better fit, but as you pointed out might not have a biological meaning. I will try and compare the equation I have with each other and see which one has the best fit AND most biological meaning. $\endgroup$ – Kay Sep 7 '15 at 12:24
  • $\begingroup$ I also have this dataset:x <- c(34.8, 32.25, 29.65, 26.75, 24.4, 21.45, 19.15, 16.2, 5.4) y <- c(0.033966667, 1.5153, 2.226266667, 2.3462, 2.096766667, 2.0986, 1.941866667, 1.5183, 0.2717) Which has two optima and gets me into trouble as the general equation does not seem to work very well on my end. $\endgroup$ – Kay Sep 7 '15 at 12:29
  • $\begingroup$ The choice of starting parameter values can be crucial. I had no problem with your second dataset when I used starting values of $z=30,a=35,b=10$, which ended up being very close to the optimum values ($z=33.1,a=35.1,b=12.3$). The instability may come from the way $z$ appears in your equation: increasing $z$ makes the first factor blow up (to infinity at $z=a$) and although it decreases the second factor it raises the exponent of the second factor. $\endgroup$ – EdM Sep 7 '15 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.