Generalized additive model, quando não sabemos que modelo usar.

Vamos olhar algo que um amigo meu me perguntou sempre sobre, a algum tempo atrás, o GAM, ou Generalized additive model.

Primeiro, vamos ver qual o interesse em outro método estatístico. Sabemos que existem muitos modelos matemáticos, que explicam partes determinísticas de processos biológicos, como crescimento populacional com o modelo do Lotka-Volterra ou o comportamento de populações no espaço com o modelo de metapopulações do Levins.

Normalmente a gente não coleta dados e ajusta um modelo de crescimento populacional neles. A gente usa coisas como análise de variância e regressão linear entrou outros modelos, porque são modelos mais simples, mais robustos que acabamos por detectar padrões nos dados, e depois vamos mais afundo nesses padrões.

Vamos pensar num exemplo aqui. Vamos gerar dados de um preditor, e duas respostas diferentes.

1
2
3
4
5
###Simulando Dados
set.seed(123)
preditor <- runif(30,0,4*pi)
resposta_pot<-5+2*(preditor)^2+rnorm(30,0,10)
resposta_seno<-5+2*sin(preditor)+rnorm(30)

Seja o preditor a quantidade de açúcar no néctar e a resposta a territorialidade na fonte desse néctar. Suponha que temos algo assim:

Bem a gente pode ajustar uma reta ali, usando uma regressão linear e vai dar um resultado legal.

E teremos uma relação significativa.

Call: lm(formula = resposta_pot ~ preditor) Residuals: Min 1Q Median 3Q Max -43.001 -24.431 -2.866 23.330 53.420 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -52.67 11.82 -4.456 0.000123 *** preditor 26.04 1.47 17.718 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 28.98 on 28 degrees of freedom Multiple R-squared: 0.9181, Adjusted R-squared: 0.9152 F-statistic: 313.9 on 1 and 28 DF, p-value: < 2.2e-16

Com um valor de R quadrado bastante alto, e fim de papo? Pode ser que sim, mas podemos fazer melhor. Veja que nesse caso, temos um modelo exponencial.

Apesar da regressão resolve ali localmente, o primeiro problema é o fato de que não podemos extrapolar baseado nesse modelo de regressão, já que ele vai falhar miseravelmente além do intervalo de dados que temos, é so pensar que a curva vai subir muito mais “rápido” que a reta da regressão para imaginar isso, outro problema é que nosso residuo não tem uma distribuição normal, temos um padrão forte nele, o padrão que vem do determinismo dos dados que não está sendo coberto.

Na figura Residual vs Fitted a gente vê do que estou falando, aqui era para nao ter padrão nenhum, no entanto temos uma curva bem fácil de identificar até no olhômetro.

Mas esse é quando a gente ta com sorte, agora imagine se os dados tem um comportamento todo estranho, um padrão, mas que não é algo com cara de reta.

Ai a regressão não vai dar em nada, no caso anterior a gente poderia resolver com uma transformação, mas nesse caso não.

E aqui a robustez não funciona

Call: lm(formula = resposta_seno ~ preditor) Residuals: Min 1Q Median 3Q Max -2.8625 -1.3596 -0.2053 1.5895 3.9992 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.46404 0.77687 7.033 1.19e-07 *** preditor -0.06799 0.09658 -0.704 0.487 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.904 on 28 degrees of freedom Multiple R-squared: 0.01739, Adjusted R-squared: -0.0177 F-statistic: 0.4956 on 1 and 28 DF, p-value: 0.4873

Mas o padrão está la.

Denovo a gente vê que não deveria acreditar na regressão, a primeira figura do Residual vs Fitted, a gente vê que estamos quebrando premissas, não há uma distribuição normal nos resíduos, tem um padrão neles, mas não conseguimos pegar no modelo.

Então é pra esse tipo de coisa que existe o GAM, a gente sabe que existe um padrão determinístico nos dados, mas não sabemos qual é ele, então ou a gente sai chutando vários modelos, ou a gente uma o Generalized additive model. De certa forma é bem simples usar, ele tem uma implementação no R no pacote mgcv, que é bem simples de usar.

a sintaxe é praticamente a mesma de um glm ou lm, so adicionando aquela função s(), então o ajuste fica assim.

> modelo_gam<- gam(resposta_pot~s(preditor,k=5,fx=T),family = gaussian) > summary(modelo_gam) Family: gaussian Link function: identity Formula: resposta_pot ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 134.639 1.618 83.22 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 908.1 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.992 Deviance explained = 99.3% GCV = 94.22 Scale est. = 78.516 n = 30

