Modelos de Regressão Não Linear

Fundamentos e Aplicações em R

1 Informação de gradiente

library(lattice)
library(latticeExtra)
library(rootSolve)

llayer <- latticeExtra::layer
#-----------------------------------------------------------------------
# Usar a informação de gradientes, em geral, acelera o procedimento de
# estimação. Para usar informação de gradiente função precisa ter um
# atributo gradiente.

data(wtloss, package = "MASS")
str(wtloss)
## 'data.frame':    52 obs. of  2 variables:
##  $ Days  : int  0 4 7 7 11 18 24 30 32 43 ...
##  $ Weight: num  184 183 180 180 178 ...
p0 <- xyplot(Weight ~ Days, data = wtloss, col = 1,
             xlab = "Dias em dieta",
             ylab = "Peso (kg)")

p0 +
    llayer(with(list(f0 = 81, f1 = 102, K = 141),
                panel.curve(f0 + f1 * 2^(-x/K), add = TRUE, col = 2)))

#-----------------------------------------------------------------------
# Modelo (é uma parametrização do monomolecular para tempo de meia vida)
#   f0 + f1 * 2^(-x/K)
# * f0: assíntota, peso final após estabilização
# * f1: total de peso à perder, f0+f1 é intercepto ou peso inicial
# * K: tempo necessário para se perder metade do passível de perder,
#   0.5*(f1-f0)

# Para passar o gradiente tem-se que conhecer as derivadas do modelo com
# relação à cada parâmetro.
model <- expression(f0 + f1 * 2^(-x/K))
pars <- c("f0","f1","K")
sapply(pars, D, expr = model)
## $f0
## [1] 1
## 
## $f1
## 2^(-x/K)
## 
## $K
## f1 * (2^(-x/K) * (log(2) * (x/K^2)))
#-----------------------------------------------------------------------
# Agora precisa-se fazer uma função que avalie essas expressões em x e
# pars e passe isso como atributo "gradient".

modelfun <- function(x, f0, f1, K){
    value <- f0 + f1 * 2^(-x/K)
    attr(value, "gradient") <-
        cbind(f0 = 1,
              f1 = 2^(-x/K),
              K = f1 * (2^(-x/K) * (log(2) * (x/K^2))))
    return(value)
}

# Avalia a função.
modelfun(x = head(wtloss$Days), f0 = 81, f1 = 102, K = 141)
## [1] 183.0000 181.0139 179.5497 179.5497 177.6308 174.3621
## attr(,"gradient")
##      f0        f1          K
## [1,]  1 1.0000000 0.00000000
## [2,]  1 0.9805283 0.01394786
## [3,]  1 0.9661738 0.02405142
## [4,]  1 0.9661738 0.02405142
## [5,]  1 0.9473608 0.03705915
## [6,]  1 0.9153151 0.05859095
#-----------------------------------------------------------------------
# Agora é só usar.

# Com informação de gradiente analítico.
n0 <- nls(Weight ~ modelfun(Days, f0, f1, K),
          data = wtloss,
          start = list(f0 = 60, f1 = 200, K = 41),
          trace = TRUE)
## 120543.9 :   60 200  41
## 5548.131 :  110.93529  70.05107  53.85617
## 383.9489 :  103.75522  78.74353  94.21701
## 57.24169 :   87.86922  96.00670 129.93709
## 39.34297 :   81.76984 102.28013 141.25126
## 39.2447 :   81.37412 102.68369 141.91071
## 39.2447 :   81.37382 102.68412 141.91036
summary(n0)
## 
## Formula: Weight ~ modelfun(Days, f0, f1, K)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## f0   81.374      2.269   35.86   <2e-16 ***
## f1  102.684      2.083   49.30   <2e-16 ***
## K   141.910      5.295   26.80   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8949 on 49 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 3.134e-08
# Sem informação de gradiente analítico.
n1 <- nls(Weight ~ f0 + f1 * 2^(-Days/K),
          data = wtloss,
          start = list(f0 = 60, f1 = 200, K = 41),
          trace = TRUE)
## 120543.9 :   60 200  41
## 5548.131 :  110.93529  70.05107  53.85617
## 383.9489 :  103.75522  78.74353  94.21701
## 57.24169 :   87.86922  96.00670 129.93709
## 39.34297 :   81.76983 102.28013 141.25126
## 39.2447 :   81.37412 102.68369 141.91071
## 39.2447 :   81.37382 102.68412 141.91037
summary(n1)
## 
## Formula: Weight ~ f0 + f1 * 2^(-Days/K)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## f0   81.374      2.269   35.86   <2e-16 ***
## f1  102.684      2.083   49.30   <2e-16 ***
## K   141.910      5.295   26.80   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8949 on 49 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 3.19e-08
# NOTA: para o presente exemplo não houve diferença em número de passos
# porque a função é muito simples.

#-----------------------------------------------------------------------
# A forma mais rápida de passar informação de gradiente é usar a função
# `deriv3()`. Tem esse nome porque faz as derivadas do modelo com
# relação aos parâmetros e retorna a avaliação de 3 funções avaliadas em
# x e theta: 1) modelo (f(x)), 2) o gradiente (df(x)/dtheta) e 3) o
# hessiano (d²f(x)/dtheta²). As derivadas são analíticas. Sendo assim, a
# `deriv3()` não sabe derivar modelos com segmentados pois não sabe
# interpretar as inequações do tipo x <= x0, por exemplo.

modelfun <- deriv3(expr = ~f0 + f1 * 2^(-x/K),
                   namevec = c("f0", "f1", "K"),
                   function.arg = function(x, f0, f1, K) NULL)

