Arquivo

Posts Tagged ‘nls’

Ilustração do método Newton-Raphson em regressão não linear

Trajetória obtida com os valores consecutivos nas iterações do método Newton-Raphson para estimação de parâmetros em um modelo de regressão não linear.

Em regressão não linear, as estimativas dos parâmetros são obtidas por emprego de algum método numérico. Um dos métodos empregados é o de Newton-Raphson. Tal método é usado para encontrar raízes de uma função. Portanto, não é um método de maximação/mínimização de funções, mas usa o fato equivalente de que a derivada no ponto crítico de uma função é zero.

Documentação sobre o método Newton-Raphson existe em abundância. O que motivou fazer esse post foi não ter encontrado nenhum que exemplificasse aplicação para estimação de parâmetros em um modelo de regressão não linear. O estimador de mínimos quadrados para um modelo de regressão é aquele que torna mínima a função

S(\theta) = \sum_{i=1}^n (y_i-f(x_i, \theta))^2.

O método Newton-Raphson vai procurar as raízes da derivada dessa função com relação à cada parâmetro. Então, se \theta é um vetor de p parâmetros, a derivada de S(\theta) será uma função p-dimensional,

S_1(\theta) = \displaystyle\frac{\partial S(\theta)}{\partial\theta}=\begin{bmatrix}  \displaystyle\frac{\partial S(\theta)}{\partial\theta_1}\\  \vdots\\  \displaystyle\frac{\partial S(\theta)}{\partial\theta_p}  \end{bmatrix}.

O método se baseia na aproximação de primeira ordem da função S_1(\theta) de forma que um valor da função para um valor \theta^{i+1} na redondeza de um valor \theta^i é

S_1(\theta^{i+1}) \approx S_1(\theta^i)+ (\theta^{i+1}-\theta^i)\displaystyle\frac{\partial S_1(\theta)}{\partial\theta}|_{\theta=\theta^i}.

O valor de S_1(\theta^{i+1}) é zero uma vez que \theta^{i+1} é a raíz da função considerando a expansão de primeira ordem, assim, reorganizando a expressão, temos que a atualização dos valores é data por

\theta^{i+1} = \theta^i- \displaystyle\frac{S_1(\theta^i)}{S_2(\theta^i)},

sendo que

S_2(\theta) = \frac{\partial S_1(\theta)}{\partial\theta}=\begin{bmatrix}  \frac{\partial S_1(\theta)}{\partial\theta_1\partial\theta_1} & \ddots\\  \ddots& \frac{\partial S_1(\theta)}{\partial\theta_p\partial\theta_p}  \end{bmatrix}

é a matriz de derivadas que tem dimensão p\times p.

Esse procedimento numérico foi implementado no script abaixo para estimação dos parâmetros do modelo Michaelis-Menten e o modelo Exponencial reparametrizado que desenvolvi na minha Dissertação de Mestrado. No início eu monto as funções separadas e no final eu faço uma função que recebe os ingredientes mínimos e faz a estimação dos parâmetros. No final do script tem um procedimento gráfico para seleção dos chutes iniciais e acompanhar a trajetória durante as iterações. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# método Newton-Raphson
# objetivo é mínimizar a soma de quadrados dos resíduos mas iremos
# encontrar as raízes da derivada dessa função pelo método de Newton-Raphson

#-----------------------------------------------------------------------------
# parâmetros, modelo e dados simulados deles

theta <- c(A=10, B=2) # parâmetros
x <- 1:12
f <- function(theta, x){
  ## função que retorna os valore preditos para theta e x fornecidos
  ## theta: vetor de parâmetros A e B do modelo Michaelis-Menten
  ## x: vetor de valores da variável independente
  haty <- theta[1]*x/(theta[2]+x)
  return(haty)
}
set.seed(222); e <- rnorm(x,0,0.1)
y <- f(theta=theta, x=x)+e
plot(y~x)
curve(f(theta, x), add=TRUE, col=2)

