R help - Aug 2020 - Change how minpack.lm:::summary.nls.lm calculates degrees of freedom

Dear R users,
I want to modify how the degrees of freedom are calculated from the
summary function of the package minpack.lm

My first thought was to replicate how minpack.lm:::summary.nls.lm
works, and modify the line of the code that is relevant for my task.
However, for some reason I am not able to use the code contained in
this function to perform this operation.

Here is a reproducible example.

First, I run a NLLS regression using minpack.lm::nls.lm

# From example 1 of the help page of ?minpack.lm::nls.lm
library(minpack.lm)
x <- seq(0,5,length=100)
getPred <- function(parS, xx) parS$a * exp(xx * parS$b) + parS$c
pp <- list(a=9,b=-1, c=6)
simDNoisy <- getPred(pp,x) + rnorm(length(x),sd=.1)
residFun <- function(p, observed, xx) observed - getPred(p,xx)
parStart <- list(a=3,b=-.001, c=1)
nls.out <- nls.lm(par=parStart, fn = residFun, observed = simDNoisy,
xx = x, control = nls.lm.control(nprint=1))
summary(nls.out)

Now, by running minpack.lm:::summary.nls.lm in the console, I get the
following
> minpack.lm:::summary.nls.lm
function (object, ...)
{
    param <- coef(object)
    pnames <- names(param)
    ibb <- chol(object$hessian)
    ih <- chol2inv(ibb)
    p <- length(param)
    rdf <- length(object$fvec) - p
    resvar <- deviance(object)/rdf
    se <- sqrt(diag(ih) * resvar)
    names(se) <- pnames
    tval <- param/se
    param <- cbind(param, se, tval, 2 * pt(abs(tval), rdf, lower.tail =
FALSE))
    dimnames(param) <- list(pnames, c("Estimate", "Std.
Error",
        "t value", "Pr(>|t|)"))
    ans <- list(residuals = object$fvec, sigma = sqrt(object$deviance/rdf),
        df = c(p, rdf), cov.unscaled = ih, info = object$info,
        niter = object$niter, stopmess = object$message, coefficients = param)
    class(ans) <- "summary.nls.lm"
    ans
}

Specifically, my task requires to modify the object p of this
function, say I want to set p <- 2.

To this purpose, I replicate what the summary function does using my
object (nls.out), however an error is returned at line
three:
> param <- coef(nls.out)
> pnames <- names(param)
> ibb <- chol(nls.out$hessian)
Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x),  :
  'data' must be of a vector type, was 'NULL'

The reason of the error is that there is no nls.out$hessian in nls.out.

I guess that I am missing something about how summary functions work,
and this is the reason why I cannot use the code from
minpack.lm:::summary.nls.lm outside the minpack.lm environment.

Does anyone have any hints on how to proceed? Any other way of dealing
with this issue will be equally appreciated.

Thank you,
Valerio

Hi Valerio,
I did a copy-paste on your reproducible example and I had no problem with
chol(nls.out$hessian).
In addition to summary() you can look at str() to display the structure of
any R object.
> str(nls.out)
List of 9
 $ par     :List of 3
  ..$ a: num 8.99
  ..$ b: num -1.01
  ..$ c: num 6.02
 $ hessian : num [1:3, 1:3] 10.3 43.2 19.9 43.2 382.2 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "a" "b" "c"
  .. ..$ : chr [1:3] "a" "b" "c"
 $ fvec    : num [1:100] 0.0232 -0.105 -0.1332 0.0388 0.1482 ...
 $ info    : int 1
 $ message : chr "Relative error in the sum of squares is at most
