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

December 7, 2012
By

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

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 language="("tr","ts","a","n")"][/c]
  ## 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 language="("f0","f1","C")"][/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.


Tags: , , ,

Comments are closed.