#-----------------------------------------------------------------------------
# função objetivo, derivada primeira e segunda com relação à theta

S0 <- function(theta, f, y, x){
  ## função que retorna a soma de quadrados dos resíduos para
  ## theta: vetor de parâmetros do modelo f
  ## f: função do modelo de regressão não linear que depende de theta e x
  ## y e x: vetores dos valores das variáveis dependente e independente
  sqr <- sum((y-f(theta, x))^2)
  return(sqr)
}

S0(theta=theta, f=f, y=y, x=x)

# derivada primeira, dimensão 1 x p
D(expression((y-A*x/(B+x))^2), "A") # derivada de S0 com relação a A
D(expression((y-A*x/(B+x))^2), "B") # derivada de S0 com relação a B

# derivada segunda, dimensão p x p
D(D(expression((y-A*x/(B+x))^2), "A"), "A") # derivada de S0 com relação a A e A
D(D(expression((y-A*x/(B+x))^2), "B"), "B") # derivada de S0 com relação a B e B
D(D(expression((y-A*x/(B+x))^2), "A"), "B") # derivada de S0 com relação a A e B
D(D(expression((y-A*x/(B+x))^2), "B"), "A") # derivada de S0 com relação a B e A

#-----------------------------------------------------------------------------
# para ter essas funções vamos usar a deriv3()

ff <- deriv3(~(y-A*x/(B+x))^2,
             c("A","B"),
             function(x, y, A, B){ NULL })

ff(x=x, y=y, theta[1], theta[2])

#-----------------------------------------------------------------------------
# mas as matrizes são 1 x p e p x p, então falta fazer a soma em n

fff <- function(x, y, A, B){
  ## função que retorna gradiente e hessiano avaliado em theta
  ## y e x: vetores dos valores das variáveis dependente e independente
  ## A e B: valores dos parâmetros do modelo Michaelis-Menten
  m <- ff(x, y, A, B)
  g <- attr(m, "gradient"); G <- apply(g, 2, sum)
  h <- attr(m, "hessian"); h <- matrix(h, ncol=2*2)
  H <- apply(h, 2, sum); H <- matrix(H, 2, 2)
  return(list(G=G, H=H))
}

fff(x=x, y=y, theta[1], theta[2])

#-----------------------------------------------------------------------------
# agora é só fazer o Newton-Raphson

th0 <- c(9,2); tol <- 1e-10; niter <- 100;
dis <- 1; i <- 0
while(dis>tol & i<=niter){
  GH <- fff(x=x, y=y, th0[1], th0[2])
  G <- GH[[1]]
  iH <- solve(GH[[2]])
  th1 <- th0-iH%*%G
  dis <- sum((th0-th1)^2)
  i <- i+1
  print(c(i, dis, th1))
  th0 <- th1
}

n0 <- nls(y~A*x/(B+x), start=c(A=9, B=2))
coef(n0)

iH/vcov(n0) # hessiano é proporcional à matriz de covariância dos parâmetros

#-----------------------------------------------------------------------------
# vamos incluir todos os elementos em uma única função!

