Ajuste de modelo de regressão não linear
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "reshape",
"car", "alr3", "nls2", "nlstools", "rootSolve",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## lattice latticeExtra reshape car alr3 nls2
## TRUE TRUE TRUE TRUE TRUE TRUE
## nlstools rootSolve plyr wzRfun
## TRUE TRUE TRUE TRUE
source("lattice_setup.R")
url <- "http://nls2.googlecode.com/svn-history/r4/trunk/R/as.lm.R"
## download.file(url, dest=basename(url))
path <- ifelse(Sys.info()["user"]=="walmes", basename(url), url)
source(path)
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()
## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 plyr_1.8.1 rootSolve_1.6.5 nlstools_1.0-0
## [5] nls2_0.2 proto_0.3-10 alr3_2.0.5 car_2.0-21
## [9] reshape_0.8.5 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 MASS_7.3-34 methods_3.1.1
## [6] nnet_7.3-8 Rcpp_0.11.2 stringr_0.6.2 tools_3.1.1
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
Ganho de peso em perus em função da metionina na dieta
##-----------------------------------------------------------------------------
## Ajuste de modelo de regressão não linear.
## turk0
str(turk0)
## 'data.frame': 35 obs. of 2 variables:
## $ A : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Gain: int 644 631 661 624 633 610 615 605 608 599 ...
xyplot(Gain~A, data=turk0, type=c("p","smooth"))
##-----------------------------------------------------------------------------
## Valores iniciais baseados na interpretação gráfica.
## Modelo: th0+th1*x/(th2+x);
start <- list(th0=625, th1=800-625, th2=0.1)
xyplot(Gain~A, data=turk0)+
layer(with(start, panel.curve(th0+th1*x/(th2+x))))
##-----------------------------------------------------------------------------
## Ajuste.
n0 <- nls(Gain~th0+th1*A/(th2+A), data=turk0,
start=start)
summary(n0)
##
## Formula: Gain ~ th0 + th1 * A/(th2 + A)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## th0 622.2417 6.2028 100.32 < 2e-16 ***
## th1 235.2169 22.5603 10.43 8.1e-12 ***
## th2 0.1548 0.0406 3.81 0.00059 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.1 on 32 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 1.18e-06
xyplot(Gain~A, data=turk0)+
layer(with(as.list(coef(n0)),
panel.curve(th0+th1*x/(th2+x), col=2)))
##-----------------------------------------------------------------------------
## Intervalos de confiança.
## Baseado na log-verossimilhança.
confint(n0)
## Waiting for profiling to be done...
## 2.5% 97.5%
## th0 609.65652 634.7144
## th1 196.53695 290.9141
## th2 0.09408 0.2665
## Baseado na aproximação quadrática da verossimilhança, conhecido como
## intervalos de Wald ou assintóticos. São simétricos por construção.
confint.default(n0)
## 2.5 % 97.5 %
## th0 610.0844 634.3990
## th1 190.9996 279.4343
## th2 0.0752 0.2343
##-----------------------------------------------------------------------------
## Valores iniciais baseados na interpretação gráfica.
## Modelo: th0+th1*(1-2^(-x/th2))
start <- list(th0=625, th1=800-625, th2=0.1)
xyplot(Gain~A, data=turk0)+
layer(with(start, panel.curve(th0+th1*(1-2^(-x/th2)))))
n1 <- nls(Gain~th0+th1*(1-2^(-A/th2)), data=turk0,
start=start)
summary(n1)
##
## Formula: Gain ~ th0 + th1 * (1 - 2^(-A/th2))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## th0 622.9581 5.9009 105.57 < 2e-16 ***
## th1 178.2519 11.6356 15.32 2.7e-16 ***
## th2 0.0973 0.0165 5.91 1.4e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.7 on 32 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 2.34e-06
xyplot(Gain~A, data=turk0)+
layer(with(as.list(coef(n1)),
panel.curve(th0+th1*(1-2^(-x/th2)), col=2)))
confint(n1)
## Waiting for profiling to be done...
## 2.5% 97.5%
## th0 610.9709 634.8552
## th1 156.3227 204.4767
## th2 0.0709 0.1396
confint.default(n1)
## 2.5 % 97.5 %
## th0 611.39253 634.5236
## th1 155.44658 201.0573
## th2 0.06505 0.1296
##-----------------------------------------------------------------------------
## Explorando o conteúdo dentro da estrutura do objeto.
## Estrutura do objeto.
str(n1)
## List of 6
## $ m :List of 16
## ..$ resid :function ()
## ..$ fitted :function ()
## ..$ formula :function ()
## ..$ deviance :function ()
## ..$ lhs :function ()
## ..$ gradient :function ()
## ..$ conv :function ()
## ..$ incr :function ()
## ..$ setVarying:function (vary = rep_len(TRUE, length(useParams)))
## ..$ setPars :function (newPars)
## ..$ getPars :function ()
## ..$ getAllPars:function ()
## ..$ getEnv :function ()
## ..$ trace :function ()
## ..$ Rmat :function ()
## ..$ predict :function (newdata = list(), qr = FALSE)
## ..- attr(*, "class")= chr "nlsModel"
## $ convInfo :List of 5
## ..$ isConv : logi TRUE
## ..$ finIter : int 4
## ..$ finTol : num 2.34e-06
## ..$ stopCode : int 0
## ..$ stopMessage: chr "converged"
## $ data : symbol turk0
## $ call : language nls(formula = Gain ~ th0 + th1 * (1 - 2^(-A/th2)), data = turk0, start = start, algorithm = "default", control = structure(list(maxiter = 50, tol = 1e-05, minFactor = 0.0009765625, ...
## $ dataClasses: Named chr "numeric"
## ..- attr(*, "names")= chr "A"
## $ control :List of 5
## ..$ maxiter : num 50
## ..$ tol : num 1e-05
## ..$ minFactor: num 0.000977
## ..$ printEval: logi FALSE
## ..$ warnOnly : logi FALSE
## - attr(*, "class")= chr "nls"
n1
## Nonlinear regression model
## model: Gain ~ th0 + th1 * (1 - 2^(-A/th2))
## data: turk0
## th0 th1 th2
## 622.9581 178.2519 0.0973
## residual sum-of-squares: 12367
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 2.34e-06
## Estrutura do summary do objeto.
sn1 <- summary(n1)
str(sn1)
## List of 11
## $ formula :Class 'formula' length 3 Gain ~ th0 + th1 * (1 - 2^(-A/th2))
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## $ residuals : num [1:35] 21.04 8.04 38.04 1.04 10.04 ...
## $ sigma : num 19.7
## $ df : int [1:2] 3 32
## $ cov.unscaled: num [1:3, 1:3] 0.090096 -0.05547 0.000102 -0.05547 0.350306 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "th0" "th1" "th2"
## .. ..$ : chr [1:3] "th0" "th1" "th2"
## $ call : language nls(formula = Gain ~ th0 + th1 * (1 - 2^(-A/th2)), data = turk0, start = start, algorithm = "default", control = structure(list(maxiter = 50, tol = 1e-05, minFactor = 0.0009765625, ...
## $ convInfo :List of 5
## ..$ isConv : logi TRUE
## ..$ finIter : int 4
## ..$ finTol : num 2.34e-06
## ..$ stopCode : int 0
## ..$ stopMessage: chr "converged"
## $ control :List of 5
## ..$ maxiter : num 50
## ..$ tol : num 1e-05
## ..$ minFactor: num 0.000977
## ..$ printEval: logi FALSE
## ..$ warnOnly : logi FALSE
## $ na.action : NULL
## $ coefficients: num [1:3, 1:4] 622.9581 178.2519 0.0973 5.9009 11.6356 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "th0" "th1" "th2"
## .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
## $ parameters : num [1:3, 1:4] 622.9581 178.2519 0.0973 5.9009 11.6356 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "th0" "th1" "th2"
## .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
## - attr(*, "class")= chr "summary.nls"
## Erro padrão residual.
sn1$sigma
## [1] 19.66
cbind(sn1$parameters, confint.default(n1))
## Estimate Std. Error t value Pr(>|t|) 2.5 % 97.5 %
## th0 622.95806 5.90089 105.57 2.856e-42 611.39253 634.5236
## th1 178.25193 11.63559 15.32 2.744e-16 155.44658 201.0573
## th2 0.09732 0.01647 5.91 1.409e-06 0.06505 0.1296
##-----------------------------------------------------------------------------
## Colocando as curvas no mesmo gráfico.
xyplot(Gain~A, data=turk0)+
layer(with(as.list(coef(n0)),
panel.curve(th0+th1*x/(th2+x), col=2)))+
layer(with(as.list(coef(n1)),
panel.curve(th0+th1*(1-2^(-x/th2)), col=4)))
## Soma de quadrados dos resíduos (menor melhor).
deviance(n0)
## [1] 12982
deviance(n1)
## [1] 12367
## R² (maior melhor).
cor(turk0$Gain, fitted(n0))^2
## [1] 0.9188
cor(turk0$Gain, fitted(n1))^2
## [1] 0.9226
## R² com quadro de anova.
R2nls(n0)
## $anova
## Df Sum Sq Mean Sq F value Pr(>F)
## 1 regression 2 146882 73441.2 181 0
## 2 residuals 32 12982 405.7 NA NA
##
## $R2
## [1] 0.9188
R2nls(n1)
## $anova
## Df Sum Sq Mean Sq F value Pr(>F)
## 1 regression 2 147497 73748.7 190.8 0
## 2 residuals 32 12367 386.5 NA NA
##
## $R2
## [1] 0.9226
##-----------------------------------------------------------------------------
## Colocar bandas de confiança.
## Modelo escrito como função dos parâmetros (theta).
f <- function(theta, xx){
with(as.list(theta), th0+th1*(1-2^(-xx/th2)))
}
## Matriz de derivadas parciais em theta (n x p).
gradient(f, x=coef(n1), xx=c(0.0, 0.2, 0.4))
## th0 th1 th2
## [1,] 1 0.0000 0.0
## [2,] 1 0.7594 -627.8
## [3,] 1 0.9421 -302.2
pred <- data.frame(A=seq(0, 0.5, l=20))
## pred$fit <- f(theta=coef(n1), xx=pred$A)
pred$fit <- predict(n1, newdata=pred)
der <- gradient(f, x=coef(n1), xx=pred$A)
str(der)
## num [1:20, 1:3] 1 1 1 1 1 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "th0" "th1" "th2"
## Etapas até o IC passando pelo erro padrão e margem de erro.
F <- der
U <- chol(vcov(n1))
pred$se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
tval <- qt(p=c(lwr=0.025, upr=0.975), df=df.residual(n1))
me <- outer(pred$se, tval, "*")
pred <- cbind(pred, sweep(me, 1, pred$fit, "+"))
str(pred)
## 'data.frame': 20 obs. of 5 variables:
## $ A : num 0 0.0263 0.0526 0.0789 0.1053 ...
## $ fit: num 623 653 679 700 717 ...
## $ se : num 5.9 4.4 4.76 5.34 5.61 ...
## $ lwr: num 611 644 669 689 706 ...
## $ upr: num 635 662 688 710 728 ...
## Equação do modelo ajustado.
coef(n1)
## th0 th1 th2
## 622.95806 178.25193 0.09732
formula(n1)
## Gain ~ th0 + th1 * (1 - 2^(-A/th2))
## Arredonda as estimativas.
theta <- mapply(round,
x=coef(n1), digits=c(2,2,4),
SIMPLIFY=FALSE)
theta
## $th0
## [1] 623
##
## $th1
## [1] 178.2
##
## $th2
## [1] 0.0973
## Equação.
formula(n1)
## Gain ~ th0 + th1 * (1 - 2^(-A/th2))
eq <- substitute(
expr=expression(Gain==th0+th1*(1-2^(-A/th2))),
env=theta)
eval(eq)
## expression(Gain == 622.96 + 178.25 * (1 - 2^(-A/0.0973)))
## Observados, preditos e a banda de confiança.
xyplot(Gain~A, data=turk0)+
as.layer(xyplot(fit~A, data=pred, type="l",
prepanel=prepanel.cbH, cty="bands",
ly=pred$lwr, uy=pred$upr,
panel=panel.cbH))+
layer(panel.key(points=FALSE, text=eval(eq),
corner=c(0.95,0.05)))
Consumo de energia em função da temperatura
##-----------------------------------------------------------------------------
## Consumo de energia (KWH/dia) em função da temperatura (F).
str(segreg)
## 'data.frame': 39 obs. of 2 variables:
## $ Temp: num 8.55 10.66 12.45 16.9 20.52 ...
## $ C : num 86.2 72.4 74.1 68.2 72.4 ...
xyplot(C~Temp, data=segreg, type=c("p","smooth"))
##-----------------------------------------------------------------------------
## Ajuste do modelo platô linear.
## f(x) = th0+th1*(x-th2)*(x>=th2)+0*(x<th2)
start <- list(th0=75, th1=0.5, th2=50)
xyplot(C~Temp, data=segreg)+
layer(with(start,
panel.curve(th0+th1*(x-th2)*(x>=th2)+0*(x<th2))))
## Ajuste.
n2 <- nls(C~th0+th1*(Temp-th2)*(Temp>=th2)+0*(Temp<th2),
data=segreg, start=start)
## Estimativas e medidas de ajuste.
summary(n2)
##
## Formula: C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) + 0 * (Temp < th2)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## th0 74.695 1.343 55.61 < 2e-16 ***
## th1 0.567 0.101 5.64 2.1e-06 ***
## th2 41.951 4.658 9.01 9.4e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.37 on 36 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 9.32e-09
## Valor de F e R².
R2nls(n2)
## $anova
## Df Sum Sq Mean Sq F value Pr(>F)
## 1 regression 2 2098 1049.07 36.34 2.308e-09
## 2 residuals 36 1039 28.87 NA NA
##
## $R2
## [1] 0.6687
## Intervalos de confiança.
confint(n2)
## Waiting for profiling to be done...
## Error: step factor 0.000488281 reduced below 'minFactor' of 0.000976562
confint.default(n2)
## 2.5 % 97.5 %
## th0 72.0626 77.3281
## th1 0.3703 0.7646
## th2 32.8212 51.0812
## Observados e preditos.
xyplot(C~Temp, data=segreg)+
layer(with(as.list(coef(n2)),
panel.curve(th0+th1*(x-th2)*(x>=th2)+0*(x<th2), col=2)))
##-----------------------------------------------------------------------------
## Análise dos resíduos.
m2 <- as.lm(n2)
par(mfrow=c(2,2)); plot(m2); layout(1)
##-----------------------------------------------------------------------------
## Colocar bandas de confiança.
f <- function(theta, xx){
with(as.list(theta), th0+th1*(xx-th2)*(xx>=th2)+0*(xx<th2))
}
## gradient(f, x=coef(n1), xx=c(0.0, 0.2, 0.4))
pred <- data.frame(Temp=sort(c(seq(10, 80, l=100),
coef(n2)["th2"]+c(-0.001,0,0.001))))
pred$fit <- predict(n2, newdata=pred)
der <- gradient(f, x=coef(n2), xx=pred$Temp)
str(der)
## num [1:103, 1:3] 1 1 1 1 1 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "th0" "th1" "th2"
F <- der
U <- chol(vcov(n2))
pred$se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
tval <- qt(p=c(lwr=0.025, upr=0.975), df=df.residual(n2))
me <- outer(pred$se, tval, "*")
pred <- cbind(pred, sweep(me, 1, pred$fit, "+"))
str(pred)
## 'data.frame': 103 obs. of 5 variables:
## $ Temp: num 10 10.7 11.4 12.1 12.8 ...
## $ fit : num 74.7 74.7 74.7 74.7 74.7 ...
## $ se : num 1.34 1.34 1.34 1.34 1.34 ...
## $ lwr : num 72 72 72 72 72 ...
## $ upr : num 77.4 77.4 77.4 77.4 77.4 ...
## Equação do modelo ajustado.
coef(n2)
## th0 th1 th2
## 74.6953 0.5674 41.9512
formula(n2)
## C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) + 0 * (Temp < th2)
## Arredonda as estimativas.
theta <- mapply(round,
x=coef(n2), digits=c(2,4,2),
SIMPLIFY=FALSE)
theta
## $th0
## [1] 74.7
##
## $th1
## [1] 0.5674
##
## $th2
## [1] 41.95
## Equação.
formula(n2)
## C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) + 0 * (Temp < th2)
eq <- substitute(
expr=c(
expression(C==th0~", se"~Temp<th2),
expression(C==th0+th1*(Temp-th2)~", se"~Temp>=th2)),
env=theta)
eval(eq)
## expression(C == 74.7 ~ ", se" ~ Temp < 41.95, C == 74.7 + 0.5674 *
## (Temp - 41.95) ~ ", se" ~ Temp >= 41.95)
## Observados, preditos e a banda de confiança.
xyplot(C~Temp, data=segreg)+
as.layer(xyplot(fit~Temp, data=pred, type="l",
prepanel=prepanel.cbH, cty="bands",
ly=pred$lwr, uy=pred$upr,
panel=panel.cbH))+
layer(panel.key(points=FALSE, text=eval(eq),
corner=c(0.05,0.95)))
Tamanho em função da idade para peixes
##-----------------------------------------------------------------------------
## Tamanho em função da idade dos peixes.
str(lakemary)
## 'data.frame': 78 obs. of 2 variables:
## $ Age : int 1 1 2 2 2 2 3 3 3 3 ...
## $ Length: int 67 62 109 83 91 88 137 131 122 122 ...
xyplot(Length~Age, data=lakemary, xlim=c(0,NA), ylim=c(0,NA))
## Ajustar o modelo von Bertalanffy.
## Original: tha*(1-exp(-thk*(x-thi)));
## Alternativo: th0+(tha-th0)*(1-2^(-x/thv));
start <- list(th0=30, tha=250, thv=3)
xyplot(Length~Age, data=lakemary)+
layer(with(start,
panel.curve(th0+(tha-th0)*(1-2^(-x/thv)), col=2)))
n3 <- nls(Length~th0+(tha-th0)*(1-2^(-Age/thv)),
data=lakemary, start=start)
summary(n3)
##
## Formula: Length ~ th0 + (tha - th0) * (1 - 2^(-Age/thv))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## th0 -6.440 20.382 -0.32 0.75
## tha 192.811 13.080 14.74 < 2e-16 ***
## thv 1.706 0.371 4.59 1.7e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11 on 75 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 6.05e-06
confint(n3)
## Waiting for profiling to be done...
## 2.5% 97.5%
## th0 -50.514 28.827
## tha 174.387 233.351
## thv 1.195 2.889
confint.default(n3)
## 2.5 % 97.5 %
## th0 -46.3877 33.508
## tha 167.1739 218.448
## thv 0.9781 2.434
xyplot(Length~Age, data=lakemary, xlim=c(0,NA), ylim=c(0,NA))+
layer(with(as.list(coef(n3)),
panel.curve(th0+(tha-th0)*(1-2^(-x/thv)), col=2)))
Abundância de peixes
##-----------------------------------------------------------------------------
str(swan96)
## 'data.frame': 27 obs. of 2 variables:
## $ Day : int 0 1 2 34 35 36 37 55 56 57 ...
## $ LCPUE: num 0.368 0.636 0 1.358 0.847 ...
xyplot(LCPUE~Day, data=swan96)
## Swan Lake, Minnesota, 1996.
## LCPUE: log da quantidade de peixes por unidade de esforço.
## Day: dia de pesça.
m0 <- lm(LCPUE~Day+I(Day^2), swan96)
c0 <- coef(m0)
names(c0) <- sprintf("b%i", 0:2)
xyplot(LCPUE~Day, data=swan96)+
layer(with(as.list(c0),
panel.curve(b0+b1*x+b2*x^2, col=2)))
## Estimar o dia correspondente ao máximo.
with(as.list(c0), -0.5*b1/b2)
## [1] 183.1
## E o intervalo de confiança? Método delta?
n4 <- nls(LCPUE~thy+thc*(Day-thx)^2, data=swan96,
start=list(thy=4, thx=180, thc=c0["b2"]))
summary(n4)
##
## Formula: LCPUE ~ thy + thc * (Day - thx)^2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## thy 4.36e+00 2.06e-01 21.14 <2e-16 ***
## thx 1.83e+02 5.96e+00 30.72 <2e-16 ***
## thc -1.27e-04 1.68e-05 -7.58 8e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.572 on 24 degrees of freedom
##
## Number of iterations to convergence: 2
## Achieved convergence tolerance: 3.3e-09
confint(n4)
## Waiting for profiling to be done...
## 2.5% 97.5%
## thy 3.936e+00 4.787e+00
## thx 1.726e+02 1.987e+02
## thc -1.619e-04 -9.262e-05
confint.default(n4)
## 2.5 % 97.5 %
## thy 3.955e+00 4.763e+00
## thx 1.714e+02 1.948e+02
## thc -1.602e-04 -9.437e-05
##-----------------------------------------------------------------------------
## Ajuste do modelo linear com dois segmentos.
## Modelo:
## f(x) = th0+ th1*x, se x<th2
## = (th0+th1*th2)+th3*(x-th2), se x>=th2
start <- list(th0=0.3, th1=3/100, th2=150, th3=-1/100)
xyplot(LCPUE~Day, data=swan96)+
layer(with(start,
panel.curve(th0+th1*x*(x<th2)+
(th1*th2+th3*(x-th2))*(x>=th2), col=2)))
n5 <- nls(LCPUE~th0+th1*Day*(Day<th2)+(th1*th2+th3*(Day-th2))*(Day>=th2),
data=swan96, start=start)
summary(n5)
##
## Formula: LCPUE ~ th0 + th1 * Day * (Day < th2) + (th1 * th2 + th3 * (Day -
## th2)) * (Day >= th2)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## th0 0.08259 0.23014 0.36 0.723
## th1 0.03789 0.00415 9.13 4.2e-09 ***
## th2 103.62090 8.70615 11.90 2.6e-11 ***
## th3 -0.00379 0.00181 -2.09 0.048 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.469 on 23 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 2.5e-07
confint(n5)
## Waiting for profiling to be done...
## Error: step factor 0.000488281 reduced below 'minFactor' of 0.000976562
confint.default(n5)
## 2.5 % 97.5 %
## th0 -0.36848 5.337e-01
## th1 0.02976 4.603e-02
## th2 86.55716 1.207e+02
## th3 -0.00734 -2.328e-04
xyplot(LCPUE~Day, data=swan96)+
layer(with(as.list(coef(n5)),
panel.curve(th0+th1*x*(x<th2)+
(th1*th2+th3*(x-th2))*(x>=th2), col=2)))
##-----------------------------------------------------------------------------
## Regiões de confiança conjunta por simulação: aceitação/rejeição.
rc <- nlsConfRegions(n5, exp=1.5, length=500)
plot(rc)
formula(n4)
rc <- nlsConfRegions(n4, exp=2.5, length=500)
plot(rc)
rc <- nlsConfRegions(n2, exp=2.5, length=500)
plot(rc)
