[R] questions about nls

[Spencer Graves]

	  There is too much here for me to parse it right now.  Did you walk 
through the code line by line looking at the output anc checking what it 
did at each step?

	  If you don't get a satisfactory reply from someone else and you can't 
figure it out from a line-by-line analysis of your code, please simplify 
the question:  Make up a toy problem that illustrates in only a few 
lines what you don't understand.  You are more likely to get the answer 
you want if it takes seconds rather than minutes to understand your 
question.

hth.  spencer graves

Mei-Chu Lo wrote:
> Dear R users,
>  
> I am new in R and I want to use the nls package to analyze some
> experimental data.  The data is in the attached file "data".  It is the
> response "Sav" measured at different "C0".  Basically, the "C0" is a
> function of C1, K2, and r, and the "Sav" is a function of C0, C1, K2,
> and r. The math equations are shown in the attached file"equations". 
> The parameters K2 and r are the physical properties I want to get from
> the non-linear regression.  The R codes I wrote is in the attached
> file"Rcode".  Basically, I wrote two functions.  First , I calculated
> the C1 for different C0 at the estimated K2 and r using the binary
> search method and implemented it in the function "fn1".  Then for the
> calculated C1 and estimated K2 and r, I used the function "fn2" to get
> the estimated Sav for different C0.  "nls" was used to minimize the
> differences between the calculated Sav and observed Sav.  When I run my
> R script, I got the error message " step factor reduced below minFactor"
> .  If I changed the minFactor to zero, the nls continued but did not
> converge (exceed the maxiteration).  I changed the tolerance to higher
> value, nls finished but the fitting is bad.  From the published results,
> the best fitted values for K2 and r are 0.2237 and 1.296*10^7,
> respectively.  I can't get these numbers using my R script.
> 
> 
> I know there are a lot math and R experts on the R mailing-list.  I
> will appreciate it if anyone can tell me what is wrong in my R script or
> in the methods I used to get these parameters.
> 
> 
> Mei-Chu Lo   
> 
> 
> ------------------------------------------------------------------------
> 
> C0      Sav
> 0.6766	6.0875
> 1.6165	6.4249
> 1.8796	6.5374
> 2.4436	7.025
> 4.8872	7.5125
> 5.3759	7.625
> 5.7518	7.8499
> 7.218	8.0749
> 9.2105	8.1125
> 12.7067	9.4624
> 12.5939	10.025
> 16.203	11.7125
> 17.2932	12.0124
> 18.1578	12.5749
> 21.5413	12.6875
> 21.5413	13.0625
> 
> 
> ------------------------------------------------------------------------
> 
> fn1<-function(C0,k2,r){
> m<-matrix(nrow=26,ncol=length(C0))
> Ct<-vector("numeric")
> C1<-vector("numeric")
> dC<-vector("numeric")
> j=1
>   repeat{
>   C1[j]<-C0[j]/2;dC[j]<-C0[j]/4
>     repeat{
>     m[1,j]=C1[j];m[2,j]=k2*C1[j]^2
>     m[26,j]=26*(k2/2)^25*r*C1[j]^26   
> 	i=3
>           repeat{
>           m[i,j]<-i*(k2/2)^(i-1)*C1[j]^i
>           i<-i+1
>           if(i>=26) break}
> 
>      Ct[j]<-sum(m[,j])
>      if(abs(Ct[j]-C0[j])<10^-6*C0[j]) break
>     
>      if(Ct[j]>C0[j]) {
>      C1[j]=C1[j]-dC[j];dC[j]<-dC[j]/2}
>      else{
>      C1[j]=C1[j]+dC[j];dC[j]<-dC[j]/2}
>     }
> 
>          
>   j=j+1
>   if (j>length(C0)) break }  
> 
> C1
> }
> 
> 
> fn2<-function(C1,C0,k2,r){
>    m<-matrix(nrow=26,ncol=length(C0))
>    m[1,]<-5.8*(1-0.018*C0)*C1;
>    m[2,]<-2^(2/3)*5.8*(1-0.018*C0)*k2*C1^2
>    m[26,]<-42*(1-0.019*C0)*26*(k2/2)^25*r*C1^26
>      
>       j=3
>          repeat{
>          m[j,]<-j^(5/3)*5.8*(1-0.018*C0)*(k2/2)^(j-1)*C1^j
>          j=j+1
>          if(j>=26)break   
>          }
>   
>   
>      p<-vector("numeric")
>       
>      k=1
>        repeat{
>        p[k]<-sum(m[,k])/C0[k]
>        k=k+1
>        if (k>length(C0))break
>        }
> p
>    }
> 
> # load data
> 
> ass<-read.table("data.txt",header=T)
> plot(ass$C0,ass$Sav,cex=1.5, xlab="C0(mg/ml)",ylab="Sav(S)")
> 
> # give initial guess of k2 and r
> k2=0.3;r=10^6
> fit<-nls(Sav~fn2(fn1(ass$C0,k2fit,rfit),ass$C0,k2fit,rfit),data=ass,start=list(k2fit=k2,rfit=r))
> 
> 
> ------------------------------------------------------------------------
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help