newton <- function(formula, data, start, tol=1e-10, niter=100, trace=TRUE){
  ## formula: expressão da resposta como função das covariáveis e parâmetros
  ## data: data.frame que contém dos dados
  ## start: vetor ou lista nomeados de valores iniciais para os parâmetros
  ## tol: norma mínima entre valores consecutivos dos parâmetros
  ## niter: número máximo de iterações do algortimo
  ## trace: se TRUE imprime o histórico de iterações
  form <- paste("~(",paste(formula[2:3], collapse="-"),")^2", sep="")
  vars <- all.vars(formula)
  nvars <- length(vars)
  pars <- 1:(nvars+1)
  l <- as.list(pars)
  names(l)[1:nvars] <- vars
  f <- as.function(l)
  body(f) <- NULL
  ff <- deriv3(as.formula(form), names(start), f)
  th0 <- unlist(start)
  i <- 0; dis <- 100; iter <<- matrix(th0, ncol=length(start))
  while(dis>tol & i<=niter){
    m <- do.call(ff, c(as.list(data), as.list(th0)))
    g <- attr(m, "gradient"); G <- apply(g, 2, sum)
    h <- attr(m, "hessian"); h <- matrix(h, ncol=length(start)^2)
    H <- apply(h, 2, sum); H <- matrix(H, nrow=length(start))
    iH <- solve(H)
    th1 <- th0-iH%*%G
    dis <- sum((th0-th1)^2)
    i <- i+1
    iter <<- rbind(iter, t(th1))
    if(trace==TRUE) print(c(i, dis, th1))
    th0 <- th1
  }
  th0 <- as.vector(th0)
  names(th0) <- names(start)
  ## retorna o último valor iterado e a uma matriz com o histórico
  return(list(theta=th0, iter=iter))
}

#-----------------------------------------------------------------------------
# vamos estudar vetor de chutes presentes nessa lista

start.l <- list(c(A=10,B=2.5), c(A=10,B=0.1), c(A=14,B=2.5),
                c(A=15.13, B=1.7), c(A=18.4,B=0.13))

data <- data.frame(x=x, y=y)
formula <- y~A*x/(B+x)
start <- start.l[[4]]
nr <- newton(formula, data, start)
nr$theta

sqe <- function(A, B, y, x){ hy <- (A*x)/(B+x); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))

A.grid <- seq(-3,15,l=100)
B.grid <- seq(-3,15,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
#png("f023.png", w=500, h=500)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
        xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
  arrows(iter[i,1], iter[i,2],
         iter[i+1,1], iter[i+1,2],
         col=2, length=0.1)
  Sys.sleep(0.1)
}
#dev.off()

#-----------------------------------------------------------------------------
# com o mouse selecione regiões da superfície para serem valores iniciais

loc <- locator()
loc <- as.data.frame(loc)
names(loc) <- c("A","B")
start.l <- split(loc, f=1:nrow(loc))

#-----------------------------------------------------------------------------
# ver como fica com o modelo Exponencial

n0 <- nls(y~A*(1-exp(-log(2)*x/B)), start=c(A=9, B=2))
coef(n0)

data <- data.frame(x=x, y=y)
formula <- y~A*(1-exp(-log(2)*x/B))
start <- list(A=5, B=0.5)
nr <- newton(formula, data, start)
nr$theta

sqe <- function(A, B, y, x){ hy <- A*(1-exp(-log(2)*x/B)); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))

A.grid <- seq(-10,20,l=100)
B.grid <- seq(-10,20,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
        xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
  arrows(iter[i,1], iter[i,2],
         iter[i+1,1], iter[i+1,2],
         col=2, length=0.1)
  Sys.sleep(0.1)
}
#-----------------------------------------------------------------------------

Anúncios

Método gráfico para valores iniciais em regressão não linear com RStudio

Obtenção de valores iniciais para o modelo duplo van Genuchten de 7 parâmetros com os recursos da função manipulate do RStudio.

No post Método gráfico interativo para valores iniciais em regressão não linear eu apresentei alguns procedimentos gráficos para obter bons valores iniciais para ajuste de modelos de regressão não linear. Ao obter valores iniciais pelo procedimento gráfico facilita-se a convergência do método de estimação e você passa a entender/interpretar melhor os parâmetros do modelo.

Nesse post vou apresentar os métodos gráficos construídos com os recursos do pacote manipulate que é parte do RStudio, um recente editor de scripts R que tem ganhado muitos adeptos por ser muito amigável. Eu já postei o uso dessas funções em RStudio e a função manipulate() onde falo sobre as opções das função.

