Archive

Posts Tagged ‘rpanel’

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.

Procedimento gráfico interativo para densidade com rpanel

Printscreen da janela gráfica e janela de controle do painel para o gráfico de densidade kernel usando funções disponíveis no pacote rpanel.

Alguns posts atrás eu apresentei procedimentos gráficos para obtenção de valores iniciais em regressão não linear. Dessa vez eu mudei o foco. Vou apresentar os procedimentos gŕaficos interativos disponíveis no pacote rpanel ilustrando o procedimento de estimação kernel de densidade.

O pacote rpanel possui funções para alterar componentes do gráfico. Os mais úteis são: deslizadores, caixas de seleção e caixas de texto. Exemplos de uso dessas funções você encontra na documentação do pacote e na página do pacote mantida pelo autor do mesmo, Adrian Bowman. Para ver o rpanel funcionando, assista o vídeo abaixo.

No CMR abaixo eu preparei um painel para controlar os argumentos kernel= e width= da função density(). Além disso, outro painel controla a posição do centro da banda. Ao alterar a opção de kernel, muda-se a função usada na estimação. Ao alterar o valor da largura de banda, muda-se a zona de influência ao redor do centro. Por último, movendo o centro sobre as observações, você verifica como a função se comporta em cada região.

No final, tem uma pequena porção de código que ilustra como obter probabilidades a partir de uma função de densidade kernel. O procedimento consiste em aproximar a função kernel e depois fazer uma integração numérica.

Para usar o rpanel, basta instalar da forma usual o pacote e suas dependências. Verifique a documentação do pacote para mais opções. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de redimento de grãos de 32 híbridos de milho

rend <- c(6388, 6166, 6047, 5889, 5823, 5513, 5202, 5172, 5166, 4975, 4778,
          4680, 4660, 5403, 5117, 5063, 4993, 4980, 4770, 4685, 4614, 4552,
          3973, 4550, 5056, 4500, 4760, 5110, 4960, 4769, 4849, 5230)

#-----------------------------------------------------------------------------
# tipos de kernel, tirado da documentação da função density()

kernels <- eval(formals(density.default)$kernel) # extrai os tipos de kernel
kernels
plot(density(0, from=-1.2, to=1.2, width=1, kernel="gaussian"), type="l",
     ylim=c(0, 2), xlab="rendimento de grãos", main="kernels")
for(i in 2:length(kernels))
  lines(density(0, width=1, kernel=kernels[i]), col=i)
legend(0.5, 2.0, legend=kernels, col=seq(kernels), lty=1, bty="n")

#-----------------------------------------------------------------------------
# painel para gráfico de densidade: largura de banda e tipo de função kernel

require(rpanel)

density.panel <- function(panel){
  ## argumento panel é uma lista, cada slot é um argumento
  ## x: vetor métrico de observações
  ## width: escalar métrico da largura de banda
  ## kernel: scalar string que informa o tipo de função kernel
  ## c0: escalar métrico que é a coodenada do centro da banda
  ## n: tamanho da amostra
  ##------------------------------------------------------------------
  ## gráfico de densidade e rug
  aux <- density(panel$x, width=panel$width, kernel=panel$kernel)
  plot(aux, main=NA)
  rug(panel$x)
  ##------------------------------------------------------------------
  ## faz o intervalo com o comprimento da banda
  arrows(panel$c0-0.5*panel$width, 0, panel$c0+0.5*panel$width, 0,
         length=0.1, code=3, angle=90, col=2)
  ##------------------------------------------------------------------
  ## faz uma seta que aponta para o valor da função
  y0 <- approx(aux$x, aux$y, xout=panel$c0)
  arrows(panel$c0, 0, panel$c0, y0$y, length=0.1, col=2)
  ##------------------------------------------------------------------
  ## desenha a função kernel acima do intevalo
  d <- density(panel$c0, width=panel$width, kernel=panel$kernel)
  lines(d$x, d$y/panel$n, col=2)
  ##------------------------------------------------------------------
  ## deve retornar a lista como última linha
  panel
  ##------------------------------------------------------------------
}

#-----------------------------------------------------------------------------
# alimenta os controladores, depois de rodar todos o painel gráfico está
# criado e então é só usar o mouse para alterar os componentes

# passar os argumentos que serão fixos, abre a janelinha
panel <- rp.control(x=rend, n=length(rend))

# controla a largura da banda
rp.slider(panel, width, 10, 1000, initval=150,
          showvalue=TRUE, action=density.panel)

# controla a coordenada do centro da banda
rp.slider(panel, c0, 3500, 6500, initval=5000,
          showvalue=TRUE, action=density.panel)

# controla o tipo de função kernel
rp.radiogroup(panel, kernel, 
              c("gaussian","epanechnikov","rectangular",
                "triangular","biweight","cosine","optcosine"),
              action=density.panel)

#-----------------------------------------------------------------------------
# obtendo valores de área abaixo da curva de densidade via aproximação da
# função densidade e integração numérica

dens <- density(rend, bw=150)
plot(dens)
rug(rend)
fdens <- approxfun(dens$x, dens$y) # cria uma função que é interpolação linear
i <- 4500; s <- 5500               # limites inferior e superior
integrate(fdens, i, s)             # faz uma integração numérica
polygon(c(i,dens$x[dens$x>i & dens$x<s],s),
        c(0,dens$y[dens$x>i & dens$x<s],0), col="gray90")
i <- 3600; s <- 6800               # limites inferior e superior
integrate(fdens, i, s)             # faz uma integração numérica
polygon(c(i,dens$x[dens$x>i & dens$x<s],s),
        c(0,dens$y[dens$x>i & dens$x<s],0), col="gray90")

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:, , ,