The R Project SVN
R
(root)/branches/R-2-6-branch/src/library/stats/tests/nls.Rout.save – Rev 44088
Rev 42766 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R version 2.7.0 Under development (unstable) (2008-01-15 r44001)Copyright (C) 2008 The R Foundation for Statistical ComputingISBN 3-900051-07-0R 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.> # File src/library/stats/tests/nls.R> # Part of the R package, http://www.R-project.org> #> # This program is free software; you can redistribute it and/or modify> # it under the terms of the GNU General Public License as published by> # the Free Software Foundation; either version 2 of the License, or> # (at your option) any later version.> #> # This program is distributed in the hope that it will be useful,> # but WITHOUT ANY WARRANTY; without even the implied warranty of> # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the> # GNU General Public License for more details.> #> # A copy of the GNU General Public License is available at> # http://www.r-project.org/Licenses/>> ## 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 modelmodel: y ~ a + b * x + c * I(x^2)data: parent.frame()a b c1.0058 0.9824 -0.0897residual sum-of-squares: 0.46Algorithm "port", convergence message: relative convergence (4)> (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 modelmodel: y ~ a + b * x + c * I(x^2)data: parent.frame()a b c1.02 0.89 0.00residual sum-of-squares: 0.468Algorithm "port", convergence message: both X-convergence and relative convergence (5)> confint(fm)Waiting for profiling to be done...2.5% 97.5%a 1.00875 1.037847b 0.84147 0.915063c 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.02413Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 0.00517 0.00764 0.68 0.52x 0.99915 0.00119 841.38 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.0132 on 7 degrees of freedomMultiple R-squared: 1, Adjusted R-squared: 1F-statistic: 7.08e+05 on 1 and 7 DF, p-value: <2e-16Correlation 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.543210.0012128 : 0.0051705 0.9991529> summary(fit, cor = TRUE)Formula: yeps ~ a + b * xParameters:Estimate Std. Error t value Pr(>|t|)a 0.00517 0.00764 0.68 0.52b 0.99915 0.00119 841.37 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.0132 on 7 degrees of freedomCorrelation of Parameter Estimates:ab -0.89Number of iterations to convergence: 1Achieved convergence tolerance: 2.74e-06> 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.070572: 0.123450 0.5432101: 6.3964587: 1.34546 0.7008402: 0.00060639084: 0.00517053 0.9991533: 0.00060639084: 0.00517052 0.999153> summary(fit2, cor = TRUE)Formula: yeps ~ a + b * xParameters:Estimate Std. Error t value Pr(>|t|)a 0.00517 0.00764 0.68 0.52b 0.99915 0.00119 841.38 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.0132 on 7 degrees of freedomCorrelation of Parameter Estimates:ab -0.89Algorithm "port", convergence message: both X-convergence and relative convergence (5)> 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 ' ' 1Residual standard error: 0.0355 on 13 degrees of freedomNumber of iterations to convergence: 2Achieved convergence tolerance: 2.43e-06>> ## 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 ' ' 1Residual standard error: 0.0355 on 13 degrees of freedomNumber of iterations to convergence: 6Achieved convergence tolerance: 1.34e-07> 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 TableModel 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.018552 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 ' ' 1Residual standard error: 0.0355 on 13 degrees of freedomNumber of iterations to convergence: 6Achieved convergence tolerance: 1.29e-07> 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 TableModel 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.018552 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 ' ' 1Residual standard error: 0.0355 on 13 degrees of freedomNumber of iterations to convergence: 5Achieved convergence tolerance: 8.48e-09> 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 TableModel 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.018552 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 ' ' 1Residual standard error: 0.0355 on 13 degrees of freedomAlgorithm "port", convergence message: relative convergence (4)> 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.5725xmid 1.28538 1.6967scal 0.97527 1.1068> (c4 <- confint(fm4, 1:2))Waiting for profiling to be done...2.5% 97.5%xmid 1.2866 1.6949scal 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 modelmodel: y ~ x^bdata: parent.frame()b0.695residual sum-of-squares: 2.39Number of iterations to convergence: 4Achieved convergence tolerance: 2.81e-07Waiting for profiling to be done...2.5% 97.5%0.68704 0.70281Nonlinear regression modelmodel: y ~ x^bdata: parent.frame()b0.695residual sum-of-squares: 2.39Algorithm "port", convergence message: relative convergence (4)Waiting for profiling to be done...2.5% 97.5%0.68704 0.70281Nonlinear regression modelmodel: y2 ~ a * x^bdata: parent.frame()a b0.502 0.724residual sum-of-squares: 2.51Number of iterations to convergence: 4Achieved convergence tolerance: 3.64e-06Waiting for profiling to be done...2.5% 97.5%a 0.49494 0.50893b 0.70019 0.74767Nonlinear regression modelmodel: y2 ~ a * x^bdata: parent.frame()a b0.502 0.724residual sum-of-squares: 2.51Algorithm "port", convergence message: relative convergence (4)Waiting for profiling to be done...2.5% 97.5%a 0.49494 0.50893b 0.70019 0.74767Nonlinear regression modelmodel: y3 ~ a * (x + exp(logc))^bdata: parent.frame()a b logc0.558 0.603 -0.176residual sum-of-squares: 2.44Number of iterations to convergence: 5Achieved convergence tolerance: 3.43e-06Waiting for profiling to be done...2.5% 97.5%a 0.35005 0.66045b 0.45122 0.91492logc -0.64620 0.40950Nonlinear regression modelmodel: y3 ~ a * (x + exp(logc))^bdata: parent.frame()a b logc0.558 0.603 -0.176residual sum-of-squares: 2.44Algorithm "port", convergence message: relative convergence (4)Waiting for profiling to be done...2.5% 97.5%a 0.35005 0.66045b 0.45122 0.91492logc -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.91489logc -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))$atau par.vals.a par.vals.b1 -4.55943 0.61825 0.000002 -3.69338 0.72777 0.000003 -1.43894 1.02459 0.000004 0.00000 1.53805 0.263315 0.70536 2.05151 0.480846 1.25688 2.65703 0.649727 1.78274 3.57029 0.816538 2.26362 5.01493 0.976389 2.70180 7.51383 1.1270710 3.08992 12.25766 1.2643411 3.40458 22.42471 1.38476$btau par.vals.a par.vals.b1 -0.74662 1.17755 0.000002 0.00000 1.53805 0.263313 0.87566 2.14030 0.556014 1.67482 3.03728 0.834065 2.43255 4.61672 1.123476 3.18042 7.90950 1.441177 3.89865 23.13399 1.794558 4.52605 166.09149 2.20381attr(,"original.fit")Nonlinear regression modelmodel: y ~ gfun(a, b, x)data: parent.frame()a b1.53805 0.26331residual sum-of-squares: 0.38938Algorithm "port", convergence message: relative convergence (4)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 ' ' 1Residual standard error: 0.22062 on 8 degrees of freedomAlgorithm "port", convergence message: relative convergence (4)attr(,"class")[1] "profile.nls" "profile"> confint(pr1)2.5% 97.5%a 0.98596 5.2072b 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)$atau par.vals.a par.vals.b1 -4.4367 0.63351 0.000002 -3.5758 0.74251 0.000003 -1.5917 1.00251 0.000004 0.0000 1.50000 0.24263$btau par.vals.a par.vals.b1 -0.74358 1.17755 0.000002 0.00000 1.50000 0.242633 1.95266 1.50000 0.538494 2.54777 1.50000 0.664525 3.23041 1.50000 0.84069attr(,"original.fit")Nonlinear regression modelmodel: y ~ gfun(a, b, x)data: parent.frame()a b1.50000 0.24263residual sum-of-squares: 0.38958Algorithm "port", convergence message: relative convergence (4)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 ' ' 1Residual standard error: 0.22068 on 8 degrees of freedomAlgorithm "port", convergence message: relative convergence (4)attr(,"class")[1] "profile.nls" "profile"> confint(m2)Waiting for profiling to be done...2.5% 97.5%a 0.90457 NAb NA 0.61016>> ## scoping problems> test <- function()+ {+ x <- seq(0,5,len=20)+ n <- 1+ y <- 2*x^2 + n + rnorm(x)+ xy <- data.frame(x=x,y=y)+ myf <- function(x,a,b,c) a*x^b+c+ nls(y ~ myf(x,a,b,n), data=xy, start=c(a=1,b=1), trace=TRUE)+ }> test()8291.9 : 1 1726.02 : 0.80544 2.42971552.85 : 1.2900 2.129070.431 : 1.9565 1.967026.555 : 1.9788 2.006426.503 : 1.9798 2.004626.503 : 1.9799 2.0046Nonlinear regression modelmodel: y ~ myf(x, a, b, n)data: xya b1.98 2.00residual sum-of-squares: 26.5Number of iterations to convergence: 6Achieved convergence tolerance: 7.45e-07> ## failed to find n in 2.2.x> ## found wrong n in 2.3.x> ## finally worked in 2.4.0>>> ## list 'start'> getExpmat <- function(theta, t)+ {+ conc <- matrix(nrow = length(t), ncol = length(theta))+ for(i in 1:length(theta)) conc[, i] <- exp(-theta[i] * t)+ conc+ }>> expsum <- as.vector(getExpmat(c(.05,.005), 1:100) %*% c(1,1))> expsumNoisy <- expsum + max(expsum) *.001 * rnorm(100)> expsum.df <-data.frame(expsumNoisy)>> ## estimate decay rates, amplitudes with default Gauss-Newton> summary (nls(expsumNoisy ~ getExpmat(k, 1:100) %*% sp, expsum.df,+ start = list(k = c(.6,.02), sp = c(1,2))))Formula: expsumNoisy ~ getExpmat(k, 1:100) %*% spParameters:Estimate Std. Error t value Pr(>|t|)k1 5.00e-02 3.07e-04 162.5 <2e-16 ***k2 4.99e-03 5.38e-05 92.8 <2e-16 ***sp1 1.00e+00 4.47e-03 224.0 <2e-16 ***sp2 9.99e-01 5.00e-03 199.8 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.00205 on 96 degrees of freedomNumber of iterations to convergence: 7Achieved convergence tolerance: 1.98e-06>> ## didn't work with port in 2.4.1> summary (nls(expsumNoisy ~ getExpmat(k, 1:100) %*% sp, expsum.df,+ start = list(k = c(.6,.02), sp = c(1,2)),+ algorithm = "port"))Formula: expsumNoisy ~ getExpmat(k, 1:100) %*% spParameters:Estimate Std. Error t value Pr(>|t|)k1 5.00e-02 3.07e-04 162.5 <2e-16 ***k2 4.99e-03 5.38e-05 92.8 <2e-16 ***sp1 1.00e+00 4.47e-03 224.0 <2e-16 ***sp2 9.99e-01 5.00e-03 199.8 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.00205 on 96 degrees of freedomAlgorithm "port", convergence message: both X-convergence and relative convergence (5)>>> cat('Time elapsed: ', proc.time() - .proctime00,'\n')Time elapsed: 7.453 0.057 7.515 0 0>