Modelos não lineares permitem o ajuste de relações entre quantidades de interesse mais complexas que relações lineares ou linearizáveis. Em muitos casos tais modelos têm a sua forma funcional específica para o problema sendo tratado, relacionada a algum mecanismo (biológico, físico, etc) inerente ao processo em questão.
Nesta seção vamos ilustrar o ajuste de modelos não lineares com dados da área de física de solos. Utilizaremos a função nls(), cujo nome é um acrônimo para non-linear least squares. Esta função é bastante flexível e incorpora diversas opções para fazer ajustes incluindo características do modelo, tipo de ajuste e controle dos algoritmos disponíveis.
Diferentemente dos modelos lineares, o ajuste de modelos não lineares não permite que as expressões dos estimadores dos parâmetros desconhecidos do modelo sejam obtidas analiticamente, sendo portanto necessário o uso de métodos numéricos. Inicialmente mostramos um ajuste feito de forma ingênua (naïve), declarando apenas a função e valores iniciais. Tal procedimento, embora simples, pode ser ineficiente para o uso de métodos numéricos. Entretanto, o ajuste com nls() pode incorporar procedimentos que tendem a aprimorar o comportamento dos métodos numéricos tais como o fornecimento de funções que informem sobre a derivada do modelo sendo ajustado, inicialização automática com valores iniciais obtidos automaticamente, linearização parcial do modelo, além da escolha e calibragem dos algoritmos. O objetivo aqui não é o de investigar todas estas opções, mas apenas fornecer os elementos iniciais para ilustrar a possibilidade de se obter tais resultados usando o R.
Inicialmente discutiremos como efetuar um único ajuste para um conjunto de dados e algumas sugestões para examinar resultados. Mais adiante 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 uma 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 Genuchten para a curva de retenção de água no solo (ou curva característica de água no solo).
Informalmente falando, a equação de van Genuchten é um dos modelos matemáticos utilizados para descrever a curva característica de água no solo. O modelo 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 de pontos com valores medidos de tensão e umidade, obtém-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 Genuchten é dado pela seguinte equação: \begin {equation} \label {eq:vanG} \theta = \theta _R + (\theta _S-\theta _R)\left [\frac {1}{1+(\alpha \Psi _m)^n}\right ]^{1-(1/n)} \end {equation} em que \(\Psi _m\) é o potencial matricial aplicado à amostra e \(\theta \) é a umidade volumétrica medida na amostra. Os parâmetros desconhecidos do 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 do conjunto de dados original que contém diversas amostras. O objetivo é determinar a curva de retenção de água no solo estimada segundo modelo de van Genuchten para cada uma das amostras. No objeto cra a primeira coluna (am) é a identificação da amostra, a segunda (pot) é o potencial aplicado e a terceira (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", head=T, 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.2628
> cra <- 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 selecionar os dados referentes a uma única amostra.
> cra30 <- subset(cra, am==30) > cra30
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.0730
No gráfico à esquerda da Figura 30.1 visualizamos os dados de umidade versus pressão aplicada na amostra.
> with(cra30, plot(u ~ pot, xlab=expression(Psi[m]), ylab=expression(theta), ylim=c(0,0.35)))
Uma melhor visualização é obtida utilizando-se no eixo horizontal o logaritmo (base 10) dos valores das pressões aplicadas conforme mostrado no gráfico à direita.
> with(cra30, plot(u ~ log10(pot), xlab=expression(log[10](Psi[m])), ylab=expression(theta), ylim=c(0,0.35)))
Portanto, os dados nas colunas u e pot do objeto de dados correspondem, respectivamente, à \(\theta \) e \(\Psi _m\) na equação 3. As demais quantidades \((\theta _S, \theta _R, \alpha , n)\) são parâmetros (coeficientes) a serem estimados a partir do ajuste do modelo teórico aos dados. Este é um modelo não linear que pode ser ajustado utilizando o método de mínimos quadrados conforme implementado em 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 cujos 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, uma lista com nomes dos elementos que devem coincidir com os utilizados na declaração do modelo. Há argumentos adicionais para controlar o comportamento do algoritmo, tal como critério de convergência. A documentação de nls() fornece mais detalhes.
A escolha dos valores iniciais é crucial e pode influenciar os 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 para curve() com o modelo de van Genuchten. Os valores foram escolhidos após uma série de tentativas.
> curve(0.05 + (0.35-0.05)/((1+(0.1*x)^1.3)^(1-1/1.3)), from=0, to=27000, add=T, 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)), data=cra30, + start=list(us=0.35,ur=0.05,alpha=0.1, n=1.3)) > summary(fit30)
Formula: u ~ ur + (us - ur)/((1 + (alpha * pot)^n)^(1 - 1/n)) Parameters: Estimate Std. Error t value us 0.324121 0.017744 18.27 ur 0.007082 0.071084 0.10 alpha 0.038780 0.026202 1.48 n 1.211816 0.105207 11.52 Residual standard error: 0.01104 on 11 degrees of freedom Number of iterations to convergence: 8 Achieved convergence tolerance: 7.748e-06
A partir do modelo ajustado pode-se calcular quantidades de interesse. Neste particular exemplo calculamos uma quantidade de interesse prático 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(coef(fit30)),abs((-n*(us-ur)*(((2*n-1)/(n-1))^(1/n-2))))) > S
[1] 0.04097127
Os valores preditos são obtidos de forma direta com fitted(fit30) ou predict(fit30). Para visualização e avaliação do modelo ajustado podemos fazer diferentes gráficos. A Figura 30.1 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 investigados em uma análise. Neste exemplo mostramos o qq-plot dos resíduos e o gráfico dos resíduos versus valores preditos.
> with(cra30, plot(log10(pot),u,xlab=expression(log[10](Psi)), + ylab=expression(theta(Umidade,g/g)))) > with(cra30, points(log10(pot),fitted(fit30), pch=3, col="red")) > legend(3, 0.30, c("observado","ajustado"),pch=c(1,3), col=c(1,2)) > with(cra30, plot(log10(pot),u,xlab=expression(log[10](Psi[m])), + ylab=expression(theta(Umidade,g/g)))) > with(cra30, points(log10(pot),fitted(fit30), type="b", pch="+", col="red")) > legend(3, 0.30, c("observado","ajustado"),pch=c(1,3), col=c(1,2)) > rs <- resid(fit30); qqnorm(rs); qqline(rs) > plot(fitted(fit30),resid(fit30)); abline(h=0, col="navy")
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 30.1 mostra à esquerda 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 gráfico da direita. As curvas com os modelos ajustados a serem desenhadas sobre o gráfico dos dados são obtidas com os comandos a seguir.
> pp <- 10^seq(1, 4.5, l=201) > lines(pp, predict(fit30, list(pot=pp))) > legend("topright", c("valores iniciais","valores ajustados"), lty=2: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 especificação do modelo e utilizar o argumento method="plinear" em nls(). Neste exemplo em particular 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 algoritmo numérico.
Vamos considerar uma situação comum na prática em que tem-se várias amostras para as quais deseja-se fazer ajustes individuais como o ilustrado anteriormente. É portanto conveniente que isto seja feito de forma automática, sem a necessidade de repetir os passos acima a cada ajuste. Neste exemplo vamos considerar duas amostras, mas o procedimento demonstrado a seguir é geral e funcionará igualmente para um maior número de amostras.
Serão mostradas duas soluções. Nesta seção o ajuste é feito para cada amostra individualmente automatizando várias chamadas à função nls() através de lapply() emulando o comportamento das várias chamadas em um loop. Na próxima seção será mostrado como obter todos os ajustes com uma única chamada à nls(). Ilustramos ambos os casos porque a forma mais adequada vai depender da situação em questão e dos objetivos da análise.
Começamos definindo uma função que contém uma chamada à nls() como acima. Estamos incluindo na função um argumento ini para passar valores iniciais. Caso não fornecidos a função assumirá os valores indicados. A seguir utilizamos a função by() para proceder o ajuste para cada amostra individualmente. Esta função retorna uma lista com dois elementos, um para cada amostra, sendo que cada um deles contém o ajuste do modelo não linear.
> fit.vG <- function(x, ini=list(us=0.30,ur=0.02,alpha=0.05,n=1.3)) + nlsfit=nls(u~ur+(us-ur)/(1+(alpha*pot)^n)^(1-1/n), data=x, + start=ini) > > allfits <- by(cra, cra$am, fit.vG) > names(allfits)
[1] "30" "41"
Neste caso, o objeto resultante allfits é uma lista de listas e portanto podemos usar funções como lapply(), sapply() ou similares para extrair resultados de interesse. Note que a primeira retorna sempre uma lista, enquanto a segunda simplifica o objeto resultante se possível. Por exemplo, quando extraindo coeficientes a função retorna uma matriz \(4 \times 2\), já que para cada uma das duas amostras são extraídos quatro coeficientes.
> lapply(allfits, summary)
$`30` Formula: u ~ ur + (us - ur)/(1 + (alpha * pot)^n)^(1 - 1/n) Parameters: Estimate Std. Error t value us 0.324120 0.017744 18.27 ur 0.007082 0.071084 0.10 alpha 0.038780 0.026202 1.48 n 1.211816 0.105207 11.52 Residual standard error: 0.01104 on 11 degrees of freedom Number of iterations to convergence: 6 Achieved convergence tolerance: 8.286e-06 $`41` Formula: u ~ ur + (us - ur)/(1 + (alpha * pot)^n)^(1 - 1/n) Parameters: Estimate Std. Error t value us 0.243148 0.009446 25.741 ur -0.122402 0.171615 -0.713 alpha 0.035928 0.022324 1.609 n 1.113320 0.079473 14.009 Residual standard error: 0.006207 on 9 degrees of freedom Number of iterations to convergence: 7 Achieved convergence tolerance: 9.37e-06
> lapply(allfits, coef)
$`30` us ur alpha n 0.324120351 0.007081962 0.038779951 1.211815750 $`41` us ur alpha n 0.24314783 -0.12240199 0.03592827 1.11332043
> sapply(allfits, coef)
30 41 us 0.324120351 0.24314783 ur 0.007081962 -0.12240199 alpha 0.038779951 0.03592827 n 1.211815750 1.11332043
Quando ajustamos o modelo apenas para uma das amostras mostramos como calcular o índice \(S\) de qualidade física do solo a partir dos coeficientes estimados. Vamos então aqui obter este índice para cada uma das amostras. Para isto simplesmente definimos uma função que recebe o modelo ajustado e usa os coeficientes para calcular o valor de \(S\). Passamos o objeto (lista) contendo todos os ajustes e a função que calcula \(S\) para sapply() que neste caso vai simplificar o resultado para formato de um vetor, já que a função calculaS retorna um escalar para cada amostra.
> calculaS <- function(fit) with(as.list(coef(fit)), + abs((-n*(us-ur)*(((2*n-1)/(n-1))^(1/n-2))))) > Sall <- sapply(allfits, calculaS) > Sall
30 41 0.04097126 0.02950320
Finalmente, para encerrar este exemplo, vamos mostrar uma possível forma de combinar a visualização dos ajustes em um único gráfico. Começamos definindo uma sequência de valores para os quais queremos visualizar os ajustes. Armazenamos os valores preditos para cada amostra no objeto allpred e optamos aqui por mostrar os ajustes para as duas amostras no mesmo gráfico.
> lpsimax <- with(cra, max(log(pot))) > pp <- 10^seq(1, lpsimax, l=501) > allpred <- lapply(allfits, predict, list(pot=pp)) > with(cra, plot(log10(pot), u, type="n", xlab=expression(log[10](Psi)), + ylab=expression(theta(Umidade,g/g)))) > with(cra, text(log10(pot), u, as.character(am))) > lapply(allpred, function(yp) lines(log10(pp), yp))
Na seção anterior obtivemos o ajuste para cada amostra separadamente fazendo várias chamadas à função nls(). Isto pode ser adequado quando deseja-se de fato ajustes individuais e se, por um lado são efetuadas várias chamadas à função, por outro o número de dados em cada uma delas é pequeno. Uma forma alternativa de obter parâmetros para cada amostra, possivelmente mais eficiente que a mostrada anteriormente, é discutida a seguir.
Nesta seção vamos considerar fazer todos os ajustes de uma só vez, isto é em uma única chamada à nls() que portanto vai utilizar todos os dados de todas as amostras. Além do aspecto computacional, isto pode ser interessante pois permite comparar e testar hipóteses para escolha entre diferentes modelos alternativos para explicar os dados. A seguir exemplificamos tal procedimento iniciando com um modelo para cada amostra e comparando com um modelo que assume que os parâmetros (\(\alpha \), \(n\)) são comuns entre as amostras. Neste caso interpreta-se que cada amostra informa sobre os respectivos valores para (\(\theta _S\), \(\theta _R\)) enquanto todas as amostras conjuntamente informam sobre (\(\alpha \), \(n\)). Após ajustar os modelos candidatos podemos fazer uma comparação formal dos ajustes através de anova(), o que não seria possível ajustando os modelos separadamente como mostrado na seção anterior.
Os dois ajustes são mostrados a seguir. O seletor \([ ]\) é usado para indicar que os dados são tratados em grupos definidos por am. No caso do modelo com parâmetros distintos informamos oito valores iniciais para os parâmetros.
> form0 <- formula(u~ur[am]+(us[am]-ur[am])*(1/(1+(alpha[am]*pot)^n[am]))^(1-1/n[am])) > mod0 <- nls(form0, cra, + start = list(us = c(0.3,0.3), ur = c(0,0), + alpha = c(0.04, 0.04), n = c(1.25, 1.25))) > mod0
Nonlinear regression model model: u ~ ur[am] + (us[am] - ur[am]) * (1/(1 + (alpha[am] * pot)^n[am]))^(1 - 1/n[am]) data: cra us1 us2 ur1 ur2 alpha1 alpha2 n1 n2 0.324120 0.243148 0.007085 -0.122402 0.038779 0.035928 1.211819 1.113321 residual sum-of-squares: 0.001688 Number of iterations to convergence: 6 Achieved convergence tolerance: 4.574e-06
Para ajuste assumindo valores comuns para os parâmetros \(\alpha \) e \(n\) não utilizamos o indicador de grupos para estes parâmetros e informamos apenas um valor inicial para cada um deles.
> form1 <- formula(u ~ ur[am] + (us[am] -ur[am])*(1/(1+(alpha*pot)^n))^(1-1/n)) > mod1 <- nls(form1, cra, + start = list(us=c(0.3,0.3), ur = c(0,0), alpha=0.04, n = 1.25)) > mod1
Nonlinear regression model model: u ~ ur[am] + (us[am] - ur[am]) * (1/(1 + (alpha * pot)^n))^(1 - 1/n) data: cra us1 us2 ur1 ur2 alpha n 0.32106 0.24870 -0.03056 -0.02759 0.03994 1.17195 residual sum-of-squares: 0.001846 Number of iterations to convergence: 5 Achieved convergence tolerance: 1.976e-06
Neste exemplo temos então um modelo inicial com oito e outro mais parcimonioso com apenas seis parâmetros e utilizamos um teste formal para orientar a escolha de modelo, que neste caso indica que o modelo mais parcimonioso com parâmetros comuns explica os dados satisfatoriamente.
> anova(mod1,mod0)
Analysis of Variance Table Model 1: u ~ ur[am] + (us[am] - ur[am]) * (1/(1 + (alpha * pot)^n))^(1 - 1/n) Model 2: u ~ ur[am] + (us[am] - ur[am]) * (1/(1 + (alpha[am] * pot)^n[am]))^(1 - 1/n[am]) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 22 0.0018462 2 20 0.0016884 2 0.00015786 0.935 0.4091
Carregando pacotes exigidos: sp