Um texto embutido e sem conjunto de caracteres especificado associado... Nome: n?o dispon?vel Url: https://stat.ethz.ch/pipermail/r-help/attachments/20071202/9577be0a/attachment.pl
OK. Since the model is linear except for A lets use brute force to repeatedly evaluate the sum of squares for values of A between -2 and 2 proceeding in steps of .01 solving the other parameters using lm. That will give us better starting values and we should be able to use nls on that.> x <- seq(-2, 2, .01) > ss <- sapply(x, function(A) sum(resid(lm(richness ~ I(area^A)))^2)) > plot(ss ~ x) > x[which.min(ss)][1] -0.45> model.lm <- lm(richness ~ I(area^-0.45)) > # use starting values based on lm and A = -0.45 > st <- c(Const = coef(model.lm)[[1]], B = coef(model.lm)[[2]], A = x[which.min(ss)]) > nls(richness ~ Const+B*(area^A), st = st)Nonlinear regression model model: richness ~ Const + B * (area^A) data: parent.frame() Const B A 33.9289 -33.4595 -0.4464 residual sum-of-squares: 8751 Number of iterations to convergence: 2 Achieved convergence tolerance: 3.368e-06 Note that our A value is suspiciously close to A = -0.5 and sqrt(area) is length so I wonder if there is an argument based on units of measurement that might support a model of the form: richness = Const + B / sqrt(area) On Dec 2, 2007 3:39 PM, Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> wrote:> > Dear Gabor, > > Thank you for your reply. > In fact I am ajusting several models at same time, like linear, log-linear, > log-log, piecewise etc. One of the models are the power model. I really need > to fit a power model because it one of the hypothesis which have been > suggested on literature. > > In addition, there are other variables which are beeing tested as > explanatory. > > Kind regards, > > miltinho > ----- Mensagem original ---- > De: Gabor Grothendieck <ggrothendieck at gmail.com> > Para: Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> > Cc: R-help <r-help at stat.math.ethz.ch> > Enviadas: Domingo, 2 de Dezembro de 2007 17:28:23 > Assunto: Re: [R] fitting "power model" in nls() > > > > Is that really the model we want? When we have problems sometimes > its just a sign that the model is not very good in the first place. > > plot(richness ~ area) > > shows most of the points crowded the left and just a few points out to > the right. This > does not seem like a very good pattern for model fitting. > > plot(richness ~ log(area)) > plot(log(richness) ~ log(area)) > > both look nicer. > > On Dec 2, 2007 2:08 PM, Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> > wrote: > > Dear all, > > I am still fighting against my "power model". > > I tryed several times to use nls() but I can?t run it. > > I am sending my variables and also the model which I would like to fit. > > As you can see, this "power model" is not the best model to be fit, but I > really need also to fit it. > > > > The model which I would like to fit is Richness = B*(Area^A) > > > > > richness<-c(44,36,31,39,38,26,37,33,34,48,25,22,44,5,9,13,17,15,21,10,16,22,13,20,9,15,14,21,23,23,32,29,20, > > > 26,31,4,20,25,24,32,23,33,34,23,28,30,10,29,40,10,8,12,13,14,56,47,44,37,27,17,32,31,26,23,31,34, > > > 37,32,26,37,28,38,35,27,34,35,32,27,22,23,13,28,13,22,45,33,46,37,21,28,38,21,18,21,18,24,18,23,22, > > 38,40,52,31,38,15,21) > > > area<-c(26.22,20.45,128.68,117.24,19.61,295.21,31.83,30.36,13.57,60.47,205.30,40.21, > > 7.99,1.18,5.40,13.37,4.51,36.61,7.56,10.30,7.29,9.54,6.93,12.60, > > 2.43,18.89,15.03,14.49,28.46,36.03,38.52,45.16,58.27,67.13,92.33,1.17, > > > 29.52,84.38,87.57,109.08,72.28,66.15,142.27,76.41,105.76,73.47,1.71,305.75, > > 325.78,3.71,6.48,19.26,3.69,6.27,1689.67,95.23,13.47,8.60,96.00,436.97, > > > 472.78,441.01,467.24,1169.11,1309.10,1905.16,135.92,438.25,526.68,88.88,31.43,21.22, > > > 640.88,14.09,28.91,103.38,178.99,120.76,161.15,137.38,158.31,179.36,214.36,187.05, > > > 140.92,258.42,85.86,47.70,44.09,18.04,127.84,1694.32,34.27,75.19,54.39,79.88, > > 63.84,82.24,88.23,202.66,148.93,641.76,20.45,145.31,27.52,30.70) > > plot(richness~area) > > > > I tryed to fit the following model: > > > > m1<-nls(richness ~ Const+B*(area^A)) > > > > Thanks a lot, > > miltinho > > Brazil. > > > > > > > > para armazenamento! > > > > [[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. > > > > > > > > ________________________________ > Abra sua conta no Yahoo! Mail, o ?nico sem limite de espa?o para > armazenamento!
Also the fitted values satisfy Const = -B = 33 (approximately) so we could try:> plot(richness ~ area) > nls(richness ~ C * (1 - 1/sqrt(area)), start = c(C = 33))Nonlinear regression model model: richness ~ C * (1 - 1/sqrt(area)) data: parent.frame() C 32.85 residual sum-of-squares: 8764 Number of iterations to convergence: 1 Achieved convergence tolerance: 5.595e-10> simple.nls <- .Last.value > points(fitted(simple.nls) ~ area, pch = "+", col = "red") >On Dec 2, 2007 4:06 PM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:> OK. Since the model is linear except for A lets use brute force to > repeatedly evaluate the sum of squares for values of A between > -2 and 2 proceeding in steps of .01 solving the other parameters using > lm. That will give us better starting values and we should be able to > use nls on that. > > x <- seq(-2, 2, .01) > > ss <- sapply(x, function(A) sum(resid(lm(richness ~ I(area^A)))^2)) > > plot(ss ~ x) > > x[which.min(ss)] > [1] -0.45 > > model.lm <- lm(richness ~ I(area^-0.45)) > > # use starting values based on lm and A = -0.45 > > st <- c(Const = coef(model.lm)[[1]], B = coef(model.lm)[[2]], A = x[which.min(ss)]) > > nls(richness ~ Const+B*(area^A), st = st) > Nonlinear regression model > model: richness ~ Const + B * (area^A) > data: parent.frame() > Const B A > 33.9289 -33.4595 -0.4464 > residual sum-of-squares: 8751 > > Number of iterations to convergence: 2 > Achieved convergence tolerance: 3.368e-06 > > Note that our A value is suspiciously close to A = -0.5 and sqrt(area) > is length so I wonder if there is an argument based on units of > measurement that might support a model of the form: > > richness = Const + B / sqrt(area) > > > > > > On Dec 2, 2007 3:39 PM, Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> wrote: > > > > Dear Gabor, > > > > Thank you for your reply. > > In fact I am ajusting several models at same time, like linear, log-linear, > > log-log, piecewise etc. One of the models are the power model. I really need > > to fit a power model because it one of the hypothesis which have been > > suggested on literature. > > > > In addition, there are other variables which are beeing tested as > > explanatory. > > > > Kind regards, > > > > miltinho > > ----- Mensagem original ---- > > De: Gabor Grothendieck <ggrothendieck at gmail.com> > > Para: Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> > > Cc: R-help <r-help at stat.math.ethz.ch> > > Enviadas: Domingo, 2 de Dezembro de 2007 17:28:23 > > Assunto: Re: [R] fitting "power model" in nls() > > > > > > > > Is that really the model we want? When we have problems sometimes > > its just a sign that the model is not very good in the first place. > > > > plot(richness ~ area) > > > > shows most of the points crowded the left and just a few points out to > > the right. This > > does not seem like a very good pattern for model fitting. > > > > plot(richness ~ log(area)) > > plot(log(richness) ~ log(area)) > > > > both look nicer. > > > > On Dec 2, 2007 2:08 PM, Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> > > wrote: > > > Dear all, > > > I am still fighting against my "power model". > > > I tryed several times to use nls() but I can?t run it. > > > I am sending my variables and also the model which I would like to fit. > > > As you can see, this "power model" is not the best model to be fit, but I > > really need also to fit it. > > > > > > The model which I would like to fit is Richness = B*(Area^A) > > > > > > > > richness<-c(44,36,31,39,38,26,37,33,34,48,25,22,44,5,9,13,17,15,21,10,16,22,13,20,9,15,14,21,23,23,32,29,20, > > > > > 26,31,4,20,25,24,32,23,33,34,23,28,30,10,29,40,10,8,12,13,14,56,47,44,37,27,17,32,31,26,23,31,34, > > > > > 37,32,26,37,28,38,35,27,34,35,32,27,22,23,13,28,13,22,45,33,46,37,21,28,38,21,18,21,18,24,18,23,22, > > > 38,40,52,31,38,15,21) > > > > > area<-c(26.22,20.45,128.68,117.24,19.61,295.21,31.83,30.36,13.57,60.47,205.30,40.21, > > > 7.99,1.18,5.40,13.37,4.51,36.61,7.56,10.30,7.29,9.54,6.93,12.60, > > > 2.43,18.89,15.03,14.49,28.46,36.03,38.52,45.16,58.27,67.13,92.33,1.17, > > > > > 29.52,84.38,87.57,109.08,72.28,66.15,142.27,76.41,105.76,73.47,1.71,305.75, > > > 325.78,3.71,6.48,19.26,3.69,6.27,1689.67,95.23,13.47,8.60,96.00,436.97, > > > > > 472.78,441.01,467.24,1169.11,1309.10,1905.16,135.92,438.25,526.68,88.88,31.43,21.22, > > > > > 640.88,14.09,28.91,103.38,178.99,120.76,161.15,137.38,158.31,179.36,214.36,187.05, > > > > > 140.92,258.42,85.86,47.70,44.09,18.04,127.84,1694.32,34.27,75.19,54.39,79.88, > > > 63.84,82.24,88.23,202.66,148.93,641.76,20.45,145.31,27.52,30.70) > > > plot(richness~area) > > > > > > I tryed to fit the following model: > > > > > > m1<-nls(richness ~ Const+B*(area^A)) > > > > > > Thanks a lot, > > > miltinho > > > Brazil. > > > > > > > > > > > > para armazenamento! > > > > > > [[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. > > > > > > > > > > > > > > ________________________________ > > Abra sua conta no Yahoo! Mail, o ?nico sem limite de espa?o para > > armazenamento! >
Um texto embutido e sem conjunto de caracteres especificado associado... Nome: n?o dispon?vel Url: https://stat.ethz.ch/pipermail/r-help/attachments/20071202/b6d8a742/attachment.pl
I played around with this a bit more and noticed that the "plinear" algorithm of nls converged using nearly every starting value I tried. In fact A = 0 was the only starting value that I could find that did not converge. Note that with "plinear" you only specify the starting values for non-linear parameters, in this case A, while the unnamed linear parameters are implied as coefficients of the columns of the matrix defined in the rhs.> nls(richness ~ cbind(1, area^A), start = c(A = 1), algorithm = "plinear")Nonlinear regression model model: richness ~ cbind(1, area^A) data: parent.frame() A .lin1 .lin2 -0.4464 33.9290 -33.4595 residual sum-of-squares: 8751 Number of iterations to convergence: 6 Achieved convergence tolerance: 4.968e-07 On Dec 2, 2007 4:06 PM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:> OK. Since the model is linear except for A lets use brute force to > repeatedly evaluate the sum of squares for values of A between > -2 and 2 proceeding in steps of .01 solving the other parameters using > lm. That will give us better starting values and we should be able to > use nls on that. > > x <- seq(-2, 2, .01) > > ss <- sapply(x, function(A) sum(resid(lm(richness ~ I(area^A)))^2)) > > plot(ss ~ x) > > x[which.min(ss)] > [1] -0.45 > > model.lm <- lm(richness ~ I(area^-0.45)) > > # use starting values based on lm and A = -0.45 > > st <- c(Const = coef(model.lm)[[1]], B = coef(model.lm)[[2]], A = x[which.min(ss)]) > > nls(richness ~ Const+B*(area^A), st = st) > Nonlinear regression model > model: richness ~ Const + B * (area^A) > data: parent.frame() > Const B A > 33.9289 -33.4595 -0.4464 > residual sum-of-squares: 8751 > > Number of iterations to convergence: 2 > Achieved convergence tolerance: 3.368e-06 > > Note that our A value is suspiciously close to A = -0.5 and sqrt(area) > is length so I wonder if there is an argument based on units of > measurement that might support a model of the form: > > richness = Const + B / sqrt(area) > > > > > > On Dec 2, 2007 3:39 PM, Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> wrote: > > > > Dear Gabor, > > > > Thank you for your reply. > > In fact I am ajusting several models at same time, like linear, log-linear, > > log-log, piecewise etc. One of the models are the power model. I really need > > to fit a power model because it one of the hypothesis which have been > > suggested on literature. > > > > In addition, there are other variables which are beeing tested as > > explanatory. > > > > Kind regards, > > > > miltinho > > ----- Mensagem original ---- > > De: Gabor Grothendieck <ggrothendieck at gmail.com> > > Para: Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> > > Cc: R-help <r-help at stat.math.ethz.ch> > > Enviadas: Domingo, 2 de Dezembro de 2007 17:28:23 > > Assunto: Re: [R] fitting "power model" in nls() > > > > > > > > Is that really the model we want? When we have problems sometimes > > its just a sign that the model is not very good in the first place. > > > > plot(richness ~ area) > > > > shows most of the points crowded the left and just a few points out to > > the right. This > > does not seem like a very good pattern for model fitting. > > > > plot(richness ~ log(area)) > > plot(log(richness) ~ log(area)) > > > > both look nicer. > > > > On Dec 2, 2007 2:08 PM, Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> > > wrote: > > > Dear all, > > > I am still fighting against my "power model". > > > I tryed several times to use nls() but I can?t run it. > > > I am sending my variables and also the model which I would like to fit. > > > As you can see, this "power model" is not the best model to be fit, but I > > really need also to fit it. > > > > > > The model which I would like to fit is Richness = B*(Area^A) > > > > > > > > richness<-c(44,36,31,39,38,26,37,33,34,48,25,22,44,5,9,13,17,15,21,10,16,22,13,20,9,15,14,21,23,23,32,29,20, > > > > > 26,31,4,20,25,24,32,23,33,34,23,28,30,10,29,40,10,8,12,13,14,56,47,44,37,27,17,32,31,26,23,31,34, > > > > > 37,32,26,37,28,38,35,27,34,35,32,27,22,23,13,28,13,22,45,33,46,37,21,28,38,21,18,21,18,24,18,23,22, > > > 38,40,52,31,38,15,21) > > > > > area<-c(26.22,20.45,128.68,117.24,19.61,295.21,31.83,30.36,13.57,60.47,205.30,40.21, > > > 7.99,1.18,5.40,13.37,4.51,36.61,7.56,10.30,7.29,9.54,6.93,12.60, > > > 2.43,18.89,15.03,14.49,28.46,36.03,38.52,45.16,58.27,67.13,92.33,1.17, > > > > > 29.52,84.38,87.57,109.08,72.28,66.15,142.27,76.41,105.76,73.47,1.71,305.75, > > > 325.78,3.71,6.48,19.26,3.69,6.27,1689.67,95.23,13.47,8.60,96.00,436.97, > > > > > 472.78,441.01,467.24,1169.11,1309.10,1905.16,135.92,438.25,526.68,88.88,31.43,21.22, > > > > > 640.88,14.09,28.91,103.38,178.99,120.76,161.15,137.38,158.31,179.36,214.36,187.05, > > > > > 140.92,258.42,85.86,47.70,44.09,18.04,127.84,1694.32,34.27,75.19,54.39,79.88, > > > 63.84,82.24,88.23,202.66,148.93,641.76,20.45,145.31,27.52,30.70) > > > plot(richness~area) > > > > > > I tryed to fit the following model: > > > > > > m1<-nls(richness ~ Const+B*(area^A)) > > > > > > Thanks a lot, > > > miltinho > > > Brazil. > > > > > > > > > > > > para armazenamento! > > > > > > [[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. > > > > > > > > > > > > > > ________________________________ > > Abra sua conta no Yahoo! Mail, o ?nico sem limite de espa?o para > > armazenamento! >