[R] Odp: [nls] singular gradient

Petr PIKAL petr.pikal at precheza.cz
Tue Oct 2 09:01:50 CEST 2007


Hi

Niner <ninerdummy at gmail.com> napsal dne 01.10.2007 18:37:17:

> Thank you very much, Petr. I define my model based on same experiment 
data. 
> The trend of (A,B) and the relationship with ev are assumed as
> 
> 1/A = a + b*ev + c*ev^2
> 1/B = d + e*exp(f*ev)
> 
> Both (1/A) and (1/B) are supposed to decrease (nonlinearly) with 
increasing 
> ev. The initial values I used are computed based on above equations. I 
knew 
> that the data is quite noisy. My purpose is to use the model such that 
when 
> given an interested ev, a curve of (N,CSR) can be built. 
> 
> Is it possible to force R to find (a,b,c,d,e,f,g) based on my equations? 
Or, 
> is there any other way I can do to improve the model? Thank you very 
much for your help.

As I said I am not an expert in  nonlinear modelling, but your data do not 
follow your expectation.

> 
> BTW, when I ran your codes in R, I got several errors like

Yes. It is an error when computing partial corelations in some groups. As 
far as some groups are computed you can use otput from nlsList as a 
starting point for nlme optimisation. For more complex explanation see 
Pinheiro&Bates, Mixed effects models with S and S+.

> "Error in sum(rr[-(1:npar)]^2) : subscript out of bounds"
> "Error in nls(formula = formula, data = data, start = start, control = 
control) : 
> singular gradient "
> "Error in nls(formula = formula, data = data, start = start, control = 
control) : 
> number of iterations exceeded maximum of 50 "
> "Error in nlsModel(formula, mf, start, wts) : 
> singular gradient matrix at initial parameter estimates "
> 
> But the plot is OK. I don't if these errors matters. Thanks again for 
your 
> response. Hope someone can help me to find a solution to my problem.

Not sure as random effects for A and B shall follouw your above mentioned 
relationships which they definitely dont. However maybe I am mistaken and 
somebody with better insight will correct me. I too would like to get 
opinion from experts.
Regards
Petr

> 
> Niner
> ----------------------------------------------------
> Niner, Seattle
> ninerdummy at gmail.com
> 
> On Oct 1, 2007, at 1:33 AM, Petr PIKAL wrote:
> 
> Althoug I am not an expert in nonlinear regression, it seems to me that 
> your data are quite noisy. I presume you have got your model from 
> elsewhere and now you try to use it for your data. Looking at it and 
using 
> nlme I suppose you will not get such definite model based on your data.
> 
> Try
> 
> data1$ev <- ordered(data1$ev)
> data1.gr <- groupedData(N~CSR|ev, data1)
> out <- nlsList(N~CSR/(A+B*CSR), data1.gr, start= c(A=1, B=1))
> out1<- nlme(out)
> ranef(out1)
>                 A             B
> 0.2  8.667265e-05 -6.267700e-04
> 0.3 -3.139092e-04  2.206547e-03
> 0.4  4.044455e-06 -2.349364e-05
> <snip>
> 
> 3.6 -7.402470e-08 -5.312590e-06
> 5   -3.844779e-06  2.942158e-05
> 
> plot(as.numeric(row.names(ranef(out1))), ranef(out1)[,1])
> plot(as.numeric(row.names(ranef(out1))), ranef(out1)[,2])
> 
> You can clearly see that A and B random efect values are oscilating 
around 
> 0 with increasing "ev" which seems to me that A and B are independent 
on 
> ev 
> 
> Regards
> Petr
> petr.pikal at precheza.cz
> 
> r-help-bounces at r-project.org napsal dne 01.10.2007 02:37:52:
> 
> Hi, I am new to R. I don't have strong background of statistics. I am 
> a student of Geotechnical Engineering. I tried to run a nonlinear 
> regression for a three-variable function, that is
> 
> N = f(CSR, ev)   # N is a function of CSR and ev, and N = CSR/(A 
> +B*CSR), wherer (A,B) are function of ev.
> 
> N, CSR and ev are observed in the experiments.
> Following is my R script.
> 
> rm(list=ls())
> library(nlme)
> 
> # assign data
> N <- c 
> 
(30.03,16.62,10.88,36.47,20.24,38.17,36.47,34.80,19.00,32.37,14.40,35.63 
> 
> , 
> 
19.00,17.79,33.98,31.58,31.58,35.63,20.24,31.58,29.27,22.18,27.77,25.60, 
> 
> 
19.00,7.05,34.80,29.27,29.27,17.79,10.42,31.58,17.79,17.20,11.36,19.00,2 
> 
> 
9.27,12.33,22.18,22.18,14.40,31.58,19.00,9.52,33.17,13.87,19.00,21.52,11 
> 
> .36,22.84,9.96,6.68,20.88,9.96,11.84,20.24,19.61,17.20,17.20)
> 
> CSR <- c 
> 
(0.25,0.42,0.12,0.438,0.49,0.42,0.47,0.46,0.24,0.45,0.37,0.46,0.337,0.36 
> 
> , 
> 
0.334,0.346,0.399,0.44,0.246,0.33,0.413,0.23,0.45,0.45,0.44,0.106,0.333, 
> 
> 
0.256,0.345,0.44,0.153,0.348,0.23,0.25,0.122,0.183,0.201,0.128,0.23,0.24 
> 
> , 
> 
0.129,0.438,0.228,0.111,0.409,0.14,0.24,0.20,0.22,0.22,0.152,0.094,0.131 
> 
> ,0.123,0.155,0.28,0.204,0.149,0.193)
> 
> ev <- c 
> 
(0.2,0.3,0.4,0.5,0.5,0.5,0.6,0.6,0.8,1,1.0,1.0,1.1,1.1,1.2,1.2,1.2,1.2,1 
> 
> . 
> 
3,1.3,1.3,1.3,1.3,1.3,1.4,1.5,1.5,1.5,1.5,1.5,1.6,1.6,1.6,1.6,1.7,1.7,1. 
> 
> 
7,1.7,1.7,1.7,1.8,1.8,1.9,1.9,1.9,1.9,1.9,2.0,2.1,2.1,2.3,2.4,2.4,2.7,2. 
> 
> 8,3.3,3.4,3.6,5.0)
> 
> # set the data frame
> data1<- data.frame(cbind(N,CSR,ev))
> # the initial values of parameters
> para1.st <- c(a=75.4,b=165.9,c=-22.4,d=0,e=0.8,f=-0.28,g=0)
> para2.st <- c(a=31.3,b=-176.9,c=53.41,d=0,e=75.2,f=-0.18,g=0)
> # call nls funciton
> try.control <- c(maxiter=100, minFactor=1/4096)
> out <- nls(N~CSR/(1/(a+b*ev+c*(ev^2))+CSR/(d+e*exp(f*ev)))-g, data1, 
> start=para1.st, control=try.control, trace=T)
> 
> The data is from experiments and the pattern of data is scatter quite 
> a bit. I want to find the best fit coefficients (a,b,c,d,e,f,g) for 
> the data.
> I have two different sets of initial values for try. But in both 
> cases I got error message of
> "singular gradint"
> 
> What can I do for this error? Is there any other nonlinear regression 
> model I can try?  This problem is kind of emergency. I really hope 
> someone can help me out. Any comment is appreciated. Thanks a lot.
> 
> Niner
> -----------------------------------------------
> Yi-Min Huang
> Civil & Environmental Engineering
> U. of Washington
> 206-6165697
> 
> ----------------------------------------------------
> Niner, Seattle
> ninerdummy at gmail.com
> 
>    [[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list