Vou apresentar o procedimento com dados reais. Preciso encontrar chutes iniciais para o modelo duplo van Genuchten (ver artigo original). Esse modelo tem 7 parâmetros na seguinte equação

\theta(\Psi) = \theta_{res}+\displaystyle \frac{\theta_{pmp}-\theta_{res}}{(1+(\alpha_{tex}\Psi)^{n_{tex}})^{1-1/n_{tex}}}+\displaystyle \frac{\theta_{sat}-\theta_{pmp}}{(1+(\alpha_{est}\Psi)^{n_{est}})^{1-1/n_{est}}},

em que \theta e \Psi são valores observados e o restante são os parâmetros da função. Em geral, quanto maior o número de parâmetros a serem estimados num modelo de regressão mais demorara será a convergência porque o espaço da solução é de alta dimensão. Além do mais, com muitos parâmetros a serem estimados, pode haver fraca de identificabilidade do modelo. Dessa forma, bons chutes iniciais podem ajudar na convergência. Uma ligeira noção do campo de variação de cada parâmetro o usuário precisa ter (limites inferior e superior). Você pode inicialmente usar intervalos amplos e ir reduzindo o comprimento do intervalo a medida que observar quais são os valores mais próximos dos ótimos.

#-----------------------------------------------------------------------------
# dados para a Curva de Retenção de Água, umidade (theta) e tensão (psi)

cra <- data.frame(psi=c(0.01,1,2,4,6,8,10,33,60,100,500,1500,2787,4727,6840,
                    7863,9030,10000,10833,13070,17360,21960,26780,44860,
                    69873,74623,87287,104757,113817,147567,162723,245310,
                    262217,298223),
                  theta=c(0.779,0.554,0.468,0.406,0.373,0.36,0.344,0.309,
                    0.298,0.292,0.254,0.241,0.236,0.233,0.223,0.202,0.172,
                    0.187,0.138,0.098,0.07,0.058,0.052,0.036,0.029,0.0213,
                    0.0178,0.0174,0.0169,0.0137,0.0126,0.0109,0.0106,0.0053))

#-----------------------------------------------------------------------------
# gráfico dos valores observados

plot(theta~log10(psi), data=cra)

#-----------------------------------------------------------------------------
# função do modelo duplo van Genuchten, 7 parâmetros

dvg <- function(x, ts, ti, tr, ae, at, ne, nt){
  tr+(ti-tr)/((1+(at*10^x)^nt)^(1-1/nt))+(ts-ti)/((1+(ae*10^x)^ne)^(1-1/ne))
}

dvg(log10(c(1,10,100,1000,10000)),         # log 10 das tensões
    0.7, 0.2, 0.05, 1.3, 0.0001, 1.5, 3.5) # valor dos parâmetros

#-----------------------------------------------------------------------------
# procedimento gráfico para obter chutes e conhecer a função dos parâmetros

require(manipulate) # carrega pacote que permite manipulação gráfica
start <- list()     # cria uma lista vazia para receber os valores finais

manipulate({
             plot(theta~log10(psi), data=cra)
             curve(dvg(x, ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt), add=TRUE)
             start <<- list(ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt)
           },
           ts=slider(0.7, 0.9, initial=0.8),
           ti=slider(0.15, 0.25, initial=0.2),
           tr=slider(0, 0.10, initial=0.05),
           ae=slider(1.01, 3, initial=1.3),
           at=slider(0, 0.0001, initial=0.00005),
           ne=slider(1.01, 3, initial=1.65),
           nt=slider(1.8, 5, initial=4.3))

#-----------------------------------------------------------------------------

start # valores salvos do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo aos dados com os chutes do procedimento gráfico

n0 <- nls(theta~tr+(ti-tr)/((1+(at*psi)^nt)^(1-1/nt))+
          (ts-ti)/((1+(ae*psi)^ne)^(1-1/ne)),
          data=cra, start=start, trace=TRUE)
