2
$\begingroup$

I am trying to define a range of uncertainty associated with each estimate of the parameters of a nonlinear model:

mod <- nls(y ~ ifelse(t <= tauLA, hEA*exp(-1/2*((t - tauEA)/muEA)^2),
                       hEA*exp(-1/2*((t - tauEA)/muEA)^2) +  hLA*exp(-1/2*(log(((t - tauLA)/muLA)^2)/sigmaLA)^2)),
           data=dat,
           start=list(hEA=0.1,
                      hLA=0.1,
                      tauEA=148,
                      tauLA=160,
                      muEA=5,
                      muLA=9,
                      sigmaLA=3.1), algorithm="port")

The result is:

Parameters:
          Estimate Std. Error t value Pr(>|t|)    
hEA      2.300e-03  6.952e-04   3.308  0.00113 ** 
hLA      8.207e-03  9.978e-03   0.823  0.41182    
tauEA    1.410e+02  2.935e+01   4.804  3.2e-06 ***
tauLA    2.458e+02  1.158e+02   2.122  0.03513 *  
muEA     2.592e+01  2.090e+01   1.240  0.21650    
muLA    -3.223e+01  1.524e+02  -0.212  0.83271    
sigmaLA -7.949e-01  3.866e+00  -0.206  0.83735    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.003139 on 186 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)

As proposed in Why is my nonlinear least squares confidence band so wide?, I used a standard bootstrapping approach that bootstraps the model residuals to calculate the range of uncertainty:

dat$res <- residuals(mod)
dat$fit <- fitted(mod)

fn_boot_p <- function(data, index) {
  data$y <- data$fit + data[index, c("res")]

  tryCatch(coef(nls(y ~ ifelse(t <= tauLA, hEA*exp(-1/2*((t - tauEA)/muEA)^2),
                                hEA*exp(-1/2*((t - tauEA)/muEA)^2) +  hLA*exp(-1/2*(log(((t - tauLA)/muLA)^2)/sigmaLA)^2)),
                    data=data,
                    start=list(hEA=0.1,
                      hLA=0.1,
                      tauEA=148,
                      tauLA=160,
                      muEA=5,
                      muLA=9,
                      sigmaLA=3.1), algorithm="port")),
           error = function(e) c(t = NA))
}

boot_param <- boot(dat, fn_boot_p, R = 1000)

mean(is.na(boot_param$t[, 1]))
    0.772

But the model fitting fails in about 77% of bootstrap replications. Is there a way to correct this ? Should I use another method to estimate uncertainty in fitted values ?

$\endgroup$
  • $\begingroup$ What do you mean when you say that the model fitting fails in 77% of the bootstrap samples (I think you mean bootstrap replications). $\endgroup$ – Michael R. Chernick May 7 '18 at 23:43
  • $\begingroup$ Yes, I mean bootstrap replications. I edited the text accordingly. $\endgroup$ – Nell May 8 '18 at 2:21
  • 2
    $\begingroup$ Then what do you mean by "model fitting failing"? $\endgroup$ – Michael R. Chernick May 8 '18 at 4:08
  • $\begingroup$ I would use the coefficients from the model fit as starting values for the bootstrap fits. Anyway, you have 7 parameters and less than 200 observations. I don't find it surprising that you have problems fitting this. Also, please answer Michael's questions. What are the errors nls is throwing? Is it "just" convergence problems? $\endgroup$ – Roland May 8 '18 at 6:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.