Exemplo: o modelo de van Genutchen
Este exemplo mostra o ajuste de um modelo não linear. Primeiro discutimos como efetuar um único ajuste para um conjunto de dados e algumas sugestões para examinar resultados. Ao final mostramos como efetuar vários ajustes de uma só vez de forma eficiente e extrair alguns resultados de particular interesse.
O exemplo mostrado aqui foi motivado por um questão levantada pelo Prof. Álvaro Pires da Silva do Departamento de Ciência do Solo da ESALQ/USP e refere-se ao ajuste da equação de van Genutchen para a curva de retenção de água no solo (ou curva de retenção de água no solo).
Informalmente falando, a equação de van Genutchen é um dos modelos matemáticos utilizados para descrever a curva característica de água no solo que caracteriza a armazenagem de água através de relação entre a umidade e o potencial matricial. Para determinação da curva característica de água o procedimento usual é o de se tomar uma amostra que é submetida a diferentes tensões em condições de laboratório. Para cada tensão aplicada a amostra perde parte do conteúdo de água e mede-se a umidade residual na amostra. A partir dos pares pontos com valores medidos de tensão e umidade, obtem-se a curva de retenção de água no solo que descreve a variação da umidade em função dos valores de tensão. O modelo de van Genutchen é dado pela seguinte equação:
\[ \theta = \theta_{R} + (\theta_{S} - \theta_{R}) \left[ \frac{1}{1 + (\alpha \Psi_{m})^{n}} \right]^{1 - (1/n)} \]
em que \(\Psi_{m}\) é o potencial matricial aplicado à amostra e \(\theta\) é a umidade volumétrica medida na amostra. O parâmetros desconhecidos do modelo modelo são \(\theta_{S}\) e \(\theta_{R}\) que correspondem à umidade volumétrica na saturação e residual, respectivamente, \(\alpha\) e \(n\) que definem o formato da curva sendo que o primeiro representa o inverso do potencial de entrada de ar e o segundo é um índice da distribuição dos tamanhos de poros. Portanto são obtidos dados para os pares de pontos \((\Psi_{m}, \theta)\) e \((\theta_{S}, \theta_{R}, \alpha, n)\) são parâmetros desconhecidos a serem estimados e que caracterizam a curva de retenção.
Para exemplificar o ajuste utilizamos dados cedidos pelo Prof. Álvaro que podem ser obtidos usando o comando mostrado a seguir. Este conjunto de dados refere-se a apenas duas amostras que são um subconjunto dos dados originais que contém diversas amostras. O objetivo é determinar da curva de retenção de água no solo estimada segundo modelo de van Genutchen para cada uma das amostras. No objeto \(\texttt{cra}\) a primeira coluna (\(\texttt{am}\)) indica o número da amostra, a segunda (\(\texttt{pot}\)) o potencial aplicado e a terceira (\(\texttt{u}\)) a umidade do solo. Vemos a seguir que dispomos de 15 pontos medidos da curva de retenção da primeira amostra e 13 para a segunda.
cra <- read.table("http://www.leg.ufpr.br/~paulojus/dados/cra.csv",
header = TRUE,
sep = ",")
head(cra) am pot u
1 30 10 0.3071
2 30 19 0.2931
3 30 30 0.2828
4 30 45 0.2753
5 30 63 0.2681
6 30 64 0.2628cra <- transform(cra,
am = as.factor(am))
summary(cra) am pot u
30:15 Min. : 10.0 Min. :0.0636
41:13 1st Qu.: 58.5 1st Qu.:0.1199
Median : 107.5 Median :0.1969
Mean : 2139.8 Mean :0.1879
3rd Qu.: 1550.0 3rd Qu.:0.2436
Max. :26300.0 Max. :0.3071 Inicialmente vamos nos concentrar na discussão do ajuste do modelo e para isto, vamos isolar os dados referentes a uma única amostra.
(cra30 <- subset(cra,
am == 30)) am pot u
1 30 10 0.3071
2 30 19 0.2931
3 30 30 0.2828
4 30 45 0.2753
5 30 63 0.2681
6 30 64 0.2628
7 30 75 0.2522
8 30 89 0.2404
9 30 105 0.2272
10 30 138 0.2120
11 30 490 0.1655
12 30 3000 0.1468
13 30 4100 0.1205
14 30 5000 0.1013
15 30 26300 0.0730No gráfico à esquerda visualizamos os dados de umidade versus pressão aplicada na amostra.
Uma melhor visualização é obtida utilizando-se no eixo horizontal o logarítmo (base 10) dos valores das pressões aplicadas conforme mostrado no gráfico à direita.
g1 <- xyplot(u ~ pot,
type = c("p", "g"),
pch = 16,
xlab = expression(Psi[m]),
ylab = list(label = expression(theta),
rot = 0),
ylim = c(0, .35),
scales = list(x = list(rot = 30)),
cra30)
g2 <- xyplot(u ~ log10(pot),
type = c("p", "g"),
pch = 16,
xlab = expression(log[10](Psi[m])),
ylab = list(label = expression(theta),
rot = 0),
ylim = c(0, .35),
scales = list(x = list(tick.number = 6,
rot = 30)),
cra30)
print(g1,
position = c(0, 0,
.5, 1),
more = TRUE)
print(g2,
position = c(.5, 0,
1, 1))
Portanto, os dados nas colunas \(\texttt{u}\) e \(\texttt{pot}\) do objeto de dados correspondem à \(\theta\) e \(\Psi_{m}\), e as demais quantidades \((\theta_{R},\theta_{S}, n, \alpha)\) são parâmetros (coeficientes) a serem estimados a partir do ajuste do modelo teórico aos dados. Este é um modelo não linear e pode ser ajustado utilizando o método de mínimos quadrados conforme implementado na \(\texttt{nls()}\).
A função possui três argumentos obrigatórios: (i) o primeiro é utilizado para declarar a expressão do modelo a ser ajustado, (ii) o segundo informa o objeto contendo o conjunto de dados cujas nomes das colunas relevantes devem ter o mesmo nome utilizado na declaração do modelo e, (iii) valores iniciais para os parâmetros a serem ajustados que devem ser passados por uma named list, isto é, uma lista com nomes dos elementos, e estes nomes também devem coincidir com os utilizados na declaração do modelo. Há argumentos adicionais para controlar o comportamento algorítimo, tal como critério de convergência. A documentação de \(\texttt{nls()}\) fornece mais detalhes.
A escolha dos valores iniciais é crucial e pode influenciar nos resultados do ajuste utilizando métodos numéricos, especialmente em exemplos como este com um pequeno número de dados. Os valores iniciais para \(\theta_{S}\) e \(\theta_{R}\) foram escolhidos inspecionando-se o gráfico e considerando a interpretação destes como valores de saturação e residual de umidade, portanto, considerando-se máximos e mínimos assintóticos para a função. A escolha de valores iniciais para os demais parâmetros é menos óbvia. Uma das formas de se obter tais valores é efetuar um ajuste aproximado, visual por tentativa e erro, traçando-se curvas sobre o gráfico dos dados. O comando a seguir ilustra como fazer tal procedimento a partir do gráfico dos dados originais mostrado anteriormente definindo uma expressão com o modelo de van Genutchen. Os valores foram escolhidos após uma séria de tentativas.
xyplot(u ~ pot,
type = c("p", "g"),
pch = 16,
xlab = expression(Psi[m]),
ylab = list(label = expression(theta),
rot = 0),
ylim = c(0, .35),
scales = list(x = list(rot = 30)),
cra30) +
layer(panel.curve(.05 + (.35 - .05) / ( (1 + (.1 * x) ^ 1.3) ^ (1 - 1 / 1.3) ),
lty = 2))
Definidos os valores iniciais prossegue-se com o ajuste do modelo conforme os comandos a seguir.
fit30 <- nls(u ~ ur + (us - ur) / ( (1 + (alpha * pot) ^ n) ^ (1 - 1 / n) ),
start = list(us = .35,
ur = .05,
alpha = .1,
n = 1.3),
data = cra30)
summary(fit30)
Formula: u ~ ur + (us - ur)/((1 + (alpha * pot)^n)^(1 - 1/n))
Parameters:
Estimate Std. Error t value Pr(>|t|)
us 0.324121 0.017744 18.27 1.41e-09 ***
ur 0.007082 0.071084 0.10 0.922
alpha 0.038780 0.026202 1.48 0.167
n 1.211816 0.105207 11.52 1.77e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01104 on 11 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 7.004e-06A partir do modelo ajustado pode-se calcular quantidades de interesse. Neste particular exemplo calculamos uma quantidade de interesse prática denotada por \(S\) que é um indicador da qualidade física do solo. Quanto maior o valor de \(S\), melhor a sua qualidade física.
(S <- with(as.list(coefficients(fit30)),
abs((- n * (us - ur) * ( ( (2 * n - 1) / (n - 1) ) ^ (1 / n - 2) ) ))))[1] 0.04097126Os valores preditos são obtidos de forma direta com \(\texttt{fitted(fit30)}\) ou \(\texttt{predict(fit30)}\). Para visualização e avaliação do modelo ajustado podemos fazer diferentes gráficos. A Figura abaixo mostra os pontos ajustados no gráfico da esquerda, e a união destes pontos no gráfico da direita. Gráficos de resíduos semelhantes aos obtidos para avaliar ajuste de modelos lineares podem e devem também ser investivados em uma análise. Neste exemplo mostramos o qq-plot dos resíduos e o gráfico dos resíduos versus valores preditos.
g1 <- xyplot(u ~ log10(pot),
type = c("p", "g"),
pch = 16,
cex = .8,
xlab = expression(log[10](Psi)),
ylab = list(label = expression(theta(Umidade, g/g)),
rot = 0),
scale = list(x = list(tick.number = 6)),
cra30,
key = list(corner = c(1, .97),
text = list(c("observado",
"ajustado")),
points = list(col = c("dodgerblue2", 2),
cex = .8,
pch = c(16, 3)))) +
as.layer(xyplot(fitted(fit30) ~ log10(pot),
col = 2,
pch = 3,
cex = .8,
cra30))
g2 <- xyplot(u ~ log10(pot),
type = c("p", "g"),
pch = 16,
cex = .8,
xlab = expression(log[10](Psi[m])),
ylab = list(label = expression(theta(Umidade, g/g)),
rot = 0),
scale = list(x = list(tick.number = 6)),
cra30,
key = list(corner = c(1, .97),
text = list(c("observado",
"ajustado")),
points = list(col = c("dodgerblue2", 2),
cex = .8,
pch = c(16, 3)))) +
as.layer(xyplot(fitted(fit30) ~ log10(pot),
type = "b",
col = 2,
pch = 3,
cex = .8,
cra30))
rs <- residuals(fit30)
g3 <- qqmath(rs,
type = c("p", "g"),
pch = 16,
xlab = "Quantis teóricos",
ylab = list(label = "Quantis \n amostrais",
rot = 0)) +
layer(panel.qqmathline(rs,
col = "dodgerblue2"))
g4 <- xyplot(rs ~ fitted(fit30),
type = c("p", "g"),
pch = 16,
cex = .8,
xlab = "Valores ajustados",
ylab = list(label = "Resíduos",
rot = 0)) +
layer(panel.abline(h = 0,
col = "navy",
lwd = 1.5))
print(g1,
position = c(0, .5,
.5, 1),
more = TRUE)
print(g2,
position = c(.5, .5,
1, 1),
more = TRUE)
print(g3,
position = c(0, 0,
.5, .5),
more = TRUE)
print(g4,
position = c(.5, 0,
1, .5))
Para obter uma melhor visualização do modelo ajustado pode-se obter valores na curva ajustada não apenas nos pontos observados, mas em uma sequência de valores ao longo do gráfico como ilustrado a seguir. A Figura à esquerda mostra o modelo definido pelos valores iniciais e o modelo ajustado na escala original. Note que neste exemplo em geral prefere-se a visualização na escala logarítmica do potencial conforme o gráfico da direita. A curva com o modelo ajustado a serem desenhadas sobre o gráfico dos dados são obtidas com comandos a seguir.
pp <- 10 ^ seq(1, 4.5,
length.out = 201)
g1 <- xyplot(u ~ pot,
type = c("p", "g"),
pch = 16,
xlab = expression(Psi[m]),
ylab = list(label = expression(theta),
rot = 0),
ylim = c(0, .35),
scales = list(x = list(rot = 30)),
key = list(corner = c(1, .97),
text = list(c("valores iniciais",
"valores ajustados")),
lines = list(col = c(1, "dodgerblue2"),
lty = 2:1)),
cra30) +
layer(panel.curve(.05 + (.35 - .05) / ( (1 + (.1 * x) ^ 1.3) ^ (1 - 1 / 1.3) ),
lty = 2)) +
layer(panel.lines(pp, predict(fit30, list(pot = pp))))
g2 <- xyplot(u ~ log10(pot),
type = c("p", "g"),
pch = 16,
cex = .8,
xlab = expression(log[10](Psi[m])),
ylab = list(label = expression(theta),
rot = 0),
scale = list(x = list(tick.number = 6,
rot = 30)),
cra30) +
layer(panel.lines(log10(pp), predict(fit30, list(pot = pp))))
print(g1,
position = c(0, 0,
.5, 1),
more = TRUE)
print(g2,
position = c(.5, 0,
1, 1))
Comentários: é importante lembrar que certos modelos não lineares são parcialmente linearizáveis e neste caso o ajuste pode ser mais preciso e numericamente estável se beneficiando disto para reduzir a dimensão do problema de otimização numérica. Para isto é necessário redefinir a especicifação do modelo e utilizar o argumento \(\texttt{method= " plinear"}\) na \(\texttt{nls()}\). Neste exemplo em particilar pode-se considerar fazer o ajuste na escala de \(log_{10}(\Psi_{m})\) já que os resultados são tipicamente visualizados desta forma. Isto reduz a escala dos valores das variáveis e também torna o problema mais estável numericamente. Por outro lado, em geral reparametrizações podem mudar a interpretação de alguns parâmetros de modelo. Finalmente cuidados usuais com ajuste de modelos utilizando métodos iterativos devem ser observados, tais como sensibilidade a valores iniciais e verificação de convergência do algorítmo numérico.
