0

I have data on the photosynthetic activity of microalgae under increasing light. PAR: is the different light exposures, (x, independent variable) ETR: is the photosynthetic activity (y, dependent variable). 

PAR ETR
1   0.409090909
46  18.81818182
126 51.54545455
234 93.6
404 173.1428571
656 246
816 272

i was able to plot the point for a sample and get the curve. But i want to fit these points to a specific non-linear model using this formula: ETR= PAR/((aPAR^2)+(bPAR)+c) I used the Nls function as part of the mosaic package that uses the nls function as shown in the code:

TABLE=read.table("sample 2.txt", header = TRUE)
TABLE
plotPoints(ETR~PAR, data= TABLE)
f<-fitModel(ETR ~ PAR/((a*PAR^2)+(b*PAR)+c),data = TABLE, start = list(a=0, b=0.004, c=1))
coef(f)
plotFun(f, add=TRUE)

My results: 

> TABLE=read.table("sample 2.txt", header = TRUE)
> TABLE
  PAR         ETR
1   1   0.4090909
2  46  18.8181818
3 126  51.5454546
4 234  93.6000000
5 404 173.1428571
6 656 246.0000000
7 816 272.0000000
> plotPoints(ETR~PAR, data= TABLE)
> f<-fitModel(ETR ~ PAR/((a*PAR^2)+(b*PAR)+c),data = TABLE, start = list(a=0, b=0.0035, c=1))
> coef(f)
            a             b             c 
 2.921397e-06 -2.027561e-03  2.718527e+00 
> plotFun(f, add=TRUE)

The problem is i do not want to have negative values for coef a,b, and c. This is because i am going to need these coefficients to get biological factors ETRmax= 1/b+2ac, Ek= c/b+2ac, alfa= 1/c which cannot be negative. I also need to use the Gauss-newton least square method and not the other algorithms where i can specify lower and upper bounds. I have tried to use the port algorithm to specify lower bounds of 0, but the problem is that b was 0, and i do not want b to be 0, i want all my variables to be greater than 0. I have tried many things but I can't solve this problem. I hope that you can help me with this, thank you.

Tamara Bou Chahine
  • 3
  • 1

1 Answers1

0

Specify lower bounds with the lower = c(0, 0, 0), alg = "port" :

fitModel(ETR ~ PAR/((a*PAR^2)+(b*PAR)+c), data = TABLE, 
  start = list(a = 0, b = 0.0035, c = 1), lower = c(0, 0, 0), alg = "port")
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
  • Thank you for your reply. I have already tried this, and b is 0, which i do not want. i want all coef to be greater than 0 – Tamara Bou Chahine Jun 05 '19 at 17:15
  • It is going to drive it to the lower bound. You can use a small number instead of 0. – G. Grothendieck Jun 05 '19 at 17:48
  • Thank you, this works. I think R is going to choose the smallest positive value possible when it needs to correct for negative coef values. But i was wondering if i could tell R that i want my variables a> 0, b>0 and c>0 without constraining it to certain value/lower bound for each. I just want my coef results to be positive. – Tamara Bou Chahine Jun 06 '19 at 09:47
  • There is no minimum on such open sets. For any proposed solution on such a set you will be able to find another proposed solution that produces a lesser objective value. For example, the minimum of x^2 on the closed set x >= 0 is 0 but there is no minimum on the open set x > 0. For any x', x'/2 will give a smaller objective. – G. Grothendieck Jun 06 '19 at 10:40
  • update: what i did was removing some data point from my data when i got negative coef. I tried to simulate the "model" curve with my data by removing the data points that were laying outside the "normal" shape of my model curve. This way i got positive coef. I know it is overfitting and modifying my data, but i think by putting constraints on my coeff i was doing that as well. – Tamara Bou Chahine Jun 06 '19 at 11:02
  • You could also try using robust nonlinear least squares. `library(robustbase); nlrob(ETR ~ PAR/((a*PAR^2)+(b*PAR)+c), data = TABLE, start = list(a = 0, b = 0.0035, c = 1), lower = c(a = 0, b = 0, c = 0), upper = c(a = 1, b = 1, c = 1), method = "MM")` – G. Grothendieck Jun 06 '19 at 11:29