# Veja que retorna 2 atributos: gradient e hessian.
modelfun
## function (x, f0, f1, K) 
## {
##     .expr3 <- 2^(-x/K)
##     .expr6 <- log(2)
##     .expr7 <- K^2
##     .expr9 <- .expr6 * (x/.expr7)
##     .expr10 <- .expr3 * .expr9
##     .value <- f0 + f1 * .expr3
##     .grad <- array(0, c(length(.value), 3L), list(NULL, c("f0", 
##         "f1", "K")))
##     .hessian <- array(0, c(length(.value), 3L, 3L), list(NULL, 
##         c("f0", "f1", "K"), c("f0", "f1", "K")))
##     .grad[, "f0"] <- 1
##     .grad[, "f1"] <- .expr3
##     .hessian[, "f1", "f1"] <- 0
##     .hessian[, "f1", "K"] <- .hessian[, "K", "f1"] <- .expr10
##     .grad[, "K"] <- f1 * .expr10
##     .hessian[, "K", "K"] <- f1 * (.expr10 * .expr9 - .expr3 * 
##         (.expr6 * (x * (2 * K)/.expr7^2)))
##     attr(.value, "gradient") <- .grad
##     attr(.value, "hessian") <- .hessian
##     .value
## }
modelfun(x = head(wtloss$Days), f0 = 81, f1 = 102, K = 141)
## [1] 183.0000 181.0139 179.5497 179.5497 177.6308 174.3621
## attr(,"gradient")
##      f0        f1          K
## [1,]  1 1.0000000 0.00000000
## [2,]  1 0.9805283 0.01394786
## [3,]  1 0.9661738 0.02405142
## [4,]  1 0.9661738 0.02405142
## [5,]  1 0.9473608 0.03705915
## [6,]  1 0.9153151 0.05859095
## attr(,"hessian")
## , , f0
## 
##      f0 f1 K
## [1,]  0  0 0
## [2,]  0  0 0
## [3,]  0  0 0
## [4,]  0  0 0
## [5,]  0  0 0
## [6,]  0  0 0
## 
## , , f1
## 
##      f0 f1            K
## [1,]  0  0 0.0000000000
## [2,]  0  0 0.0001367437
## [3,]  0  0 0.0002357982
## [4,]  0  0 0.0002357982
## [5,]  0  0 0.0003633250
## [6,]  0  0 0.0005744211
## 
## , , K
## 
##      f0           f1             K
## [1,]  0 0.0000000000  0.0000000000
## [2,]  0 0.0001367437 -0.0001958968
## [3,]  0 0.0002357982 -0.0003352850
## [4,]  0 0.0002357982 -0.0003352850
## [5,]  0 0.0003633250 -0.0005114491
## [6,]  0 0.0005744211 -0.0007943076
#-----------------------------------------------------------------------
# Alguns modelos não lineares possuem operadores/funções que não estão
# na tabela de derivadas considerada pela `deriv3()`. Um exemplo é o
# modelo linear-platô que possui os operadores < e >= na função. Dessa
# forma não tem como usar informação de derivada analítica por meio da
# `deriv3()`.  Pode-se fazer à mão ou pode-se empregar derivadas
# numéricas. A função `rootSolve::gradiente()` pode ser usada.

# Dá mensagem de erro!
# mlp <- deriv3(~f0 + tx * x * (x < xde) + tx * xde * (x >= xde) ,
#               c("f0", "tx", "xde") ,
#               function(x, f0, tx, xde) NULL)

f <- function(theta, xi) {
    with(as.list(theta),
         f0 + tx * xi * (xi < xde) + tx * xde * (xi >= xde))
}

f(c(f0 = 10, tx = 1, xde = 3), xi = 1:4)
## [1] 11 12 13 13
# É um gradiente numérico (rootSolve).
gradient(f,
         x = c(f0 = 3, tx = 0.9, xde = 3.5),
         xi = seq(0, 5, by = 0.5))
##       f0  tx xde
##  [1,]  1 0.0 0.0
##  [2,]  1 0.5 0.0
##  [3,]  1 1.0 0.0
##  [4,]  1 1.5 0.0
##  [5,]  1 2.0 0.0
##  [6,]  1 2.5 0.0
##  [7,]  1 3.0 0.0
##  [8,]  1 3.5 0.0
##  [9,]  1 3.5 0.9
## [10,]  1 3.5 0.9
## [11,]  1 3.5 0.9
#-----------------------------------------------------------------------
# Criando uma função cujo gradiente é obtido numericamente.  Existe
# custo em obter gradiente numérico, melhor seria analítico.

modelfun.n <- function(x, f0, tx, xde) {
    f <- function(theta, xi) {
        with(as.list(theta),
             f0 + tx * xi * (xi < xde) + tx * xde * (xi >= xde))
    }
    value <- f(theta = c(f0 = f0, tx = tx, xde = xde), xi = x)
    attr(value, "gradient") <-
        gradient(f,
                 x = c(f0 = f0, tx = tx, xde = xde),
                 xi = x)
    return(value)
}

modelfun.n(x = seq(0, 5, by = 0.5), f0 = 3, tx = 0.9, xde = 3.5)
##  [1] 3.00 3.45 3.90 4.35 4.80 5.25 5.70 6.15 6.15 6.15 6.15
## attr(,"gradient")
##       f0  tx xde
##  [1,]  1 0.0 0.0
##  [2,]  1 0.5 0.0
##  [3,]  1 1.0 0.0
##  [4,]  1 1.5 0.0
##  [5,]  1 2.0 0.0
##  [6,]  1 2.5 0.0
##  [7,]  1 3.0 0.0
##  [8,]  1 3.5 0.0
##  [9,]  1 3.5 0.9
## [10,]  1 3.5 0.9
## [11,]  1 3.5 0.9
#-----------------------------------------------------------------------
# Criando uma função cujo gradiente é obtido analiticamente.

modelfun.a <- function(x, f0, tx, xde){
    f <- function(theta, xi) {
        with(as.list(theta),
             f0 + tx * xi * (xi < xde) + tx * xde * (xi >= xde))
    }
    xo <- x < xde
    value <- f(theta = c(f0 = f0, tx = tx, xde = xde), xi = x)
    attr(value, "gradient") <-
        cbind(1,
              x * (xo) + xde * (!xo),
              0 * (xo) + tx * (!xo))
    return(value)
}

modelfun.n(x = seq(0, 5, by = 0.5), f0 = 3, tx = 0.9, xde = 3.5)
##  [1] 3.00 3.45 3.90 4.35 4.80 5.25 5.70 6.15 6.15 6.15 6.15
## attr(,"gradient")
##       f0  tx xde
##  [1,]  1 0.0 0.0
##  [2,]  1 0.5 0.0
##  [3,]  1 1.0 0.0
##  [4,]  1 1.5 0.0
##  [5,]  1 2.0 0.0
##  [6,]  1 2.5 0.0
##  [7,]  1 3.0 0.0
##  [8,]  1 3.5 0.0
##  [9,]  1 3.5 0.9
## [10,]  1 3.5 0.9
## [11,]  1 3.5 0.9
modelfun.a(x = seq(0, 5, by = 0.5), f0 = 3, tx = 0.9, xde = 3.5)
##  [1] 3.00 3.45 3.90 4.35 4.80 5.25 5.70 6.15 6.15 6.15 6.15
## attr(,"gradient")
##       [,1] [,2] [,3]
##  [1,]    1  0.0  0.0
##  [2,]    1  0.5  0.0
##  [3,]    1  1.0  0.0
##  [4,]    1  1.5  0.0
##  [5,]    1  2.0  0.0
##  [6,]    1  2.5  0.0
##  [7,]    1  3.0  0.0
##  [8,]    1  3.5  0.9
##  [9,]    1  3.5  0.9
## [10,]    1  3.5  0.9
## [11,]    1  3.5  0.9
#-----------------------------------------------------------------------
# Avaliando os custos computacionais.

