Arquivo

Posts Tagged ‘não linear’

Interface gráfica para ajuste de modelos de regressão não linear

Ajuste do modelo van Genuchten à curva de retenção de água do solo com interface gráfica desenvolvida com as ferramentas do pacote rpanel.

Modelos de regressão não linear são considerados quando alguma informação a priori existe sobre o fenômeno. Essa informação pode ser, por exemplo, que a curva seja sempre crescente, típico para curvas de crescimento/acúmulo. Em geral, esses modelos têm parâmetros com interpretação física/química/biológica e alguns parâmetros não têm, mas estão presentes para conferir flexibilidade.

Talvez um dos grandes gargalos no ajuste de modelos não lineares seja atribuição de valores iniciais para o método numérico de estimação. Para parâmetros com interpretação é possível obter valores ao ver gráficos de dispersão das variáveis ou pelo menos se tem idéia do intervalo em que o parâmetro se encontra. Para parâmetros sem interpretação a obtenção de valores iniciais é realmente trabalhosa.

Um procedimento muito útil para obter valores iniciais é fazer um grid de valores iniciais e plotar a curva corresponde para cada conjunto sobre o diagrama de dispersão dos dados. Você usa como valores iniciais os correspondentes à curva que melhor se aproxima do comportamento dos dados. Por outro lado, criar um grid de valores, sobrepor cada curva ao gráfico, escolher à curva que melhor se aproxima aos dados, fornecer os valores iniciais e estimar os parâmetros é algo de demanda tempo se o número de curvas à ajustar for grande. Seria de grande ajuda se esse processo pudesse ser automatizado ou facilitado, por exemplo, por meio de uma interface gráfica com botões de ação e controle de parâmetros por meio do mouse.

O pacote rpanel possui um conjunto de funções básicas que permitem construir interfaces gráficas para o usuário interagir com o R. O conjunto de funções compreende desde caixas de texto, caixas de seleção, deslizadores, até botões para disparo de alguma função. Usando essas ferramentas foi possível construir uma interface que permite selecionar o conjunto de dados, pré-ajustar a curva aos dados por meio de deslizadores e ajustar o modelo de regressão não linear. O código abaixo fornece as ferramentas para obter o que se vê nos gifs animados dessa matéria. Modificações devem ser feitas para outros dados/modelos. Minha intenção no futuro é aprimorar essas funções de forma que o usuário forneça o modelo e os limites dos intervalos para os parâmetros diretamente pela interface. Em caso de sucesso, isso pode virar um pacote ou complementar algum dos pacotes existentes dedicados à regressão não linear.

O primeiro exemplo usa gráficos do pacote graphics, básico do R. Os dados são de 10 curvas de retenção de água no solo onde ajustou-se o modelo van Genuchten (1980). O segundo conjunto de dados é de um experimento para avaliar o impacto da desfolha na produtividade do algodão. Para esse caso adotou-se o modelo potência e gráficos do pacote lattice. Clique na caixa abaixo para ver o código. Até à próxima ridícula.

#-----------------------------------------------------------------------------
# definições da sessão

require(rpanel)
require(lattice)
require(latticeExtra)
require(grid)

#-----------------------------------------------------------------------------
# dados de curva de retenção de água no solo

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

#-----------------------------------------------------------------------------
# ver

cras <- subset(cra, condi=="LVA3,5")
cras <- with(cras, aggregate(cbind(umid),
                             list(posi=posi, tens=tens, prof=prof), mean))
cras$caso <- with(cras, interaction(posi, prof))

xyplot(umid~log10(tens)|posi, groups=prof, data=cras, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste do modelo van Genuchten de modo iterativo com rpanel

# nlsajust é excutada quando clicar no botão "Ajustar"
nlsajust <- function(panel){
  ## ajuste do modelo não linear
  n0 <- try(nls(umid~tr+(ts-tr)/(1+(a*tens)^n)^(1-1/n), # modelo
                data=da,                                # dados
                start=start))                           # valores iniciais
  ## em caso de não convergência imprime mensagem de erro
  if(class(n0)=="try-error"){
    par(usr=c(0, 1, 0, 1))
    text(0.5, 0.5, "Não convergiu!\nAproxime mais.", col="red", cex=2)
  } else {
    ## coloca à curva ajusta sobre os pontos
    with(as.list(coef(n0)), curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n),
                                  add=TRUE, col=2))
    ## salva o ajuste numa lista
    aju[[i]] <<- n0
  }
  panel
}

