R : Copyright 2006, The R Foundation for Statistical Computing Version 2.3.0 beta (2006-04-13 r37760) ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## tests of nls, especially of weighted fits > > .proctime00 <- proc.time() > library(stats) > options(digits=5) # to avoid trivial printed differences > > postscript("nls-test.ps") > > ## selfStart.default() w/ no parameters: > logist <- deriv( ~Asym/(1+exp(-(x-xmid)/scal)), c("Asym", "xmid", "scal"), + function(x, Asym, xmid, scal){} ) > logistInit <- function(mCall, LHS, data) { + xy <- sortedXyData(mCall[["x"]], LHS, data) + if(nrow(xy) < 3) stop("Too few distinct input values to fit a logistic") + Asym <- max(abs(xy[,"y"])) + if (Asym != max(xy[,"y"])) Asym <- -Asym # negative asymptote + xmid <- NLSstClosestX(xy, 0.5 * Asym) + scal <- NLSstClosestX(xy, 0.75 * Asym) - xmid + value <- c(Asym, xmid, scal) + names(value) <- mCall[c("Asym", "xmid", "scal")] + value + } > logist <- selfStart(logist, initial = logistInit) ##-> Error in R 1.5.0 > str(logist) function (x, Asym, xmid, scal) - attr(*, "initial")=function (mCall, LHS, data) - attr(*, "class")= chr "selfStart" > > ## lower and upper in algorithm="port" > set.seed(123) > x <- runif(200) > a <- b <- 1; c <- -0.1 > y <- a+b*x+c*x^2+rnorm(200, sd=0.05) > plot(x,y) > curve(a+b*x+c*x^2, add = TRUE) > nls(y ~ a+b*x+c*I(x^2), start = c(a=1, b=1, c=0.1), algorithm = "port") Nonlinear regression model model: y ~ a + b * x + c * I(x^2) data: parent.frame() a b c 1.00584 0.98239 -0.08972 residual sum-of-squares: 0.46041 > (fm <- nls(y ~ a+b*x+c*I(x^2), start = c(a=1, b=1, c=0.1), + algorithm = "port", lower = c(0, 0, 0))) Nonlinear regression model model: y ~ a + b * x + c * I(x^2) data: parent.frame() a b c 1.02306 0.88971 0.00000 residual sum-of-squares: 0.46824 > confint(fm) Waiting for profiling to be done... 2.5% 97.5% a 1.00875 1.037847 b 0.84147 0.915063 c NA 0.042693 > > > ## weighted nls fit: unsupported < 2.3.0 > set.seed(123) > y <- x <- 1:10 > yeps <- y + rnorm(length(y), sd = 0.01) > wts <- rep(c(1, 2), length = 10); wts[5] <- 0 > fit0 <- lm(yeps ~ x, weights = wts) > summary(fit0, cor = TRUE) Call: lm(formula = yeps ~ x, weights = wts) Residuals: Min 1Q Median 3Q Max -0.01562 -0.00723 -0.00158 0.00403 0.02413 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.00517 0.00764 0.68 0.52 x 0.99915 0.00119 841.38 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0132 on 7 degrees of freedom Multiple R-Squared: 1, Adjusted R-squared: 1 F-statistic: 7.08e+05 on 1 and 7 DF, p-value: <2e-16 Correlation of Coefficients: (Intercept) x -0.89 > cf0 <- coef(summary(fit0))[, 1:2] > fit <- nls(yeps ~ a + b*x, start = list(a = 0.12345, b = 0.54321), + weights = wts, trace = TRUE) 112.14 : 0.12345 0.54321 0.0012128 : 0.0051705 0.9991529 > summary(fit, cor = TRUE) Formula: yeps ~ a + b * x Parameters: Estimate Std. Error t value Pr(>|t|) a 0.00517 0.00764 0.68 0.52 b 0.99915 0.00119 841.37 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0132 on 7 degrees of freedom Correlation of Parameter Estimates: a b -0.89 > stopifnot(all.equal(residuals(fit), residuals(fit0), 1e5, + check.attributes = FALSE)) > stopifnot(df.residual(fit) == df.residual(fit0)) > cf1 <- coef(summary(fit))[, 1:2] > fit2 <- nls(yeps ~ a + b*x, start = list(a = 0.12345, b = 0.54321), + weights = wts, trace = TRUE, algorithm = "port") 0 56.0706: 0.123450 0.543210 1 6.39646: 1.34546 0.700840 2 0.000606391: 0.00517053 0.999153 3 0.000606391: 0.00517051 0.999153 > summary(fit2, cor = TRUE) Formula: yeps ~ a + b * x Parameters: Estimate Std. Error t value Pr(>|t|) a 0.00517 0.00764 0.68 0.52 b 0.99915 0.00119 841.38 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0132 on 7 degrees of freedom Correlation of Parameter Estimates: a b -0.89 > cf2 <- coef(summary(fit2))[, 1:2] > rownames(cf0) <- c("a", "b") > # expect relative errors ca 2e-08 > stopifnot(all.equal(cf1, cf0, 1e-6), all.equal(cf1, cf0, 1e-6)) > stopifnot(all.equal(residuals(fit2), residuals(fit0), 1e5, + check.attributes = FALSE)) > > > DNase1 <- subset(DNase, Run == 1) > DNase1$wts <- rep(8:1, each = 2) > fm1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), + data = DNase1, weights = wts) > summary(fm1) Formula: density ~ SSlogis(log(conc), Asym, xmid, scal) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > > ## directly > fm2 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), + data = DNase1, weights = wts, + start = list(Asym = 3, xmid = 0, scal = 1)) > summary(fm2) Formula: density ~ Asym/(1 + exp((xmid - log(conc))/scal)) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > stopifnot(all.equal(coef(summary(fm2)), coef(summary(fm1)), tol = 1e-6)) > stopifnot(all.equal(residuals(fm2), residuals(fm1), tol = 1e-5)) > stopifnot(all.equal(fitted(fm2), fitted(fm1), tol = 1e-6)) > fm2a <- nls(density ~ Asym/(1 + exp((xmid - log(conc)))), + data = DNase1, weights = wts, + start = list(Asym = 3, xmid = 0)) > anova(fm2a, fm2) Analysis of Variance Table Model 1: density ~ Asym/(1 + exp((xmid - log(conc)))) Model 2: density ~ Asym/(1 + exp((xmid - log(conc))/scal)) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 14 0.01855 2 13 0.01643 1 0.00212 1.68 0.22 > > ## and without using weights > fm3 <- nls(~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))/scal))), + data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1)) > summary(fm3) Formula: 0 ~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))/scal))) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > stopifnot(all.equal(coef(summary(fm3)), coef(summary(fm1)), tol = 1e-6)) > ft <- with(DNase1, density - fitted(fm3)/sqrt(wts)) > stopifnot(all.equal(ft, fitted(fm1), tol = 1e-6)) > # sign of residuals is reversed > r <- with(DNase1, -residuals(fm3)/sqrt(wts)) > all.equal(r, residuals(fm1), tol = 1e05) [1] TRUE > fm3a <- nls(~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))))), + data = DNase1, start = list(Asym = 3, xmid = 0)) > anova(fm3a, fm3) Analysis of Variance Table Model 1: 0 ~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))))) Model 2: 0 ~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))/scal))) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 14 0.01855 2 13 0.01643 1 0.00212 1.68 0.22 > > ## using conditional linearity > fm4 <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)), + data = DNase1, weights = wts, + start = list(xmid = 0, scal = 1), algorithm = "plinear") > summary(fm4) Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal)) Parameters: Estimate Std. Error t value Pr(>|t|) xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** .lin 2.3350 0.0966 24.2 3.5e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > cf <- coef(summary(fm4))[c(3,1,2), ] > rownames(cf)[2] <- "Asym" > stopifnot(all.equal(cf, coef(summary(fm1)), tol = 1e-6, + check.attributes = FALSE)) > stopifnot(all.equal(residuals(fm4), residuals(fm1), tol = 1e-5)) > stopifnot(all.equal(fitted(fm4), fitted(fm1), tol = 1e-6)) > fm4a <- nls(density ~ 1/(1 + exp((xmid - log(conc)))), + data = DNase1, weights = wts, + start = list(xmid = 0), algorithm = "plinear") > anova(fm4a, fm4) Analysis of Variance Table Model 1: density ~ 1/(1 + exp((xmid - log(conc)))) Model 2: density ~ 1/(1 + exp((xmid - log(conc))/scal)) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 14 0.01855 2 13 0.01643 1 0.00212 1.68 0.22 > > ## using 'port' > fm5 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), + data = DNase1, weights = wts, + start = list(Asym = 3, xmid = 0, scal = 1), + algorithm = "port") > summary(fm5) Formula: density ~ Asym/(1 + exp((xmid - log(conc))/scal)) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > stopifnot(all.equal(coef(summary(fm5)), coef(summary(fm1)), tol = 1e-6)) > stopifnot(all.equal(residuals(fm5), residuals(fm1), tol = 1e-5)) > stopifnot(all.equal(fitted(fm5), fitted(fm1), tol = 1e-6)) > > ## check profiling > pfm1 <- profile(fm1) > pfm3 <- profile(fm3) > for(m in names(pfm1)) stopifnot(all.equal(pfm1[[m]], pfm3[[m]], tol=1e-5)) > pfm5 <- profile(fm5) > for(m in names(pfm1)) stopifnot(all.equal(pfm1[[m]], pfm5[[m]], tol=1e-5)) > (c1 <- confint(fm1)) Waiting for profiling to be done... 2.5% 97.5% Asym 2.14936 2.5725 xmid 1.28538 1.6967 scal 0.97527 1.1068 > (c4 <- confint(fm4, 1:2)) Waiting for profiling to be done... 2.5% 97.5% xmid 1.2866 1.6949 scal 0.9757 1.1063 > stopifnot(all.equal(c1[2:3, ], c4, tol = 1e-3)) > > ## some low-dimensional examples > npts <- 1000 > set.seed(1001) > x <- runif(npts) > b <- 0.7 > y <- x^b+rnorm(npts, sd=0.05) > a <- 0.5 > y2 <- a*x^b+rnorm(npts, sd=0.05) > c <- 1.0 > y3 <- a*(x+c)^b+rnorm(npts, sd=0.05) > d <- 0.5 > y4 <- a*(x^d+c)^b+rnorm(npts, sd=0.05) > m1 <- c(y ~ x^b, y2 ~ a*x^b, y3 ~ a*(x+exp(logc))^b) > s1 <- list(c(b=1), c(a=1,b=1), c(a=1,b=1,logc=0)) > for(p in 1:3) { + fm <- nls(m1[[p]], start = s1[[p]]) + print(fm) + print(confint(fm)) + fm <- nls(m1[[p]], start = s1[[p]], algorithm="port") + print(fm) + print(confint(fm)) + } Nonlinear regression model model: y ~ x^b data: parent.frame() b 0.69489 residual sum-of-squares: 2.3866 Waiting for profiling to be done... 2.5% 97.5% 0.68704 0.70281 Nonlinear regression model model: y ~ x^b data: parent.frame() b 0.69489 residual sum-of-squares: 2.3866 Waiting for profiling to be done... 2.5% 97.5% 0.68704 0.70281 Nonlinear regression model model: y2 ~ a * x^b data: parent.frame() a b 0.50191 0.72370 residual sum-of-squares: 2.5143 Waiting for profiling to be done... 2.5% 97.5% a 0.49494 0.50893 b 0.70019 0.74767 Nonlinear regression model model: y2 ~ a * x^b data: parent.frame() a b 0.50191 0.72370 residual sum-of-squares: 2.5143 Waiting for profiling to be done... 2.5% 97.5% a 0.49494 0.50893 b 0.70019 0.74767 Nonlinear regression model model: y3 ~ a * (x + exp(logc))^b data: parent.frame() a b logc 0.55767 0.60255 -0.17648 residual sum-of-squares: 2.4382 Waiting for profiling to be done... 2.5% 97.5% a 0.35005 0.66045 b 0.45122 0.91492 logc -0.64620 0.40950 Nonlinear regression model model: y3 ~ a * (x + exp(logc))^b data: parent.frame() a b logc 0.55766 0.60256 -0.17645 residual sum-of-squares: 2.4382 Waiting for profiling to be done... 2.5% 97.5% a 0.35005 0.66045 b 0.45122 0.91492 logc -0.64620 0.40950 > > fm <- nls(y2~x^b, start=c(b=1), algorithm="plinear") > confint(profile(fm)) 2.5% 97.5% 0.70019 0.74767 > fm <- nls(y3 ~ (x+exp(logc))^b, start=c(b=1, logc=0), algorithm="plinear") > confint(profile(fm)) 2.5% 97.5% b 0.45121 0.91489 logc -0.64628 0.40953 > > > ## more profiling with bounds > npts <- 10 > set.seed(1001) > a <- 2 > b <- 0.5 > x <- runif(npts) > y <- a*x/(1+a*b*x) + rnorm(npts, sd=0.2) > gfun <- function(a,b,x) { + if(a < 0 || b < 0) stop("bounds violated") + a*x/(1+a*b*x) + } > m1 <- nls(y ~ gfun(a,b,x), algorithm = "port", + lower = c(0,0), start = c(a=1, b=1)) > (pr1 <- profile(m1)) $a tau par.vals.a par.vals.b 1 -4.55943 0.61825 0.00000 2 -3.69338 0.72777 0.00000 3 -1.43894 1.02459 0.00000 4 0.00000 1.53805 0.26331 5 0.70536 2.05151 0.48084 6 1.25688 2.65703 0.64972 7 1.78274 3.57029 0.81653 8 2.26362 5.01493 0.97638 9 2.70180 7.51383 1.12707 10 3.08992 12.25766 1.26434 11 3.40458 22.42470 1.38476 $b tau par.vals.a par.vals.b 1 -0.74662 1.17755 0.00000 2 0.00000 1.53805 0.26331 3 0.87566 2.14030 0.55601 4 1.67482 3.03728 0.83406 5 2.43255 4.61672 1.12347 6 3.18042 7.90950 1.44117 7 3.89865 23.13401 1.79455 8 4.52605 166.09148 2.20381 attr(,"original.fit") Nonlinear regression model model: y ~ gfun(a, b, x) data: parent.frame() a b 1.53805 0.26331 residual sum-of-squares: 0.38938 attr(,"summary") Formula: y ~ gfun(a, b, x) Parameters: Estimate Std. Error t value Pr(>|t|) a 1.53805 0.61727 2.4917 0.03742 * b 0.26331 0.35188 0.7483 0.47569 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.22062 on 8 degrees of freedom attr(,"class") [1] "profile.nls" "profile" > confint(pr1) 2.5% 97.5% a 0.98596 5.2072 b NA 1.0733 > > gfun <- function(a,b,x) { + if(a < 0 || b < 0 || a > 1.5 || b > 1) stop("bounds violated") + a*x/(1+a*b*x) + } > m2 <- nls(y ~ gfun(a,b,x), algorithm = "port", + lower = c(0, 0), upper=c(1.5, 1), start = c(a=1, b=1)) > profile(m2) $a tau par.vals.a par.vals.b 1 -4.4367 0.63351 0.00000 2 -3.5758 0.74251 0.00000 3 -1.5917 1.00251 0.00000 4 0.0000 1.50000 0.24263 $b tau par.vals.a par.vals.b 1 -0.74358 1.17755 0.00000 2 0.00000 1.50000 0.24263 3 1.95266 1.50000 0.53849 4 2.54777 1.50000 0.66452 5 3.23041 1.50000 0.84069 attr(,"original.fit") Nonlinear regression model model: y ~ gfun(a, b, x) data: parent.frame() a b 1.50000 0.24263 residual sum-of-squares: 0.38958 attr(,"summary") Formula: y ~ gfun(a, b, x) Parameters: Estimate Std. Error t value Pr(>|t|) a 1.50000 0.59807 2.5080 0.03648 * b 0.24263 0.35567 0.6822 0.51438 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.22068 on 8 degrees of freedom attr(,"class") [1] "profile.nls" "profile" > confint(m2) Waiting for profiling to be done... 2.5% 97.5% a 0.90457 NA b NA 0.61016 > > cat('Time elapsed: ', proc.time() - .proctime00,'\n') Time elapsed: 4.744 0.021 4.912 0 0 >