How to get p value and confidence intervals for nls functions?

Question

I have 2 questions.

1) How can I have p.value for my 2 functions? My hypothesis is that I have a correlation between my function and my data.

2) How can I have a confidence intervals for my 2 functions?

library(ggplot2)
g <- function (x, a,b,c) a * (1-exp(-(x-c)/abs(b)))
X1 <- c(129.08,109.92,85.83,37.72)
Y1 <- c(0.7,0.5,0.39,-1.36)
dt1 <- data.frame(x1=X1,y1=Y1)
model1 <- nls(Y1 ~ g(X1, a, b, c), 
          start = list(a=0.5, b=60, c=50),control=nls.control(maxiter = 200))

ggplot(data = dt1,aes(x = x1, y = y1)) + 
     theme_bw() + geom_point() + 
     geom_smooth(data=dt1, method="nls", formula=y~g(x, a, b, c),
       se=F, start=list(a=0.5, b=60, c=50))


f <- function (x, a, b, c) a*(x^2)+b*x+c   
X2 <- c(589.62,457.92,370.16,295.98,243.99,199.07,159.91,142.63,
124.15, 101.98, 87.93, 83.16, 82.2, 74.48, 47.68, 37.51, 31,
27.9, 21.24,18.28)
Y2 <- c(0.22,0.37,0.49,0.65,0.81,0.83,1,0.81,0.65,0.44,0.55,0.63,
0.65,0.55,0.37,0.32,0.27,0.22,0.17,0.14)
dt2 <- data.frame(x2=X2,y2=Y2)
model2 <- nls(Y2 ~ f(X2, a, b, c), 
           start = list(a=-1, b=3, c=0),control=nls.control(maxiter = 200))
ggplot(data = dt2,aes(x = x2, y = y2)) + 
      theme_bw() + geom_point() + 
      geom_smooth(data=dt2, method="nls", formula=y~f(x, a, b, c),
       se=F, start=list(a=-1, b=3, c=0))

Thank you in advance

Answer 1

1. - You could try (this is an approximation)

library(nls2)  
summary(as.lm(model))  

• You can obtain a p-value for all parameters used in your model using

  summary(model)

• You can get p values for a model by comparing it to another ("nested") model using

  anova(model1, model2)

  where model 2 is a simplified version of model 1 (it is your null hypothesis)

• You can use methods such a bootstrapping, to get a measure of the probability of fit of your complete model.

  2.

• You can possibly get full model confidence interval using (this is an approximation)

  library(nls2) predict(as.lm(model2), interval = "confidence")

• You can obtain the confidence interval of the parameters using

  confint(model)

• You can get more information about these parameter intervals using

  profile(model)

  plot(profile(model))

• You can obtain the pair wise confidence interval for two of your parameters (for both plotting and to get the matrix) using

  ellipse.nls(model)

Answer 2

Regarding confidence intervals, other answers here seem to have issues with the use of functions (as.lm.nls, as.proto.list) that for some reason are not defined for some users (like me). After some surfing, I found an answer that works for me, requiring only the MASS package. At the urging of @etov, I am posting the answer I found here. It is originally from https://www.r-bloggers.com/predictnls-part-1-monte-carlo-simulation-confidence-intervals-for-nls-models/ and appears to be by someone named Andrej who uses the handle anspiess. This function by Andrej, in his words, "takes an nls object, extracts the variables/parameter values/parameter variance-covariance matrix, creates an “augmented” covariance matrix (with the variance/covariance values from the parameters and predictor variables included, the latter often being zero), simulates from a multivariate normal distribution (using mvrnorm of the ‘MASS’ package), evaluates the function (object$call$formula) on the values and finally collects statistics". So it is a Monte-Carlo-based method of getting confidence intervals for an nls model. His code:

predictNLS <- function(
object, 
newdata,
level = 0.95, 
nsim = 10000,
...
)
{
  require(MASS, quietly = TRUE)

  ## get right-hand side of formula
  RHS <- as.list(object$call$formula)[[3]]
  EXPR <- as.expression(RHS)

  ## all variables in model
  VARS <- all.vars(EXPR)

  ## coefficients
  COEF <- coef(object)

  ## extract predictor variable    
  predNAME <- setdiff(VARS, names(COEF))  

  ## take fitted values, if 'newdata' is missing
  if (missing(newdata)) {
    newdata <- eval(object$data)[predNAME]
    colnames(newdata) <- predNAME
  }

  ## check that 'newdata' has same name as predVAR
  if (names(newdata)[1] != predNAME) stop("newdata should have name '", predNAME, "'!")

  ## get parameter coefficients
  COEF <- coef(object)

  ## get variance-covariance matrix
  VCOV <- vcov(object)

  ## augment variance-covariance matrix for 'mvrnorm' 
  ## by adding a column/row for 'error in x'
  NCOL <- ncol(VCOV)
  ADD1 <- c(rep(0, NCOL))
  ADD1 <- matrix(ADD1, ncol = 1)
  colnames(ADD1) <- predNAME
  VCOV <- cbind(VCOV, ADD1)
  ADD2 <- c(rep(0, NCOL + 1))
  ADD2 <- matrix(ADD2, nrow = 1)
  rownames(ADD2) <- predNAME
  VCOV <- rbind(VCOV, ADD2) 

  ## iterate over all entries in 'newdata' as in usual 'predict.' functions
  NR <- nrow(newdata)
  respVEC <- numeric(NR)
  seVEC <- numeric(NR)
  varPLACE <- ncol(VCOV)   

  ## define counter function
  counter <- function (i) 
  {
    if (i%%10 == 0) 
      cat(i)
    else cat(".")
    if (i%%50 == 0) 
      cat("\n")
    flush.console()
  }

  outMAT <- NULL 

  for (i in 1:NR) {
    counter(i)

    ## get predictor values and optional errors
    predVAL <- newdata[i, 1]
    if (ncol(newdata) == 2) predERROR <- newdata[i, 2] else predERROR <- 0
    names(predVAL) <- predNAME  
    names(predERROR) <- predNAME  

    ## create mean vector for 'mvrnorm'
    MU <- c(COEF, predVAL)

    ## create variance-covariance matrix for 'mvrnorm'
    ## by putting error^2 in lower-right position of VCOV
    newVCOV <- VCOV
    newVCOV[varPLACE, varPLACE] <- predERROR^2

    ## create MC simulation matrix
    simMAT <- mvrnorm(n = nsim, mu = MU, Sigma = newVCOV, empirical = TRUE)

    ## evaluate expression on rows of simMAT
    EVAL <- try(eval(EXPR, envir = as.data.frame(simMAT)), silent = TRUE)
    if (inherits(EVAL, "try-error")) stop("There was an error evaluating the simulations!")

    ## collect statistics
    PRED <- data.frame(predVAL)
    colnames(PRED) <- predNAME   
    FITTED <- predict(object, newdata = data.frame(PRED))
    MEAN.sim <- mean(EVAL, na.rm = TRUE)
    SD.sim <- sd(EVAL, na.rm = TRUE)
    MEDIAN.sim <- median(EVAL, na.rm = TRUE)
    MAD.sim <- mad(EVAL, na.rm = TRUE)
    QUANT <- quantile(EVAL, c((1 - level)/2, level + (1 - level)/2))
    RES <- c(FITTED, MEAN.sim, SD.sim, MEDIAN.sim, MAD.sim, QUANT[1], QUANT[2])
    outMAT <- rbind(outMAT, RES)
  }

  colnames(outMAT) <- c("fit", "mean", "sd", "median", "mad", names(QUANT[1]), names(QUANT[2]))
  rownames(outMAT) <- NULL

  cat("\n")

  return(outMAT)  
}

And then he writes: "The input is an ‘nls’ object, a data.frame ‘newdata’ of values to be predicted with the value x_new in the first column and (optional) “errors-in-x” (as sigma) in the second column. The number of simulations can be tweaked with nsim as well as the alpha-level for the confidence interval. The output is f(x_new, beta) (fitted value), mu(y_n) (mean of simulation), sigma(y_n) (s.d. of simulation), median(y_n) (median of simulation), mad(y_n) (mad of simulation) and the lower/upper confidence interval."

He has some additional text explaining this further and giving a usage example, but I don't feel like it's really appropriate for me to copy his entire blog post into this answer, so please visit his page, if it still exists, for further details. Anyway it's pretty simple and self-explanatory, and worked for me right out of the box, on the first try. Thanks Andrej!

Answer 3

A note regarding confidence intervals (2 above), and the answer by @Etienne Low-Décarie:

Even after attaching nls2, the as.lm functions is sometimes unavailable. Based on this (now stale) reference (originally authored by delichon), here's the function's source:

as.lm.nls <- function(object, ...) {
    if (!inherits(object, "nls")) {
        w <- paste("expected object of class nls but got object of class:",
        paste(class(object), collapse = " "))
        warning(w)
    }

    gradient <- object$m$gradient()
    if (is.null(colnames(gradient))) {
        colnames(gradient) <- names(object$m$getPars())
    }

    response.name <- if (length(formula(object)) == 2) "0" else
        as.character(formula(object)[[2]])

    lhs <- object$m$lhs()
    L <- data.frame(lhs, gradient)
    names(L)[1] <- response.name

    fo <- sprintf("%s ~ %s - 1", response.name,
    paste(colnames(gradient), collapse = "+"))
    fo <- as.formula(fo, env = as.proto.list(L))

    do.call("lmst(fo, offset = substitute(fitted(object)))")
}

Then use predict the standard way:

predCI <- predict(as.lm.nls(fittednls), interval = “confidence”, level = 0.95)

Thanks @waybackmachine

Answer 4

I was also banging my head on this one and eventually found predictNLS() function in the propagate package.

For example:

library(propagate)
Y    <- c(282, 314, 581, 846, 1320, 2014, 2798, 4593, 6065, 7818, 9826)
temp <- data.frame(y = Y, x = seq(1:11))
mod  <- nls(y ~ exp(a + b * x), data = temp, start = list(a = 0, b = 1))

(PROP1 <- predictNLS(mod, newdata = data.frame(x = c(12,13)), interval = "prediction"))

Hope this helps.

Link to R documentation