microbenchmark::microbenchmark(
                    a = modelfun.a(x = seq(0, 5, by = 0.5),
                                   f0 = 3, tx = 0.9, xde = 3.5),
                    n = modelfun.n(x = seq(0, 5, by = 0.5),
                                   f0 = 3, tx = 0.9, xde = 3.5),
                    times = 500)
## Unit: microseconds
##  expr     min       lq      mean  median       uq      max neval cld
##     a  63.735  70.5305  94.79558  73.770  93.1915  554.345   500  a 
##     n 216.808 238.5580 319.48861 259.743 311.1210 1404.394   500   b
#-----------------------------------------------------------------------
# Com isso podemos usar modelos não lineares que possuem funções
# densidade ou acumuladas de distribuições de probabilidade (orgiva na
# normal, beta, gama, etc), funções que possuem integral dentro, funções
# compostas (como linear-platô), etc. O gradiente numérico pode ser
# usado não só para `nls()` mas para `optim()` ou outras funções de
# otmização.

2 Modelos parcialmente lineares

#-----------------------------------------------------------------------
# Um modelo é parcilamente linear quando fixando o valor de alguns
# parâmetros o modelo se torna um modelo linear. Os modelos dessa sessão
# são parcialmente lineares. O modelo para perda de peso, se conhecido
# K, ele se torna a equação da reta. O modelo quadrático centrado, se
# conhecido xc, se torna o modelo quadrático sem termo linear. Pode-se
# tirar vantagem do modelo parcilamente linear para estimação que passa
# ter 2 passos: 1) analítico para porção linear e 2) númerico para
# porção não linear. O custo computacional de um modelo não linear está
# na parte não linear que requer otimização numérica. Deixando parte dos
# parâmetros para serem estimados analíticamente reduz-se o número de
# dimensões da otimização numérica. Quando mais facilitarmos isso
# melhor.

#-----------------------------------------------------------------------
# Ajustando com `algorthm = "plinear"`
#   f0 + f1 * 2^(-x/K)
# se K é conhecido então temos
#   f0 + f1 * g(x)
# que é a equação da reta. O trabalho de otimização se reduz à estimar K
# pois f0 e f1 é calculado por (X'X)^{-1} X'y quando g(x) for conhecido.
# De 3 dimensões passa-se para 1.

data(wtloss, package = "MASS")

p0 <- xyplot(Weight ~ Days, data = wtloss,
             xlab = "Dias em dieta",
             ylab = "Peso (kg)")
p0 +
    llayer(with(list(f0 = 81, f1 = 102, K = 141),
                panel.curve(f0 + f1 * 2^(-x/K),
                            add = TRUE, col = 2)))

n0 <- nls(Weight ~ f0 + f1 * 2^(-Days/K),
          data = wtloss,
          start = list(f0 = 81, f1 = 102, K = 141),
          trace = TRUE)
## 87.49479 :   81 102 141
## 39.24471 :   81.3749 102.6829 141.9110
## 39.2447 :   81.37382 102.68411 141.91036
summary(n0)
## 
## Formula: Weight ~ f0 + f1 * 2^(-Days/K)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## f0   81.374      2.269   35.86   <2e-16 ***
## f1  102.684      2.083   49.30   <2e-16 ***
## K   141.910      5.295   26.80   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8949 on 49 degrees of freedom
## 
## Number of iterations to convergence: 2 
## Achieved convergence tolerance: 1.47e-07
# f0: df(x)/df0 = 1
# f1: df(x)/df1 = 2^(-Days/K)
# 2 parâmetros na parte analítica, 1 parâmetro na parte numérica.

n1 <- nls(Weight ~ cbind(f0 = 1, f1 = 2^(-Days/K)),
          data = wtloss,
          start = list(K = 141),
          trace = TRUE,
          algorithm = "plinear")
## 39.26875 :  141.00000  81.75839 102.34194
## 39.2447 :  141.90792  81.37485 102.68320
## 39.2447 :  141.91037  81.37381 102.68412
summary(n1)
## 
## Formula: Weight ~ cbind(f0 = 1, f1 = 2^(-Days/K))
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## K        141.910      5.295   26.80   <2e-16 ***
## .lin.f0   81.374      2.269   35.86   <2e-16 ***
## .lin.f1  102.684      2.083   49.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8949 on 49 degrees of freedom
## 
## Number of iterations to convergence: 2 
## Achieved convergence tolerance: 8.722e-07
# Perfil não disponível para plinear.
confint.default(n1)
##             2.5 %    97.5 %
## K       131.53331 152.28744
## .lin.f0  76.92664  85.82098
## .lin.f1  98.60198 106.76626
p0 +
    llayer(with(as.list(coef(n0)),
                panel.curve(f0 + f1 * 2^(-x/K),
                            add = TRUE, col = 2)))

#-----------------------------------------------------------------------
# Usando para os dados de ganho de peso de peru.

data(turk0, package = "alr4")
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 ...
p1 <- xyplot(Gain ~ A, data = turk0, col = 1,
             xlab = "Metionina",
             ylab = "Ganho de peso")
p1 +
    llayer(with(list(f0 = 620, f1 = 160, k = 9),
                panel.curve(f0 + f1 * (1 - exp(-k * x)),
                            add = TRUE, col = 2)))

# Modelo monomolecular, na verdade é uma parametrização do anterior.
n0 <- nls(Gain ~ f0 + f1 * (1 - exp(-k * A)),
          data = turk0,
          start = list(f0 = 620, f1 = 160, k = 9))
summary(n0)
## 
## Formula: Gain ~ f0 + f1 * (1 - exp(-k * A))
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## f0  622.958      5.901  105.57  < 2e-16 ***
## f1  178.252     11.636   15.32 2.74e-16 ***
## k     7.122      1.205    5.91 1.41e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.66 on 32 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 3.853e-06
# f0: df(x)/df0 = 1
# f1: df(x)/df1 = 1 - exp(-k * A)
# 2 parâmetros na parte analítica, 1 parâmetro na parte numérica.
D(expression(f0 + f1 * (1 - exp(-k * A))), "k")
## f1 * (exp(-k * A) * A)
n1 <- nls(Gain ~ cbind(1, 1 - exp(-k * A)),
          data = turk0,
          start = list(k = 9),
          trace = TRUE,
          algorithm = "plinear")
