Colleagues, I am fitting an Emax model using nls. The code is: START <- list(EMAX=INITEMAX, EFFECT=INITEFFECT, C50=INITC50) CONTROL <- list(maxiter=1000, warnOnly=T) #FORMULA <- as.formula(YVAR ~ EMAX - EFFECT * XVAR^GAMMA / (XVAR^GAMMA + C50^GAMMA)) ## alternate version of formula FORMULA <- as.formula(YVAR ~ EMAX - EFFECT / (1 + (C50/XVAR)^GAMMA)) FIT <- nls(FORMULA, start=START, control=CONTROL, trace=T) If GAMMA equals 10-80, nls converges successfully and the fit tracks the fit from a smoother (Supersmoother). However, if I attempt to estimate GAMMA using: START <- list(EMAX=INITEMAX, EFFECT=INITEFFECT, C50=INITC50, GAMMA=INITGAMMA) GAMMA increases rapidly to > 500 and nls terminates with: Error in chol2inv(object$m$Rmat()) : element (4, 4) is zero, so the inverse cannot be computed In addition: Warning message: In nls(FORMULA, start = START, control = CONTROL, trace = T) : singular gradient I also tried fixing GAMMA to > 1000 and I get a similar error message: Error in chol2inv(object$m$Rmat()) : element (2, 2) is zero, so the inverse cannot be computed In addition: Warning message: In nls(FORMULA, start = START, control = CONTROL, trace = T) : singular gradient The data do not suggest a very large value for GAMMA so I am surprised that the estimate is increasing so rapidly. I attempted to use the port algorithm with an upper bound on GAMMA but the upper bound is reached rapidly, suggesting that the data support a large value for GAMMA. A subset of the data (with added noise) is shown below. A GAMMA value of 1280 triggers the error with this subset XVAR <- c(26, 31.3, 20.9, 24.8, 22.9, 4.79, 19.6, 18, 19.6, 9.69, 21.7, 26.6, 27.8, 9.12, 10.5, 20.1, 16.7, 14.1, 10.2, 19.2, 24.7, 34.6, 26.6, 25.1, 5.98, 13.4, 15.7, 9.59, 7.39, 21.5, 15.7, 12.4, 19.2, 17.8, 19.7, 27.1, 25.6, 36.4, 22.9, 8.68, 27, 25.9, 33.3, 24.2, 21.4, 31, 19.1, 18.7, 23.5, 19.4, 10.3, 12.8, 13.9, 18.5, 21, 15.2, 18.9, 9.12, 16.9, 12.9, 29.5, 15.5, 7.34, 8.97, 8.04, 23.7, 16.3, 37.6, 35.2, 13.7, 28.1, 29.5, 15.1, 26, 6.52) YVAR <- c(-34.2, -84.2, -71.1, -91.9, -104.1, -23.2, -27.2, -13.4, -143.2, 24.7, -72.1, -38, 25.2, -8, -34.1, -15.1, -112.6, -93.5, -130.9, -127.8, -118.7, -53.5, -29.8, 98, 0, -37.6, -99.4, 57.9, 0.2, -62.2, -27.3, 8.3, -51.6, -111.6, -25.6, -51.7, -106.4, -85.1, -63.1, -60.8, -27.7, -20.7, 22.9, -49.4, -85.7, -90.9, -107, -20.6, -36.3, -40.2, 39.8, -55, -54.5, -103.9, -53.1, -2.3, -72.3, -65.6, -57.8, -64.4, -129.1, 10.4, -9.9, -29.6, -40.8, 52, -94, 8.8, -98.8, 28, -16.3, -99.2, -48.5, -111.9, -15.4) I suspect that I am making a conceptual error in the use of nls. Any help would be appreciated. If a different function to fit nonlinear regression would work better, please direct me. Dennis Dennis Fisher MD P < (The "P Less Than" Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-866-PLessThan (1-866-753-7784) www.PLessThan.com
Dennis Fisher wrote on 10/11/2011 07:20:35 PM:> > Colleagues, > > I am fitting an Emax model using nls. The code is: > START <- list(EMAX=INITEMAX, EFFECT=INITEFFECT, C50=INITC50) > CONTROL <- list(maxiter=1000, warnOnly=T) > #FORMULA <- as.formula(YVAR ~ EMAX - EFFECT * XVAR^GAMMA / > (XVAR^GAMMA + C50^GAMMA)) ## alternate version of formula > FORMULA <- as.formula(YVAR ~ EMAX - EFFECT / (1 +(C50/XVAR)^GAMMA))> FIT <- nls(FORMULA, start=START, control=CONTROL, trace=T) > > If GAMMA equals 10-80, nls converges successfully and the fit tracks > the fit from a smoother (Supersmoother). However, if I attempt to > estimate GAMMA using: > START <- list(EMAX=INITEMAX, EFFECT=INITEFFECT, > C50=INITC50, GAMMA=INITGAMMA) > GAMMA increases rapidly to > 500 and nls terminates with: > Error in chol2inv(object$m$Rmat()) : > element (4, 4) is zero, so the inverse cannot be computed > In addition: Warning message: > In nls(FORMULA, start = START, control = CONTROL, trace = T) : > singular gradient > > I also tried fixing GAMMA to > 1000 and I get a similar error message: > Error in chol2inv(object$m$Rmat()) > : element (2, 2) is zero, so the inverse cannot be computed > In addition: Warning message: > In nls(FORMULA, start = START, control = CONTROL, trace = T) > : singular gradient > > The data do not suggest a very large value for GAMMA so I am > surprised that the estimate is increasing so rapidly. I attempted > to use the port algorithm with an upper bound on GAMMA but the upper > bound is reached rapidly, suggesting that the data support a large > value for GAMMA. > > A subset of the data (with added noise) is shown below. A GAMMA > value of 1280 triggers the error with this subset > > XVAR <- c(26, 31.3, 20.9, 24.8, 22.9, 4.79, 19.6, 18, 19.6, 9.69, > 21.7, 26.6, 27.8, 9.12, 10.5, 20.1, 16.7, 14.1, 10.2, 19.2, 24.7, 34.6, > 26.6, 25.1, 5.98, 13.4, 15.7, 9.59, 7.39, 21.5, 15.7, 12.4, 19.2, > 17.8, 19.7, 27.1, 25.6, 36.4, 22.9, 8.68, 27, 25.9, 33.3, 24.2, > 21.4, 31, 19.1, 18.7, 23.5, 19.4, 10.3, 12.8, 13.9, 18.5, 21, 15.2, > 18.9, 9.12, 16.9, 12.9, 29.5, 15.5, 7.34, 8.97, 8.04, 23.7, > 16.3, 37.6, 35.2, 13.7, 28.1, 29.5, 15.1, 26, 6.52) > > > YVAR <- c(-34.2, -84.2, -71.1, -91.9, -104.1, -23.2, -27.2, -13.4, > -143.2, 24.7, -72.1, -38, 25.2, -8, -34.1, -15.1, -112.6, -93.5,-130.9,> -127.8, -118.7, -53.5, -29.8, 98, 0, -37.6, -99.4, 57.9, 0.2, -62.2, > -27.3, 8.3, -51.6, -111.6, -25.6, -51.7, -106.4, -85.1, > -63.1, -60.8, -27.7, -20.7, 22.9, -49.4, -85.7, -90.9, -107, -20.6, > -36.3, -40.2, 39.8, -55, -54.5, -103.9, -53.1, -2.3, -72.3, > -65.6, -57.8, -64.4, -129.1, 10.4, -9.9, -29.6, -40.8, 52, -94, 8.8, > -98.8, 28, -16.3, -99.2, -48.5, -111.9, -15.4) > > I suspect that I am making a conceptual error in the use of nls. > Any help would be appreciated. If a different function to fit > nonlinear regression would work better, please direct me. > > Dennis > > Dennis Fisher MD > P < (The "P Less Than" Company) > Phone: 1-866-PLessThan (1-866-753-7784) > Fax: 1-866-PLessThan (1-866-753-7784) > www.PLessThan.comIt's great that you provided code, data, and error messages for clarity. But, without knowing what the starting values are (INITEMAX, INITEFFECT, INITC50, and INITGAMMA), I can't reproduce what you're doing. I'm also a bit confused about GAMMA. In your first fit you don't provide a starting value for GAMMA, so it is being treated as an independent variable. But in your second fit, you do provide starting values for GAMMA, so it is being treated as a parameter to be estimated. Which is correct? Finally, I was able to look at the plot of XVAR vs. YVAR, and there is so much noise that it's not surprising you're having difficulty fitting a four-parameter non-linear model. Jean [[alternative HTML version deleted]]
Dennis Fisher <fisher@plessthan.com> wrote on 10/12/2011 08:06:12 AM:> > jean > > initial values: > INITEMAX <- -25 > INITEFFECT <- 25 > INITC50 <- 14 > GAMMA <- INITGAMMA <- 500 > > see below for other issues. > > dennis > > Dennis Fisher MD > P < (The "P Less Than" Company) > Phone: 1-866-PLessThan (1-866-753-7784) > Fax: 1-866-PLessThan (1-866-753-7784) > www.PLessThan.com > > On Oct 12, 2011, at 5:40 AM, Jean V Adams wrote: > > > It's great that you provided code, data, and error messages for > > clarity. But, without knowing what the starting values are > > (INITEMAX, INITEFFECT, INITC50, and INITGAMMA), I can't reproduce > > what you're doing. > > > > I'm also a bit confused about GAMMA. In your first fit you don't > > provide a starting value for GAMMA, so it is being treated as an > > independent variable. But in your second fit, you do provide > > starting values for GAMMA, so it is being treated as a parameter to > > be estimated. Which is correct? > > your assessment is correct. one approach applies the initial value > for GAMMA; the other estimates a value. > > > Finally, I was able to look at the plot of XVAR vs. YVAR, and there > > is so much noise that it's not surprising you're having difficulty > > fitting a four-parameter non-linear model. > > i sent only a subset of the data and i added noise (to preserve > confidentiality). the entire dataset looks much better, although > somewhat noisy. in that context, i expected GAMMA to estimated at a > small value, i.e., a large value for GAMMA would imply a sharp > distinction, something not supported by the data. > > > > JeanI can't get convergence either if I try to estimate GAMMA as a parameter in the model. If you have a look at the predicted relations from nls() for the different fixed values of GAMMA, you can see that orders of magnitude changes in GAMMA have pretty minor effects on the curve. This may or may not be the case with the full non-noised data set, but for this subset it certainly seems that a four-parameter sigmoid is overkill. A three-parameter sigmoid does just as good a job (see FIT2 below). Jean P.S. It's a good idea to cc: the list when replying so that the thread of conversation is maintained for the benefit of others. XVAR <- c(26, 31.3, 20.9, 24.8, 22.9, 4.79, 19.6, 18, 19.6, 9.69, 21.7, 26.6, 27.8, 9.12, 10.5, 20.1, 16.7, 14.1, 10.2, 19.2, 24.7, 34.6, 26.6, 25.1, 5.98, 13.4, 15.7, 9.59, 7.39, 21.5, 15.7, 12.4, 19.2, 17.8, 19.7, 27.1, 25.6, 36.4, 22.9, 8.68, 27, 25.9, 33.3, 24.2, 21.4, 31, 19.1, 18.7, 23.5, 19.4, 10.3, 12.8, 13.9, 18.5, 21, 15.2, 18.9, 9.12, 16.9, 12.9, 29.5, 15.5, 7.34, 8.97, 8.04, 23.7, 16.3, 37.6, 35.2, 13.7, 28.1, 29.5, 15.1, 26, 6.52) YVAR <- c(-34.2, -84.2, -71.1, -91.9, -104.1, -23.2, -27.2, -13.4, -143.2, 24.7, -72.1, -38, 25.2, -8, -34.1, -15.1, -112.6, -93.5, -130.9, -127.8, -118.7, -53.5, -29.8, 98, 0, -37.6, -99.4, 57.9, 0.2, -62.2, -27.3, 8.3, -51.6, -111.6, -25.6, -51.7, -106.4, -85.1, -63.1, -60.8, -27.7, -20.7, 22.9, -49.4, -85.7, -90.9, -107, -20.6, -36.3, -40.2, 39.8, -55, -54.5, -103.9, -53.1, -2.3, -72.3, -65.6, -57.8, -64.4, -129.1, 10.4, -9.9, -29.6, -40.8, 52, -94, 8.8, -98.8, 28, -16.3, -99.2, -48.5, -111.9, -15.4) CONTROL <- list(maxiter=1000, warnOnly=T) FORMULA <- as.formula(YVAR ~ EMAX - EFFECT / (1 + (C50/XVAR)^GAMMA)) START <- list(EMAX=-25, EFFECT=25, C50=14) GAMMAi <- c(1, 10, 100, 800) plot(XVAR, YVAR) for(i in seq(GAMMAi)) { GAMMA <- GAMMAi[i] FIT <- nls(FORMULA, start=START, control=CONTROL, trace=T) points(XVAR, predict(FIT), pch=i, col=i+1) } FIT2 <- nls(YVAR ~ b + (a-b)/(1 + 10^(XVAR-c)), start=list(a=-20, b=-60, c=10), control=CONTROL, trace=T) points(XVAR, predict(FIT2), pch=16, col="purple") [[alternative HTML version deleted]]