summary(n0) # quadro de estimativas

#-----------------------------------------------------------------------------
# gráfico dos dados com a curva estimada

lis <- c(list(x=NULL), as.list(coef(n0)), body(dvg))
plot(theta~log10(psi), cra,
     ylab=expression(Conteúdo~de~água~no~solo~(theta*","~g~g^{-1})),
     xlab=expression(Tensão~matricial~(Psi*","~kPa)), xaxt="n")
tmp <- as.function(lis)
curve(tmp, add=TRUE, col=2, lwd=1.5)
axis(1, at=-2:5, label=as.character(10^(-2:5)), lwd.ticks=2)
s <- log10(sort(sapply(1:9, function(x) x*10^(-3:6))))
axis(1, at=s, label=NA, lwd=0.5)
abline(v=-2:6, h=seq(0,1,0.05), lty=3, col="gray50")

#-----------------------------------------------------------------------------

O procedimento permite encontrar chutes iniciais para duas ou mais curvas. Basta escrever adequadamente as funções dentro da manipulate(). Abaixo estão os códigos para obter os valores inicias para o modelo logístico usado no contexto de epidemiologia. Em duas variedades de uma espécie vegetal foram observadas a severidade da doença ao longo do tempo. O objetivo é encontrar estimativas dos parâmetros para a função (que no caso é a logística) que associa severidade ao tempo e comparar as funções entre as variedades. A comparação das funções pode ser feita fazendo-se testes de hipóteses para funções dos parâmetros (em geral contrastes) ou testes envolvendo modelos aninhados, como você verá a seguir.

#-----------------------------------------------------------------------------
# dados de severidade em função do tempo para duas variedades

da <- read.table("http://www.leg.ufpr.br/~walmes/ridiculas/epidem.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# gráficos exploratórios

require(lattice)
xyplot(V1+V2~DAI, data=da, type=c("p","a"))

# ajustar modelos. separado ou junto?
# qual modelo? quais os chutes?

#-----------------------------------------------------------------------------
# mudando a disposição dos dados

require(reshape)
db <- melt(da, id=c("DAI","rep"))
names(db)[3:4] <- c("vari","sever")
str(db)

#-----------------------------------------------------------------------------
# carrega pacote e faz a função para encontrar via deslizadores valores ótimos

require(manipulate)
start <- list()

manipulate({
  plot(sever~DAI, db, col=ifelse(db$vari=="V1",1,2))
  A1 <- a1; S1 <- s1; X01 <- x01
  A2 <- a2; S2 <- s2; X02 <- x02
  curve(A1/(1+exp(-(x-X01)/S1)), add=TRUE, col=1)
  curve(A2/(1+exp(-(x-X02)/S2)), add=TRUE, col=2)
  start <<- list(A=c(A1,A2), S=c(S1,S2), X0=c(X01,X02))
  },
  a1=slider(90, 110, initial=95, step=1),
  a2=slider(90, 110, initial=95, step=1),
  s1=slider(2, 8, initial=6, step=0.1),
  s2=slider(2, 8, initial=6, step=0.1),
  x01=slider(10, 40, initial=20, step=1),
  x02=slider(10, 40, initial=20, step=1)
  )

start # valores do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo não restrito: cada variedade com os seus parâmetros

n0 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S[vari])),
          data=db, start=start)
summary(n0)

#-----------------------------------------------------------------------------
# ajuste do modelo restrito: as variades possuem mesmo parâmetro S

start$S <- mean(start$S)
start

n1 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S)),
          data=db, start=start)
summary(n1)

#-----------------------------------------------------------------------------
# testa a hipótese do modelo nulo n1 tendo como alternativa o modelo n0

anova(n1, n0)

#-----------------------------------------------------------------------------

Obtenção de valores iniciais para o ajuste de duas curvas.

