Arquivo

Posts Tagged ‘curve’

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

Obtenção de valores iniciais para o modelo duplo van Genuchten de 7 parâmetros com os recursos da função manipulate do RStudio.

No post Método gráfico interativo para valores iniciais em regressão não linear eu apresentei alguns procedimentos gráficos para obter bons valores iniciais para ajuste de modelos de regressão não linear. Ao obter valores iniciais pelo procedimento gráfico facilita-se a convergência do método de estimação e você passa a entender/interpretar melhor os parâmetros do modelo.

Nesse post vou apresentar os métodos gráficos construídos com os recursos do pacote manipulate que é parte do RStudio, um recente editor de scripts R que tem ganhado muitos adeptos por ser muito amigável. Eu já postei o uso dessas funções em RStudio e a função manipulate() onde falo sobre as opções das função.

Vou apresentar o procedimento com dados reais. Preciso encontrar chutes iniciais para o modelo duplo van Genuchten (ver artigo original). Esse modelo tem 7 parâmetros na seguinte equação

\theta(\Psi) = \theta_{res}+\displaystyle \frac{\theta_{pmp}-\theta_{res}}{(1+(\alpha_{tex}\Psi)^{n_{tex}})^{1-1/n_{tex}}}+\displaystyle \frac{\theta_{sat}-\theta_{pmp}}{(1+(\alpha_{est}\Psi)^{n_{est}})^{1-1/n_{est}}},

em que \theta e \Psi são valores observados e o restante são os parâmetros da função. Em geral, quanto maior o número de parâmetros a serem estimados num modelo de regressão mais demorara será a convergência porque o espaço da solução é de alta dimensão. Além do mais, com muitos parâmetros a serem estimados, pode haver fraca de identificabilidade do modelo. Dessa forma, bons chutes iniciais podem ajudar na convergência. Uma ligeira noção do campo de variação de cada parâmetro o usuário precisa ter (limites inferior e superior). Você pode inicialmente usar intervalos amplos e ir reduzindo o comprimento do intervalo a medida que observar quais são os valores mais próximos dos ótimos.

#-----------------------------------------------------------------------------
# dados para a Curva de Retenção de Água, umidade (theta) e tensão (psi)

cra <- data.frame(psi=c(0.01,1,2,4,6,8,10,33,60,100,500,1500,2787,4727,6840,
                    7863,9030,10000,10833,13070,17360,21960,26780,44860,
                    69873,74623,87287,104757,113817,147567,162723,245310,
                    262217,298223),
                  theta=c(0.779,0.554,0.468,0.406,0.373,0.36,0.344,0.309,
                    0.298,0.292,0.254,0.241,0.236,0.233,0.223,0.202,0.172,
                    0.187,0.138,0.098,0.07,0.058,0.052,0.036,0.029,0.0213,
                    0.0178,0.0174,0.0169,0.0137,0.0126,0.0109,0.0106,0.0053))

#-----------------------------------------------------------------------------
# gráfico dos valores observados

plot(theta~log10(psi), data=cra)

#-----------------------------------------------------------------------------
# função do modelo duplo van Genuchten, 7 parâmetros

dvg <- function(x, ts, ti, tr, ae, at, ne, nt){
  tr+(ti-tr)/((1+(at*10^x)^nt)^(1-1/nt))+(ts-ti)/((1+(ae*10^x)^ne)^(1-1/ne))
}

dvg(log10(c(1,10,100,1000,10000)),         # log 10 das tensões
    0.7, 0.2, 0.05, 1.3, 0.0001, 1.5, 3.5) # valor dos parâmetros

#-----------------------------------------------------------------------------
# procedimento gráfico para obter chutes e conhecer a função dos parâmetros

require(manipulate) # carrega pacote que permite manipulação gráfica
start <- list()     # cria uma lista vazia para receber os valores finais

manipulate({
             plot(theta~log10(psi), data=cra)
             curve(dvg(x, ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt), add=TRUE)
             start <<- list(ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt)
           },
           ts=slider(0.7, 0.9, initial=0.8),
           ti=slider(0.15, 0.25, initial=0.2),
           tr=slider(0, 0.10, initial=0.05),
           ae=slider(1.01, 3, initial=1.3),
           at=slider(0, 0.0001, initial=0.00005),
           ne=slider(1.01, 3, initial=1.65),
           nt=slider(1.8, 5, initial=4.3))

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

start # valores salvos do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo aos dados com os chutes do procedimento gráfico

n0 <- nls(theta~tr+(ti-tr)/((1+(at*psi)^nt)^(1-1/nt))+
          (ts-ti)/((1+(ae*psi)^ne)^(1-1/ne)),
          data=cra, start=start, trace=TRUE)
summary(n0) # quadro de estimativas

#-----------------------------------------------------------------------------
# gráfico dos dados com a curva estimada

lis <- c(list(x=NULL), as.list(coef(n0)), body(dvg))
plot(theta~log10(psi), cra,
     ylab=expression(Conteúdo~de~água~no~solo~(theta*","~g~g^{-1})),
     xlab=expression(Tensão~matricial~(Psi*","~kPa)), xaxt="n")
tmp <- as.function(lis)
curve(tmp, add=TRUE, col=2, lwd=1.5)
axis(1, at=-2:5, label=as.character(10^(-2:5)), lwd.ticks=2)
s <- log10(sort(sapply(1:9, function(x) x*10^(-3:6))))
axis(1, at=s, label=NA, lwd=0.5)
abline(v=-2:6, h=seq(0,1,0.05), lty=3, col="gray50")

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