E temos um sumario bem parecido com o da regressão linear, o modelo significativo e tal, e podemos plotar o modelo.

E ele consegue capitar bem o modelo original não? Apesar de o intercepto não ter ficado muito correto, bem, existe um padrão nos dados, nos simulamos ele assim, e ele respondeu que existe esse padrão, e deu a forma mais ou menos correta.

Agora vamos ver para o segundo caso, onde a regressão linear não era de muita serventia.

> modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian) > summary(modelo_gam) Family: gaussian Link function: identity Formula: resposta_seno ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.9750 0.2667 18.66 3.46e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 5.861 0.00181 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.401 Deviance explained = 48.4% GCV = 2.5598 Scale est. = 2.1332 n = 30

E de novo conseguimos detectar que algo está acontecendo, esses dados não apenas acaso puro, e vamos ver como é o modelo que ele está propondo gráficamente.

Agora sim não? Mas que magia é essa? Bem o GAM funciona a partir de regressões locais como se medíssemos como os dados mudam ao longo da nossa amostra. Algo como o modelos de Broken-Stick em regressão segmentada ou em piecewise.

Então a é algo como a gente pega uma janelinha nos dados e faz uma regressão.

Depois outra janelinha e faz outra regressão.

E a gente vai ficando com um monte dessas regressões.

Depois é só fazer uma função, tipo piecewise bem grandona e juntar todas elas.

No nosso modelo clássico de regressão, a gente tem

Y_i = \alpha + \beta \cdot X_i + \varepsilon_i

onde o \varepsilon_i=N(0,\sigma^2),

já que falamos que estamos um erro com distrubição normal dos resíduos.

A mudança que o GAM faz é

Y_i = \alpha + f (X_i) + \varepsilon_i

onde o \varepsilon_i=N(0,\sigma^2)

e o  f (X_i) é uma função que descreve os dados, e não mais um coeficiente como era na regressão.

So que você pode pensar, ei, então esse negocio vai ajustar qualquer coisa, qualquer conjunto de dados ele vai achar uma uma função que explica!

Não é tão simples assim, podemos fazer um teste, gerar dados sem padrão nenhum e ver o que ele encontra.

E para um conjunto de dados sem padrão.

Family: gaussian Link function: identity Formula: resposta_aleatoria ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.2277 0.1509 -1.509 0.144 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 1.592 0.208 R-sq.(adj) = 0.0754 Deviance explained = 20.3% GCV = 0.81982 Scale est. = 0.68318 n = 30

Ele nos diz, não tem padrão nenhum aqui uai. O modelo simplesmente cobre tudo, não existe previsibilidade nenhuma.

Agora nem tudo são flores, e temos que tomar algumas decisões, sempre existe um preço a se pagar por algo tão legal assim, que é como definir a função que descreve os dados, veja que quando estavamo vendo como é feito ela, que é definida naquela função s(), separamos os dados em quadradinhos, fazendo pequenas regressões dentro deles certo, mas qual o tamanho ideal de quadradinho a usar? Isso é o que esta sendo definido naquele argumento k dentro da função s, e esta ligado ao nossos graus de liberdade, veja que quanto menor o quadradinho, menos dados por quadradinho então é mais facil definir errado a inclinação da reta, se tiver somente um ponto num quadradinho, nem tem como definir a inclinação da reta, agora se o quadradinho for muito grande, bem a regressão linear é como se fosse so um quadrado, então nao tem porque desse método, além de que podemos peder mudanças nos dados, então isso pode ser um problema, veja que no caso dos dados oscilando na função seno, se a gente usar um k muito baixo

modelo_gam<- gam(resposta_seno~s(preditor,k=2,fx=T),family = gaussian)

Avaliando os resíduos ainda vemos que algo não está correto.

Agora se aumentarmos o k
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)

E podemos ir aumentando até ficar aceitavel a distribuição dos residuos, a partir do momento que estiver aceitavel, aumentar mais é disperdicio pois estamos diminuindo os graus de confiança sem obter nenhum ganho em ter um modelo com melhor ajuste.

E é claro, alguém automatizou isso, então se você coloca -1 no k, a gente usa um algoritimo que define esse valor automaticamente.

