Arquivo

Archive for dezembro \07\UTC 2012

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.

Fatorial duplo com dois tratamentos adicionais

Número de sementes em 40 infectadas por Fusarium (tranformado) em função da aplicação de fungicida em combinação com polímero.

Experimentos fatoriais são aqueles em que dois ou mais fatores estão presentes. Nos fatoriais incompletos nem todas as combinações possíveis estão presentes. Algumas das combinações podem não ser de interesse e por isso não são realizadas, alguma combinação é acidentalmente perdida durante o experimento ou certas combinações não identificam novos níveis, como por exemplo aplicação do produto A e B na dose 0 que, apesar dos rótulos diferentes, são o mesmo tratamento.

Um caso comum de fatorial incompleto são os fatoriais (completos) com a presença de tratamentos adicionais. Experimentos fatorais são realizados para verificar se os fatores interagem e então escolher combinação de níveis mais promissora para o processo estudado. Caso não exista interação, entende-se que os fatores têm efeito aditivos. Em qualquer caso, adicionam-se tratamentos adicionais ao experimento quando se deseja comparar o promissor aos de referência ou extremos.

Esse é um tipo de ensaio típico na avaliação do efeito de insumos agrícolas (fungicidas, insecticidas e herbicidas). Os ensaios testam produtos e suas doses e precisam ter as chamadas testemunhas para avaliar a eficiência das combinações geradas. Devido não ser um delineamento de níveis completamente cruzados a análise do experimento fatorial com tratamentos adicionais têm alguns detalhes que exigem mais atenção.

Nessa matéria faço análise de um fatorial com dois tratamentos adicionais. Os dados analisados foram retirados do Zimmermann (2004) de um experimento que estudou três fungicidas (benlate, captam, derosal) em três formas de aplicação (puro, antes de polímero, misturado com polímero). Adicionou-se dois tratamentos que foram a testemunha e aplicação do polímero puro. A resposta
avaliada foi o número de sementes infectadas com Fusarium em um total de 40 sementes.

O código está comentado para informar o leitor. São realizados a importação e reordenação dos níveis, a análise de variância com decomposição das somas de quadrados, contrastes entre efeitos, e testes sobre as médias ajustadas. Foi dada prioridade para os procedimentos matriciais de análise que podem ser adaptados para situações de desbalanceamento. Até a próxima ridícula.

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

require(lattice)
require(latticeExtra)
require(grid)
require(MASS)
require(doBy)
require(gmodels)
require(multcomp)

#-----------------------------------------------------------------------------
# leitura

fu <- read.table(
  "http://www.leg.ufpr.br/~walmes/data/zimmermann_fusarium.txt",
                 header=TRUE, sep="\t")
str(fu)

#-----------------------------------------------------------------------------
# modificar ordem dos níveis, primeiro níveis fatoriais, depois adicionais

levels(fu$fungicida)               # ordem correta
levels(fu$aplicacao)               # precisa alterar
levels(fu$aplicacao)[c(4,1,2,3,5)] # para ficar assim

fu$aplicacao <- factor(fu$aplicacao,
                       levels=levels(fu$aplicacao)[c(4,1,2,3,5)])
levels(fu$aplicacao) # agora ordem correta

#-----------------------------------------------------------------------------
# análise exploratória

xtabs(~fungicida+aplicacao, fu)

# os níveis não se cruzam completamente, portanto é um fatorial incompleto
# no caso 3x3+2, fatorial com 2 tratamentos adicionais

xyplot(infec40~tratamentos, fu) # resposta do tipo contagem de sucessos
xyplot(infec40~tratamentos, fu, jitter.x=TRUE)

#-----------------------------------------------------------------------------
# Vamos primeito analisar fora da estrutura fatorial, ou seja, considerando
# 3x3+2 = 11 níveis de um único fator.
# Vamos considerar distribuição para reproduzir os resultados do livro
# do Zimmermann (2004). Todavia, não é a análise mais correta devido os dados
# serem discretos. A distribuição binomial parece ser mais indicada.

m0 <- lm(infec40~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m0); layout(1) # não tão ruim, tentar boxcox?

m0 <- lm(infec40+1~tratamentos, data=fu) # soma 1 para eliminar os 0
boxcox(m0) # aponta log

m1 <- lm(log(infec40+1)~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m1); layout(1) # alguma melhoria?

# Zimmermann usa asin(sqrt(x/40)).
m2 <- lm(asin(sqrt(infec40/40))~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m2); layout(1) # alguma melhoria?

#-----------------------------------------------------------------------------
# Visualmente todos os resíduos são semelhantes, apresentam fugas.