## 13227.65 :    9.0000 619.8563 169.2483
## 12393.93 :    6.822887 623.571298 180.124056
## 12367.55 :    7.143394 622.916003 178.125359
## 12367.42 :    7.120436 622.961583 178.262539
## 12367.42 :    7.12237 622.95774 178.25095
## 12367.42 :    7.122209 622.958057 178.251915
summary(n1)
## 
## Formula: Gain ~ cbind(1, 1 - exp(-k * A))
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## k        7.122      1.205    5.91 1.41e-06 ***
## .lin1  622.958      5.901  105.57  < 2e-16 ***
## .lin2  178.252     11.636   15.32 2.74e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.66 on 32 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 6.438e-06
p1 +
    llayer(with(as.list(coef(n0)),
                panel.curve(f0 + f1 * (1 - exp(-k * x)),
                            add = TRUE, col = 2)))

#-----------------------------------------------------------------------
# Modelo linear platô.

# b0: df(x)/db0 = 1
# b1: df(x)/db1 = x * (x < xb) + xb * (x >= xb)
# 2 parâmetros na parte analítica, 1 parâmetro na parte numérica.

n0 <- nls(Gain ~ b0 + b1 * A * (A < xb) + b1 * xb * (A >= xb),
          data = turk0,
          start = list(b0 = 600, b1 = 700, xb = 0.2),
          trace = TRUE)
## 60557 :  600.0 700.0   0.2
## 12736.14 :  629.6000000 713.3333333   0.2249048
## 12735.03 :  629.6000002 713.3333292   0.2244393
n1 <- nls(Gain ~ cbind(b0 = 1, b1 = A * (A < xb) + xb * (A >= xb)),
          data = turk0,
          start = list(xb = 0.2),
          trace = TRUE,
          algorithm = "plinear")
## 13618.1 :    0.2000 627.1939 787.2034
## 12742.28 :    0.222013 629.343549 720.398441
## 12735.03 :    0.2244153 629.5974606 713.4026302
## 12735.03 :    0.2244393 629.5999998 713.3333401
summary(n1)
## 
## Formula: Gain ~ cbind(b0 = 1, b1 = A * (A < xb) + xb * (A >= xb))
## 
## Parameters:
##          Estimate Std. Error t value Pr(>|t|)    
## xb        0.22444    0.01816   12.36 1.00e-13 ***
## .lin.b0 629.60000    5.55362  113.37  < 2e-16 ***
## .lin.b1 713.33334   64.38570   11.08 1.74e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.95 on 32 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 7.385e-08
p1 +
    llayer(with(as.list(coef(n0)),
                panel.curve(b0 + b1 * x * (x < xb) +
                            b1 * xb * (x >= xb),
                            add = TRUE, col = 2)))

#-----------------------------------------------------------------------
# Modelo van Genuchten para retenção de água no solo.

cra <- read.table("http://www.leg.ufpr.br/~walmes/data/cra_manejo.txt",
                  header = TRUE, sep = "\t")

# Troca tensão 0 por 0.1.
cra$tens[cra$tens == 0] <- 0.1
str(cra)
## 'data.frame':    1464 obs. of  6 variables:
##  $ condi: chr  "LVA3,5" "LVA3,5" "LVA3,5" "LVA3,5" ...
##  $ prof : int  20 20 20 20 20 20 20 20 20 20 ...
##  $ posi : chr  "EL" "EL" "EL" "EL" ...
##  $ rep  : int  1 2 3 1 2 3 1 2 3 1 ...
##  $ tens : num  0.1 0.1 0.1 1 1 1 2 2 2 4 ...
##  $ umid : num  0.61 0.54 0.57 0.56 0.52 0.56 0.54 0.51 0.52 0.48 ...
#-----------------------------------------------------------------------
# Selecionando um subconjunto para a presente demonstração.

cras <- droplevels(subset(cra, condi == "B"))

p2 <- xyplot(umid ~ tens, data = cras, col = 1,
             ylab = expression("Conteúdo de água do solo" ~ (g ~ g^{-1})),
             xlab = expression("Tensão matricial (kPa)"),
             scales = list(x = list(log = 10)))
p2

start <- list(tr = 0.3, ts = 0.6, al = 0.5, n = 1.3)
p2 +
    llayer(panel.curve(tr + (ts - tr)/(1 + (al * 10^x)^n)^(1 - 1/n),
                       add = TRUE, col = 2),
           data = start)

#-----------------------------------------------------------------------
# Modelo van Genuchten original
#   f(x) = tr + (ts - tr)/(1 + (al * tens)^n)^(1 - 1/n)
# tr: assintota inferior;
# ts: assintota superior;
# al e n: parâmetros de forma e escala;

D(expression(tr + (ts - tr)/(1 + (al * tens)^n)^(1 - 1/n)), "tr")
## 1 - 1/(1 + (al * tens)^n)^(1 - 1/n)
D(expression(tr + (ts - tr)/(1 + (al * tens)^n)^(1 - 1/n)), "ts")
## 1/(1 + (al * tens)^n)^(1 - 1/n)
#-----------------------------------------------------------------------

n0 <- nls(umid ~ tr + (ts - tr)/(1 + (al * tens)^n)^(1 - 1/n),
          data = subset(cras),
          start = start,
          trace = TRUE)
## 0.01577304 :  0.3 0.6 0.5 1.3
## 0.002424601 :  0.3021933 0.6282174 0.4502268 1.2946416
## 0.00241469 :  0.3021613 0.6282819 0.4555884 1.2952971
## 0.00241469 :  0.3021867 0.6282759 0.4554458 1.2953747
## 0.00241469 :  0.3021857 0.6282761 0.4554512 1.2953720
summary(n0)
## 
## Formula: umid ~ tr + (ts - tr)/(1 + (al * tens)^n)^(1 - 1/n)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## tr 0.302186   0.020365  14.839 2.87e-13 ***
## ts 0.628276   0.006019 104.384  < 2e-16 ***
## al 0.455451   0.089933   5.064 3.98e-05 ***
## n  1.295372   0.049501  26.169  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01025 on 23 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 4.704e-07
n1 <- nls(umid ~ cbind(tr = 1 - 1/(1 + (al * tens)^n)^(1 - 1/n),
                       ts = 1/(1 + (al * tens)^n)^(1 - 1/n)),
          data = cras,
          start = list(al = 0.5, n = 1.3),
          trace = TRUE,
          algorithm = "plinear")