Ainda é necessário fazer a checagem das pressuposições envolvidas no modelo de regressão clássico: homocedasticidade (a variância dos desvios de regressão é constante ao longo a curva), normalidade (os erros aleatórios têm distribuição normal) e independência entre observações (normalmente a independência é garantida por um plano de amostragem adequado). Consulte Análise de resíduos em regressão não linear para saber como obter os 4 gráficos de análise de resíduos.

Esse post foi um aprimoramento de uma resposta na R-help e estimulado por citações na lista. Até a próxima ridícula.

Análise de resíduos em regressão não linear

Conjunto de gráficos para análise dos resíduos de um modelo de regressão não linear.

A análise dos resíduos de um modelo é feita para verificar a plausividade das pressuposições envolvidas. Os modelos linerares de regressão clássico, ou seja, aqueles em que as observações são realizações independentes (independência) e apresentam a mesma dispersão (homogeneidade de variância), podem ser ajustados no R com a função lm() e aov(). Para objetos dessas classes existe um método plot.lm() que apresenta os gráficos de análise de resíduo. Porém, modelos de regressão não linear podem ser ajustados com a função nls() que não possui um método para análise de resíduos.

Para aplicação de inferência (testes de hipótese, intervalos de confiança, etc), o modelo não linear é aproximado linearmente. Com isso eu quero dizer que podemos usar as mesmas técnicas de análise de resíduos para modelos de regressão não linear. A matriz gradiente do modelo (de derivadas da função em relação ao vetor de parâmetros) pode ser usada dentro da função lm() e com isso podemos obter facilmente os gráficos de análise de resíduos.

Esse procedimento ainda permite obter a análise de variância para o modelo de regressão não linear, fazer a partição ortogonal da soma de quadrados total em devido ao modelo de regressão e devido aos desvios de regressão. No entanto, esse não deveria ser o quadro apresentado se compararmos o que obtemos com modelo de regressão linear. No modelo linear particionamos a soma de quadrados total corrigida para o intercepto. O nosso modelo de regressão não linear se torna um modelo de intercepto se V e D forem zero (ver modelo do CMR). Portanto, o modelo nulo seria aquele em que estimamos apenas A e que podemos comparar com o modelo completo por meio da função/método anova.nls(). Esse é o teste de hipótese para o modelo de regressão não linear. Vale lembrar que alguns modelos não lineares não podem ser reduzidos a modelos de intercepto via restrição paramétrica. Para esses modelos deve-se ter cuidados ao usar o R², inclusive.

Abaixo apresento o procedimento para ajuste de um modelo não linear, extração da matriz gradiente e obtenção dos gráficos de resíduos. No final eu fiz uma função (R2()) para facilitar a obtenção do quadro de análise de variância e R². Não é porque você tem uma função para isso que vai sair aplicando ela em qualquer modelo não linear, ok? Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de postássio liberado em função do tempo

klib <- data.frame(k=c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67,
                     91.28, 105.22, 72.74, 81.88, 97.62, 90.14,
                     89.88, 113.22, 90.91, 115.39, 112.63, 87.51,
                     104.69, 120.58, 114.32, 130.07, 117.65, 111.69,
                     128.54, 126.88, 127.00, 134.17, 149.66, 118.25,
                     132.67, 154.48, 129.11, 151.83, 147.66, 127.30),
                   t=rep(c(15, 30, 45, 60, 75, 90,
                     120, 150, 180, 210, 240, 270), each=3))

#-----------------------------------------------------------------------------
# ajustando o modelo de regressão não linear aos dados

n0 <- nls(k~A*t/(V+t)+D*t, data=klib, start=list(A=90, V=15, D=0.21))
summary(n0)

#-----------------------------------------------------------------------------
# extraindo a matriz gradiente avaliada nas estimativas dos parâmetros

F <- attr(n0$m$fitted(), "gradient")
F

#-----------------------------------------------------------------------------
# passando a matriz gradiente para a lm(), importante remover intercepto

