How to fit regression to custom model in R

Question

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.

Answer 1

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.

If you search for "nonlinear least squares" or "nls" you should turn up a number of posts, some with useful advice. I agree with the comments you've already received about seeking advice.

Answer 2

@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.

Answer 3

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

Install SPSS and go to regression >> non linear