## 0.002509097 :  0.5000000 1.3000000 0.3070695 0.6322750
## 0.002414817 :  0.4541042 1.2950590 0.3019482 0.6281356
## 0.00241469 :  0.4554840 1.2953564 0.3021801 0.6282775
## 0.00241469 :  0.4554498 1.2953727 0.3021859 0.6282760
summary(n1)
## 
## Formula: umid ~ cbind(tr = 1 - 1/(1 + (al * tens)^n)^(1 - 1/n), ts = 1/(1 + 
##     (al * tens)^n)^(1 - 1/n))
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## al      0.455450   0.089933   5.064 3.98e-05 ***
## n       1.295373   0.049501  26.169  < 2e-16 ***
## .lin.tr 0.302186   0.020365  14.839 2.87e-13 ***
## .lin.ts 0.628276   0.006019 104.384  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01025 on 23 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 6.964e-06
p2 +
    llayer(panel.curve(tr + (ts - tr)/(1 + (al * 10^x)^n)^(1 - 1/n),
                       add = TRUE, col = 2),
           data = as.list(coef(n0)))

#-----------------------------------------------------------------------
# Em todos os exemplos pode-se até pensar que não existe vantagem em
# usar o algortimo plinear porque o número de passos para convergência
# foi igual.  Vale ressaltar que os modelos considerados são bem
# regulares/comportados e os valores iniciais foram bem próximos das
# estimativas. O gargalo da estimação são os valores iniciais. Se os
# valores forncecidos são bons espera-se pouca contribuição do método de
# estimação.

3 Funções self-start

#-----------------------------------------------------------------------
# Em algumas situações, é possível obter valores iniciais para os
# parâmetros baseados em estatísticas da amostra. O primeiro passo é
# definir quais os melhores "pré-estimadores" para cada parâmetro. O
# segundo é escrever uma função que passe esses valores para a `nls()`
# ou outra função de estimação. Modelos que auto providênciam valores
# iniciais são chamadas de self-start.

#-----------------------------------------------------------------------
# Self-start para o modelo de perda de peso (wtloss).

data(wtloss, package = "MASS")

p0 <- xyplot(Weight ~ Days, data = wtloss, col = 1,
             xlab = "Dias em dieta",
             ylab = "Peso (kg)")

p0 +
    llayer(with(list(f0 = 81, f1 = 102, K = 141),
                panel.curve(f0 + f1 * 2^(-x/K),
                            add = TRUE, col = 2)))

#-----------------------------------------------------------------------
# Pré-estimativas (uma idéia realmente simples!).

pe <- with(wtloss, {
    f1 <- diff(range(Weight))           # Amplitude.
    f0 <- min(Weight)                   # Mínimo.
    m <- f0 + 0.5 * f1                  # Semiamplitude.
    mm <- which.min(abs(Weight - m))[1] # y mais perto da semiamplitude.
    K <- Days[mm]                       # x correspondente.
    list(f0 = f0, f1 = f1, K = K)
})

p0 +
    llayer(with(pe,
                panel.abline(v = K,
                             h = c(f0 + c(0, 0.5, 1) * f1),
                             lty = 2)))

#-----------------------------------------------------------------------
# Criando função que retorna as estimativas.

model.init <- function(x, y){
    f1 <- diff(range(y))           # Amplitude.
    f0 <- min(y)                   # Mínimo.
    m <- f0 + 0.5 * f1             # Semiamplitude.
    mm <- which.min(abs(y - m))[1] # y mais perto da semiamplitude.
    K <- x[mm]                     # x correspondente.
    list(f0 = f0, f1 = f1, K = K)
}

model.init(x = wtloss$Days, y = wtloss$Weight)
## $f0
## [1] 112.6
## 
## $f1
## [1] 71.75
## 
## $K
## [1] 88
#-----------------------------------------------------------------------
# Colocando isso dentro de uma função self-start.

model.init <- function(mCall, data, LHS){
    xy <- data.frame(sortedXyData(mCall[["input"]], LHS, data))
    xx <- xy[["x"]]
    yy <- xy[["y"]]
    f0 <- min(yy)                   # Mínimo da amostra.
    f1 <- diff(range(yy))           # Amplitude.
    m <- f0 + 0.5 * f1              # Semiamplitude.
    mm <- which.min(abs(yy - m))[1] # y mais perto da semiamplitude.
    K <- xx[mm]                     # x correspondente.
    value <- c(f0, f1, K)
    names(value) <- mCall[c("f0","f1","K")]
    value
}

# Saí o gradiente atributo gradiente junto.
SSmodel <- selfStart(model = ~f0 + f1 * 2^(-input/K),
                     parameters = c("f0", "f1", "K"),
                     initial = model.init,
                     template = function(input, f0, f1, K) NULL)
SSmodel
## function (input, f0, f1, K) 
## {
##     .expr3 <- 2^(-input/K)
##     .value <- f0 + f1 * .expr3
##     .grad <- array(0, c(length(.value), 3L), list(NULL, c("f0", 
##         "f1", "K")))
##     .grad[, "f0"] <- 1
##     .grad[, "f1"] <- .expr3
##     .grad[, "K"] <- f1 * (.expr3 * (log(2) * (input/K^2)))
##     attr(.value, "gradient") <- .grad
##     .value
## }
## attr(,"initial")
## function(mCall, data, LHS){
##     xy <- data.frame(sortedXyData(mCall[["input"]], LHS, data))
##     xx <- xy[["x"]]
##     yy <- xy[["y"]]
##     f0 <- min(yy)                   # Mínimo da amostra.
##     f1 <- diff(range(yy))           # Amplitude.
##     m <- f0 + 0.5 * f1              # Semiamplitude.
##     mm <- which.min(abs(yy - m))[1] # y mais perto da semiamplitude.
##     K <- xx[mm]                     # x correspondente.
##     value <- c(f0, f1, K)
##     names(value) <- mCall[c("f0","f1","K")]
##     value
## }
## attr(,"pnames")
## [1] "f0" "f1" "K" 
## attr(,"class")
## [1] "selfStart"
# As pré-estimativas avaliadas.
getInitial(Weight ~ SSmodel(Days, f0, f1, K),
           data = wtloss)
##     f0     f1      K 
## 112.60  71.75  88.00
#-----------------------------------------------------------------------
# Usando o modelo self-start.

n0 <- nls(Weight ~ SSmodel(Days, f0, f1, K),
          data = wtloss,
          trace = TRUE)
## 1159.454 :  112.60  71.75  88.00
## 100.393 :   89.68663  94.11074 128.73510
## 39.59629 :   81.85524 102.19312 141.26936
## 39.2447 :   81.37408 102.68373 141.91136
## 39.2447 :   81.37382 102.68411 141.91036
summary(n0)
## 
## Formula: Weight ~ SSmodel(Days, f0, f1, K)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## f0   81.374      2.269   35.86   <2e-16 ***
## f1  102.684      2.083   49.30   <2e-16 ***
## K   141.910      5.295   26.80   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8949 on 49 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 9.199e-08
#-----------------------------------------------------------------------
# Modelo quadrático com ponto crítico (QCC) desenvolvido na minha Tese.
# Taxa de crescimento do disco de micélio para de 16 isolados de fungo
# em função da temperatura.
# Dados de Paulo Lichtemberg (Doutorando em Produção Vegetal UFPR).

