5
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ Does "summary(model1)" deliver what you want? $\endgroup$ – Cyan Mar 27 '12 at 2:01
  • 1
    $\begingroup$ it doesn't, summary(model) gives pvalue for lm functions $\endgroup$ – Kristina Mar 27 '12 at 5:07
  • $\begingroup$ @Kristina does the method described in the answer below for linearizing your models so summary can produce the values you want work for you? $\endgroup$ – Etienne Low-Décarie Mar 27 '12 at 19:27
4
$\begingroup$

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)

$\endgroup$
4
  • $\begingroup$ I've installed and attached the nls2 library, but don't seem to have the as.lm function available. Any ideas? $\endgroup$ – Daniel Kessler Jul 23 '13 at 16:03
  • 1
    $\begingroup$ @Daniel Kessler, try quantitativeconservationbiology.wordpress.com/2013/07/02/… $\endgroup$ – etov Aug 19 '14 at 12:12
  • 1
    $\begingroup$ 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()? $\endgroup$ – bhaller Dec 9 '16 at 11:54
  • $\begingroup$ @bhaller - right. The content was a bit too long for a comment; posted now as an answer. $\endgroup$ – etov Dec 11 '16 at 7:36
3
$\begingroup$

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!

$\endgroup$
2
$\begingroup$

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

$\endgroup$
5
  • $\begingroup$ 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/… $\endgroup$ – bhaller Dec 11 '16 at 11:38
  • $\begingroup$ So, would you post it inline rather than linking to an external page? :) $\endgroup$ – etov Dec 12 '16 at 7:17
  • $\begingroup$ 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? $\endgroup$ – bhaller Dec 12 '16 at 9:14
  • $\begingroup$ @bhaller, I think you should. Besides the fact links might become stale, this seems a more statistically sound solution. $\endgroup$ – etov Dec 13 '16 at 6:57
  • $\begingroup$ OK, done. See my answer. $\endgroup$ – bhaller Dec 14 '16 at 8:06
1
$\begingroup$

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

$\endgroup$

Your Answer

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

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