[R] Determining Starting Values for Model Parameters in Nonlinear Regression
Paul Bernal
p@u|bern@|07 @end|ng |rom gm@||@com
Sun Aug 20 04:58:13 CEST 2023
Dear Dr. Nash,
Thank you for your kind feedback. I was able to successfully find the
starting values for the nonlinear regression model parameters, here is my
code:
#Specifying Nonlinear Regression Model
mymod <- y ~ 1/(Beta1*x1 + Beta2*x2 + Beta3*x3)
strt <- c(Beta1 = 1, Beta2 = 2, Beta3 = 3)
trysol<-nlxb(formula=mymod, data=mod14data_random, start=strt, trace=TRUE)
trysol
#Examining structure of trysol model
str(trysol)
#Extracting Coefficients for Beta1, Beta2, and Beta3
trysol$coefficients[[1]]
trysol$coefficients[[2]]
trysol$coefficients[[3]]
The dataset is the same I shared previously, but for safety, I will share
it again below:
>dput(mod14data_random)
structure(list(Mixture = c(17, 14, 5, 1, 11, 2, 16, 7, 19, 23,
20, 6, 13, 21, 3, 18, 15, 26, 8, 22), x1 = c(69.98, 72.5, 77.6,
79.98, 74.98, 80.06, 69.98, 77.34, 69.99, 67.49, 67.51, 77.63,
72.5, 67.5, 80.1, 69.99, 72.49, 64.99, 75.02, 67.48), x2 = c(29,
25.48, 21.38, 19.85, 22, 18.91, 29.99, 19.65, 26.99, 29.49, 32.47,
20.35, 26.48, 31.47, 16.87, 27.99, 24.49, 31.99, 24.96, 30.5),
x3 = c(1, 2, 1, 0, 3, 1, 0, 2.99, 3, 3, 0, 2, 1, 1, 3, 2,
3, 3, 0, 2), y = c(1.4287, 1.4426, 1.4677, 1.4774, 1.4565,
1.4807, 1.4279, 1.4684, 1.4301, 1.4188, 1.4157, 1.4686, 1.4414,
1.4172, 1.4829, 1.4291, 1.4438, 1.4068, 1.4524, 1.4183)), row.names =
c(NA,
-20L), class = "data.frame")
I was able to get the results you shared:
> trysol
residual sumsquares = 1.5412e-05 on 20 observations
after 29 Jacobian and 43 function evaluations
name coeff SE tstat pval gradient
JSingval
Beta1 0.00629212 5.997e-06 1049 2.425e-42 4.049e-08
721.8
Beta2 0.00867741 1.608e-05 539.7 1.963e-37 -2.715e-08
56.05
Beta3 0.00801948 8.809e-05 91.03 2.664e-24 1.497e-08
10.81
Now my question is, how could I perform adequacy checking on this nonlinear
regression model? I was trying to use the deviance metric as follows:
#We will fit several models
nlregmod2 <- nls(y ~ 1/(Beta1*x1),
start =
list(Beta1 = trysol$coefficients[[1]]),
data=mod14data_random)
nlregmod3 <- nls(y ~ 1/(Beta2*x2),
start =
list(Beta2 = trysol$coefficients[[2]]),
data=mod14data_random)
nlregmod4 <- nls(y ~ 1/(Beta3*x3),
start =
list(Beta3 = trysol$coefficients[[3]]),
data=mod14data_random)
nlregmod5 <- nls(y ~ 1/(Beta1*x1+Beta2*x2),
start =
list(Beta1 = trysol$coefficients[[1]],
Beta2 = trysol$coefficients[[2]]),
data=mod14data_random)
nlregmod6 <- nls(y ~ 1/(Beta1*x1+Beta3*x3),
start =
list(Beta1 = trysol$coefficients[[1]],
Beta3 = trysol$coefficients[[3]]),
data=mod14data_random)
nlregmod7 <- nls(y ~ 1/(Beta2*x2+Beta3*x3),
start =
list(Beta2 = trysol$coefficients[[2]],
Beta3 = trysol$coefficients[[3]]),
data=mod14data_random)
deviance(nlregmod)
deviance(nlregmod2)
deviance(nlregmod3)
deviance(nlregmod4)
deviance(nlregmod5)
deviance(nlregmod6)
deviance(nlregmod7)
However, for some models I get the following error:
nlregmod4 <- nls(y ~ 1/(Beta3*x3),
start =
list(Beta3 = trysol$coefficients[[3]]),
data=mod14data_random)
Error in numericDeriv(form[[3L]], names(ind), env, central = nDcentral) :
Missing value or an infinity produced when evaluating the model
deviance(nlregmod5)
Error: object 'nlregmod5' not found
deviance(nlregmod6)
Error: object 'nlregmod6' not found
deviance(nlregmod7)
Error: object 'nlregmod7' not found
I spent a significant amount of time looking for a way to assess the
nonlinear regression model goodness of fit and to perform adequacy checking
using R (and in general) but was not able to find anything.
I would have thought that calculating the deviance for different models
would be a good way to perform adequacy checking, but got some errors. The
best I could find was this link
https://cran.r-project.org/web/packages/nlstools/vignettes/vignetteJSS.pdf.
Any help and/or guidance will be greatly appreciated.
Kind regards,
Paul
El El sáb, 19 de ago. de 2023 a la(s) 9:30 p. m., J C Nash <
profjcnash using gmail.com> escribió:
> The cautions people have given about starting values are worth heeding.
> That nlxb() does well in many cases is useful,
> but not foolproof. And John Fox has shown that the problem can be tackled
> very simply too.
>
> Best, JN
>
>
> On 2023-08-19 18:42, Paul Bernal wrote:
> > Thank you so much Dr. Nash, I truly appreciate your kind and valuable
> contribution.
> >
> > Cheers,
> > Paul
> >
> > El El sáb, 19 de ago. de 2023 a la(s) 3:35 p. m., J C Nash <
> profjcnash using gmail.com <mailto:profjcnash using gmail.com>> escribió:
> >
> > Why bother. nlsr can find a solution from very crude start.
> >
> > Mixture <- c(17, 14, 5, 1, 11, 2, 16, 7, 19, 23, 20, 6, 13, 21, 3,
> 18, 15, 26, 8, 22)
> > x1 <- c(69.98, 72.5, 77.6, 79.98, 74.98, 80.06, 69.98, 77.34, 69.99,
> 67.49, 67.51, 77.63,
> > 72.5, 67.5, 80.1, 69.99, 72.49, 64.99, 75.02, 67.48)
> > x2 <- c(29, 25.48, 21.38, 19.85, 22, 18.91, 29.99, 19.65, 26.99,
> 29.49, 32.47,
> > 20.35, 26.48, 31.47, 16.87, 27.99, 24.49, 31.99, 24.96,
> 30.5)
> > x3 <- c(1, 2, 1, 0, 3, 1, 0, 2.99, 3, 3, 0, 2, 1, 1, 3, 2,
> > 3, 3, 0, 2)
> > y <- c(1.4287, 1.4426, 1.4677, 1.4774, 1.4565,
> > 1.4807, 1.4279, 1.4684, 1.4301, 1.4188, 1.4157, 1.4686,
> 1.4414,
> > 1.4172, 1.4829, 1.4291, 1.4438, 1.4068, 1.4524, 1.4183)
> > mydata<-data.frame(Mixture, x1, x2, x3, y)
> > mydata
> > mymod <- y ~ 1/(Beta1*x1 + Beta2*x2 + Beta3*x3)
> > library(nlsr)
> > strt<-c(Beta1=1, Beta2=2, Beta3=3)
> > trysol<-nlxb(formula=mymod, data=mydata, start=strt, trace=TRUE)
> > trysol
> > # or pshort(trysol)
> >
> >
> > Output is
> >
> > residual sumsquares = 1.5412e-05 on 20 observations
> > after 29 Jacobian and 43 function evaluations
> > name coeff SE tstat pval
> gradient JSingval
> > Beta1 0.00629212 5.997e-06 1049 2.425e-42
> 4.049e-08 721.8
> > Beta2 0.00867741 1.608e-05 539.7 1.963e-37
> -2.715e-08 56.05
> > Beta3 0.00801948 8.809e-05 91.03 2.664e-24
> 1.497e-08 10.81
> >
> > J Nash
> >
> >
> > On 2023-08-19 16:19, Paul Bernal wrote:
> > > Dear friends,
> > >
> > > Hope you are all doing well and having a great weekend. I have
> data that
> > > was collected on specific gravity and spectrophotometer analysis
> for 26
> > > mixtures of NG (nitroglycerine), TA (triacetin), and 2 NDPA (2 -
> > > nitrodiphenylamine).
> > >
> > > In the dataset, x1 = %NG, x2 = %TA, and x3 = %2 NDPA.
> > >
> > > The response variable is the specific gravity, and the rest of the
> > > variables are the predictors.
> > >
> > > This is the dataset:
> > > dput(mod14data_random)
> > > structure(list(Mixture = c(17, 14, 5, 1, 11, 2, 16, 7, 19, 23,
> > > 20, 6, 13, 21, 3, 18, 15, 26, 8, 22), x1 = c(69.98, 72.5, 77.6,
> > > 79.98, 74.98, 80.06, 69.98, 77.34, 69.99, 67.49, 67.51, 77.63,
> > > 72.5, 67.5, 80.1, 69.99, 72.49, 64.99, 75.02, 67.48), x2 = c(29,
> > > 25.48, 21.38, 19.85, 22, 18.91, 29.99, 19.65, 26.99, 29.49, 32.47,
> > > 20.35, 26.48, 31.47, 16.87, 27.99, 24.49, 31.99, 24.96, 30.5),
> > > x3 = c(1, 2, 1, 0, 3, 1, 0, 2.99, 3, 3, 0, 2, 1, 1, 3, 2,
> > > 3, 3, 0, 2), y = c(1.4287, 1.4426, 1.4677, 1.4774, 1.4565,
> > > 1.4807, 1.4279, 1.4684, 1.4301, 1.4188, 1.4157, 1.4686,
> 1.4414,
> > > 1.4172, 1.4829, 1.4291, 1.4438, 1.4068, 1.4524, 1.4183)),
> row.names =
> > > c(NA,
> > > -20L), class = "data.frame")
> > >
> > > The model is the following:
> > > y = 1/(Beta1x1 + Beta2x2 + Beta3x3)
> > >
> > > I need to determine starting (initial) values for the model
> parameters for
> > > this nonlinear regression model, any ideas on how to accomplish
> this using
> > > R?
> > >
> > > Cheers,
> > > Paul
> > >
> > > [[alternative HTML version deleted]]
> > >
> > > ______________________________________________
> > > R-help using r-project.org <mailto:R-help using r-project.org> mailing list
> -- To UNSUBSCRIBE and more, see
> > > https://stat.ethz.ch/mailman/listinfo/r-help <
> https://stat.ethz.ch/mailman/listinfo/r-help>
> > > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > <http://www.R-project.org/posting-guide.html>
> > > and provide commented, minimal, self-contained, reproducible code.
> >
>
[[alternative HTML version deleted]]
More information about the R-help
mailing list