8
$\begingroup$

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:

Plot

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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ 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.) $\endgroup$ Commented Oct 17, 2014 at 22:06
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Oct 17, 2014 at 22:44
  • $\begingroup$ Hello all. I edited the question, sorry if I wasn't clear the first time around. Maybe have another look? Thanks. $\endgroup$ Commented Oct 18, 2014 at 1:51
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Oct 18, 2014 at 4:32

3 Answers 3

Reset to default
10
$\begingroup$

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:

enter image description here

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Commented Oct 18, 2014 at 15:19
  • $\begingroup$ Actually, I figured it out! It took me a while but I got there eventually :) $\endgroup$ Commented Oct 20, 2014 at 18:41
2
$\begingroup$

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

Improve this answer
$\endgroup$
0
-1
$\begingroup$

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

Install SPSS and go to regression >> non linear

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Commented Oct 18, 2014 at 3:35
  • 10
    $\begingroup$ I don't think 'Install SPSS' counts as a simple solution. $\endgroup$ Commented Oct 18, 2014 at 5:28

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.