vg <- function(panel){
  ## seleciona o subconjunto dos dados
  i <<- panel$caso
  da <<- subset(cras, caso==i)
  ## lista com valores iniciais vindos dos deslizadores
  start <<- panel[c("tr","ts","a","n")]
  ## diagrama de dispersão
  plot(umid~log10(tens), data=da,
       ylab=expression(Conteúdo~de~água~(g~g^{-1})),
       xlab=expression(log[10]~Tensão~matricial~(kPa)),
       main=expression(f(x)==tr+(ts-tr)/(1+(a*x)^n)^{1-1/n}))
  ## sobrepõe a curva controlada pelos deslizadores
  with(start, curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n),
                    add=TRUE, col=2, lty=2))
  panel
}

par(mar=c(4.1,4.2,3.1,1))
# cria objetos vazios que serão preenchidos durante processo
da <- c(); start <- list(); aju <- list(); i <- c()
# abre interface gráfica
panel <- rp.control()
# controla parâmetros do modelo
rp.slider(panel, ts, 0.4, 0.8, initval=0.6, showvalue=TRUE, action=vg)
rp.slider(panel, tr, 0, 0.5, initval=0.3, showvalue=TRUE, action=vg)
rp.slider(panel, a, 0.1, 5, initval=1.3, showvalue=TRUE, action=vg)
rp.slider(panel, n, 1, 5, initval=1.6, showvalue=TRUE, action=vg)
# seleciona o conjunto de dados
rp.listbox(panel, caso, vals=levels(cras$caso),
           title="Subconjunto", action=vg)
# cria botão "Ajustar"
rp.button(panel, action=nlsajust, title="Ajustar")

sapply(aju, coef)    # estimativas dos coeficientes
lapply(aju, summary) # summary dos modelos ajustados

#-----------------------------------------------------------------------------
# dados de peso de capulhos 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")
str(cap)
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")

#-----------------------------------------------------------------------------
# o exemplo a seguir mostra como usar com gráficos do pacote lattice

# função adapatada para outro modelo de regressão não linear
nlsajust <- function(panel){
  n0 <- try(nls(pcapu~f0-f1*desf^exp(C), data=da, start=start))
  if(class(n0)=="try-error"){
    trellis.focus("panel", nivel, 1, highlight=FALSE)
    grid.text(x=0.5, y=0.5, label="Não convergiu!\nAproxime mais.",
              gp=gpar(col="red"))
    trellis.unfocus()
  }
  trellis.focus("panel", nivel, 1, highlight=FALSE)
  with(as.list(coef(n0)), panel.curve(f0-f1*x^exp(C), add=TRUE, col=2))
  trellis.unfocus()
  print(coef(n0))
  aju[[i]] <<- n0
  panel
}

# função adapatada para outro modelo de regressão não linear
ptn <- function(panel){
  nivel <<- as.numeric(panel$nivel)
  i <<- levels(cap$estag)[nivel]
  da <<- subset(cap, estag==i)
  start <<- panel[c("f0","f1","C")]
  print({
    xyplot(pcapu~desf|estag, data=cap, layout=c(5,1), col=1,
           xlab="Nível de desfolha artificial", ylab="Peso de capulhos",
           main=expression(f(x)==f[0]-f[1]*x^exp(C)),
           strip=strip.custom(bg="gray90"))
  })
    trellis.focus("panel", nivel, 1, highlight=FALSE)
    with(start, panel.curve(f0-f1*x^exp(C), add=TRUE, col=2, lty=2))
    trellis.unfocus()
  panel
}

da <- c(); start <- list(); aju <- list(); i <- c(); nivel <- c()
panel <- rp.control()
rp.slider(panel, f0, 10, 50, initval=30, showvalue=TRUE, action=ptn)
rp.slider(panel, f1, 0, 30, initval=10, showvalue=TRUE, action=ptn)
rp.slider(panel, C, -5, 5, initval=0, showvalue=TRUE, action=ptn)
rp.listbox(panel, nivel, vals=1:nlevels(cap$estag),
           title="Painel", action=ptn)
rp.button(panel, action=nlsajust, title="Ajustar")

sapply(aju, coef)
lapply(aju, summary)

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

Ajuste do modelo potência à dados de redução na produção de algodão devido à desfolha com interface gráfica desenvolvida com as ferramentas do pacote rpanel.

Reparametrização do modelo quadrático para ponto crítico

Perfil de log-verossimilhaça representada pelo módulo da raíz da função deviance para cada um dos parâmetros do modelo quadrático reparametrizado. Linhas tracejadas indicam os limites dos intervalos de confiança de 99, 95, 90, 80 e 50%. Padrões de perfil diferentes de um V indicam intervalos assimétricos.

