# How to fit regression to custom model in R

This is for my honors thesis. I have a large data set, of which I'm sharing only what I call the "Low phosphorus" series:

> P0
R   N P D.weight
1  r1   0 0     63.8
2  r2   0 0     34.2
3  r3   0 0     24.9
4  r4   0 0     30.4
5  r5   0 0     33.3
6  r1  45 0     24.5
7  r2  45 0     20.1
8  r3  45 0     23.7
9  r4  45 0     20.0
10 r5  45 0     66.8
11 r1  90 0     27.8
12 r2  90 0     17.2
13 r3  90 0     36.4
14 r4  90 0     33.5
15 r5  90 0     14.0
16 r1 180 0     20.6
17 r2 180 0      9.7
18 r3 180 0      8.8
19 r4 180 0     14.4
20 r5 180 0     21.6
21 r1 360 0     18.4
22 r2 360 0      8.9
23 r3 360 0     31.4
24 r4 360 0     13.3
25 r5 360 0     21.9

• R is rep
• N is nitrogen applied to the soil
• P is phosphorus applied to the soil
• D.weight is average dry weight of plants in grams

The way to visualize these data is to put N on the x-axis and dry weight on the y axis: I have to make a nonlinear regression of these data, but I don't want to fit it to a quadratic model; instead, I wanna fit it to the equation below (an alternative to the Mitscherlich equation):

$Y = a - b \times\exp(-cx)$

• Y is dry weight
• a is a fitted parameter representing the maximum biomass
• b is a fitted parameter representing the initial level of the added nutrient in the soil
• c is a fitted parameter representing the rate of increase in biomass with increasing nutrient amendment
• x is, in this case, nitrogen level

The problem is, I just don't know how to code for this. I have been going crazy trying to find out how to "tell" R that I wanna use that equation for my regression, and not $Y = ax + b$ (like in a linear regression), or $Y = ax^2 + bx + c$ like in a quadratic regression, etc.

• Looks like you are seeking academic statistical consultation on how to do nonlinear modeling with repeated measures. That's not part of the goals of SO. (It's also unclear how you plan to use the N and P values.) – DWin Oct 17 '14 at 22:06
• You need to consult with a statistician, not a programmer. Try Cross Validated instead; that's the proper venue for such questions. What type of model you fit is mostly dependent on the hypothesis you want to test. – MrFlick Oct 17 '14 at 22:44
• Hello all. I edited the question, sorry if I wasn't clear the first time around. Maybe have another look? Thanks. – XGF Oct 18 '14 at 1:51
• This is illustrated in many places here on CV. Among my answers I can point to stats.stackexchange.com/a/70184, which uses a similar function (easily modified to yours) and makes the fit with a half dozen lines of code; and stats.stackexchange.com/a/64039, which fits a generalized linear model. There are also quite a few threads dealing with exactly your exponential model (but they are hard to find in a search). But as previous comments suggest, there's much more to this analysis than just fitting a curve to data. Look for statistical support at your institution if you can. – whuber Oct 18 '14 at 4:32

You could fit a nonlinear regression. This would be suitable if - aside from the nonlinear relationship between $E(Y)$ and $x$ - your model assumptions were similar to ordinary regression (independence and constant variance, for example).

In R, see ?nls.

The difficulty with this particular model is finding suitable starting values. With a reparameterization, however, you might be able to get it into the form of one of the available self-start functions and save some trouble there (specifically, I believe SSasymp is the self start function for the reparameterized model you need). However, I managed to find reasonable enough start values and got convergence:

 nlsfit <- nls(D.weight ~ a - b * exp(-c*N) ,P0,start=list(a=10,b=-20,c=.05))

summary(nlsfit)

Formula: D.weight ~ a - b * exp(-c * N)

Parameters:
Estimate Std. Error t value Pr(>|t|)
a  16.208572   6.222312   2.605  0.01617 *
b -22.000400   7.552922  -2.913  0.00806 **
c   0.011082   0.009454   1.172  0.25364
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.63 on 22 degrees of freedom

Number of iterations to convergence: 11
Achieved convergence tolerance: 8.554e-06


It seems to fit tolerably well, though the constant variance assumption may be suspect: There's also some suggestion of the possibility of right skewness.

• Thank you, @Glen_b! nls() is exactly what I was looking for. Anybody: ideas for a code to get good initial values? Using selfStart, like @Glen_b suggested, or any other method. I'm trying by myself as well, but any help is welcome. – XGF Oct 18 '14 at 15:19
• Actually, I figured it out! It took me a while but I got there eventually :) – XGF Oct 20 '14 at 18:41

@Glen_b provides the right path you should pursue. Here are some issues I would like to point out for this nonlinear modeling scenario.

First, the nonlinear fit @Glen_b showed has an underlying assumption that each observation, either at the same N or not, are independent. I noticed that you have a variable "R" which is the rep. I am not sure what it means, but it seems to me that there are 5 R's, i.e., r1, r2, r3,r4, and r5. So I am concerned that the data from each rep might be correlated (you might want to provide more information about your variables). If they're correlated, you might want to use the pooled mean at each nitrogen level N to fit the nonlinear least squares regression model. Otherwise, you should get similar point estimates, but the estimates of the standard errors are different.

Second, as for the initial values, here are the tricks you can use. Note that $a$ is the weight value when $x$ goes to positive infinity. So you can use the pooled mean of weight data from the largest $x$ (i.e., nitrogen level) at hand for $a$, since this would be the closest guess. Similarly, when $x=0, y = a - b$, thus you can then use the pooled mean of weight at $x=0$ to obtain the initial estimate for $a-b$. Thus you can figure out the initial guess for $b$. After that, you can choose data at an arbitrary $x$, and substitute $a, b$ to solve for the initial estimate $c$.

Third, if the data from each rep are truly correlated. You might want to first build a nonlinear random-effects model, assuming that at different level of R, the nonlinear curves are different, i.e., the $a, b, c$ are different. After you build the model, you can test whether each of the three random coefficients are significant or not. For example, random $a$ is essentially random intercept. If there is just random intercept, then that means all the nonlinear curves are the same shape, but just shift up and down by random $a$. I've not used R to do nonlinear mixed-effects model in R. But in SAS, the PROC NLMIXED exactly does this.

I assumed that since you are not statistician I will give you a simple solution.

Install SPSS and go to regression >> non linear • I tried since your data is not fully attached, so that the model is not really non linier, if you assumed the graph relationship between N vs DW is non linear such mitscherlich you can go by using steps above, $b=max$ $a=min$ $c=any coefficient$ so this like try and error but the final result is the iteration of the approximation parameter. – valerie Oct 18 '14 at 3:35
• I don't think 'Install SPSS' counts as a simple solution. – Glen Oct 18 '14 at 5:28