The R Project SVN R

Rev

Rev 37035 | Rev 39452 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed 


R : Copyright 2006, The R Foundation for Statistical Computing
Version 2.3.0 Under development (unstable) (2006-01-24 r37156)
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 Paraneter 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.00517052 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 Paraneter 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.91488
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.42471    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.13399    1.79455
8  4.52605  166.09149    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:  5.15 0.01 5.23 NA NA 
>