cm <- read.table("http://www.leg.ufpr.br/~walmes/data/cresmicelial.txt",
                 header = TRUE,
                 sep = "\t",
                 colClasses = c("factor", "numeric", "numeric"))
str(cm)
## 'data.frame':    192 obs. of  3 variables:
##  $ isolado: Factor w/ 16 levels "1","10","12",..: 1 1 1 8 8 8 11 11 11 12 ...
##  $ temp   : num  20 20 20 20 20 20 20 20 20 20 ...
##  $ mmdia  : num  14.25 13.25 13.25 9.75 12.75 ...
#-----------------------------------------------------------------------
# Ver os dados.

p1 <- xyplot(mmdia ~ temp | isolado,
             data = cm,
             type = c("p", "a"),
             col = 1,
             ylab = "Taxa de crescimento em diâmetro do disco micelial",
             xlab = "Temperatuda da câmara de crescimento")
print(p1)

#-----------------------------------------------------------------------
# Ajustar para as 16 situações requer 16 conjutos de 3 valores iniciais
# (48 no total). Esse modelo tem interpretação cartesiana dos parâmetros,
# é fácil conseguí-los. Dada relação entre os parâmetros nas duas
# parametrizações do quadrático, podemos ajustar com a lm() e passar os
# valores para nls(). Vamos incorporar isso na self-start.

#-----------------------------------------------------------------------
# Modelo quadrático original (MQO)
#                               f(x) = f0+tx0*x+cv*x^2
# f0: valor da função na origem
# tx0: taxa na origem
# cv: curvatura
#
# Modelo quadrático centrado no ponto crítico (MQC)
#                               f(x) = fc+cv*(x-xc)^2
# fc: valor da função no ponto crítico (mínimo ou máximo)
# xc: ponto crítico
# cv: curvatura
#
# xc é calculado por xc=-tx0/(2*cv);
# fc é clauclado por fc=f0+tx0*xc+cv*xc^2;

#-----------------------------------------------------------------------
# Fazendo a self-start para o modelo quadrático.

mqc.init <- function(mCall, data, LHS) {
    xy <- data.frame(sortedXyData(mCall[["input"]], LHS, data))
    m0 <- lm(y ~ poly(x, 2, raw = TRUE), data = xy)
    c0 <- coef(m0)
    cv <- c0[3]
    xc <- -0.5 * c0[2]/cv
    fc <- c0[1] + c0[2] * xc + cv * xc
    value <- c(fc, xc, cv)
    names(value) <- mCall[c("fc", "xc", "cv")]
    value
}

SSmqc <- selfStart(model = ~fc + cv * (input - xc)^2,
                   parameters = c("fc", "xc", "cv"),
                   initial = mqc.init,
                   template = function(input, fc, xc, cv) NULL)

# São mais que pré-estimativas, são as próprias estimativas!
getInitial(mmdia ~ SSmqc(temp, fc, xc, cv),
           data = subset(cm, isolado == "1"))
##          fc          xc          cv 
## 144.8294997  22.6213375  -0.2639147
#-----------------------------------------------------------------------
# Ajustando o modelo para todos os isolados, usar lapply.

cms <- split(cm, f = cm$isolado)

# São as estimativas.
sapply(cms,
       FUN = getInitial,
       object = mmdia ~ SSmqc(temp, fc, xc, cv))
##              1         10          12          13          15          18          19
## fc 144.8294997 98.4113423 103.4151912 193.5447032 36.29401324 141.2557967 128.8486643
## xc  22.6213375 22.3112078  22.7514927  23.8436386 22.81546454  23.1069169  22.4207221
## cv  -0.2639147 -0.1829921  -0.1834489  -0.3264447 -0.06088136  -0.2497444  -0.2396385
##              2         24          25           4           5         6           7
## fc 222.8228983 71.0619101 50.68873740 167.7179993 169.7166482 217.37310 147.3360725
## xc  23.4325194 22.5735045 21.28915437  24.3359842  23.0839096  23.41014  22.5886086
## cv  -0.3908379 -0.1290247 -0.09392596  -0.2720432  -0.3013447  -0.38257  -0.2698738
##             8           9
## fc 92.4982155 159.9247939
## xc 21.8205085  23.3071114
## cv -0.1769139  -0.2771935
# Ajuste em lote.
n0 <- lapply(cms,
             FUN = nls,
             formula = mmdia ~ SSmqc(temp, fc, xc, cv))

