2

I'm trying to calculate the confidence interval from my nls model. And I tried the same code as this checked answer: How to calculate confidence intervals for Nonlinear Least Squares in r?

But I get a strange error:

Error in eval(form[[3]]) : object 'a' not found
4.
eval(form[[3]])
3.
eval(form[[3]])
2.
predFit.nls(gloss.nls, newdata = data.frame(stimulus = seq(0, 
1, by = 0.1)), interval = "confidence", level = 0.9)
1.
predFit(gloss.nls, newdata = data.frame(stimulus = seq(0, 1, 
by = 0.1)), interval = "confidence", level = 0.9)

I nearly use the same code as the answer above, only differing in data:

gloss.nls <- nls(
                normP ~ a[1]*stimulus^3+a[2]*stimulus,
                data = data.mlds %>% filter(overall == TRUE),
                start = list(a=c(0.4,0.6))
                )

predFit(gloss.nls, newdata = data.frame(stimulus=seq(0, 1, by = 0.1)), interval = "confidence", level= 0.9)

And Here's my data:

id  rank  stimulus  pscale  normP  overall
0   1   0.000   0.0000000   0.00000000  TRUE
0   2   0.125   0.3151757   0.05889716  TRUE
0   3   0.250   0.9225827   0.17240385  TRUE
0   4   0.375   1.4164383   0.26469110  TRUE
0   5   0.500   1.7400011   0.32515557  TRUE
0   6   0.625   2.3531344   0.43973235  TRUE
0   7   0.750   3.1662257   0.59167546  TRUE
0   8   0.875   4.3538122   0.81360082  TRUE
0   9   1.000   5.3512879   1.00000000  TRUE
1   1   0.000   0.0000000   0.00000000  FALSE
3
  • If I copy the same code from that reference answer, it works well on my PC.
    – T3-Tan
    Dec 16, 2021 at 9:29
  • What does a do? You've referenced an external vector in your formula but I'm not sure what that adds? Dec 16, 2021 at 9:45
  • @AndyBaxter It's just a start point for the nls calculation." start = list(a=c(0.4,0.6)"
    – T3-Tan
    Dec 17, 2021 at 0:20

1 Answer 1

2

Short answer

Try

a <- c(0.4,0.6)
predFit(gloss.nls, newdata = data.frame(stimulus=seq(0, 1, by = 0.1)), interval = "confidence", level= 0.9)

Long answer

First, note that your model is linear in parameters, so you can just estimate the model in plain ols, where confidence intervals are straightforward.

library(tidyverse)
gloss.lm  <-  lm(normP ~ I(stimulus^3)+stimulus,
                data = data.mlds %>% filter(overall == TRUE)  )
predict(gloss.lm, newdata = data.frame(stimulus=seq(0, 1, by = 0.1)), interval = "confidence", level= 0.9)
           fit         lwr        upr
1  0.005554547 -0.02791979 0.03902889
2  0.061136954  0.03572392 0.08654999
3  0.119410056  0.09931945 0.13950067
4  0.183064551  0.16435972 0.20176938
5  0.254791132  0.23459593 0.27498634
6  0.337280497  0.31518226 0.35937873
7  0.433223342  0.41047149 0.45597519
8  0.545310361  0.52354420 0.56707652
9  0.676232250  0.65548488 0.69697962
10 0.828679707  0.80399326 0.85336616
11 1.005343426  0.96738940 1.04329745

If you insist on estimating the model using nonlinear least squares, then

gloss.nls <-  nls(normP ~ a[1]*stimulus^3+a[2]*stimulus,
                data = data.mlds %>% filter(overall == TRUE) ,
                start=list(a=c(.5, .5)) )

Annoyingly, predict.nls does not seem have confidence interval calculation, so this does not produce confidence intervals.

predict(gloss.nls, newdata = data.frame(stimulus=seq(0, 1, by = 0.1)), interval = "confidence", level= 0.9)
 [1] 0.00000000 0.05704647 0.11672200 0.18165566 0.25447650 0.33781360
 [7] 0.43429601 0.54655279 0.67721301 0.82890574 1.00426003

Luckily, investr::predFit has an implementation for confidence interval calculation.

library(investr)
predFit(gloss.nls, interval='prediction', newdata = data.frame(stimulus=seq(0, 1, by = 0.1), confidence=.9))

... but this returns an error (which you encountered in your question).

I did not dig too deep into predFit.nls code but it seems that it predFit silently runs gloss.nls$call in the background, and if it does not find everything it needs, it returns a weird error. It is enough to create an object into the namespace with the same shape as a to resolve the error.

a <- coef(gloss.nls)
investr::predFit(gloss.nls, interval='prediction', newdata = data.frame(stimulus=seq(0, 1, by = 0.1), confidence=.9))
             fit          lwr        upr
 [1,] 0.00000000 -0.050130071 0.05013007
 [2,] 0.05704647  0.006916398 0.10717654
 [3,] 0.11672200  0.066591929 0.16685207
 [4,] 0.18165566  0.131525585 0.23178573
 [5,] 0.25447650  0.204346430 0.30460657
 [6,] 0.33781360  0.287683526 0.38794367
 [7,] 0.43429601  0.384165935 0.48442608
 [8,] 0.54655279  0.496422720 0.59668286
 [9,] 0.67721301  0.627082944 0.72734309
[10,] 0.82890574  0.778775670 0.87903581
[11,] 1.00426003  0.954129960 1.05439010

Interestingly, the values in a do not make any difference. Try, e.g. a <- c(7500,-100) and you will get the same results. This might be a bug in investr?

a <- c(7500,-100)
predFit(gloss.nls, interval='prediction', newdata = data.frame(stimulus=seq(0, 1, by = 0.1), confidence=.9))
             fit          lwr        upr
 [1,] 0.00000000 -0.050130071 0.05013007
 [2,] 0.05704647  0.006916398 0.10717654
 [3,] 0.11672200  0.066591929 0.16685207
 [4,] 0.18165566  0.131525585 0.23178573
 [5,] 0.25447650  0.204346430 0.30460657
 [6,] 0.33781360  0.287683526 0.38794367
 [7,] 0.43429601  0.384165935 0.48442608
 [8,] 0.54655279  0.496422720 0.59668286
 [9,] 0.67721301  0.627082944 0.72734309
[10,] 0.82890574  0.778775670 0.87903581
[11,] 1.00426003  0.954129960 1.05439010

Data:

data.mlds <- structure(list(id = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L),
    rank = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 1L), stimulus = c(0,
    0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1, 0), pscale = c(0,
    0.3151757, 0.9225827, 1.4164383, 1.7400011, 2.3531344, 3.1662257,
    4.3538122, 5.3512879, 0), normP = c(0, 0.05889716, 0.17240385,
    0.2646911, 0.32515557, 0.43973235, 0.59167546, 0.81360082,
    1, 0), overall = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,
    TRUE, TRUE, FALSE)), row.names = c(NA, -10L), class = "data.frame")
2
  • Thanks for your time doing these test. It works! It might to be a little bug. Maybe I should ask the author.
    – T3-Tan
    Dec 17, 2021 at 0:13
  • And sorry I'm new here. I don't have enough reputation to cast a vote for your great answer.😥
    – T3-Tan
    Dec 17, 2021 at 0:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.