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

April 9, 2011

This post was kindly contributed by Ridículas - go there to comment and to read the full post.

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.