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.
