|
Extensões de Modelos de Regressão
|
Tópicos em Modelos Não Lineares
Prof. Walmes M. Zeviani & Prof. Paulo Justiniano R. Jr
2017-09-04
1 Ajuste do modelo
# Lista de pacotes.
pkg <- c("lattice",
"latticeExtra",
"reshape",
"car",
"alr3",
"nlme",
"nls2",
"nlstools",
"rootSolve",
"plyr",
"wzRfun")
sapply(pkg,
FUN = library,
character.only = TRUE,
logical.return = TRUE)#-----------------------------------------------------------------------
# Ajuste de modelo de regressão não linear.
# Estrutura.
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(panel.curve(th0 + th1 * x/(th2 + x)),
data = start)
#-----------------------------------------------------------------------
# Ajuste.
n0 <- nls(Gain ~ th0 + th1 * A/(th2 + A),
data = turk0,
start = start,
trace = FALSE)
# Resumo do ajuste.
summary(n0)##
## Formula: Gain ~ th0 + th1 * A/(th2 + A)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## th0 622.24167 6.20283 100.316 < 2e-16 ***
## th1 235.21690 22.56028 10.426 8.09e-12 ***
## th2 0.15476 0.04059 3.812 0.000591 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.14 on 32 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 1.151e-06# Verificando.
xyplot(Gain ~ A, data = turk0) +
layer(panel.curve(th0 + th1 * x/(th2 + x)),
data = as.list(coef(n0)))
1.1 Intervalos de confiança para os parâmetros
# Intervalos de Wald: são baseado na aproximação quadrática da
# verossimilhança, conhecido como intervalos de Wald ou
# assintóticos. São simétricos por construção.
wci <- confint.default(n0)
wci## 2.5 % 97.5 %
## th0 610.08434930 634.3989826
## th1 190.99955017 279.4342411
## th2 0.07519662 0.2343179# Intervalos de confiança baseados no perfil da log-verossimilhança.
pci <- confint(n0)
pci## 2.5% 97.5%
## th0 609.656524 634.7143666
## th1 196.536953 290.9141161
## th2 0.094081 0.2665297# Perfil.
par(mfrow = c(2, 2))
plot(profile(n0))
layout(1)
# Intervalo de confiança bootstrap (bootstrap não paramétrico).
bci <- nlsBoot(nls = n0, niter = 999)
summary(bci)##
## ------
## Bootstrap statistics
## Estimate Std. error
## th0 622.2094327 5.7903030
## th1 237.8503806 21.8630652
## th2 0.1610904 0.0413194
##
## ------
## Median of bootstrap estimates and percentile confidence intervals
## Median 2.5% 97.5%
## th0 622.1983058 610.5366995 633.4165784
## th1 235.8907183 201.2129137 284.8400970
## th2 0.1535691 0.0975758 0.25524631.2 Regiões de confiança
#-----------------------------------------------------------------------
# Região de confiança para os parâmetros.
# Contornos baseados na função residual sum of squares.
rss <- nlsContourRSS(n0)## 0%33%66%100%
## RSS contour surface array returnedplot(rss)
# Regiões baseadas em aceitação-rejeição.
cr <- nlsConfRegions(n0, exp = 2)## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%100%
## Confidence regions array returnedplot(cr)
1.3 Inferência para funções dos parâmetros
Método Delta é um procedimento para aproximar a distribuição de uma função de um vetor aleatório através de uma expansão de Taylor. É simples, mas útil, para deduzir a distribuição assintótica de variáveis aleatórias. O Teorema do Limite Central pode ser considerado como um caso particular do Método Delta. Para aplicar o método delta precisamos da função \(f(\theta)\), sua deriavada de primeira order \(f'(\theta)\), das estimativas \(\hat{\theta}\) e matriz de covariância das estimativas \(\text{Var}(\hat{\theta})\).
Vamos considerar para demostração que desejamos inferir sobre o tempo de \(q = 0.9\) de vida, ou seja, o valor em \(x = \theta_q\) no qual a função vale \(\theta_0 + 0.9\theta_1\).
\[ \begin{align*} \theta_0 + q\theta_1 &= \theta_0 + \frac{\theta_1 x}{\theta_2 + x}\\ q \theta_1 &= \frac{\theta_1 x}{\theta_2 + x}\\ q &= \frac{x}{\theta_2 + x}\\ x &= q\theta_2 + qx\\ x - qx &= q\theta_2\\ x &= \frac{q\theta_2}{1 - q} \end{align*} \]
Portanto, as funções necessárias são \[ \begin{align*} g(\theta) &= \left(\frac{q}{1 - q}\right)\theta_2\\ g'(\theta) &= \left(\frac{q}{1 - q}\right). \end{align*} \]
#-----------------------------------------------------------------------
# Método delta.
# Usando função pronta do pacote `car`.
dm <- deltaMethod(n0,
g = "(0.9/0.1) * th2")
dm## Estimate SE 2.5 % 97.5 %
## (0.9/0.1) * th2 1.392815 0.3653363 0.6767695 2.108861# Fazendo passo a passo.
v <- vcov(n0) # Matriz de covariâncias.
gval <- cbind((0.9/0.1) * coef(n0)["th2"]) # Valor de g().
gr <- cbind(th0 = 0, th1 = 0, th2 = (0.9/0.1)) # Valor de g'().
# Erro padrão de g(theta).
sqrt(gr %*% v %*% t(gr))## [,1]
## [1,] 0.3653363#-----------------------------------------------------------------------
# Usando bootstrap.
str(bci)## List of 4
## $ coefboot: num [1:999, 1:3] 626 623 616 617 625 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:3] "th0" "th1" "th2"
## $ rse : num [1:999] 20.9 18.6 18.5 17 21.1 ...
## $ bootCI : num [1:3, 1:3] 622.198 235.891 0.154 610.537 201.213 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "th0" "th1" "th2"
## .. ..$ : chr [1:3] "Median" "2.5%" "97.5%"
## $ estiboot: num [1:3, 1:2] 622.209 237.85 0.161 5.79 21.863 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "th0" "th1" "th2"
## .. ..$ : chr [1:2] "Estimate" "Std. error"
## - attr(*, "class")= chr "nlsBoot"# Função avaliada nas estimativas bootstrap.
bg <- bci$coefboot[, "th2"] * (0.9/0.1)
# Erro padrão.
sd(bg)## [1] 0.3718746# IC percentil
quantile(bg, probs = c(0.025, 0.975))## 2.5% 97.5%
## 0.8781822 2.2972166#-----------------------------------------------------------------------
# Usando os valores internos da região de confiança conjunta.
# DANGER: existe preocupação com relação ao uso dessa abordagem. A
# primeira refere-se ao nível nominal de cobertura desse IC. Erro padrão
# não deve ser calculado pois a distribuição dos valores é uniforme
# dentro da RC e não representam, portanto, a distribuição do estimador.
range(cr$cr[, "th2"] * (0.9/0.1))## [1] 0.7043424 2.87103071.4 Bandas de confiança
Bandas de confiança são um caso particular de função das estimativas dos parâmetros. A forma de obter uma região de confiança é via aproximação linear, ou seja, método delta.
# Modelo escrito como função dos parâmetros (theta).
formula(n0)## Gain ~ th0 + th1 * A/(th2 + A)f <- function(theta, A) {
with(as.list(theta), th0 + th1 * A/(th2 + A))
}
# Matriz de derivadas parciais em \hat{\theta} (n x p).
gradient(f, x = coef(n0), A = c(0, 0.2, 0.4))## th0 th1 th2
## [1,] 1 0.0000000 0.000
## [2,] 1 0.5637657 -373.797
## [3,] 1 0.7210361 -305.719# IMPORTANT: esse gradiente é obtido numericamente. Você pode definir
# uma função que retorne o gradiante usando as derivadas
# analíticas. Isse deve ser mais seguro e mais barato
# computacionalmente.
pred <- data.frame(A = seq(0, 0.5, length.out = 31))
# Valores ajustados (g(\theta)).
# pred$fit <- f(theta = coef(n1), xx = pred$A)
pred$fit <- predict(n0, newdata = pred)
# Matriz gradiente (um gradiente para cada observação, g'(\theta)).
grad <- gradient(f, x = coef(n0), A = pred$A)
str(grad)## num [1:31, 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.
U <- chol(vcov(n0))
pred$se <- sqrt(apply(grad %*% t(U), 1, function(x) sum(x^2)))
tval <- qt(p = c(lwr = 0.025, upr = 0.975), df = df.residual(n0))
me <- outer(pred$se, tval, "*")
pred <- cbind(pred, sweep(me, 1, pred$fit, "+"))
str(pred)## 'data.frame': 31 obs. of 5 variables:
## $ A : num 0 0.0167 0.0333 0.05 0.0667 ...
## $ fit: num 622 645 664 680 693 ...
## $ se : num 6.2 4.74 4.99 5.45 5.73 ...
## $ lwr: num 610 635 654 669 681 ...
## $ upr: num 635 655 674 691 705 ...# Equação do modelo ajustado.
coef(n0)## th0 th1 th2
## 622.2416659 235.2168957 0.1547573formula(n0)## Gain ~ th0 + th1 * A/(th2 + A)# Arredonda as estimativas.
theta <- mapply(round,
x = coef(n0),
digits = c(1, 1, 2),
SIMPLIFY = FALSE)
theta## $th0
## [1] 622.2
##
## $th1
## [1] 235.2
##
## $th2
## [1] 0.15# Equação.
eq <- substitute(
expr = expression(Gain == th0 + th1 * A/(th2 + A)),
env = theta)
eval(eq)## expression(Gain == 622.2 + 235.2 * A/(0.15 + A))# 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)))
# Adiciona as bandas do bootstrap.
xyplot(Gain ~ A, data = turk0) +
layer({
th <- bci$coefboot
for (i in 1:nrow(th)) {
panel.curve(th[i, "th0"] +
th[i, "th1"] * x/(th[i, "th2"] + x),
col = rgb(0, 1, 0, 0.05),
from = 0,
to = 0.5)
}
}) +
as.layer(xyplot(fit ~ A,
data = pred,
type = "l",
cty = "bands",
ly = pred$lwr,
uy = pred$upr,
panel = panel.cbH))
2 Ajuste interativo com a rp.nls
# help(rp.nls, help_type = "html")
library(rpanel)
data(turk0, package = "alr3")
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 ...plot(Gain ~ A,
data = turk0,
xlab = "Metionine",
ylab = "Weight gain")
turk.fit <- rp.nls(model = Gain ~ Int + (Asy - Int) * A/(Half + A),
data = turk0,
start = list(Int = c(600, 650),
Asy = c(750, 850),
Half = c(0, 0.2)),
extra = function(Int, Asy, Half) {
abline(h = c(Asy, Int),
v = Half,
col = "green")
},
xlab = "Metionine",
ylab = "Weight gain")
summary(turk.fit)
confint(turk.fit)
f <- formula(turk.fit)
plot(Gain ~ A,
data = turk0,
xlab = "Metionine",
ylab = "Weight gain")
with(as.list(coef(turk.fit)), {
eval(call("curve", f[[3]],
add = TRUE,
xname = intersect(all.vars(f[-2]),
names(turk0))))
})#-----------------------------------------------------------------------
# Um exemplo com grupos.
xyplot(rate ~ conc,
groups = state,
data = Puromycin,
type = c("p", "smooth"),
span = 1,
auto.key = TRUE)
Puro.fit <- rp.nls(model = rate ~
Int + (Top - Int) * conc/(Half + conc),
data = Puromycin,
start = list(Int = c(20, 70),
Top = c(100, 200),
Half = c(0, 0.6)),
subset = "state",
gpar = list(try = list(col = 2, lty = 2),
fit = list(col = "blue", lwd = 1.5)),
xlab = "Concentration", ylab = "Rate",
xlim = c(0, 1.2), ylim = c(40, 220))
length(Puro.fit)
sapply(Puro.fit, coef)
sapply(Puro.fit, logLik)
sapply(Puro.fit, deviance)#-----------------------------------------------------------------------
# Ajuste de curvas de retenção de água do solo.
cra <- read.table("http://www.leg.ufpr.br/~walmes/data/cra_manejo.txt",
header = TRUE, sep = "\t")
cra$tens[cra$tens == 0] <- 0.1
cras <- subset(cra, condi == "LVA3,5")
cras <- aggregate(umid ~ posi + prof + tens, data = cras, FUN = mean)
cras$caso <- with(cras, interaction(posi, prof))
cras$ltens <- log(cras$tens)
xyplot(umid ~ ltens | posi,
groups = prof,
data = cras,
type = c("p", "a"))
# modelo: van Genuchten com retrição de Mualem.
# x: representado por ltens (log da tensão matricial, psi).
# y: representado por umid, conteúdo de água do solo (%).
# th1: assíntota inferior, mínimo da função, quando x -> +infinito.
# th2: assíntota superior, máximo da função, quando x -> -infinito.
# th3: locação, translada o ponto de inflexão.
# th4: forma, altera a taxa ao redor da inflexão.
model <- umid ~ Ur + (Us - Ur)/(1 + exp(n * (al + ltens)))^(1 - 1/n)
start <- list(Ur = c(init = 0.2,
from = 0,
to = 0.5),
Us = c(init = 0.6,
from = 0.4,
to = 0.8),
al = c(1, -2, 4),
n = c(1.5, 1, 4))
cra.fit <- rp.nls(model = model,
data = cras,
start = start,
subset = "caso")
t(sapply(cra.fit, coef))#-----------------------------------------------------------------------
# Curva de produção em função da desfolha do algodão.
cap <- read.table("http://www.leg.ufpr.br/~walmes/data/algodão.txt",
header = TRUE, sep = "\t", encoding = "latin1")
cap$desf <- cap$desf/100
cap <- subset(cap, select = c(estag, desf, pcapu))
cap$estag <- factor(cap$estag,
labels = c("vegetativo", "botão floral",
"florescimento", "maçã", "capulho"))
str(cap)## 'data.frame': 125 obs. of 3 variables:
## $ estag: Factor w/ 5 levels "vegetativo","botão floral",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ desf : num 0 0 0 0 0 0.25 0.25 0.25 0.25 0.25 ...
## $ pcapu: num 33.2 28.7 31.5 28.9 36.4 ...xyplot(pcapu ~ desf | estag,
data = cap,
layout = c(5, 1),
xlab = "Nível de desfolha artificial",
ylab = "Peso de capulhos")
# modelo: potência.
# x: representado por desf (nível de desfolha artifical).
# y: representado por pcapu (peso de capulhos), produto do algodão.
# th1: intercepto, valor da função quando x=0 (situação ótima).
# th2: diferença no valor da função para x=0 e x=1 (situação extrema).
# th3: forma, indica como a função decresce, se th3=0 então função linear.
model <- pcapu ~ f0 - delta * desf^exp(curv)
start <- list(f0 = c(30, 25, 35),
delta = c(8, 0, 16),
curv = c(0, -2, 4))
cap.fit <- rp.nls(model = model,
data = cap,
start = start,
subset = "estag")
model <- pcapu ~ f0 - f1 * desf^((log(5) - log(f1))/log(xde))
start <- list(f0 = c(30, 25, 35),
f1 = c(8, 0, 16),
xde = c(0.5, 0, 1))
cap.fit <- rp.nls(model = model,
data = cap,
start = start,
subset = "estag",
assignTo = "cap.fit",
extraplot = function(f0, f1, xde) {
abline(v = xde,
h = c(f0, f0 - f1),
lty = 2, col = 2)
})
length(cap.fit)
sapply(cap.fit, coef)
lapply(cap.fit, summary)3 Ajuste do modelo em lote
# dput(rowMeans(sapply(cra.fit, coef)))
start <- structure(c(0.259625743051275,
0.574758577621488,
1.29703804811817,
-0.187699896144559),
.Names = c("Ur", "Us", "n", "al"))
n00 <- nlsList(umid ~ Ur + (Us - Ur)/
(1 + exp(n * (al + ltens)))^(1 - 1/n) | caso,
data = cras,
start = start)
coef(n00)## Ur Us n al
## EL.20 0.3605169 0.5765759 1.439755 -0.26933232
## L.20 0.3137652 0.6585131 1.442087 0.99563732
## EL.40 0.3480770 0.5740054 1.385218 -0.35662783
## L.40 0.3281452 0.6057806 1.523760 0.18444816
## EL.80 0.2668528 0.5847261 1.150290 0.07666282
## L.80 0.3584887 0.5391944 1.472546 -0.51452771
## EL.120 0.2539915 0.5445529 1.159476 -1.05656438
## L.120 0.2960853 0.5450071 1.145878 -0.45324686
## EL.160 -0.3297885 0.5730737 1.032742 -0.17350316
## L.160 0.4001424 0.5461564 1.218629 -0.30996410summary(n00)## Call:
## Model: umid ~ Ur + (Us - Ur)/(1 + exp(n * (al + ltens)))^(1 - 1/n) | caso
## Data: cras
##
## Coefficients:
## Ur
## Estimate Std. Error t value Pr(>|t|)
## EL.20 0.3605169 0.010057241 35.8465048 3.186134e-09
## L.20 0.3137652 0.008152311 38.4878823 2.074393e-08
## EL.40 0.3480770 0.012896217 26.9906290 4.238990e-08
## L.40 0.3281452 0.007102429 46.2018356 8.240905e-09
## EL.80 0.2668528 0.073303673 3.6403747 8.533414e-03
## L.80 0.3584887 0.009593639 37.3673290 2.499936e-10
## EL.120 0.2539915 0.093451255 2.7179039 1.026747e-01
## L.120 0.2960853 0.091527727 3.2349246 6.939301e-03
## EL.160 -0.3297885 1.812500736 -0.1819522 9.147074e-01
## L.160 0.4001424 0.038014115 10.5261540 3.811227e-06
## Us
## Estimate Std. Error t value Pr(>|t|)
## EL.20 0.5765759 0.007623101 75.63535 3.626719e-11
## L.20 0.6585131 0.010474804 62.86640 1.099543e-09
## EL.40 0.5740054 0.007535098 76.17757 8.486112e-11
## L.40 0.6057806 0.008105641 74.73568 4.622320e-10
## EL.80 0.5847261 0.008380935 69.76860 4.219218e-10
## L.80 0.5391944 0.007279014 74.07519 4.126993e-12
## EL.120 0.5445529 0.006691655 81.37791 1.837328e-09
## L.120 0.5450071 0.007563875 72.05395 1.302356e-10
## EL.160 0.5730737 0.008121304 70.56425 1.013153e-08
## L.160 0.5461564 0.007710834 70.82974 4.357888e-11
## n
## Estimate Std. Error t value Pr(>|t|)
## EL.20 1.439755 0.08495393 16.947481 2.797433e-07
## L.20 1.442087 0.05604556 25.730621 2.293179e-07
## EL.40 1.385218 0.08048934 17.209956 6.185113e-07
## L.40 1.523760 0.06798641 22.412710 6.165847e-07
## EL.80 1.150290 0.06528229 17.620246 1.558735e-06
## L.80 1.472546 0.10747671 13.701067 1.012889e-07
## EL.120 1.159476 0.09092132 12.752524 1.034334e-04
## L.120 1.145878 0.09518905 12.037919 5.598647e-06
## EL.160 1.032742 0.07580939 13.622876 1.589559e-04
## L.160 1.218629 0.13158926 9.260857 8.075290e-06
## al
## Estimate Std. Error t value Pr(>|t|)
## EL.20 -0.26933232 0.2720137 -0.9901424 0.196816363
## L.20 0.99563732 0.2473118 4.0258385 0.006963389
## EL.40 -0.35662783 0.2901567 -1.2290869 0.171502396
## L.40 0.18444816 0.1919174 0.9610809 0.386954741
## EL.80 0.07666282 0.5051586 0.1517599 0.877989257
## L.80 -0.51452771 0.3027482 -1.6995237 0.010642627
## EL.120 -1.05656438 0.5091227 -2.0752645 0.192171874
## L.120 -0.45324686 0.6486830 -0.6987185 0.418336670
## EL.160 -0.17350316 0.7963848 -0.2178635 0.897980505
## L.160 -0.30996410 0.7519629 -0.4122067 0.554456588
##
## Residual standard error: 0.007388536 on 60 degrees of freedom#-----------------------------------------------------------------------
n00 <- nlsList(umid ~ Ur + (Us - Ur)/
(1 + exp(n * (al + ltens)))^(1 - 1/n) | caso,
data = cras,
start = start)4 Ajuste de curvas com estrutura experimental
4.1 Desfolha do algodão
#-----------------------------------------------------------------------
# Ajuste dos modelos separados por um fator categórico.
library(latticeExtra)
library(nlme)
xyplot(pcapu ~ desf | estag,
data = cap,
layout = c(5, 1),
xlab = "Nível de desfolha artificial",
ylab = "Peso de capulhos")
# ~f0 - f1 * desf^((log(5) - log(f1))/log(xde))
# ~f0 - f1 * desf^exp(curv)
# A nlsList() é basicamente a nls() dentro de um lapply().
n00 <- nlsList(pcapu ~ f0 - f1 *
desf^((log(5) - log(f1))/log(xde)) | estag,
data = cap,
start = c(f0 = 30, f1 = 8, xde = 0.8))
coef(n00)## f0 f1 xde
## vegetativo 31.14929 10.166278 0.8538319
## botão floral 29.26289 6.117394 0.9440160
## florescimento 28.50115 13.014666 0.2021585
## maçã 28.26343 21.608571 0.7243578
## capulho 34.15315 8.277792 0.4954646# Como prover valores iniciais a partir deste ajuste.
L <- cbind(1, contrasts(cap$estag))
L## botão floral florescimento maçã capulho
## vegetativo 1 0 0 0 0
## botão floral 1 1 0 0 0
## florescimento 1 0 1 0 0
## maçã 1 0 0 1 0
## capulho 1 0 0 0 1# \mu = L \beta, então \beta = L^{-1} \mu.
start <- solve(L) %*% as.matrix(coef(n00))
c(start)## [1] 31.14929457 -1.88640219 -2.64814939 -2.88586070 3.00385297
## [6] 10.16627822 -4.04888413 2.84838777 11.44229229 -1.88848629
## [11] 0.85383191 0.09018408 -0.65167345 -0.12947410 -0.35836736#-----------------------------------------------------------------------
# Ajuste conjunto.
n0 <- gnls(model = pcapu ~ f0 - f1 * desf^((log(5) - log(f1))/log(xde)),
params = f0 + f1 + xde ~ estag,
# start = c(30, 0, 0, 0, 0,
# 8, 0, 0, 0, 0,
# 0.1, 0, 0, 0, 0),
start = c(start),
data = cap)
summary(n0)## Generalized nonlinear least squares fit
## Model: pcapu ~ f0 - f1 * desf^((log(5) - log(f1))/log(xde))
## Data: cap
## AIC BIC logLik
## 687.3049 732.5579 -327.6524
##
## Coefficients:
## Value Std.Error t-value p-value
## f0.(Intercept) 31.149266 1.0413031 29.913735 0.0000
## f0.estagbotão floral -1.886383 1.5260306 -1.236137 0.2190
## f0.estagflorescimento -2.648075 1.8939030 -1.398211 0.1649
## f0.estagmaçã -2.885844 1.4707311 -1.962183 0.0523
## f0.estagcapulho 3.003895 1.8900426 1.589327 0.1149
## f1.(Intercept) 10.166268 1.8711175 5.433260 0.0000
## f1.estagbotão floral -4.048876 2.6572896 -1.523687 0.1305
## f1.estagflorescimento 2.848397 2.7639130 1.030567 0.3050
## f1.estagmaçã 11.442298 2.6457253 4.324825 0.0000
## f1.estagcapulho -1.888474 2.7631564 -0.683448 0.4958
## xde.(Intercept) 0.853835 0.0809136 10.552434 0.0000
## xde.estagbotão floral 0.090181 0.1268542 0.710905 0.4786
## xde.estagflorescimento -0.651684 0.1549329 -4.206236 0.0001
## xde.estagmaçã -0.129476 0.1005280 -1.287964 0.2005
## xde.estagcapulho -0.358374 0.2776476 -1.290752 0.1995
##
## Correlation:
## f0.(I) f0stgf f0.stgf f0stgç f0.stgc f1.(I) f1stgf
## f0.estagbotão floral -0.682
## f0.estagflorescimento -0.550 0.375
## f0.estagmaçã -0.708 0.483 0.389
## f0.estagcapulho -0.551 0.376 0.303 0.390
## f1.(Intercept) 0.533 -0.363 -0.293 -0.377 -0.293
## f1.estagbotão floral -0.375 0.541 0.206 0.266 0.207 -0.704
## f1.estagflorescimento -0.361 0.246 0.661 0.255 0.199 -0.677 0.477
## f1.estagmaçã -0.377 0.257 0.207 0.532 0.208 -0.707 0.498
## f1.estagcapulho -0.361 0.246 0.198 0.255 0.654 -0.677 0.477
## xde.(Intercept) -0.708 0.483 0.389 0.501 0.390 -0.530 0.373
## xde.estagbotão floral 0.452 -0.721 -0.248 -0.320 -0.249 0.338 -0.728
## xde.estagflorescimento 0.370 -0.252 -0.730 -0.262 -0.204 0.277 -0.195
## xde.estagmaçã 0.570 -0.389 -0.313 -0.677 -0.314 0.427 -0.300
## xde.estagcapulho 0.206 -0.141 -0.113 -0.146 -0.806 0.154 -0.109
## f1.stgf f1stgç f1.stgc xd.(I) xdstgf xd.stgf xdstgç
## f0.estagbotão floral
## f0.estagflorescimento
## f0.estagmaçã
## f0.estagcapulho
## f1.(Intercept)
## f1.estagbotão floral
## f1.estagflorescimento
## f1.estagmaçã 0.479
## f1.estagcapulho 0.458 0.479
## xde.(Intercept) 0.359 0.375 0.359
## xde.estagbotão floral -0.229 -0.239 -0.229 -0.638
## xde.estagflorescimento -0.434 -0.196 -0.187 -0.522 0.333
## xde.estagmaçã -0.289 -0.445 -0.289 -0.805 0.513 0.420
## xde.estagcapulho -0.105 -0.109 -0.576 -0.291 0.186 0.152 0.235
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.43631478 -0.59093061 -0.03534145 0.60230115 2.92708061
##
## Residual standard error: 3.547339
## Degrees of freedom: 125 total; 110 residual# Estimativas dos parâmetros.
kronecker(diag(3), L) %*% cbind(coef(n0))## [,1]
## [1,] 31.1492661
## [2,] 29.2628832
## [3,] 28.5011907
## [4,] 28.2634222
## [5,] 34.1531615
## [6,] 10.1662678
## [7,] 6.1173915
## [8,] 13.0146649
## [9,] 21.6085661
## [10,] 8.2777934
## [11,] 0.8538355
## [12,] 0.9440167
## [13,] 0.2021511
## [14,] 0.7243590
## [15,] 0.4954611coef(n00)## f0 f1 xde
## vegetativo 31.14929 10.166278 0.8538319
## botão floral 29.26289 6.117394 0.9440160
## florescimento 28.50115 13.014666 0.2021585
## maçã 28.26343 21.608571 0.7243578
## capulho 34.15315 8.277792 0.4954646# Quadros de testes de Wald.
anova(n0, type = "sequential")## Denom. DF: 110
## numDF F-value p-value
## f0.(Intercept) 1 13064.520 <.0001
## f0.estag 4 124.019 <.0001
## f1.(Intercept) 1 1355.491 <.0001
## f1.estag 4 23.604 <.0001
## xde.(Intercept) 1 379.711 <.0001
## xde.estag 4 6.769 1e-04anova(n0, type = "marginal")## Denom. DF: 110
## numDF F-value p-value
## f0.(Intercept) 1 1016.8540 <.0001
## f0.estag 4 3.4885 0.0101
## f1.(Intercept) 1 33.5458 <.0001
## f1.estag 4 11.5520 <.0001
## xde.(Intercept) 1 126.5385 <.0001
## xde.estag 4 6.7691 0.0001#-----------------------------------------------------------------------
# Vendo o ajuste.
pred <- expand.grid(desf = seq(0, 1, length.out = 30),
estag = levels(cap$estag))
pred$y <- predict(n0, newdata = pred)
xyplot(pcapu ~ desf | estag,
data = cap,
layout = c(5, 1),
xlab = "Nível de desfolha artificial",
ylab = "Peso de capulhos") +
as.layer(xyplot(y ~ desf | estag,
data = pred,
type = "l",
col = 1))
#-----------------------------------------------------------------------
# Bandas de confiança.
# Ajuste do modelo com a parametrização de estimativas por nível.
n1 <- gnls(model = pcapu ~ f0 - f1 * desf^((log(5) - log(f1))/log(xde)),
params = f0 + f1 + xde ~ 0 + estag,
start = unlist(coef(n00)),
data = cap)
coef(n1)## f0.estagvegetativo f0.estagbotão floral f0.estagflorescimento
## 31.1492661 29.2628832 28.5011909
## f0.estagmaçã f0.estagcapulho f1.estagvegetativo
## 28.2634222 34.1531614 10.1662677
## f1.estagbotão floral f1.estagflorescimento f1.estagmaçã
## 6.1173915 13.0146650 21.6085661
## f1.estagcapulho xde.estagvegetativo xde.estagbotão floral
## 8.2777934 0.8538355 0.9440167
## xde.estagflorescimento xde.estagmaçã xde.estagcapulho
## 0.2021511 0.7243590 0.4954611# Verifica que a medida de ajuste é a mesma.
logLik(n0)## 'log Lik.' -327.6524 (df=16)logLik(n1)## 'log Lik.' -327.6524 (df=16)# Retorna f(x, theta), f'(x, theta) e f''(x, theta).
model <- deriv3(expr = ~f0 - f1 * x^((log(5) - log(f1))/log(xde)),
namevec = c("f0", "f1", "xde"),
function.arg = function(x, f0, f1, xde) {
NULL
})
# Obtendo a matriz do modelo: derivadas parciais em theta avaliados em
# cada x_i.
# Obtendo os componentes para um vetor x e valores dos parâmetros.
model(x = c(0.1, 0.5, 0.9), f0 = 30, f1 = 8, xde = 0.5)## [1] 28.32113 25.00000 22.55160
## attr(,"gradient")
## f0 f1 xde
## [1,] 1 0.4872781 7.563348
## [2,] 1 0.0000000 6.780719
## [3,] 1 -0.7895278 1.535399
## attr(,"hessian")
## , , f0
##
## f0 f1 xde
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
##
## , , f1
##
## f0 f1 xde
## [1,] 0 -0.20233784 -0.1836806
## [2,] 0 0.00000000 1.8033688
## [3,] 0 0.01500133 0.5710994
##
## , , xde
##
## f0 f1 xde
## [1,] 0 -0.1836806 -5.553252
## [2,] 0 1.8033688 16.372971
## [3,] 0 0.5710994 5.473148# Grid de valores para predição com bandas.
desf.seq <- seq(0.01, 0.99, length.out = 41)
pred <- expand.grid(desf = desf.seq,
estag = levels(cap$estag),
KEEP.OUT.ATTRS = FALSE)
pred$y <- predict(n1, newdata = pred)
# Cria lista com as estimativas de cada nível separadas.
ths <- split(coef(n1), rep(1:5, times = 3))
# Acerta os nomes.
ths <- lapply(ths, "names<-", c("f0", "f1", "xde"))
ths## $`1`
## f0 f1 xde
## 31.1492661 10.1662677 0.8538355
##
## $`2`
## f0 f1 xde
## 29.2628832 6.1173915 0.9440167
##
## $`3`
## f0 f1 xde
## 28.5011909 13.0146650 0.2021511
##
## $`4`
## f0 f1 xde
## 28.263422 21.608566 0.724359
##
## $`5`
## f0 f1 xde
## 34.1531614 8.2777934 0.4954611mm <- lapply(ths,
FUN = function(i) {
d3 <- do.call(what = model,
args = c(list(x = desf.seq),
as.list(i)))
# Retorna apenas o gradiente.
attr(d3, "gradient")
})
str(mm)## List of 5
## $ 1: num [1:41, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:3] "f0" "f1" "xde"
## $ 2: num [1:41, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:3] "f0" "f1" "xde"
## $ 3: num [1:41, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:3] "f0" "f1" "xde"
## $ 4: num [1:41, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:3] "f0" "f1" "xde"
## $ 5: num [1:41, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:3] "f0" "f1" "xde"# Para criar uma matriz bloco diagonal a partir de lista de matrizes.
X <- as.matrix(Matrix::bdiag(mm))
# DANGER: a ordem das colunas de X deve ser a mesma dos elementos de
# coef(n1).
# Para colocar as colunas da matriz na ordem correta.
j <- c(matrix(seq_along(coef(n1)),
byrow = TRUE,
ncol = length(n1$pmap)))
X <- X[, j]
coef(n1)## f0.estagvegetativo f0.estagbotão floral f0.estagflorescimento
## 31.1492661 29.2628832 28.5011909
## f0.estagmaçã f0.estagcapulho f1.estagvegetativo
## 28.2634222 34.1531614 10.1662677
## f1.estagbotão floral f1.estagflorescimento f1.estagmaçã
## 6.1173915 13.0146650 21.6085661
## f1.estagcapulho xde.estagvegetativo xde.estagbotão floral
## 8.2777934 0.8538355 0.9440167
## xde.estagflorescimento xde.estagmaçã xde.estagcapulho
## 0.2021511 0.7243590 0.4954611# Erro padrão da predição.
U <- chol(vcov(n1))
pred$se <- sqrt(apply(X %*% t(U),
MARGIN = 1,
FUN = function(x) {
sum(x^2)
}))
# Quantil para obter os limites.
tval <- qt(p = c(lwr = 0.025, fit = 0.5, upr = 0.975),
df = nrow(cap) - ncol(X))
# Margem de erro = quantil * se.
me <- outer(pred$se, tval, FUN = "*")
# Acrescenta os limites da banda de confiança.
pred <- cbind(pred,
sweep(me,
MARGIN = 1,
STATS = pred$y,
FUN = "+"))
str(pred)## 'data.frame': 205 obs. of 7 variables:
## $ desf : num 0.01 0.0345 0.059 0.0835 0.108 ...
## $ estag: Factor w/ 5 levels "vegetativo","botão floral",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ y : atomic 31.1 31.1 31.1 31.1 31.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ se : num 1.04 1.04 1.04 1.04 1.04 ...
## $ lwr : num 29.1 29.1 29.1 29.1 29.1 ...
## $ fit : num 31.1 31.1 31.1 31.1 31.1 ...
## $ upr : num 33.2 33.2 33.2 33.2 33.2 ...xyplot(pcapu ~ desf | estag,
data = cap,
col = 1,
layout = c(5, 1),
xlim = c(-0.2, 1.2),
xlab = "Desfolha artifical",
ylab = "Peso de capulhos (g)") +
as.layer(xyplot(fit + lwr + upr ~ desf | estag,
data = pred,
type = "l",
col = 1,
lty = c(1, 2, 2))) +
layer({
xde <- coef(n00)$xde[which.packet()]
f0 <- coef(n00)$f0[which.packet()]
f1 <- coef(n00)$f1[which.packet()]
panel.abline(v = xde,
lty = 2,
col = "orange")
panel.abline(h = c(f0, f0 - f1),
lty = 2,
col = "gray50")
})
4.2 Produção da soja
| Extensões de Modelos de Regressão |
Prof. Paulo Justiniano R. Jr & Prof. Walmes M. Zeviani |