O emprego de funções polinomiais é muito frequente em ciências agrárias. O polinômio de segundo grau ou modelo quadrático é a escolha mais natural e imediata de modelo para quando a variável dependente apresentam uma tendência/relação não linear com a variável independente. Não raramente, o emprego desse modelo motiva também a estimação do ponto crítico (mínimo ou máximo). Para isso, o que se faz é estimar os parâmetros do modelo e via regras do cálculo obter a estimativa pontual do ponto crítico. O problema é que normalmente se necessita e não se faz, por desconhecimento de métodos, associação de incerteza à essa estimativa. É sobre isso que vou tratar agora: como associar incerteza ao ponto crítico de um modelo quadrático.

O modelo quadrático tem essa equação

m(x) = a+b\cdot x + c\cdot x^2,

em que a é o intercepto, ou seja, E(Y|x=0) = ab é valor da derivada do modelo na origem, ou seja, d(m_0(x))/dx = b, e c mede a curvatura. Derivando o modelo, igualando a zero e resolvendo, obtemos que o ponto crítico é

x_m = -b/(2\cdot c),

portanto, x_m é uma função das estimativas dos parâmetros, e como tal, deve apresentar incerteza associada.

Uma forma de associar incerteza é empregar o método delta para encontrar erros padrões. No R, o mesmo está disponível em msm::deltamethod(). Outra forma é usar um simulação boostratp e aplicar nas amostras bootstrap a relação acima. E o terceiro, e que eu mais gosto pois tem uma série de vantagens, é usar o perfil de log-verossimilhança.

Para usar o perfil de verossimilhança eu tive que escrever o modelo de forma que o ponto crítico fosse um parâmetro dele. Reparametrizar um modelo nem sempre é uma tarefa fácil. Nesse caso o procedimento escrito é curto e simples de entender. Pórem eu demorei cerca de uma hora com tentativa-erro até usar a seguinte abordagem: imagine uma função que tenha ponto de máximo. Podemos aproximar essa função ao redor de um valor x_0 por série de Taylor. Para uma aproximação de segundo grau temos

f(x) \approx f(x_0)+f'(x_0)(x-x_0)+0.5\cdot f''(x_0)(x-x_0)^2.

Agora imagine que x_0 é um ponto crítico da função. O termo linear da aproximação é 0 no ponto crítico por definição. Se a função que estamos aproximando for perfeitamente quadrática, esse procedimento nos leva a reparametrização que desejamos. Assim, o modelo quadrático reparametrizado é

n(x) = y_m+c\cdot (x-x_m)^2.

em que x_m é o nosso parâmetro de interesse, o ponto crítico, y_m é o valor da função no ponto crítico, ou seja, E(Y|x=x_m) = y_m, c é a métade da segunda derivada da função no ponto crítico, ou seja, c=0.5\cdot f''(x_m), é uma curvatura. Então nos reparametrizamos o model quadrático para ter o ponto crítico como parâmetro. Nessa reparametrização obtvemos 3 parâmetros com interpretação simples. Em termos de interpretação, se c=0 não temos curvatura, portanto, não temos ponto crítico e devemos ajustar um modelo linear da reta. Se temos curvatura podemos estimar o ponto crítico e o valor da função no ponto crítico. Então imagine que você estuda doses de nutrientes para uma cultura, não é fundamental saber qual a dose que confere o máximo e qual a produção esperada nessa dose? Melhor que isso é poder associar incerteza à essas quantidades.

O modelo reparametrizado, diferente do original, é um modelo não linear. Por isso, para estimar os parâmetros no R usaremos a função nls() e não a lm(). (Visite os posts passados para ver mais coisas relacionadas à modelos não lineares.) Para obter os intervalos de perfil de log-verossimilhança usaremos os métodos confint() e profile(). Consulte a documentação dessas funções para saber o que elas de fato fazem.

No CMR abaixo eu ilustro a aplicação desses métodos. É recomendando ao leitor que compare os intervalos de confiança, que repita os procedimento simulando os dados com outros parâmetros, tamanhos e espaços, e se for o caso, que use dados reais. Futuros artigos do Rídiculas vão incluir comparação de modelos quadráticos (teste de hipótese) e o modelo quadrático para duas regressoras (bidimensional). Até a próxima rídicula.

#-----------------------------------------------------------------------------
# dados simulados

set.seed(23456)
x <- 1:12
y <- 3+0.17*x-0.01*x^2+rnorm(x,0,0.05)
plot(x, y)                         # dados observados
curve(3+0.17*x-0.01*x^2, add=TRUE) # curva que gerou os dados