`ftol'."
 $ diag    :List of 3
  ..$ a: num 9.98
  ..$ b: num 88.6
  ..$ c: num 10
 $ niter   : int 8
 $ rsstrace: num [1:9] 1966.2 327.2 104.8 53.9 33.2 ...
 $ deviance: num 1.06
 - attr(*, "class")= chr "nls.lm"

Also
> nls.out$hessian
         a         b         c
a 10.26361  43.17086  19.89616
b 43.17086 382.17773 166.43747
c 19.89616 166.43747 100.00000

HTH,
Eric



On Wed, Aug 19, 2020 at 12:17 PM Valerio Leone Sciabolazza <
sciabolazza at gmail.com> wrote:
> Dear R users,
> I want to modify how the degrees of freedom are calculated from the
> summary function of the package minpack.lm
>
> My first thought was to replicate how minpack.lm:::summary.nls.lm
> works, and modify the line of the code that is relevant for my task.
> However, for some reason I am not able to use the code contained in
> this function to perform this operation.
>
> Here is a reproducible example.
>
> First, I run a NLLS regression using minpack.lm::nls.lm
>
> # From example 1 of the help page of ?minpack.lm::nls.lm
> library(minpack.lm)
> x <- seq(0,5,length=100)
> getPred <- function(parS, xx) parS$a * exp(xx * parS$b) + parS$c
> pp <- list(a=9,b=-1, c=6)
> simDNoisy <- getPred(pp,x) + rnorm(length(x),sd=.1)
> residFun <- function(p, observed, xx) observed - getPred(p,xx)
> parStart <- list(a=3,b=-.001, c=1)
> nls.out <- nls.lm(par=parStart, fn = residFun, observed = simDNoisy,
> xx = x, control = nls.lm.control(nprint=1))
> summary(nls.out)
>
> Now, by running minpack.lm:::summary.nls.lm in the console, I get the
> following
> > minpack.lm:::summary.nls.lm
> function (object, ...)
> {
>     param <- coef(object)
>     pnames <- names(param)
>     ibb <- chol(object$hessian)
>     ih <- chol2inv(ibb)
>     p <- length(param)
>     rdf <- length(object$fvec) - p
>     resvar <- deviance(object)/rdf
>     se <- sqrt(diag(ih) * resvar)
>     names(se) <- pnames
>     tval <- param/se
>     param <- cbind(param, se, tval, 2 * pt(abs(tval), rdf, lower.tail
> FALSE))
>     dimnames(param) <- list(pnames, c("Estimate", "Std.
Error",
>         "t value", "Pr(>|t|)"))
>     ans <- list(residuals = object$fvec, sigma =
sqrt(object$deviance/rdf),
>         df = c(p, rdf), cov.unscaled = ih, info = object$info,
>         niter = object$niter, stopmess = object$message, coefficients >
param)
>     class(ans) <- "summary.nls.lm"
>     ans
> }
>
> Specifically, my task requires to modify the object p of this
> function, say I want to set p <- 2.
>
> To this purpose, I replicate what the summary function does using my
> object (nls.out), however an error is returned at line three:
> > param <- coef(nls.out)
> > pnames <- names(param)
> > ibb <- chol(nls.out$hessian)
> Error in array(x, c(length(x), 1L), if (!is.null(names(x)))
> list(names(x),  :
>   'data' must be of a vector type, was 'NULL'
>
> The reason of the error is that there is no nls.out$hessian in nls.out.
>
> I guess that I am missing something about how summary functions work,
> and this is the reason why I cannot use the code from
> minpack.lm:::summary.nls.lm outside the minpack.lm environment.
>
> Does anyone have any hints on how to proceed? Any other way of dealing
> with this issue will be equally appreciated.
>
> Thank you,
> Valerio
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>
	[[alternative HTML version deleted]]

Eric,
you are right! If I run the code on a different instance, it works.
There must be something in my old R environment that messed this up.
Thank you very much for checking this out.
Best,
Valerio

On Wed, Aug 19, 2020 at 11:36 AM Eric Berger <ericjberger at gmail.com>
wrote:
>
> Hi Valerio,
> I did a copy-paste on your reproducible example and I had no problem with
chol(nls.out$hessian).
> In addition to summary() you can look at str() to display the structure of
any R object.
>
> > str(nls.out)
>
> List of 9
>  $ par     :List of 3
>   ..$ a: num 8.99
>   ..$ b: num -1.01
>   ..$ c: num 6.02
>  $ hessian : num [1:3, 1:3] 10.3 43.2 19.9 43.2 382.2 ...
>   ..- attr(*, "dimnames")=List of 2
>   .. ..$ : chr [1:3] "a" "b" "c"
>   .. ..$ : chr [1:3] "a" "b" "c"
>  $ fvec    : num [1:100] 0.0232 -0.105 -0.1332 0.0388 0.1482 ...
>  $ info    : int 1
>  $ message : chr "Relative error in the sum of squares is at most
`ftol'."
>  $ diag    :List of 3
>   ..$ a: num 9.98
>   ..$ b: num 88.6
>   ..$ c: num 10
>  $ niter   : int 8
>  $ rsstrace: num [1:9] 1966.2 327.2 104.8 53.9 33.2 ...
>  $ deviance: num 1.06
>  - attr(*, "class")= chr "nls.lm"
>
> Also
> > nls.out$hessian
>          a         b         c
> a 10.26361  43.17086  19.89616
> b 43.17086 382.17773 166.43747
> c 19.89616 166.43747 100.00000
>
> HTH,
> Eric
>
>
>
> On Wed, Aug 19, 2020 at 12:17 PM Valerio Leone Sciabolazza <sciabolazza
at gmail.com> wrote:
>>
>> Dear R users,
>> I want to modify how the degrees of freedom are calculated from the
>> summary function of the package minpack.lm
>>
>> My first thought was to replicate how minpack.lm:::summary.nls.lm
>> works, and modify the line of the code that is relevant for my task.
>> However, for some reason I am not able to use the code contained in
>> this function to perform this operation.
>>
>> Here is a reproducible example.
>>
>> First, I run a NLLS regression using minpack.lm::nls.lm
>>
>> # From example 1 of the help page of ?minpack.lm::nls.lm
>> library(minpack.lm)
>> x <- seq(0,5,length=100)
>> getPred <- function(parS, xx) parS$a * exp(xx * parS$b) + parS$c
>> pp <- list(a=9,b=-1, c=6)
>> simDNoisy <- getPred(pp,x) + rnorm(length(x),sd=.1)
>> residFun <- function(p, observed, xx) observed - getPred(p,xx)
>> parStart <- list(a=3,b=-.001, c=1)
>> nls.out <- nls.lm(par=parStart, fn = residFun, observed = simDNoisy,
>> xx = x, control = nls.lm.control(nprint=1))
>> summary(nls.out)
>>
>> Now, by running minpack.lm:::summary.nls.lm in the console, I get the
following
>> > minpack.lm:::summary.nls.lm
>> function (object, ...)
>> {
>>     param <- coef(object)
>>     pnames <- names(param)
>>     ibb <- chol(object$hessian)
>>     ih <- chol2inv(ibb)
>>     p <- length(param)
>>     rdf <- length(object$fvec) - p
>>     resvar <- deviance(object)/rdf
>>     se <- sqrt(diag(ih) * resvar)
>>     names(se) <- pnames
>>     tval <- param/se
>>     param <- cbind(param, se, tval, 2 * pt(abs(tval), rdf,
lower.tail = FALSE))
>>     dimnames(param) <- list(pnames, c("Estimate",
"Std. Error",
>>         "t value", "Pr(>|t|)"))
>>     ans <- list(residuals = object$fvec, sigma =
sqrt(object$deviance/rdf),
>>         df = c(p, rdf), cov.unscaled = ih, info = object$info,
>>         niter = object$niter, stopmess = object$message, coefficients =
param)
>>     class(ans) <- "summary.nls.lm"
>>     ans
>> }
>>
>> Specifically, my task requires to modify the object p of this
>> function, say I want to set p <- 2.
>>
>> To this purpose, I replicate what the summary function does using my
>> object (nls.out), however an error is returned at line three:
>> > param <- coef(nls.out)
>> > pnames <- names(param)
>> > ibb <- chol(nls.out$hessian)
>> Error in array(x, c(length(x), 1L), if (!is.null(names(x)))
list(names(x),  :
>>   'data' must be of a vector type, was 'NULL'
>>
>> The reason of the error is that there is no nls.out$hessian in nls.out.
>>
>> I guess that I am missing something about how summary functions work,
>> and this is the reason why I cannot use the code from
>> minpack.lm:::summary.nls.lm outside the minpack.lm environment.
>>
>> Does anyone have any hints on how to proceed? Any other way of dealing
>> with this issue will be equally appreciated.
>>
>> Thank you,
>> Valerio
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.