m0 <- lm(k~-1+F, data=klib)

#-----------------------------------------------------------------------------
# gráfico de análise dos resíduos

#png("f007.png", w=500, h=500);
par(mfrow=c(2,2), mar=c(5.1,4.1,4.1,2.1))
plot(m0)
mtext("Análise de resíduos para modelo de regressão não linear",
      outer=TRUE, line=-2, cex=1.4) 
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# veja que o ajuste é o mesmo pelas medidas abaixo

cbind(fitted(m0), fitted(n0))       # valores ajustados
cbind(residuals(m0), residuals(n0)) # valores preditos
c(summary(m0)$sig, summary(n0)$sig) # estimativa do desvio padrão residual
vcov(m0); vcov(n0)                  # matriz de covariância das estimativas

#-----------------------------------------------------------------------------
# quadro de anova com SQ de regressão e SQ de resíduo

anova(m0) # partição da soma de quadrados total
anova(n0, lm(k~1, klib)) # SQ do modelo não linear corrigido para intercepto

#-----------------------------------------------------------------------------
# função que retorna a anova e R2 para modelos de regressão não linear

R2 <- function(nls.obj){
  da <- eval(nls.obj$data)
  resp.name <- all.vars(summary(nls.obj)$formula)[1]
  form <- paste(resp.name, "~1", sep="")
  m0 <- lm(form, da)
  an <- anova(nls.obj, m0)
  sqn <- deviance(nls.obj)
  sqe <- deviance(m0)
  r2 <- 1-(sqn/sqe)
  aov <- data.frame(fv=c("regression","residuals"),
                    gl=c(-an$Df[2],an$Res.Df[1]),
                    sq=c(-an$Sum[2],an$Res.Sum[1]))
  aov$qm <- aov$sq/aov$gl
  aov$F <- c(aov$qm[1]/aov$qm[2], NA)
  aov$"Pr(>F)" <- c(1-pf(aov$F[1], df1=aov$gl[1], df2=aov$gl[2]), NA)
  names(aov) <- c(" ","Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")
  return(list(anova=aov, R2=r2))
}

R2(n0)

#-----------------------------------------------------------------------------

Bandas de confiança para modelo de regressão não linear

Potássio liberado em função do tempo de incubação.

Obter o intervalo de confiança para a resposta predita em um modelo de regressão linear é uma tarefa simples. A matriz do modelo é conhecida, ou seja, a matriz de derivadas parciais do modelo com relação a cada um dos parâmetros possui valores conhecidos e que não dependem dos valores dos parâmetros. Para um modelo de regressão não linear, a matriz de derivadas (gradiente) é desconhecida e dependente dos valores dos parâmetros. É por isso que se usa um procedimento iterativo para encontrar as estimativas que minimizam a soma de quadrados dos resíduos.

Uma vez que temos valores para os parâmetros (estimativas), podemos avaliar o gradiente com esses valores. Com isso somos capazes de calcular o intervalo de confiança para os valores preditos por um modelo de regressão não linear. Essa é uma tarefa meramente matricial para qual o R possui muitos recursos.

No R, obtemos esse intervalo para objetos de classe lm com a opção predict(..., interval="confidence"). Para modelos de classe nls, o método predict não possui essa opção. Sendo assim, temos que fazer as coisas no “braço”, o que é um exercício interessante em termos de programação e teoria de modelos de regressão. Você vai se sentir seguro ao saber como as coisas funcionam. O mesmo raciocínio pode ser empregado para modelos lineares generalizados.

No CMR abaixo eu forneço os procedimentos para obter intervalos de confiança para valores preditos de modelos de regressão não linear. Os aspectos teóricos sobre a obtenção dos intervalos você encontra nos bons livros de regressão e/ou modelos (não) lineares. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de postássio liberado em função do tempo