lapply(list(m0, m1, m2), function(i) shapiro.test(residuals(i)))
lapply(list(m0, m1, m2), function(i) ks.test(scale(residuals(i)), "pnorm"))

# Pelos p-valores, a tranformação Box-Cox foi melhor, embora os resíduos ainda
# não apresentem um comportamento satísfatório.
# O foco não é testar transformações, então vamos seguir o que Zimmermann fez
# e trabalhar com asin(sqrt(x/40)).

#-----------------------------------------------------------------------------
# análise exploratória com a variável transformada

fu$y <- asin(sqrt(fu$infec40/40)) # cria a transformada

xyplot(y~tratamentos, fu, jitter.x=TRUE)

xyplot(y~fungicida, groups=aplicacao, fu, jitter.x=TRUE,
       type=c("p","a"), auto.key=TRUE)

# Apresenta indicativos de efeito de interação e que os adicionais
# são piores que o os fatoriais.

#-----------------------------------------------------------------------------
# análise na estrutura fatorial 3x3+2

m3 <- lm(y~fungicida*aplicacao, data=fu)
anova(m3)

# Os graus de liberdade são o que são pois:
# fungicida 4: são 5 niveis (3 fatoriais, 2 adicionais), 5-1=4;
# aplicacao 4: são 5, mas 2 estão confundidos (os adicionais), 5-2-1=2;
# fun:apli  4: (3-1)*(3-1)=4;

# O que queremos é a anova com somas de quadrados particionadas da seguinte
# forma
# fonte de variação     gl
# fungicidas             2
# aplicacao              2
# fun:apli               4
# fatorial vs adicionais 1
# polimero vs testemunha 1

#-----------------------------------------------------------------------------
# Na minha opinião, obter esse quadro de análise de variância é desperdício
# de tempo pois:
# * não tem sentido prático comparar o efeito médio dos fatoriais contra o
#   efeito médio dos adicionais, principalmente sabendo que a interação
#   fungicida:aplicacao foi significativa;
# * a interação fungicida:aplicacao significativa torna não interpretável
#   o teste para os efeitos principais;
# * polímero vs testemunha tem só um grau de liberdade, posso comprar os
#   efeitos por um testes t;

coef(m3) # tem NA para os efeitos não estimados (não presentes)
efstm <- colnames(vcov(m3)) # nomes dos efeitos que foram estimados
b <- coef(m3)[efstm] # só os estimados
b

# matriz de coeficientes para estimação das médias ajustadas
X <- popMatrix(m3, effect=c("fungicida", "aplicacao"))
str(X)

Xn <- attr(X, "grid") # 25 linhas, 5x5=25 (faz pensando que é completo)
Xo <- unique(fu[,c("fungicida","aplicacao")]) # que ocorrem (3x3+2)
Xo

#-----------------------------------------------------------------------------
# calculando as médias ajustadas matricialmente

rownames(X) <- apply(Xn, 1, paste, collapse=":")
on <- apply(Xo, 1, paste, collapse=":")
X <- X[on,efstm] # só linhas e colunas que correspondem à efeitos que ocorrem

X%*%b                                             # são as médias ajustadas
with(fu,
     tapply(y, list(fungicida, aplicacao), mean)) # são médias amostrais

# elas coincidem por ser um experimento regular (balanceado/ortogonal)

#-----------------------------------------------------------------------------
# fazendo contrastes
# não funcionam com presenta de NA
# contrast::contrast()
# multcomp::glht()
# para isso usar a gmodels::estimable()

m4 <- aov(formula(m3), data=fu) # aov não NA em coef

# polímero contra a testemunha
p.vs.t <- matrix(X[10,]-X[11,], nrow=1)
rownames(p.vs.t) <- "p.vs.t"
estimable(m4, cm=p.vs.t)

# fatorial contra adicional
f.vs.a <- matrix(apply(X[1:9,], 2, mean)-apply(X[10:11,], 2, mean), 1)
rownames(f.vs.a) <- "f.vs.a"
estimable(m4, cm=f.vs.a)

estimable(m4, cm=rbind(p.vs.t, f.vs.a)) # os dois de uma só vez

# E aí, complicado fazer tudo isso, né? Qual a sensação em saber de que
# pouco disso vai ser útil? Se deu interação teremos que estudar níveis
# de um fator fixando os níveis do outro.

#-----------------------------------------------------------------------------
# fazendo a decomposição das somas de quadrados dentro da anova

cf <- contrasts(C(fu$fungicida, helmert))
cf # o segredo está em montar os contrastes corretamente
cf[,3] <- c(-1,-1,-1,1.5,1.5) # fatorial contra adicional
cf[,4] <- c(0,0,0,-1,1)       # polímero vs testemunha
cf
contrasts(fu$fungicida) <- cf

