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

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))


• Does "summary(model1)" deliver what you want?
– Cyan
Mar 27, 2012 at 2:01
• it doesn't, summary(model) gives pvalue for lm functions Mar 27, 2012 at 5:07
• @Kristina does the method described in the answer below for linearizing your models so summary can produce the values you want work for you? Mar 27, 2012 at 19:27

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)

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)

• I've installed and attached the nls2 library, but don't seem to have the as.lm function available. Any ideas? Jul 23, 2013 at 16:03
• @Daniel Kessler, try quantitativeconservationbiology.wordpress.com/2013/07/02/…
– etov
Aug 19, 2014 at 12:12
• I'm having the same problem; the as.lm function does not exist. The link posted by @etov is stale. Please post solutions inline, rather than links to external pages that may become stale. So, what is the solution to the missing as.lm()? Dec 9, 2016 at 11:54
• @bhaller - right. The content was a bit too long for a comment; posted now as an answer.
– etov
Dec 11, 2016 at 7:36

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! 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()
names(L)[1] <- response.name

fo <- sprintf("%s ~ %s - 1", response.name,
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

• Yes, I actually found that through Google, but as.proto.list() was also not available, so I wasn't able to use that. I found a completely different solution that seems to work well, at r-bloggers.com/… Dec 11, 2016 at 11:38
• So, would you post it inline rather than linking to an external page? :)
– etov
Dec 12, 2016 at 7:17
• Hahaha, touché. I would love to, except that it doesn't fit in a comment. :-> Do you think I should post it as a solution? Dec 12, 2016 at 9:14
• @bhaller, I think you should. Besides the fact links might become stale, this seems a more statistically sound solution.
– etov
Dec 13, 2016 at 6:57
• OK, done. See my answer. Dec 14, 2016 at 8:06

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.