Bem é isso ai, uma pequena introdução ao GAM, veja que aquela similaridade com a ideia de regressão não é so para ser bonitinho, isso permite trabalharmos com mais de um termo, interação entre termos, variáveis aleatória, bem tudo que da pra fazer com regressão da para fazer com GAM, sendo que regressão linear pode ser considerado um caso especial de GAM, que é quando a gente cai na função da reta. Eu acho que ando muito desanimado em casa, mas acho que preciso voltar a fazer posts com análise de dados que é muito legal pra ver se animo. O script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Zuur, AF 2012 Beginner’s Guide to Generalized Additive Models with R

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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
###Simulando Dados
set.seed(123)
preditor <- runif(30,0,4*pi)
resposta_pot<-5+2*(preditor)^2+rnorm(30,0,10)
resposta_seno<-5+2*sin(preditor)+rnorm(30)
 
##
jpeg("01.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
dev.off()
 
##
jpeg("02.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
modelo_linear <- lm(resposta_pot~preditor)
abline(modelo_linear, col="red")
dev.off()
 
##
summary(modelo_linear)
 
##
jpeg("03.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
curve(5+2*x^2,add=T,lwd=2,col="blue")
legend("topleft",col=c("blue","red"),lwd=c(2,1),legend=c("Realidade","Modelo"),bty="n")
 
modelo_linear <- lm(resposta_pot~preditor)
abline(modelo_linear, col="red")
novos_dados <- seq(0,4*pi,0.25)
pred_interval <- predict(modelo_linear, newdata=data.frame(preditor=novos_dados), interval="prediction",level = 0.95)
lines(novos_dados, pred_interval[,2], col="orange", lty=2)
lines(novos_dados, pred_interval[,3], col="orange", lty=2)
dev.off()
 
##
jpeg("04.jpg")
par(mfrow=c(2,2))
plot(modelo_linear)
dev.off()
 
##
jpeg("05.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12),xlab="Preditor",ylab="Resposta")
curve(5+2*sin(x),add=T,lwd=2,col="blue")
legend("topleft",col=c("blue","red"),lwd=c(2,1),legend=c("Realidade","Modelo"),bty="n")
 
modelo_linear <- lm(resposta_seno~preditor)
abline(modelo_linear, col="red")
novos_dados <- seq(0,4*pi,0.25)
pred_interval <- predict(modelo_linear, newdata=data.frame(preditor=novos_dados), interval="prediction",level = 0.95)
lines(novos_dados, pred_interval[,2], col="orange", lty=2)
lines(novos_dados, pred_interval[,3], col="orange", lty=2)
dev.off()
 
##
summary(modelo_linear)
 
##
jpeg("06.jpg")
par(mfrow=c(2,2))
plot(modelo_linear)
dev.off()
 
##Pacote utilizado  para o GAM
##install.packages("mgcv")
library(mgcv)
 
## Ajuste de modelo
modelo_gam<- gam(resposta_pot~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
## Visualização do modelo
jpeg("07.jpg")
plot(modelo_gam)
dev.off()
 
##
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
##
jpeg("08.jpg")
plot(modelo_gam)
dev.off()
 
##
jpeg("09.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4),y=c(2,4,4,2,2),type="l",lty=2)
selecao <- preditor>=4 & preditor<=6
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],4,6,add=T,lwd=2,col="red")
dev.off()
 
##
jpeg("10.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4)+3,y=c(2,4,4,2,2)*2+1,type="l",lty=2)
selecao <- preditor>=7 & preditor<=9
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],7,9,add=T,lwd=2,col="red")
dev.off()
 
##
jpeg("11.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4),y=c(2,4,4,2,2),type="l",lty=2)
selecao <- preditor>=4 & preditor<=6
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],4,6,add=T,lwd=2,col="red")
 
points(x=c(4,4,6,6,4)+3,y=c(2,4,4,2,2)*2+1,type="l",lty=2)
selecao <- preditor>=7 & preditor<=9
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],7,9,add=T,lwd=2,col="red")
dev.off()
 
##Testando gam contra dados aleatórios
resposta_aleatoria <- rnorm(30)
 
##
jpeg("12.jpg")
plot(resposta_aleatoria~preditor,pch=19,frame=F)
dev.off()
 
modelo_gam<- gam(resposta_aleatoria~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
##
jpeg("13.jpg")
plot(modelo_gam)
dev.off()
 
##Definindo um valor de K
modelo_gam<- gam(resposta_seno~s(preditor,k=2,fx=T),family = gaussian)
jpeg("14.jpg")
hist(resid(modelo_gam))
dev.off()
 
##
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)
jpeg("15.jpg")
hist(resid(modelo_gam))
dev.off()
 
##Auto ajuste de k
modelo_gam<- gam(resposta_seno~s(preditor,k=-1,fx=T),family = gaussian)
jpeg("16.jpg")
hist(resid(modelo_gam))
dev.off()