m5 <- aov(y~fungicida*aplicacao, data=fu)
summary(m5, expand=FALSE,
        split=list("fungicida"=list(
                     "fungicida"=1:2,
                     "fat.vs.adi"=3,
                     "pol.vs.tes"=4)))

# Esse procedimento depende de balanceamento pois as somas de quadrados são
# sequênciais.

#-----------------------------------------------------------------------------
# Uma vez que deu interação é de prática fazer o desdobramento que consiste
# em estudar o efeito de aplicação dentro dos níveis de fungicida e
# vice-versa. Além do mais, por ter tratamentos adicionais, pode-se comparar
# com os níveis dos fatoriais. Vamos listas o número de hipóteses que teremos
# que fazer:
# * 3 níveis de aplicação em cada fungicida: 3*choose(3,2)=9
# * 3 níveis de fungicida em cada aplicação: 3*choose(3,2)=9
# * testemunha contra cada um dos 9 níveis do 3x3: 9
# * polimero contra cada um dos 9 níveis do 3x3: 9
# Com isso são 4*9=36 comparações. Já que é assim, vou fazer todos
# contra todos, choose(11, 2)=55, que o usuário pode olhar a que lhe
# interessa, afinal de contas, de 36 para 55 não têm diferença

mm <- model.matrix(m4)
dim(mm)
mm <- mm[,efstm]

# modelo de classe lm equivalente à m5, sem NA
m6 <- lm(y~-1+mm, data=fu)
anova(m6)

cbind(coef(m4), coef(m6)) # são o mesmo modelo

#-----------------------------------------------------------------------------
# médias ajustadas (ma)

cima <- confint(glht(m6, linfct=X)) # intervalos de confiança para ma
cima # tem cobertura *global* de 95% (e não indivídual)

str(cima)
cima$confint

result <- cbind(Xo, cima$confint)
str(result)

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm")

#-----------------------------------------------------------------------------
# carrega algumas funções que serão úteis para comparar médias e fazer
# gráficos com intervalos de confiança

objs <- ls()
objs <- c(objs, "apc", "prepanel.cbH", "panel.cbH")
source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functionsi.R",
       encoding="latin1")

rm(list=ls()[!(ls()%in%objs)])
ls()

#-----------------------------------------------------------------------------
# comparando as médias ajustadas, choose(11, 2), contrastes de Tukey

str(X)
Xc <- apc(X)
str(Xc)

#-----------------------------------------------------------------------------
# fazendo as comparações com a glht, método para correção de p-valor
# fdr, o single-step demora muito

c0 <- summary(glht(m6, linfct=Xc), test=adjusted(type="fdr"))
c0

c0$focus <- "comparacoes" # para poder usar a cld()
cld(c0)

result$cld <- cld(c0)$mcletters$Letters

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

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm",
        panel=function(x, y, z, centers, subscripts, ...){
          panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...)
          panel.text(centers, z,
                     paste(format(centers, digits=2),
                           result$cld[subscripts], sep=" "), pos=3)
        })

#-----------------------------------------------------------------------------
# gráfico final

result$desloc <- 0.05*c(0,-1,1,0,0)[match(result$aplicacao,
                                          levels(fu$aplicacao))]

#png("f048.png", 500, 400)
xyplot(y~fungicida, groups=aplicacao, data=fu, jitter.x=TRUE, amount=0.05,
       xlab="Fungicida", ylab=expression(asin(sqrt(x/40))),
       prepanel=prepanel.cbH, ly=result$lwr, uy=result$upr,
       auto.key=list(title="Aplicação", cex.title=1.2, columns=2,
         type="o", divide=1, lines=TRUE, points=FALSE))+
  as.layer(xyplot(Estimate~fungicida, groups=aplicacao, data=result, type="l",
                  ly=result$lwr, uy=result$upr,
                  desloc=result$desloc,
                  cty="bars",
                  panel.groups=panel.cbH,
                  panel=panel.superpose))
trellis.focus("panel", 1, 1, h=FALSE)
x <- as.numeric(result$fungicida)+result$desloc
y <- result$Estimate
lab <- paste(format(y, digits=2), result$cld, sep=" ")
a <- paste(rep("0", max(nchar(lab))), collapse="")
grid.rect(x=unit(x, "native"), y=unit(y, "native"),
          width=unit(1, data=a, "strwidth"),
          height=unit(1.2, data=a, "strheight"),
          hjust=-0.1, gp=gpar(fill="gray90", col=NA, fontsize=10))
grid.text(lab, x=unit(x, "native"), y=unit(y, "native"),
          hjust=-0.2, gp=gpar(col="black", fontsize=10))
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------
Categorias:delineamento Tags:, ,