#-----------------------------------------------------------------------------
# ajuste do y = a+b*x+c*x², é um modelo linear

m0 <- lm(y~x+I(x^2))
summary(m0)

#-----------------------------------------------------------------------------
# método delta

require(msm)
coef <- coef(m0)[2:3]; vcov <- vcov(m0)[2:3,2:3]
xm <- -0.5*coef[1]/coef[2]; xm                         # estimativa
sd.xm <- deltamethod(~(-0.5*x1/x2), coef, vcov); sd.xm # erro padrão
xm+c(-1,1)*qt(0.975, df=df.residual(m0))*sd.xm         # IC 95%

#-----------------------------------------------------------------------------
# simulação (bootstrap paramétrico)

require(MASS)
aa <- mvrnorm(1000, mu=coef, Sigma=vcov)
xm.boot <- apply(aa, 1, function(i) -0.5*i[1]/i[2])
mean(xm.boot)                      # média bootstrap
sd(xm.boot)                        # desvio-padrão bootstrap
quantile(xm.boot, c(0.025, 0.975)) # IC boostrap

#-----------------------------------------------------------------------------
# ajuste do y = ym+c*(x-xm)², é um modelo não linear

plot(x, y)
abline(h=3.7, v=8.5, lty=3) # aproximação visual do chute
start <- list(ym=3.7, xm=8.5, c=coef(m0)[3])
max(fitted(m0))             # boa aproximação de chute
x[which.max(fitted(m0))]    # boa aproximação de chute
n0 <- nls(y~ym+c*(x-xm)^2, start=start)
summary(n0) # estimativa do xm que é parâmetro do modelo, estimado diretamente

#-----------------------------------------------------------------------------
# gráfico dos observados com a curva ajustada

plot(x, y)
with(as.list(coef(n0)),{
  curve(ym+c*(x-xm)^2, add=TRUE)
  abline(h=ym, v=xm, lty=3)
})

#-----------------------------------------------------------------------------
# veja a igualdade dos modelos

cbind(coef(m0),
      with(as.list(coef(n0)), c(ym+c*xm^2, -2*c*xm, c)))
deviance(m0)
deviance(n0)
plot(fitted(m0), fitted(n0)); abline(0,1)

#-----------------------------------------------------------------------------
# intervalos de confiança

confint.default(n0) # IC assintóticos, são simétricos por construção
confint(n0)         # IC baseado em perfil de log-verossimilhança

#-----------------------------------------------------------------------------
# gráficos de perfil de verossimilhança

#png(file="f030.png", width=500, height=175)
par(mfrow=c(1,3))
plot(profile(n0))
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# representação do ajuste

pred <- data.frame(x=seq(0,13,l=30))
pred <- cbind(as.data.frame(predict(m0, newdata=pred, interval="confidence")),
              pred)
str(pred)
ic <- confint(n0); ic

plot(x, y, type="n")              # janela vazia
with(pred,                        # bandas de confiança
     polygon(c(x,rev(x)), c(lwr,rev(upr)), col="gray90", border=NA))
points(x, y)                      # valores observados
with(pred, lines(x, fit))         # valores preditor
abline(v=ic[2,], h=ic[1,], lty=2) # IC para xm e ym

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

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á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()

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

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

Quem alguma vez já ajustou modelo de regressão não linear à dados, concorda comigo que o gargalo da tarefa é encontrar bons valores iniciais (ou c.h.u.t.e., cálculo hipotético universal técnico estimativo) para fornecer ao método iterativo de obtenção de estimativas. Nessa ridícula vou fornecer um procedimento interativo para obter bons valores. Basicamente você vai alterar valores dos parâmetros com o mouse e ver a curva se mover até representar os dados.

 

Deslizadores para alterar os valores dos parâmetros do modelo de regressão não linear.

Por definição, um modelo é não linear quando a primeira derivada da função com relação a algum dos parâmetros, ainda depende de algum dos parâmetros. Para ilustrar, vejamos como fica as derivadas para o modelo Michaelis-Mentem:

f(x) = \displaystyle \frac{\beta_0 \times x}{\beta_1 + x}
\displaystyle \frac{\partial f(x)}{\partial \beta_0} = \displaystyle \frac{x}{\beta_1+x}
\displaystyle \frac{\partial f(x)}{\partial \beta_1} = \displaystyle \frac{\beta_0 x}{(\beta_1+x)^2} ,

#-----------------------------------------------------------------------------
# dados para a curva característica de água do solo