O procedimento permite encontrar chutes iniciais para duas ou mais curvas. Basta escrever adequadamente as funções dentro da manipulate(). Abaixo estão os códigos para obter os valores inicias para o modelo logístico usado no contexto de epidemiologia. Em duas variedades de uma espécie vegetal foram observadas a severidade da doença ao longo do tempo. O objetivo é encontrar estimativas dos parâmetros para a função (que no caso é a logística) que associa severidade ao tempo e comparar as funções entre as variedades. A comparação das funções pode ser feita fazendo-se testes de hipóteses para funções dos parâmetros (em geral contrastes) ou testes envolvendo modelos aninhados, como você verá a seguir.

#-----------------------------------------------------------------------------
# dados de severidade em função do tempo para duas variedades

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

#-----------------------------------------------------------------------------
# gráficos exploratórios

require(lattice)
xyplot(V1+V2~DAI, data=da, type=c("p","a"))

# ajustar modelos. separado ou junto?
# qual modelo? quais os chutes?

#-----------------------------------------------------------------------------
# mudando a disposição dos dados

require(reshape)
db <- melt(da, id=c("DAI","rep"))
names(db)[3:4] <- c("vari","sever")
str(db)

#-----------------------------------------------------------------------------
# carrega pacote e faz a função para encontrar via deslizadores valores ótimos

require(manipulate)
start <- list()

manipulate({
  plot(sever~DAI, db, col=ifelse(db$vari=="V1",1,2))
  A1 <- a1; S1 <- s1; X01 <- x01
  A2 <- a2; S2 <- s2; X02 <- x02
  curve(A1/(1+exp(-(x-X01)/S1)), add=TRUE, col=1)
  curve(A2/(1+exp(-(x-X02)/S2)), add=TRUE, col=2)
  start <<- list(A=c(A1,A2), S=c(S1,S2), X0=c(X01,X02))
  },
  a1=slider(90, 110, initial=95, step=1),
  a2=slider(90, 110, initial=95, step=1),
  s1=slider(2, 8, initial=6, step=0.1),
  s2=slider(2, 8, initial=6, step=0.1),
  x01=slider(10, 40, initial=20, step=1),
  x02=slider(10, 40, initial=20, step=1)
  )

start # valores do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo não restrito: cada variedade com os seus parâmetros

n0 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S[vari])),
          data=db, start=start)
summary(n0)

#-----------------------------------------------------------------------------
# ajuste do modelo restrito: as variades possuem mesmo parâmetro S

start$S <- mean(start$S)
start

n1 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S)),
          data=db, start=start)
summary(n1)

#-----------------------------------------------------------------------------
# testa a hipótese do modelo nulo n1 tendo como alternativa o modelo n0

anova(n1, n0)

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

Obtenção de valores iniciais para o ajuste de duas curvas.

Ainda é necessário fazer a checagem das pressuposições envolvidas no modelo de regressão clássico: homocedasticidade (a variância dos desvios de regressão é constante ao longo a curva), normalidade (os erros aleatórios têm distribuição normal) e independência entre observações (normalmente a independência é garantida por um plano de amostragem adequado). Consulte Análise de resíduos em regressão não linear para saber como obter os 4 gráficos de análise de resíduos.

Esse post foi um aprimoramento de uma resposta na R-help e estimulado por citações na lista. Até a próxima ridícula.

Anúncios

Como fazer legendas em gráficos

Gráficos com legendas e assinatura.

 

Para os gráficos do pacote graphics a adição de legendas é feita por meio da função legend(). Essa função possuí diversas opções para confecção de diversos tipos de legenda (veja os exemplos da documentação).

Nessa ridícula eu fiz alguns exemplos de como colocar legendas em gráficos de funções, gráfico com pontos e retas, gráfico do ajuste de um modelo de regressão, gráfico de barras, legenda fora da área gráfica, escolha da posição da legenda com o mouse e como colocar assinatura nos gráficos. A documentação da função esclarece o significado dos argumentos usados. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# legenda para uma função

curve(x+x^2, 0, 3, lty=1, col=2, lwd=2)
legend("topleft", legend="f(x)", lty=1, col=2, lwd=2, bty="n")

#-----------------------------------------------------------------------------
# legenda para gráfico de pontos ligados

x <- 1:10
y <- x+rnorm(x,0.5)
plot(y~x, type="o", pch=19)
legend("bottomright", legend="valores observados",
       lty=1, col=1, bty="n", pch=19)

#-----------------------------------------------------------------------------
# legenda para pontos e reta de regressão

m0 <- lm(y~x)
plot(y~x)
abline(m0, col=2, lty=2, lwd=2)
legend("top", legend=c("valores observados", "valores ajustados"),
       lty=c(NA,2), col=c(1,2), lwd=1:2, bty="n", pch=c(1,NA))

#-----------------------------------------------------------------------------
# legenda para gráficos de barras

z <- matrix(1:4, 2, 2)
colnames(z) <- c("A","B")
barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend("topleft", legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")

#-----------------------------------------------------------------------------
# legenda fora da área gráfica disposto em colunas

barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend(x=1, y=4.5, xpd=TRUE, ncol=2, legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")

#-----------------------------------------------------------------------------
# escolhendo a posição da legenda com o mouse

barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend(locator(1), xpd=TRUE, ncol=2, legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")

#-----------------------------------------------------------------------------
# nota de rodapé no gráfico com assinatura, data e versão

barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend("topleft", xpd=TRUE, ncol=2, legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")
mtext(side=4, at=par("usr")[3], adj=0,
      cex=0.8, col="gray40", line=-1,
      text=paste("Walmes Zeviani --",
        format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
        R.version.string))

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

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

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.