klib <- data.frame(k=c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67,
                     91.28, 105.22, 72.74, 81.88, 97.62, 90.14,
                     89.88, 113.22, 90.91, 115.39, 112.63, 87.51,
                     104.69, 120.58, 114.32, 130.07, 117.65, 111.69,
                     128.54, 126.88, 127.00, 134.17, 149.66, 118.25,
                     132.67, 154.48, 129.11, 151.83, 147.66, 127.30),
                   t=rep(c(15, 30, 45, 60, 75, 90,
                     120, 150, 180, 210, 240, 270), each=3))

#-----------------------------------------------------------------------------
# criando função que representa o modelo, retorna gradidente e hessiano

expo.der <- deriv3(~A*(1-exp(-log(2)*t/V))+D*t,
                   c("A", "V", "D"),
                   function(t, A, V, D) NULL)
str(expo.der)

#-----------------------------------------------------------------------------
# diagnose gráfica e primeiro chute

plot(k~t, data=klib, xlab="Período de incubacão (dias)",
       ylab="Potássio liberado acumulado (mg/kg de solo)")
A <- 90; V <- 20; D <- 0.2
curve(expo.der(x, A, V, D), add=TRUE, col=2)
start <- list(A=A, V=V, D=D)

#-----------------------------------------------------------------------------
# ajustando o modelo aos dados a partir dos valores iniciais via gráfico

n0 <- nls(k~expo.der(t, A, V, D), data=klib, start=start)
summary(n0)
confint(n0)

#-----------------------------------------------------------------------------
# valores preditos, gradiente e hessiano avaliado nos valores estimados

str(n0)
str(n0$m$fitted())
c(n0$m$fitted())
attr(n0$m$fitted(), "gradient")
attr(n0$m$fitted(), "hessian")

#-----------------------------------------------------------------------------
# obtenção dos valores preditos

pred <- data.frame(t=seq(0,300,l=100))
der <- do.call(expo.der, args=c(list(t=pred$t), as.list(coef(n0))))

F <- attr(der, "gradient") # gradiente avaliado no novo t
U <- chol(vcov(n0))
se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2))) # erro padrão

#-----------------------------------------------------------------------------
# gráficos dos observados, preditos com IC, legenda e equações

#png("f004.png", w=500, h=400); par(mar=c(5.1,4.1,2.1,2.1))
plot(k~t, data=klib, xlab="Período de incubacão (dias)",
     ylab="Potássio liberado acumulado (mg/kg de solo)",
     xlim=c(0,300), ylim=c(0,160))
matlines(pred$t, c(der)+
         outer(se, qt(c(.5, .025,.975), df=df.residual(n0))),
         type="l", col=c(1,2,2), lty=c(1,2,2))
legend("bottomright",
       legend=c("valores observados", "valores preditos",
         "intervalo de confiança (95%)"),
       lty=c(NA,1,2), col=c(1,1,2), pch=c(1,NA,NA), bty="n")
cf <- format(coef(n0), digits=3)
text(par("usr")[1], par("usr")[4], adj=c(-0.05,1.5),
     label=substitute(hat(k)[total]==a%.%(1-e^{-ln(2)%.%t/v})+d%.%t,
       list(a=cf[1], v=cf[2], d=cf[3])))
abline(v=coef(n0)["V"], h=coef(n0)["A"], col="gray70")
curve(expo.der(x, coef(n0)["A"], coef(n0)["V"], 0), add=TRUE, col=3)
curve(expo.der(x, 0, 0, coef(n0)["D"]), add=TRUE, col=3)
text(225, coef(n0)["A"], pos=3,
     label=substitute(hat(k)[fácil]==a%.%(1-e^{-ln(2)%.%t/v}),
       list(a=cf[1], v=cf[2])))
text(173, 38, pos=3, srt=18,
     label=substitute(hat(k)[difícil]==d%.%t, list(d=cf[3])))
#dev.off()

#-----------------------------------------------------------------------------