cra <- data.frame(th=c(0.3071, 0.2931, 0.2828, 0.2753, 0.2681,
                    0.2628, 0.2522, 0.2404, 0.2272, 0.212,
                    0.1655, 0.1468, 0.1205, 0.1013, 0.073),
                  psi=c(10, 19, 30, 45, 63, 64, 75, 89, 105,
                    138, 490, 3000, 4100, 5000, 26300))
str(cra)

#-----------------------------------------------------------------------------
# criando função que representa o modelo

vanG <- function(x, thmin, thmax, alpha, n){
  thmin+(thmax-thmin)/((1+(alpha*10^x)^n)^(1-1/n))
}

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

plot(th~log10(psi), cra)
thmin <- 0.05; thmax <- 0.4; alpha <- 0.03; n <- 1.5     # primeiro chute
curve(vanG(x, thmin, thmax, alpha, n), add=TRUE, col=2)

#-----------------------------------------------------------------------------
# usando os recursos gráficos do R para mudar os valores dos parâmetros

library(gWidgetsRGtk2)     # carrega o pacote que cria janelas interativas

#-----------------------------------------------------------------------------
# lista com os chutes para intervalo, slots devem ter os nomes dos parâmetros

limits <- list(thmin=c(0,0.1),        # chute para o intevalo de thmin
               thmax=c(0.25,0.440),   # chute para o intevalo de thmax
               alpha=c(0.01, 0.1),    # chute para o intevalo de alpha
               n=c(1.01,2))           # chute para o intevalo de n
start <- list()                       # lista com os valores para chute

#-----------------------------------------------------------------------------
# função que será atualizada a cada movimento do deslizador
# parâmetros dentro de svalue() são controlados, nomes igual aos da lista

plot.chute <- function(...){
  ## faz o gráfico de dispersão
  plot(th~log10(psi), cra)
  ## sobrepõe a curva com os valores dos deslizadores
  curve(vanG(x, svalue(thmin), svalue(thmax), svalue(alpha), svalue(n)),
        add=TRUE, col=2)
  ## reescreve o start com os valores dos delizadores, para usar na nls()
  start <<- list(thmin=svalue(thmin), thmax=svalue(thmax),
                 alpha=svalue(alpha), n=svalue(n))
}

#-----------------------------------------------------------------------------
# criação da janela com deslizadores
# na primeira chamada escolher uma das opções (sempre escolho a 1)
##  Select a GUI toolkit 
##   1: gWidgetsRGtk2
##   2: gWidgetstcltk
# essa função pode estar num arquivo fn.R e carregada com source("fn.R")

w <- gwindow("Caixa com deslizadores para controlar parâmetros")
tbl <- glayout(cont=w)
for(i in 1:length(limits)){
  tbl[i,1] <- paste("Controle", names(limits)[i])
  tbl[i,2, expand=TRUE] <- (assign(names(limits)[i],
             gslider(from=limits[[i]][1],
                     to=limits[[i]][2],
                     by=diff(limits[[i]])/20,
                     value=mean(limits[[i]]),
                     container=tbl, handler=plot.chute)))
}

#-----------------------------------------------------------------------------
# agora com a caixa criada, basta chamar a função e mover os deslizadores

plot.chute()
start

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

n0 <- nls(th~thmin+(thmax-thmin)/((1+(alpha*psi)^n)^(1-1/n)),
          data=cra, start=start)
summary(n0)

#-----------------------------------------------------------------------------
# fazendo o gráfico com os dados e a curva estimada

plot(th~log10(psi), cra)
for(i in 1:length(coef(n0))){ assign(names(coef(n0))[i], coef(n0)[i]) }
curve(vanG(x, thmin, thmax, alpha, n), add=TRUE, col=2)

#-----------------------------------------------------------------------------
# outra forma de adicionar a curva (ainda pode-se usar predict())

lis <- list(x=NULL, body(vanG))
lis <- append(lis, coef(n0), after=1)
plot(th~log10(psi), cra)
tmp <- as.function(lis)
curve(tmp, add=TRUE, col=3)

#-----------------------------------------------------------------------------
# ainda outra maneira

lis <- c(list(x=NULL), as.list(coef(n0)), body(vanG))
plot(th~log10(psi), cra)
tmp <- as.function(lis)
curve(tmp, add=TRUE, col=4)

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

Para outros modelos/dados basta fazer as modificações pertinentes. As funcionalidades contidas no pacote gWidgetsRGtk2 são verdadeiras ferramentas no ensino de estatística. Sempre que possível vou postar novidades usando esse pacote. Até a próxima ridícula.