c0 <- t(sapply(n0, coef))
c0
##           fc       xc          cv
## 1  15.747872 22.62134 -0.26391471
## 10 11.402460 22.31121 -0.18299214
## 12 12.630189 22.75149 -0.18344892
## 13 15.738307 23.84364 -0.32644466
## 15  5.991539 22.81546 -0.06088136
## 18 13.680697 23.10692 -0.24974439
## 19 13.757928 22.42072 -0.23963854
## 2  17.378770 23.43252 -0.39083792
## 24  8.228228 22.57350 -0.12902469
## 25 10.118459 21.28915 -0.09392596
## 4  13.223529 24.33598 -0.27204322
## 5  16.096225 23.08391 -0.30134475
## 6  16.667437 23.41014 -0.38257005
## 7  15.730306 22.58861 -0.26987384
## 8  12.123749 21.82051 -0.17691388
## 9  15.807922 23.30711 -0.27719350
ic0 <- lapply(n0, confint.default)
ic0
## $`1`
##         2.5 %    97.5 %
## fc 14.4326624 17.063082
## xc 22.1462456 23.096429
## cv -0.3211214 -0.206708
## 
## $`10`
##         2.5 %     97.5 %
## fc 10.9737562 11.8311634
## xc 22.0551298 22.5672857
## cv -0.2022811 -0.1637032
## 
## $`12`
##         2.5 %     97.5 %
## fc 10.7184582 14.5419191
## xc 21.8126286 23.6903568
## cv -0.2655552 -0.1013426
## 
## $`13`
##         2.5 %     97.5 %
## fc 14.4660573 17.0105557
## xc 23.6042484 24.0830288
## cv -0.3774317 -0.2754576
## 
## $`15`
##         2.5 %       97.5 %
## fc  4.7740561  7.209021508
## xc 21.0626370 24.568292125
## cv -0.1128637 -0.008898969
## 
## $`18`
##         2.5 %     97.5 %
## fc 11.1510401 16.2103547
## xc 22.3202574 23.8935763
## cv -0.3551974 -0.1442914
## 
## $`19`
##         2.5 %     97.5 %
## fc 11.7527082 15.7631470
## xc 21.5493807 23.2920636
## cv -0.3287375 -0.1505396
## 
## $`2`
##         2.5 %     97.5 %
## fc 15.9461134 18.8114258
## xc 23.1803316 23.6847072
## cv -0.4493119 -0.3323639
## 
## $`24`
##         2.5 %      97.5 %
## fc  5.9016375 10.55481760
## xc 20.8180431 24.32896597
## cv -0.2307185 -0.02733085
## 
## $`25`
##         2.5 %      97.5 %
## fc  7.8178177 12.41910097
## xc 17.1062047 25.47210407
## cv -0.2130253  0.02517337
## 
## $`4`
##       2.5 %     97.5 %
## fc 11.49458 14.9524756
## xc 23.96148 24.7104893
## cv -0.34066 -0.2034265
## 
## $`5`
##         2.5 %    97.5 %
## fc 14.0025901 18.189861
## xc 22.5393259 23.628493
## cv -0.3887714 -0.213918
## 
## $`6`
##         2.5 %     97.5 %
## fc 12.8654097 20.4694640
## xc 22.7211688 24.0991153
## cv -0.5379462 -0.2271939
## 
## $`7`
##         2.5 %     97.5 %
## fc 13.8599378 17.6006746
## xc 21.9183626 23.2588546
## cv -0.3514993 -0.1882483
## 
## $`8`
##         2.5 %      97.5 %
## fc  8.8379606 15.40953712
## xc 19.2978930 24.34312409
## cv -0.3343823 -0.01944545
## 
## $`9`
##         2.5 %     97.5 %
## fc 14.3175646 17.2982790
## xc 22.9204540 23.6937689
## cv -0.3384766 -0.2159104
#-----------------------------------------------------------------------
# Análise dos resíduos em série.

res <- do.call(c, lapply(n0, residuals))
qqmath(~res | cm$isolado)

#-----------------------------------------------------------------------
# Colocando os resultados no gráfico.

# Formato de lista com vetores nomeados é mais adequado agora.
c0 <- lapply(n0, coef)

xyplot(mmdia ~ temp | isolado,
       data = cm,
       col = 1,
       ylab = "Taxa de crescimento em diâmetro do disco micelial",
       xlab = "Temperatuda da câmara de crescimento",
       panel = function(...) {
           panel.xyplot(...)
           wp <- which.packet()
           with(as.list(c0[[wp]]),
                panel.curve(fc + cv * (x - xc)^2,
                            col = 2, lwd = 2))
           panel.abline(v = c0[[wp]]["xc"], lty = 2, col = 2)
           panel.abline(v = ic0[[wp]]["xc", ], lty = 3, col = 2)
       })

#-----------------------------------------------------------------------
# Para essas situações é mais indicado usar a `nls::nlsList()`, pois ela
# trata melhor as saídas colocando todos os ajustes no mesmo output.
# Outra opção de análise é tratar o efeito de isolado como aleatório
# (pois são uma amostra) e verificar os componentes de variância em
# relação à cada parâmetro. Para inferir sobre a temperatura de ótimo
# pode-se usar o modelo quadrático original e obter IC por meio do
# método delta. Eles serão exatamente iguais aos IC assintóticos feitos
# a partir do modelo quadrático centrado.

#-----------------------------------------------------------------------
# Inferência para xc via método delta.

m0 <- lm(mmdia ~ temp + I(temp^2), cms[[1]])
n0 <- nls(mmdia ~ SSmqc(temp, fc, xc, cv), cms[[1]])

# mesma estimativa e erro padrão, logo mesmo IC
car::deltaMethod(m0, "-0.5 * b1/b2",
                 parameterNames = c("b0", "b1", "b2"))
##              Estimate      SE   2.5 % 97.5 %
## -0.5 * b1/b2  22.6213  0.2424 22.1462 23.096
summary(n0)$coeff["xc", 1:2]
##   Estimate Std. Error 
## 22.6213375  0.2423983
#-----------------------------------------------------------------------
# Mais funções self-start você encontra procurando aquelas que começam
# com SS.

apropos("^SS", where = TRUE)
##            33            33            33            33            33            33 
##     "SSasymp"  "SSasympOff" "SSasympOrig"     "SSbiexp"         "SSD"       "SSfol" 
##            33            33            33            33             1             1 
##       "SSfpl"  "SSgompertz"     "SSlogis"    "SSmicmen"     "SSmodel"       "SSmqc" 
##            33 
##   "SSweibull"
# Mais delas no pacote drc, HydroMe e outros.

4 Busca em grid

#-----------------------------------------------------------------------
# Não raramente tem-se modelos que não apresentam interpretação
# cartesiana para os parâmetros dificultando assim o usuário obter
# valores iniciais pela inspeção de gráficos ou medidas descritivas. Em
# geral, gasta-se muito tempo testando valores iniciais e uma forma de
# evitar isso é criar um grid de valores e avaliar a função para
# comninações dos parâmetros. Aquela que apresentar menor soma de
# quadrados pode ser usada como valor inicial de onde espera-se sucesso
# na convergência.  Para essa tarefa usa-se funções disponíveis no
# pacote `nls2`, em especial a função `nls()` que tem o algortimo
# "grid-search" ou "brute-force".

#-----------------------------------------------------------------------
# Curvas de lactação.

milk <- read.table("http://www.leg.ufpr.br/~walmes/data/saxton_lactacao1.txt",
                   header = TRUE, sep = "\t", colClasses = c("factor", NA, NA))
str(milk)
## 'data.frame':    29 obs. of  3 variables:
##  $ vaca : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 2 ...
##  $ leite: num  27.4 30.5 29.9 30.2 29.1 28.6 27.4 25.4 22.8 35.8 ...
##  $ dia  : int  23 54 83 114 146 175 205 266 297 17 ...
xyplot(leite ~ dia, groups = vaca, data = milk, type = "o")

milks <- subset(milk, vaca == "3")
p4 <- xyplot(leite ~ dia,
             groups = vaca,
             data = milks,
             col = 1,
             type = "o")
p4

#-----------------------------------------------------------------------
# Ajustar o modelo de Wood f(x) = A * x^B * exp(-C * x), como obter
# valores inicais para esses parâmetros? O modelo de Wood usa o núcleo
# da distribuição Gama.

