Página Inicial > gráficos, regressão > Método gráfico para valores iniciais em regressão não linear com 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.

About these ads
  1. Nenhum comentário ainda.
  1. No trackbacks yet.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s

Seguir

Obtenha todo post novo entregue na sua caixa de entrada.

Junte-se a 52 outros seguidores

%d blogueiros gostam disto: