Regressão segmentada ou piecewise.

Depois de um tempo sumido, espero voltar finalmente a postar aqui. Estou devagar para estudar ultimamente, é de dezembro meu último post aqui, mas agora que acabou minha faculdade, vou voltar a estudar novas coisas e postar regularmente, existe muito coisa que pensei em postar e parei.

Mas para retomar as atividades, vamos falar sobre a piecewise regression, ou regressão segmentada. A gente ta acostumado a modelos de diferenças de médias, tipo anova, e modelos de regressão linear, ou a combinação desses dois. Eles resolvem bem quase todos os casos, mas as vezes a gente já sabe um pouco mais sobre o sistema que estamos estudando, e podemos fazer melhor, utilizando esse conhecimento.

Por exemplo, sabemos que plantas muitas vezes são limitadas, no seu crescimento, pelo mínimo de alguns nutrientes, a lei de Liebig. Podemos expandir esse conceito e pensar assim, enquanto não há um mínimo de nutrientes no solo, o crescimento é estagnado, mas a partir de algum ponto começamos a ter um incremento no crescimento de uma espécie dado um acréscimo nos nutrientes.

Podemos simular alguns dados para essa situação, então imaginando que temos uma medida dos nutrientes num vaso e do tamanho da planta após um determinado tempo fixo, podemos visualizar esses dados da seguinte forma.

Veja então que no eixo x temos a concentração de nutrientes, e no eixo y temos o tamanho das plantulas, então enquanto não temos uma certa concentração parece que todo mundo tem o mesmo tamanho, mas a partir de uma concentração, quanto maior a concentração, maior ficam as plantulas. Se a gente pensar, isso pode ocorrer em um bocado de situações, o que interessa é que, a partir de um ponto, as coisas mudam, algo começa a ocorrer, ou para, aqui temos um ponto crítico, algo entre o 4 e o 6 em que abaixo dele, a concentração de nutrientes é irrelevante.

Podemos representar isso matematicamente como uma equação do tipo piecewise, nesse caso, o que usamos para simular os dados foi:

  f(x) = \left\{          \begin{array}{ll}              8 & \quad x\ \textless \  5 \\              1+2x & \quad x \geq 5          \end{array}      \right.

Sendo x o nosso preditor, concentração de nutriente e f(x) o tamanho, para a parte determinística do modelo.
Esses são os valores que usamos para gerar essas dados de exemplo, que no caso é partido no 5, antes do 5 temos uma reta f(x)=8+0x(eu omiti o 0x, porque tudo que é multiplicado por zero da zero, mas temos uma reta, só que sem inclinação, ou melhor, com inclinação zero) e depois do 5 temos uma reta f(x)=1+2x.

Nesse caso, o partição, o valor 5 tem um significa biológico importante, ele é o limite mínimo de nutrientes, para começarmos a observar alguma coisa, antes dele nada acontece, mas como podemos recuperar esse valor dos dados? A primeira opção é usar um tipo de força bruta, testar o modelo acima.

Para aquela equação podemos ajustar ela usando uma variável tipo dummy usando a função lm, ou usar a função nls e ajustar diretamente a equação piecewise usando um ifelse, aqui vamos usar diretamente a equação com nls que é mais direto a interpretação dos parâmetros(uma reta é a+bx, então vamos ter dois a, o a1 e a2 e dois b, o b1 e b2 para as duas retas), ajustamos duas retas e um valor de quebra, no caso ajustamos o modelo para muitas valores de quebra e vemos qual da o melhor ajuste usando logLikelihood

Então criamos um vetor de possíveis valores de quebra.

1
2
quebras <- seq(1,9,0.5)
minimo_quadrado <- vector(length = length(quebras))

E testamos todos eles

1
2
3
4
for(i in 1:length(quebras)){
 piecewise <- nls(tamanho ~ ifelse(nutriente < quebras[i],a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
 minimo_quadrado[i] <- logLik(piecewise)
}

E podemos avaliar graficamente o resultado.

Certo, como esperávamos, o valor da simulação é recuperado, veja que o pico de ajusta está bem no 5.

Podemos então ajustar o modelo usando 5 como ponto de quebra e ver se recuperamos os coeficientes corretamente

> modelo<-nls(tamanho ~ ifelse(nutriente < 5,a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1)) > summary(modelo) Formula: tamanho ~ ifelse(nutriente < 5, a1 + b1 * nutriente, a2 + b2 * nutriente) Parameters: Estimate Std. Error t value Pr(>|t|) a1 7.7307 0.3410 22.668 <2e-16 *** a2 3.2449 1.1431 2.839 0.0063 ** b1 0.2193 0.1249 1.756 0.0846 . b2 1.7430 0.1369 12.734 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9221 on 56 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 1.57e-08

E temos as inclinações, zero para a primeira equação esta no intervalo do b1, que é o nome que demos para o parâmetro, e assim como o 2 está no intervalo do b2, os interceptos também são recuperados. Vamos ver ele graficamente, e além disso, vamos comparar com um modelo linear.

Veja que ele fica bem mais certinho, mas o modelo linear também detectaria uma tendência aqui, porém não nos forneceria essa informação, de a partir de onde, de qual concentração é relevante, seriamos enganados pensando que a concentração de nutrientes é relevante em qualquer intensidade quando ela so importa em uma parte do preditor, isso mesmo avaliando os resíduos dos modelos.

Temos um ajuste menor (lembre-se que estamos falando de números negativos, então -78 é maior que -111), mas não temos uma tendência clara nos resíduos para esse caso, rejeitar o modelo linear, somente observando os resíduos, então aqui interessa muito o que pensamos a priori da construção do modelo.

Como podemos ter casos onde a busca precisa ser mais precisa para o ponto de quebra, temos funções prontas que vão cuidar dessas atividades, por exemplo o pacote segmented, que ajusta esse tipo de modelo diretamente, fazendo um update do modelo linear.

Basicamente temos o seguinte resultado:

> modelo_pieciwise<- segmented(modelo_linear, seg.Z = ~nutriente, psi=1) > summary(modelo_pieciwise) ***Regression Model with Segmented Relationship(s)*** Call: segmented.lm(obj = modelo_linear, seg.Z = ~nutriente, psi = 1) Estimated Break-Point(s): Est. St.Err 3.799 0.304 Meaningful coefficients of the linear terms: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.8733 0.4533 17.369 <2e-16 *** nutriente 0.1274 0.1967 0.648 0.52 U1.nutriente 1.9158 0.2222 8.621 NA --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.045 on 56 degrees of freedom Multiple R-Squared: 0.9576, Adjusted R-squared: 0.9554 Convergence attained in 7 iterations with relative change 3.485019e-16

Veja que temos o Estimated Break-Point, que da 3.799, um pouco inferior ao que realmente usamos, e temos apenas um intercepto, mas a primeira inclinação fica em zero, e a segunda em 0.12+1.91, diferença de inclinação para depois do ponto de quebra, que é 2, então ele estima bem os parâmetros, mas é um modelo um pouco diferente do que usamos (só um intercepto), eu não li muito sobre esse pacote, para saber ajustar exatamente como fizemos ali em cima, mas é um pacote a ser explorado, para quem quer utilizar esse tipo de modelo e o sistema de plot dele é bem simples.

Bem é isso ai, depois de um tempo enroscado, de volta a atividade, o script vai estar lá no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Michael J. Crawley 2012 The R Book, 2nd ed. Wiley 1076pp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
set.seed(31415)
 
nutriente <- runif(60,0,10)
 
fn_bs <- function(x){
    saida<-vector(length=length(x))
    for(i in 1:length(x)){       
        if(x[i]<5){
            saida[i] <- 8+0*x[i]
        }else{
            saida[i] <- 1+2*x[i]
        }
    }
    saida <- saida+rnorm(length(x),0,1)
    return(saida)
}
 
tamanho <- fn_bs(nutriente)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
 
#####################################
## Força bruta
#####################################
quebras <- seq(1,9,0.5)
minimo_quadrado <- vector(length = length(quebras))
 
for(i in 1:length(quebras)){
 piecewise <- nls(tamanho ~ ifelse(nutriente < quebras[i],a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
 minimo_quadrado[i] <- logLik(piecewise)
}
 
minimo_quadrado <- unlist(minimo_quadrado)
 
plot(quebras,minimo_quadrado,type="b",pch=19,xlab="Valor de quebra",ylab="LogLikelihood",frame=F,xlim=c(0,10))
 
modelo<-nls(tamanho ~ ifelse(nutriente < 5,a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
summary(modelo)
 
modelo_linear<- lm(tamanho~nutriente)
summary(modelo_linear)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
 
abline(modelo_linear,lwd=3,lty=3,col="blue")
 
curve(coef(modelo)[1]+coef(modelo)[3]*x,0,5,add=T,lwd=3,lty=3,col="red")
curve(coef(modelo)[2]+coef(modelo)[4]*x,5,12,add=T,lwd=3,lty=3,col="red")
abline(v=5,lwd=1,lty=2,col="red")
legend("topleft",lwd=3,lty=3,col=c("blue","red"),legend=c("Modelo Linear","Modelo piecewise"),bty="n")
 
par(mfrow=c(2,1))
plot(resid(modelo_linear),main=paste("Loglikelihood do modelo liner =",round(logLik(modelo_linear),3)),frame=F,xlab="Amostras",ylab="Resíduo",ylim=c(-5,5))
abline(h=0,lty=2,lwd=2)
plot(resid(modelo),main=paste("Loglikelihood do modelo piecewise =",round(logLik(modelo),3)),frame=F,xlab="Amostras",ylab="Resíduo",ylim=c(-5,5))
abline(h=0,lty=2,lwd=2)
 
#####################################
## Pacote Segmented
#####################################
##install.packages("segmented")
library(segmented)
 
modelo_pieciwise<- segmented(modelo_linear, seg.Z = ~nutriente, psi=1)
summary(modelo_pieciwise)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
plot(modelo_pieciwise,add=T,lwd=2,lty=2,col="red")

Leave a Reply

Your email address will not be published. Required fields are marked *