--- title: "An introduction to nlsr" author: "John C. Nash" date: "`r Sys.Date()`" output: pdf_document bibliography: ImproveNLS.bib vignette: > %\VignetteIndexEntry{nlsr Introduction} %\usepackage[utf8]{inputenc} %\VignetteEngine{knitr::rmarkdown} --- **R** has several tools for estimating nonlinear models and minimizing sums of squares functions. The base R distribution includes a function `nls()` that is feature-rich but has some annoying weaknesses (@JNABRefactoringNLS22). Some of these are addressed in the package **nlsr** (@nlsr2019manual), for which this vignette provides an introduction with examples. ## What does nlsr do? **nlsr** has functions to find a set of parameters that minimizes the sum of squares (or deviance) of differences between a model specified either as a formula using function `nlxb()` or as both residual and Jacobian R functions using function `nlfb()`. As described below, the Jacobian may be approximated by particular control settings. Existing R tools that are similar to `nlxb()` are `nls()` and `minpack.lm::nlsLM()`; `minpack.lm::nls.lm()` is similar to `nlfb()`. There are several important differences in approach and results, some of which are detailed in this document. ### Output object The output object of `nlxb()` is smaller than that of `nls()` and is focused on the solution and its properties rather than the **modeling** task. That is, package `nlsr` emphasizes the solution of the nonlinear least squares problem rather than the estimation of a nonlinear model that fits or explains the data. `nls()` and `nlsLM` return an object of class `nls`, which allows for a number of specialized modeling and diagnostic extensions. For compatibility, the `nlsr` package has function `wrapnlsr()`, for which `nlsr()` is an alias. This uses `nlxb()` to find good parameters, then calls `nls()` to return the class `nls` object. Unless particular modeling features are needed, the use of `wrapnlsr()` is unnecessary and wasteful of resources. ### Jacobian calculation Internally, `nlxb()` uses symbolic and automatic differentiation tools to create the residual and Jacobian functions needed for `nlfb()` then uses it to solve the relevant nonlinear least squares minimization problem. As with other nonlinear least squares packages, we provide the Jacobian code as the "gradient" attribute of the Jacobian function. This lets us embed the code for the Jacobian as this attribute of the **residual** function so that the call to `nlfb()` can be made with the same name used for both residual and Jacobian function arguments. This programming trick saves a lot of trouble for the package developer, but it can be a nuisance for users trying to understand the code. Generally `nls()` and `nlsLM()` use a fairly simple forward difference approximation of derivatives for the Jacobian, though a central approximation can be specified in control parameters. Package `nlsr` provides four approximation options, with a further control `ndstep` for the size of the step used in the approximation. ### Stabilization of Gauss-Newton computations The iteration in all the major tools mentioned is the solution of the Gauss-Newton equations. `nls()` uses a variant of an approach suggested by @Hartley61, while `nlsr` and `minpack.lm` use variants of the method published by @Marquardt1963, though it was suggested earlier by @Levenberg1944. The details of the method in `nlsr` are discussed in @JNABRefactoringNLS22. Though it is unlikely to be useful in general, control settings for `nlxb()` or `nlfb()` allow for those routines to perform a variety of Hartley and Marquardt algorithms. This has been useful for learning about the algorithms, but tangential to simply finding parameters of nonlinear models. In general, the Levenberg-Marquardt stabilization is important in obtaining solutions. The default method of `nls()` terminates with `singular gradient` unnecessarily in too many instances. ### Programming language An important choice made in developing `nlsr` was to code entirely within the R programming language. `nls()` uses a mix of R, C and Fortran, as does `minpack.lm`. Generally, I believe that keeping to a single programming language allows for easier maintenance and upgrades. On the other hand, R is usually somewhat slower in performing computations because it keeps track of names and structures and because it is usually interpreted rather than compiled. In recent years the performance penalty for using code entirely in R has been much reduced with the just-in-time compiler and other improvements, so that using an all-R package offers acceptable performance. Indeed, in `nlsr` the use of R may be less of a performance cost than the attempt to ensure the solution has been found, which forces more iterations to be used. ## An illustrative example The Hobbs weed infestation problem (@cnm79 [page 120]) is a growth curve modeling task which is seemingly straightforward but turns out to be quite nasty. The data and a graph of it is given below. ```{r ex01, echo=TRUE} weed <- c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443, 38.558, 50.156, 62.948, 75.995, 91.972) tt <- 1:12 weeddf <- data.frame(tt, weed) plot(weeddf, main="Hobbs weed infestation data") ``` Most nonlinear least squares problems encountered by R users have data, but it is not strictly required, as we shall see below in the section on functional specification of nonlinear least squares problems. What is required is a definition of the (residual) functions that are to be squared and then their squares summed to provide the **loss function** or **objective function** to be minimized by adjusting the **parameters** that define the set of functions. In the case of the Hobbs problem, I was told the scientists who collected the data believe that the growth of weeds followed a 3-parameter logistic sigmoid growth curve where $y$ (the `weed` data) changes with time ($t$) (the `tt` index) according to the function **Logistic3U**: $$ y \approx b_1 / (1 + b_2 * exp(- b_3 * t)) $$ The problem (using the 3 parameter logistic form above) has very bad scaling and regions in the parameter space where the sum of squares objective is nearly "flat", so that its Hessian, or matrix of second derivatives, is almost singular. When this occurs, methods cannot easily make progress towards the optimal solution. Solution methods also need an initial guess for the parameters to be adjusted. There are often some ways to provide such initial parameters, though in my experience they take quite a lot of work to set up reliably, but R does have some worthwhile possibilities in some, though not all, of the **selfStart** modelling functions, especially those in the package `nlraa` (@MiguezNLRAA2021). There are, however, alternatives to the model above. For example, a simple scaling can make the solution to the problem easier to find, such as **Logistic3S**: $$ y \approx 100 * c_1 / (1 + 10 * c_2 * exp(- 0.1 * c_3 * t)) $$ Another form is **Logistic3T**: $$ y \approx Asym / (1 + exp((xmid - t)/scal)) $$ The functions above are equivalent. Their parameters are related as follows: $$ Asym = b_1 = 100 * c_1 $$ $$ exp(xmid/scal) = b_2 = 10 * c_2 $$ $$ 1/scal = b_3 $$ `nlsr` tries to find solutions from very poor initial parameter values. While this is not always possible, `nslr` has generally been able to succeed in finding solutions, even when all parameters are started at 1. Let us compare its results with those of `nls()` and `nlsLM()`. ### Problem setup ```{r ex02set, echo=TRUE} # model formulas frmu <- weed ~ b1/(1+b2*exp(-b3*tt)) frms <- weed ~ 100*c1/(1+10*c2*exp(-0.1*c3*tt)) frmt <- weed ~ Asym /(1 + exp((xmid-tt)/scal)) # # Starting parameter sets stu1<-c(b1=1, b2=1, b3=1) sts1<-c(c1=1, c2=1, c3=1) stt1<-c(Asym=1, xmid=1, scal=1) ``` ### Solution attempts with nls() ```{r ex02nls} unls1<-try(nls(formula=frmu, start=stu1, data=weeddf)) summary(unls1) snls1<-try(nls(formula=frms, start=sts1, data=weeddf)) summary(snls1) tnls1<-try(nls(formula=frmt, start=stt1, data=weeddf)) summary(tnls1) cat("\n") ``` ### Solution attempts with nlsr tools ```{r ex02nlsr} library(nlsr) unlx1<-try(nlxb(formula=frmu, start=stu1, data=weeddf)) print(unlx1) pshort(unlx1) # A short form output snlx1<-try(nlxb(formula=frms, start=sts1, data=weeddf)) # pshort(snlx1) print(snlx1) tnlx1<-try(nlxb(formula=frmt, start=stt1, data=weeddf)) # pshort(tnlx1) print(tnlx1) ct<-as.list(tnlx1$coefficients) # Need a list to use ct$... in next line cat("exp(xmid/scal)=",exp(ct$xmid/ct$scal),"\n") cat("\n") # explicit residuals (weighted if there are weights) rtnlx1<-residuals(tnlx1) print(rtnlx1) cat("explicit sumsquares =", sum(rtnlx1^2),"\n") ``` ### Solution attempts with minpack.lm NOTE: The original version of this vignette passed all package checks for Linux, Windows and (Intel-based) Macintosh computers, but failed on Macintosh machines using the M1 chip. ```{r ex02minpack} library(minpack.lm) unlm1<-try(nlsLM(formula=frmu, start=stu1, data=weeddf)) summary(unlm1) # Use summary() to get display unlm1 # Note the difference. Use this form to get sum of squares # pnls(unlm1) # Short form of output snlm1<-try(nlsLM(formula=frms, start=sts1, data=weeddf)) # pnls(snlm1) summary(snlm1) tnlm1<-try(nlsLM(formula=frmt, start=stt1, data=weeddf)) if (inherits(tnlm1, "try-error")) { cat("Failure to compute solution -- likely singular Jacobian\n") } else { pnls(tnlm1) # short form to give sum of squares summary(tnlm1) } ``` ### Solution attempts with wrapnlsr() wrapper ```{r ex02wrapnlsr} ## Try the wrapper. Calling wrapnlsr() instead of nlsr() is equivalent ## unlw1<-try(nlsr(formula=frmu, start=stu1, data=weeddf)) print(unlw1) # using 'unlx1' gives name of object as 'x' only snlw1<-try(nlsr(formula=frms, start=sts1, data=weeddf)) # pshort(snlx1) print(snlw1) tnlw1<-try(nlsr(formula=frmt, start=stt1, data=weeddf)) print(tnlw1) ``` ### Commentary There are several details that may be important for users. - `nlsr` is set up to use `print()` to output standard errors and singular values of the Jacobian (for diagnostic purposes). By contrast, `minpack.lm` and `nls()` use `summary()`, which does NOT display the sum of squares, while `print()` gives the sum of squares, but not the standard errors. - The singular values displayed by `print.nlsr()` (the internal name for the adaptation of the generic `print()`) are displayed in a column to the right of the coefficient and standard error display, but are NOT specific to the parameters. Their position is purely for efficient use of page space. - The most common use of the singular values is to gauge how "nearly singular" the Jacobian is at the solution, and the ratio of the largest to smallest of the singular values is a simple but effective measure. In the above example, we note that the scaled logistic has the smallest ratio. - The result from `nlsLM` for the transformed model has a very large sum of squares, which may suggest that the program has failed. Since neither `nls()` nor `nlsLM()` offer the singular values, we can "cheat" and use `nlxb()`, though the Jacobian used will be the analytic one used by this last program rather than the forward difference approximation that is generally used by the others. ```{r ex03, echo=TRUE} stspecial<- c(Asym = 35.532, xmid = 43376, scal = -2935.4) # force to exit at once by setting femax to 1 (maximum number of sum of squares evaluations) getsvs<-nlxb(formula=frmt, start=stspecial, data=weeddf, control=list(femax=1)) print(getsvs) ``` - The trick used above to get singular values is convenient, but it is actually quite easy to get the Jacobian from a class `nls` object such as `tnlm1` as follows. ```{r ex04, echo=TRUE, eval=FALSE} if (inherits(tnlm1, "try-error")) { cat("Cannot compute solution -- likely singular Jacobian\n") } else { Jtm <- tnlm1$m$gradient() svd(Jtm)$d # Singular values } ``` We see that there are differences in detail, but the more important result is that two out of three singular values are essentially 0. Our Jacobian is singular, and no method of the Gauss-Newton type should be able to continue. Indeed, from this set of parameters, `nlxb` also stops. ```{r ex05, echo=TRUE} stspecial<- c(Asym = 35.532, xmid = 43376, scal = -2935.4) # force to exit at once by setting femax to 1 (maximum number of sum of squares evaluations) badstart<-nlxb(formula=frmt, start=stspecial, data=weeddf) print(badstart) ``` ## selfStart options The form of the logistic given by `frmt` is available as a **selfStart** model, needing no starting parameters with `nls()`. The code is in base R in location ` src/library/stats/R/zzModels.R`. ```{r ex06, echo=TRUE} frmtss <- weed ~ SSlogis(tt, Asym, xmid, scal) ssts1<-nls(formula=frmtss, data=weeddf) summary(ssts1) require(minpack.lm) sstm1<-nlsLM(formula=frmtss, data=weeddf) summary(sstm1) ``` The essence of selfStart models is to provide starting parameters for the iterative methods used to find the final parameters. However, an examination of the code for the `SSlogis()` function shows that it makes use of a DIFFERENT solver and returns the final results of that solver. The computational effort in doing this is, I believe, greater than that expended by `nlxb()` from an extremely crude start (all 1's). An alternative selfStart for this version of the logistic is included in package `nlsr` as `SSlogisJN.R`. Users should be aware that the programming of selfStart models involves quite esoteric aspects of R which I find prone to errors and difficult to grasp easily. Developers who have built them for us deserve our thanks. A particularly important collection of such tools is in the package @MiguezNLRAA2021. ### SSlogisJN.R for the 3-parameter logistic `SSlogisJN.R` uses ideas for selfStart from @Ratkowsky1983 for the Logistic3T form of the model. It produces only an approximate set of starting parameters, unlike `SSlogis.R`. The core of the idea is as follows, where we use the abbreviated output function to save page space. ```{r ex07, echo=TRUE} Asymt <- 2*max(weed) Asymt pw <- Asymt/weed pw1<-pw-1 lpw1<-log(pw1) clin <- coef(lm(lpw1~tt)) clin<-as.list(clin) scalt <- -1/clin$tt scalt xmidt <- scalt*as.numeric(clin[1]) xmidt library(nlsr) try1<-nlxb(frmt, data=weeddf, start=c(Asym=1, xmid=1, scal=1)) pshort(try1) try2<-nlxb(frmt, data=weeddf, start=c(Asym=Asymt, xmid=xmidt, scal=scalt)) pshort(try2) ``` We see a dramatic reduction in the computational effort as measured by residual and Jacobian evaluations. ### Running selfStart models with nlxb() There are two important steps in using selfStart models with `nlxb()`: - we must provide the starting parameters explicitly to `nlxb()`, which seems counter to the goal of selfStart models. However, we use function `getInitial()` to exploit particular selfStart modeling functions. - we must point to a mechanism to compute the Jacobian, since the selfStart function is unlikely to be in the derivatives table of analytic or automatic derivatives. Note that Jacobian ("gradient") code is generally part of selfStart functions. If code is included for the analytic Jacobian, it is likely worth using. However, because it is easy to make errors in coding derivatives, it may also be worth checking results of such code, for example, using tools from package `numDeriv`. I have made more such errors over my career than I like to admit. These steps can be automated, and this has been carried out in function `nlsrSS()`. The automation is built into the functions `nls()` and `minpack.lm::nlsLM()`, but I have preferred to make the process explicit and give it a separate function. Let us illustrate how to use the `getInitial()` function (part of base R) to do this. ```{r ex08, eval=TRUE} frmtssJ <- weed ~ SSlogisJN(tt, Asym, xmid, scal) # The model formula lstrt <- getInitial(frmtssJ, weeddf) # Here we get the starting parameters cat("starting parameters:\n") print(lstrt) # Next line uses the start. We indicate that the "approximate" Jacobian code is in # the selfStart code, though in fact it is not an approximation in this case. sstx1<-nlxb(formula=frmtssJ, start=lstrt, data=weeddf, control=list(japprox="SSJac")) print(sstx1) # If no Jacobian code is included in selfStart module, we can actually use an # approximate Jacobian. See "gradient" elements for differences in result. sstx1a<-nlxb(formula=frmtssJ, start=lstrt, data=weeddf, control=list(japprox="jacentral")) print(sstx1a) ## We can automate the above in function nlsrSS() sstSS<-nlsrSS(formula=frmtssJ, data=weeddf) print(sstSS) ``` ### Starting parameters for the Logistic3U model A similar approach for the Logistic3U model can be used. ```{r ex09, echo=TRUE} b1t <- 2*max(weed) b1t lpw1<-log(b1t/weed-1) clin <- as.list(coef(lm(lpw1~tt))) b3t <- -clin$tt b3t b2t <- exp(as.numeric(clin[1])) b2t library(nlsr) try1<-nlxb(frmu, data=weeddf, start=c(b1=1, b2=1, b3=1)) pshort(try1) try2<-nlxb(frmu, data=weeddf, start=c(b1=b1t, b2=b2t, b3=b3t)) pshort(try2) ``` Once again, there is a significant saving in the number of iterations. We have not coded a selfStart function for Logistic3U, preferring to use the scaled form Logistic3S for which we do not need good starting values. ## Jacobian computation `nlxb()` offers more options for how the Jacobian is computed than either `nls()` or `minpack.lm::nlsLM()`. The choice is made using the `control` list element `japprox`. If this is NOT specified, `nlxb()` attempts to build Jacobian code using analytic or automatic derivative codes. This is not, of course, always possible, and sometimes we will get an "error" message that required information is not in the derivatives table. This will be the case when we use a function like `SSlogisJN()` unless we indicate that the code is available from a selfStart module by using a value of `SSJac` for `japprox`. The example also showed how we can specify an approximation method, which is also useful when a formula does not have analytic derivatives available. In such cases `control` element `japprox` can be set to one of `jafwd`, `jaback`, `jacentral` or `jand`, giving forward, backward, central and package `numDeriv` approximations. I recommend `jacentral` as a reasonable compromise between accuracy of approximation and amount of computational work. ## Bounds constraints on parameters In some cases, we know that parameters cannot be bigger or smaller than some externally known limits. Cash reserves cannot be negative. Animals have a minimum need for water. Airplanes cannot carry more than a known or legislated cargo weight. Such limits can be built into models, but there are some important details for using the tools in R. - `nls()` can only impose bounds if the `algorithm="port"` argument is used in the call. Unfortunately, the documentation warns us: *The algorithm = "port" code appears unfinished, and does not even check that the starting value is within the bounds. Use with caution, especially where bounds are supplied.* - `nlsLM()` includes bounds in the standard call, but I have observed cases where it fails to get the correct answer. From my examination of the code, I believe the authors have not taken into account all possibilities, though it should be noted that I regard all programs as having some weakness in regard to constrained optimization. Software creators have to work with assumptions on the extremity of scale that they are willing to countenance, and sometimes problems will be outside the scope envisaged. ```{r ex10, echo=TRUE} # Start MUST be feasible i.e. on or within bounds anlshob1b <- nls(frms, start=sts1, data=weeddf, lower=c(0,0,0), upper=c(2,6,3), algorithm='port') pnls(anlshob1b) # check the answer (short form) # nlsLM seems NOT to work with bounds in this example anlsLM1b <- nlsLM(frms, start=sts1, data=weeddf, lower=c(0,0,0), upper=c(2,6,3)) pnls(anlsLM1b) # also no warning if starting out of bounds, but gets good answer!! st4<-c(c1=4, c2=4, c3=4) anlsLMob <- nlsLM(frms, start=st4, data=weeddf, lower=c(0,0,0), upper=c(2,6,3)) pnls(anlsLMob) # Try nlsr::nlxb() anlx1b <- nlxb(frms, start=sts1, data=weeddf, lower=c(0,0,0), upper=c(2,6,3)) pshort(anlx1b) ``` ## Functional specification of nonlinear least squares problems We illustrate how to solve nonlinear least squares problems where the set of functions to be squared is provided by an R function, as is the Jacobian matrix. Note that `nlsr::nlfb()` gets this information from the object returned by the `jacfn` argument in the "gradient" attribute of that object. This is a programming artifice to allow a simplification of the `nlxb()` code. The example is the well-known Banana shaped valley of @Rosenbrock1960. ```{r exrosen, eval=TRUE} frres<-function(x) { ## Rosenbrock as residuals x1 <- x[1] x2 <- x[2] res<-c( 10 * (x2 - x1 * x1), (1 - x1) ) } frjac<-function(x) { ## Rosenbrock residual derivatives matrix x1 <- x[1] x2 <- x[2] J<-matrix(NA, nrow=2, ncol=2) J[1,1]<- -20*x1 J[1,2]<- 10 J[2,1]<- -1 J[2,2]<- 0 attr(J,"gradient")<-J # NEEDED for nlfb() J } require(minpack.lm) require(nlsr) strt<-c(-1.2,1) rosnlf<-nlfb(strt, resfn=frres, jacfn=frjac) print(rosnlf) rosnlm<-nls.lm(par=strt, fn=frres, jac=frjac) summary(rosnlm) ``` I find it awkward to solve problems specified this way with `nls()` and do not believe it is worth pursuing that topic. ## Weighted nonlinear regression Regression attempts to explain a *dependent* (or *predicted*) variable as a function of *independent* (or *explanatory* or *predictor*) variables and estimated parameters. The estimation method commonly used is minimization of the sum of squared residuals. However, these residuals may be weighted. In the tools available in **R** for nonlinear least squares, weights can be specified which multiply the **squared** residuals, and we need a weight for each term in the sum. Each residual is therefore multiplied by the square root of its respective weight. Typically, the weights are given as a vector of numbers. A popular choice is the reciprocal of a measure of the variance of each data observation. Replicated data allows for sample variances to be used; otherwise we need some proxy. The concept is that the weights are proportional to our belief that the observations are correct. ### Static weights Let us demonstrate with weights that are simply the reciprocal of the independent variable. This may not make sense for statistical modeling, but it lets us see that the `residuals()` function returns weighted residuals. ```{r staticwts} wts <- 1/weeddf$tt # wts are reciprocal of time value tnlx1w<-try(nlxb(formula=frmt, start=stt1, data=weeddf, weights=wts)) # pshort(tnlx1) print(tnlx1w) ct<-as.list(tnlx1w$coefficients) # Need a list to use ct$... in next line cat("exp(xmid/scal)=",exp(ct$xmid/ct$scal),"\n") cat("\n") rtnlx1w<-residuals(tnlx1w) print(rtnlx1w) cat("explicit sumsquares =", sum(rtnlx1w^2),"\n") ``` ### Weights that are functions of the model parameters Possible choices for weights include the reciprocal of the squared residuals or squared fitted values. Unfortunately, such quantities depend on the parameter values. Consider our Logistic example. We want to predict $y$ with $$ fitted(b_1, b_2, b_3, t) = b_1 / (1 + b_2 * exp(- b_3 * t)) $$ Thus the raw residuals are $$ rawres(b_1, b_2, b_3, t) = y - fitted(b_1, b_2, b_3, t) $$ But if the weights are $(1/ fitted(b_1, b_2, b_3, t)$, then the residuals for which the sum of squares is to be calculated are $$ resid(b_1, b_2, b_3, t) = rawres(b_1, b_2, b_3, t) * (1/fitted(b_1, b_2, b_3, t) $$ $$ = y/(b_1 / (1 + b_2 * exp(- b_3 * t))) - 1 $$ However, we still need to keep in mind that the user needs to use the formula for the fitted values to get predictions. While the software needs to minimize the weighted (hence modified) sum of squares, it also must keep track of the fitted/predicted values. `minpack.lm` offers a function `wfct()` that can be used within the argument supplied as `weights` in the call to `nlsLM()` or `nls()`, though it usually fails with `nlsr::nlxb()`. The action of this function was surprising to me, as `wfct()` takes an expression as its argument (and may itself be part of an expression in the argument after the `=` sign). If the argument to `wfct` contains `fitted`, `resid` or `error`, then the call to `nlsLM` or `nls` will trigger a further call to one of these functions, evaluate the `weights` and then run again to compute the solution. This can be seen by setting `trace = TRUE`. The process can be verified explicitly for the example used in the `wfct` documentation. Unfortunately, `wfct()` appears to crash `knitr` or `rmarkdown`, so we have evaluated the example outside of Rstudio and included the output below. ```{r puromycinwfct, echo=TRUE, eval=FALSE} library(minpack.lm) library(nlsr) # for pnls Treated <- Puromycin[Puromycin$state == "treated", ] # First get raw estimates wtt3nlm0<-nlsLM(rate ~ Vm * conc/(K + conc), data = Treated, trace=TRUE, start = c(Vm = 200, K = 0.05)) fit0<-fitted(wtt3nlm0) # Static wts using fit0 wtt3nlm0s<-nlsLM(rate ~ Vm * conc/(K + conc), data = Treated, trace=TRUE, start = c(Vm = 200, K = 0.05), weights = wfct(1/fit0^2)) pnls(wtt3nlm0s) # and run directly, noting the 2 phase operation wtt3nlm<-nlsLM(rate ~ Vm * conc/(K + conc), data = Treated, trace=TRUE, start = c(Vm = 200, K = 0.05), weights = wfct(1/fitted^2)) pnls(wtt3nlm) cat("weights from wtt3nlm\n") as.numeric(wtt3nlm$weights) ``` ``` It. 0, RSS = 1636.59, Par. = 200 0.05 It. 1, RSS = 1205.62, Par. = 211.157 0.0616271 It. 2, RSS = 1195.57, Par. = 212.511 0.0638418 It. 3, RSS = 1195.45, Par. = 212.666 0.0640939 It. 4, RSS = 1195.45, Par. = 212.682 0.0641186 It. 5, RSS = 1195.45, Par. = 212.684 0.064121 It. 0, RSS = 0.22811, Par. = 200 0.05 It. 1, RSS = 0.227222, Par. = 200.913 0.0494747 It. 2, RSS = 0.227221, Par. = 200.847 0.0494303 It. 3, RSS = 0.227221, Par. = 200.842 0.049426 It. 4, RSS = 0.227221, Par. = 200.841 0.0494256 wtt3nlm0s -- ss= 0.2272213 : Vm = 200.841 K = 0.04942558; 4 itns It. 0, RSS = 1636.59, Par. = 200 0.05 It. 1, RSS = 1205.62, Par. = 211.157 0.0616271 It. 2, RSS = 1195.57, Par. = 212.511 0.0638418 It. 3, RSS = 1195.45, Par. = 212.666 0.0640939 It. 4, RSS = 1195.45, Par. = 212.682 0.0641186 It. 5, RSS = 1195.45, Par. = 212.684 0.064121 It. 0, RSS = 0.22811, Par. = 200 0.05 It. 1, RSS = 0.227222, Par. = 200.913 0.0494747 It. 2, RSS = 0.227221, Par. = 200.847 0.0494303 It. 3, RSS = 0.227221, Par. = 200.842 0.049426 It. 4, RSS = 0.227221, Par. = 200.841 0.0494256 wtt3nlm -- ss= 0.2272213 : Vm = 200.841 K = 0.04942558; 4 itns Show in New Window weights from wtt3nlm [1] 3.910941e-04 3.910941e-04 9.460635e-05 9.460635e-05 5.539227e-05 5.539227e-05 [7] 3.687173e-05 3.687173e-05 2.745956e-05 2.745956e-05 2.475956e-05 2.475956e-05 ``` The documentation of `wfct` suggests that it can also use the name of the response (dependent) variable or the name of the predictor (independent) variable. However, some models have more than one independent variable, and I have not explored what happens in such cases. It would be helpful if `nlsr` had the capability to include functional weights. Duncan Murdoch suggested a patch that allows this. We modify the `nlfb()` routine to detect and act on functional weights. ```{r formulawts} library(nlsr) Treated <- Puromycin[Puromycin$state == "treated", ] # "regression" formula w1frm <- rate ~ Vm * conc/(K + conc) # Explicit dynamic weighted nlxb w3frm <- rate/(Vm * conc/(K + conc)) ~ 1 dynweighted <- nlxb(w3frm, data = Treated, start = c(Vm = 201.003, K = 0.04696), trace=TRUE, control=list(prtlvl=0)) pshort(dynweighted) # Compute the model from the ORIGINAL model (w1frm) using parameters from dynweighted dyn0 <- nlxb(w1frm, start=coefficients(dynweighted), data=Treated, control=list(jemax=0)) wfdyn0 <- 1/fitted(dyn0)^2 # weights print(as.numeric(wfdyn0)) pshort(dyn0) # Shows sumsquares without weights, but computes weights we used formulaweighted <- nlxb(w1frm, data = Treated, start = c(Vm = 201.003, K = 0.04696), weights = ~ 1/fitted^2, trace=TRUE, control=list(prtlvl=0)) pshort(formulaweighted) wtsfromfits<-1/fitted(formulaweighted)^2 print(as.numeric(wtsfromfits)) print(as.numeric(formulaweighted$weights)) ``` There are differences in the (weighted) sums of squares. The **dynweighted** case minimizes the residuals weighted by the actual CURRENT fits as expressed by the formula. The **formulaweighted** case computes the weights at each iteration from the previous (or original) fit. It is a form of iteratively weighted least squares. The case **wtt3nlm** above uses `nlsLM` with weights specified in `wfct()` and computes the weights once from an unweighted model, then runs a second, statically weighted calculation. Keeping track of these differences requires a lot of care, and I advise documenting what is done carefully. ## Ongoing efforts The examples above show how to use package `nlsr` as it exists at the end of December 2022. The story, however, is not finished. It would be very nice to be able to identify and work with models that are partially linear. Solving such models is possible with the `algorithm="plinear"` option of `nls()`, but identifying which parameters enter linearly is not easy for most users. Moreover, the specification of the model to the `VARPRO` solver that is called by the `plinear` option is inconsistent with the general specification used by `nls()` and other nonlinear least squares tools for R. There are also a number of possible tweaks to details in `nlsr` and efforts to include such improvements are a continuing interest. I welcome communication and collaboration to continue the improvement of the package and R in general. ## References