# Grid de valores iniciais (as amplitudes vem do conhecimento prévio)
vi.grid <- expand.grid(A = seq(5, 30, l = 10),
                       B = seq(0, 1, l = 10),
                       C = seq(0, 0.1, l = 50))
nrow(vi.grid)
## [1] 5000
# Testa toda tupla de parâmetros e fica com a melhor.
n0 <- nls2::nls2(leite ~ A * dia^B * exp(-C * dia),
                 data = milks,
                 start = vi.grid,
                 algorithm = "grid-search")
coef(n0)
##            A            B            C 
## 21.666666667  0.222222222  0.004081633
# Usa a melhor tupla como valores iniciais.
n1 <- nls(formula(n0), data = milks, start = coef(n0))
summary(n1)
## 
## Formula: leite ~ A * dia^B * exp(-C * dia)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## A 2.023e+01  1.465e+00   13.81 2.46e-06 ***
## B 2.529e-01  2.135e-02   11.85 6.93e-06 ***
## C 4.803e-03  2.252e-04   21.33 1.26e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7134 on 7 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 1.984e-06
p4 +
    llayer(with(as.list(coef(n1)),
                panel.curve(A * x^B * exp(-C * x),
                            add = TRUE, col = 2)))

#-----------------------------------------------------------------------
# No caso de modelos onde é dificil montar uma função self-start pode-se
# usar a `nls2()` com o "grid-search" dentro e os valores inicais seriam
# aqueles selecionados.

5 Ajuste com interface gráfica

library(rpanel)

#-----------------------------------------------------------------------
# Exemplos da documentação.

# library(wzRfun)
# help(rp.nls, h = "html")

# Se não possuir o pacote `wzRfun`, carregue a função da fonte.
source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/rp.nls.R")

#-----------------------------------------------------------------------
# Ganho de peso de perus.

library(rpanel)

data(turk0, package = "alr4")
str(turk0)

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)

# Guarde isso para reusar.
dput(as.list(round(coef(turk_fit), 4)))
# list(Int = 622.2417, Asy = 857.4585, Half = 0.1548)

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))))
})

#-----------------------------------------------------------------------
# 1. 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.
# ltens  : representado por ltens (log da tensão matricial, psi).
# umid   : representado por umid, conteúdo de água do solo ( % ).
# Us     : assíntota inferior, mínimo da função, quando x -> +infinito.
# Ur     : assíntota superior, máximo da função, quando x -> -infinito.
# al     : locação, translada o ponto de inflexão.
# n      : 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")

sapply(cra_fit, coef)

#-----------------------------------------------------------------------
# 2. 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)

xyplot(pcapu ~ desf | estag, data = cap, layout = c(5, 1),
       xlab = "Nível de desfolha artificial", ylab = "Peso de capulhos")

# modelo: potência.
# desf  : representado por desf (nível de desfolha artifical).
# pcapu : representado por pcapu (peso de capulhos), produto do algodão.
# f0    : intercepto, valor da função quando x=0 (situação ótima).
# delta : diferença no valor da função para x=0 e x=1 (situação
#         extrema).
# curv  : 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")

length(cap_fit)
sapply(cap_fit, coef)
lapply(cap_fit, summary)

#-----------------------------------------------------------------------
# 3. Curva de produção em função do nível de potássio no solo.

soja <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
                   header = TRUE, sep = "\t", encoding = "latin1",
                   dec = ",")

soja$agua <- factor(soja$agua)
str(soja)

xyplot(rengrao ~ potassio|agua, data=soja)

# modelo: linear-plato.
# x: representado por potássio, conteúdo de potássio do solo.
# y: representado por rengrao, redimento de grãos por parcela.
# f0: intercepto, valor da função quando x=0.
# tx: taxa de incremento no rendimento por unidade de x.
# brk: valor acima do qual a função é constante.

model <- rengrao ~ f0 + tx * potassio * (potassio < brk) +
    tx * brk * (potassio >= brk)
start <- list(f0 = c(15, 5, 25),
              tx = c(0.2, 0, 1),
              brk = c(50, 0, 180))

pot_fit <- rp.nls(model = model, data = soja, start = start,
                  subset = "agua")

sapply(pot_fit, coef)

#-----------------------------------------------------------------------
# 4. Curva de lactação.

lac <- read.table("http://www.leg.ufpr.br/~walmes/data/saxton_lactacao1.txt",
                  header = TRUE, sep = "\t", encoding = "latin1")

lac$vaca <- factor(lac$vaca)
str(lac)

xyplot(leite ~ dia | vaca, data = lac)

# modelo: de Wood (nucleo da densidade gama).
# x: representado por dia, dia após parto.
# y: representado por leite, quantidade produzida.
# th1: escala, desprovido de interpretação direta.
# th2: forma, desprovido de interpretação direta.
# th3: forma, desprovido de interpretação direta.

model <- leite ~ th1 * dia^th2 * exp(-th3 * dia)
start <- list(th1 = c(15, 10, 20),
              th2 = c(0.2, 0.05, 0.5),
              th3 = c(0.0025, 0.001, 0.008))

lac_fit <- rp.nls(model = model,
                  data = lac,
                  start = start,
                  subset = "vaca",
                  xlim = c(0, 310))

# dput(sapply(lac_fit, coef))

# structure(c(15.753532986063, 0.194320322994364, 0.00237895534802468,
# 20.6944991160729, 0.220198842519765, 0.00391704656088934, 20.2268920495124,
# 0.252907554946877, 0.00480308843437403), .Dim = c(3L, 3L), .Dimnames = list(
#     c("th1", "th2", "th3"), c("1", "2", "3")))

# save(lac_fit, "ajustes.RData")
# load("ajustes.RData")

#-----------------------------------------------------------------------
# 5. Curva de secagem do solo em microondas.

sec <- read.table("http://www.leg.ufpr.br/~walmes/data/emr11.txt",
                  header = TRUE, sep = "\t", encoding = "latin1")
sec$nome <- factor(sec$nome)
str(sec)

xyplot(umrel ~ tempo | nome, data = sec)

# modelo: logístico.
# x: representado por tempo, período da amostra dentro do microondas.
# y: representado por umrel, umidade relativa o conteúdo total de água.
# th1: assíntota superior.
# th2: tempo para evaporar metade do conteúdo total de água.
# th3: proporcional à taxa máxima do processo.

model <- umrel ~ th1/(1 + exp(-(tempo - th2)/th3))
start <- list(th1 = c(1, 0.8, 1.2),
              th2 = c(15, 0, 40),
              th3 = c(8, 2, 14))

sec_fit <- rp.nls(model = model,
                  data = sec,
                  start = start,
                  subset = "nome")

sapply(sec_fit, coef)
lapply(sec_fit, summary)