Arquivo

Posts Tagged ‘RStudio’

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

RStudio e a função manipulate()

Manipulando ângulo de observação e gradiente de cores de um gráfico tridimensional com as ferramentas do RStudio.

O RStudio é um dos vários disponíveis editores para script R. Uma das suas peculiaridades é o pacote manipulate. Esse pacote é parte integrante do RStudio e fica disponível apenas com a sua instalação.

A função manipulate::manipulate() permite que você use os gráficos do R de forma interativa. Ou seja, através de deslizadores (slider), seletores (picker) e opções (checkbox) você é capaz de mudar aspectos gráficos. A documentação da função traz exemplos interessantes. Nessa ridícula eu implementei dois exemplos. O primeiro é para escolher o ângulo de observação de uma superfície de resposta (inclui escolha do gradiente de cores). O segundo permite que você entenda como funciona a estimação kernel de densidade.

Com a manipulate() fica muito mais simples construir funções para obter valores iniciais em modelos de regressão não linear. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados para um gráfico tridimensional

da <- expand.grid(x=seq(0,10,l=30), z=seq(0,10,l=30))
da$y <- with(da, x+z+0.2*x*z) # gera dados

#-----------------------------------------------------------------------------
# escolher o ângulo de observação e o esquema de cores
 
manipulate({
             ## faz o gráfico tridimensional
             colr <- colorRampPalette(c(c1, c2, c3), space="rgb")
             arrows <- arr
             wireframe(y~x+z, data=da, scales=list(arrows=arrows),
                       col="gray30", col.contour="gray30",
                       col.regions=colr(100),  drape=TRUE,
                       screen=list(z=z.angle, x=x.angle)
                       )
           },
           ## controla o valor dos angulos e das cores
           z.angle=slider(0, 360, step=10, initial=40),
           x.angle=slider(-180, 0, step=5, initial=-60),
           arr=checkbox(FALSE, "show.arrows"),
           c1=picker("red","yellow","orange","green","blue","pink","violet"),
           c2=picker("red","yellow","orange","green","blue","pink","violet"),
           c3=picker("red","yellow","orange","green","blue","pink","violet")
           )

#-----------------------------------------------------------------------------
# gráfico de densidade controlando o bandwidth e tipo de função kernel

x <- rgamma(300, 3, 7)

manipulate({
             plot(density(x, bw=bw, kernel=kernel))
             if(show.rug==TRUE) rug(x)
           },
           kernel=picker("gaussian", "epanechnikov", "rectangular",
             "triangular", "biweight","cosine",
             "optcosine"),
           bw=slider(0.01, 0.15, step=0.003, initial=0.05),
           show.rug=checkbox(TRUE, "show rug")
           )

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

Realçador de código R para blogs

Quem trabalha com programação reconhece que a organização do código facilita muito a releitura e checagem. O R, uma vez que é uma linguagem, possui editores que apresentam recursos para organização do código, como Emacs, RStudio, Tinn-R, Rkward entre diversos. A maioria dos editores apresenta um padrão de identação do código e realçador (highlighter) de texto. A necessidade de divulgação de código R, de forma organizada e atraente, fez com que os realçadores aparecessem nas páginas de blogs, wikis, etc. O wordpress.com possui um realçador de código R nativo. Para quem não tem wordpress.com, pode usar o Pretty R syntax highlighter que possui, além dos realces usuais, links para as páginas de documentação das funções. Abaixo está código R com o realçador do wordpress.com e com o pretty-R. Perceba que as funções (negrito azul) possuem links para página de ajuda.

# importa dados
soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)
 
# ajusta um modelo e pede anova
m1 <- aov(rengrao~bloco+agua*potassio, soja)
anova(m1)
 
# cria uma lista com as variáveis resposta
respostas <- do.call(c, apply(soja[,4:7], 2, list))
do.call(c, respostas)
 
# faz o ajuste para todas as respostas
ajustes <- lapply(respostas,
                  function(r){
                    m0 <- aov(r~bloco+agua*potassio, data=soja)
                    m0
                  })
 
# pede todas as anovas
lapply(ajustes, anova)
# importa dados
soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)
 
# ajusta um modelo e pede anova
m1 <- aov(rengrao~bloco+agua*potassio, soja)
anova(m1)
 
# cria uma lista com as variáveis resposta
respostas <- do.call(c, apply(soja[,4:7], 2, list))
do.call(c, respostas)
 
# faz o ajuste para todas as respostas
ajustes <- lapply(respostas,
                  function(r){
                    m0 <- aov(r~bloco+agua*potassio, data=soja)
                    m0
                  })
 
# pede todas as anovas
lapply(ajustes, anova)

Created by Pretty R at inside-R.org

Além de no R você fazer praticamente tudo, você pode apresentar de forma elegante. Até a próxima ridícula.