---
title: "Aula TP 9 -- Resolução FT 7"
author: "João Malato e Tiago A. Marques"
date: "Novembro de 2019"
output:
  prettydoc::html_pretty:
    theme: cayman
    highlight: github
    toc: true
    toc_depth: 1
---

<!-- Para ter o documendo justificado dos dois lados -->
<style>
body {
text-align: justify}
</style>


# Disclaimer

Este documento fornece algumas propostas de resolução aos exercícios da folha de trabalho `TP7.pdf`, resolvida durante a aula prática número 9. Como sempre, é necessário notar que pode haver (quase) infinitas maneiras de solucionar cada um dos exercícios. Se esta resolução não for *exactamente igual* à que foi apresentada na aula, ou à que cada um de vós conseguiu criar, não tem mal. O importante aqui é que sejam capazes de compreender o que fizeram e entender o que é apresentado neste documento para que possam mais à frente reproduzi-lo de forma independente.

Neste documento assumimos que os ficheiros de dados se encontram na mesma directoria do .Rmd. Caso não seja o caso, será necessário primeiro identificar essa directoria.

Os gráficos aqui apresentados são previamente trabalhados para a melhor compreensão dos resultados. Não é obrigatório representarem sempre gráficos assim, embora seja aconselhado.

Além do HTML em si, iremos disponibilizar o .Rmd para que possam também ver como se faz este relatório e quem sabe colher alguma inspiração para implementar nas aulas.  
Nota: caso queiram usar este template precisam instalar o package `prettydoc` para isso acontecer, Podem dar uma vista de olhos aos diversos temas e como instalar [aqui](https://prettydoc.statr.me/index.html).

Feito este preliminar, passemos à resolução.

<br>

***

# Exercício 1.

**Estimou a mortalidade de aves em alguns parques eólicos, em função de variáveis ambientais (`DataTP7parqueseolicos.csv`). Verifique se a variável `altitude` é um bom preditor da mortalidade, através de um modelo de regressão simples. Caracterize a relação entre as variáveis. Elabore uma função que gere estimativas de mortalidade em função da altitude. Avalie o desempenho do modelo de regressão.**

Como sempre, lemos os dados e fazemos uma rápida análise antes de partirmos para a resolução do exercício.

```{r}
dpeol <- read.csv2("DataTP7parqueseolicos.csv")
```

```{r}
str(dpeol)
summary(dpeol)
```

Todas as colunas apresentam valores numéricos.  
Podemos agora ver como se comporta a variável `mortalidade` em função da altitude.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
plot(dpeol$altitude, dpeol$mortalidade, xlab="Altitude (m)", ylab="Mortalidade", ylim=c(0,140), xlim=c(0,2000), pch=20, col="cornflowerblue", las=1, cex=2)
box(lwd=2)
```

Parece haver uma dependência das mortalidade em função da altitude. À medida que a altitude aumenta, também a mortalidade aumenta. Podemos implementar um modelo linear simples para tentar compreender melhor esta relação.

```{r}
mod.alt <- lm(mortalidade ~ altitude, data = dpeol)
```

Podemos ver o nosso modelo graficamente adicionando uma linha (função `abline()`) ao gráfico anterior.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
plot(dpeol$altitude, dpeol$mortalidade, xlab="Altitude (m)", ylab="Mortalidade", ylim=c(0,140), xlim=c(0,2000), pch=20, col="cornflowerblue", las=1, cex=2)
abline(mod.alt, lwd = 2, lty = 2, col = "gray35")
box(lwd=2)
```

Este resultado vem corroborar a clara evidência de que à medida que a altitude aumenta, a mortalidade também aumenta.  
Podemos então ver o output do modelo (função `summary()`).

```{r}
summary(mod.alt)
```

> **Pro tip:** Quando se olha para o "sumário" de um modelo há várias coisas que podemos ter em atenção. Enumerando algumas: 1) os valores dos coeficientes do modelo (`Estimate`), 2) o erro associado (`Std. Error`) e 3) o p-value para o teste-t do coeficiente (`Pr(>|t|)`), que nos diz indica se a variável é considerada significante no modelo (p-value < $\alpha$). Podemos também ver 4) o coeficiente de determinação (`Multiple R-squared`) que nos indica o quão próximos os dados estão da linha ajustada do modelo, ou seja explica a variabilidade do modelo (podendo ser interpretado em percentagem) e 5) o p-value associado ao teste-F (`p-value`) para a significância global do modelo. Dá para estar atendo e ver mais resultados, mas se conseguirem idenificar esses já é bom para conseguirem interpretar qualquer modelo linear.

Tendo o modelo ajustado, vamos agora olhar para os residuos do modelo.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
res.alt <- residuals(mod.alt)
plot(res.alt, pch=20, col="coral3", las=1, cex=2, ylab = "Resíduos")
abline(h = 0, lty = 2, col="gray35")
box(lwd=2)
```

Seriam Gaussianos?

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
hist(res.alt, main="Resíduos do modelo Mortalidade~Altitude", xlab="Resíduos", col="coral3", las=1)
box(lwd=2)
```

Sabendo o valor do `Intercept` (`r coef(mod.alt)[[1]]`) e da variável `altitude` (`r coef(mod.alt)[2]`) no modelo,
$$Y\ \tilde\ \beta_0 + \beta_1 X_1\\ $$ ou $$\\ \text{Mortalidade}\ \tilde\ -15.9787257 + 0.0504822\cdot\text{Altitude}$$
podemos prever a mortalidade. A maneira mais clara de perceber o que se passa é criar uma função personalizada.

```{r}
getMort <- function(altitude) {
  # tal como vimos na fórmula
  mortalidade <- -15.9787257 + 0.0504822 * altitude
  return(mortalidade)
}
```

Tendo esta função, podemos prever a mortalidade para a altitude de, por exemplo, 1500 metros. 

```{r}
getMort(1500)
```

Fazendo esta previsão visualmente.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
# gráfico mortalidade~altitude
plot(dpeol$altitude, dpeol$mortalidade, xlab = "Altitude (m)", ylab = "Mortalidade", ylim = c(0, 140), xlim = c(0, 2000), pch = 20, col = "cornflowerblue", las = 1, cex = 2)
# linha do modelo
abline(mod.alt, lwd = 2, lty = 2, col = "gray35")
# linha dos 1500 metros
abline(v = 1500, lty = 3, col = "gray50")
# mortalidade a 1500 metros
abline(h = getMort(1500), lty = 3, col = "gray50")
# cruzamento das linhas
points(x = 1500, y = getMort(1500), pch = 4, col = "firebrick", cex = 2, lwd = 3)
box(lwd=2)
```

Se quisermos ver qual a predição para cada um dos pontos observados, que corresponde à projecção desses valores na linha preditiva do modelo podemos fazê-lo aravés da função `predict()`. Esta função não faz mais do que calcular (através de uma função semelhante à que foi criada, `getMort()`) o valor da mortalidade associada aos diferentes pontos de altitude observados.

```{r}
# valores observados e previstos de mortalidade para cada registo de altitude
cbind(altitude  = dpeol$altitude,
      mort.obs  = dpeol$mortalidade,
      mort.prev = predict(mod.alt))

# podemos ver o nossa função lado-a-lado com a função predict()
cbind(altitude     = dpeol$altitude,
      mort.getMort = getMort(dpeol$altitude),
      mort.prev    = predict(mod.alt))
```

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
# gráfico mortalidade~altitude
plot(dpeol$altitude, dpeol$mortalidade, xlab="Altitude (m)", ylab="Mortalidade", ylim=c(0,140), xlim=c(0,2000), pch=20, col="cornflowerblue", las=1, cex=2)
# linha do modelo
abline(mod.alt, lwd = 2, lty = 2, col = "gray35")
# predição dos pontos na linha
points(dpeol$altitude, predict(mod.alt), pch = 4, cex = 2, col = 2, lwd = 3)
# podemos adicionar um segmento de recta que mostre o ajuste dos pontos observados à recta do modelo
segments(x0=dpeol$altitude, y0=dpeol$mortalidade, x1=dpeol$altitude, y1=predict(mod.alt), lty=3, col="gray50")
box(lwd=2)
```

Se quisermos prever a mortalidade para outros valores da altitude podemos usar a função `predict()` ou a `predict.lm()`, indicando o modelo e o valor da variável explicativa que pretendemos estudar.

```{r}
# usando a função
predict.lm(mod.alt, newdata = data.frame(altitude = 1000))

# ou podemos usar a nossa função...
getMort(1000)
```

<br>

***

# Exercício 2.

**Realize uma regressão múltipla aos dados `DataTP7parqueseolicos.csv` e interprete os resultados. Pretende obter um sub-modelo para utilizar noutras regiões com parques eólicos. Obtenha um sub-modelo para esse efeito.**

Vamos agora criar um modelo linear com mais covariáveis. Como logo à partida não sabemos as relações entre as variáveis podemos começar por criar um modelo com todas as variáveis.

> **Pro tip:** Quando não temos grande conhecimento de potenciais relações entre variáveis podemos usar a função `pairs()` que cria gráficos com todas as variáveis presentes nos dados para que possamos visualmente estudar as relações de dependência. No exemplo em baixo é possível ver que o segundo gráfico da primeira linha represente a relação da mortalidade em função da altitude, como já estudámos no exercício anterior. Podemos talmbém ver que a variável `proximidade.aero` também poderá ter um efeito positivo e linear na mortalidade.

```{r, fig.width=7, fig.height=7, fig.align='center'}
pairs(dpeol, pch=16)
```

```{r, tidy=T}
# modelo completo, com todas as covariáveis presentes nos dados
mod.multi <- lm(mortalidade ~ ., data = dpeol)
# o ponto "." é uma maneira mais rápida de incluir todas as vars

# mod.multi <- lm(mortalidade ~ densidade.aero + altitude + proximidade.aero + estradas, data = dpeol)

summary(mod.multi)
```

Vendo o sumário do modelo podemos criar um modelo mais simples apenas com as covariáveis `altitude` e `proximidade.aero`.

```{r, tidy=T}
mod.multi2 <- lm(mortalidade ~ altitude + proximidade.aero, data = dpeol)
summary(mod.multi2)
```

<br>

***

# Exercício 3.

**Efectue uma análise da regressão aos dados `DataTP7densidade.csv`, os quais são relativos à densidade de uma espécie em função de variáveis ambientais. Seja a variável dependente a Densidade e a independente a Temperatura. Explore os resultados e efectue uma avaliação do cumprimento dos pressupostos.**

Como sempre, lemos os dados. **Ter sempre atenção** à estrutura dos dados. Abrir os dados antes de os importar para ver como foram recolhidos, qual a separação das colunas e qual o tipo de registo quanto aos números decimais.

```{r, tidy=T}
# ler os dados, ter em atenção o tipo de dados e como se separam os numeros decimais (virgulas ou pontos)
ddens <- read.csv2("DataTP7densidade.csv", dec=".")
```

Ver a estrutura e sumário dos dados.

```{r}
str(ddens) # tudo valores numéricos

summary(ddens)
```

Temos `r nrow(ddens)` observações e `r ncol(ddens)` variáveis recolhidas. Vamos usar a `Densidade` como variável dependente e a `Temperatura` como variável independente. Antes de partirmos para a construção do modelo podemos estudar visualmente a relação entre estas variáveis.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
plot(ddens$Temperatura, ddens$Densidade, xlim = c(5, 30), xlab = "Temperatura", ylab = "Densidade", pch = 20, cex = 2, col = "brown", las = 1)
box(lwd=2)
```

Podemos ajustar o modelo linear simples.

```{r}
mod.temp <- lm(Densidade ~ Temperatura, data = ddens)
summary(mod.temp)
```

Pelo sumário do modelo é possível dizer que a temperatura tem um efeito negativo na densidade da espécie. Quando a temperatura aumenta, o densidade vai reduzir. Mais correctamente podemos dizer que por grau de temperatura incrementado o modelo estima que a densidade vai decrescer em aproximadamente 1.45.  
Vemos a linha do modelo ajustado no gráfico anterior.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
plot(ddens$Temperatura, ddens$Densidade, xlim = c(5, 30), xlab = "Temperatura", ylab = "Densidade", pch = 20, cex = 2, col = "brown", las = 1)
abline(mod.temp, lwd = 2, lty = 2, col = "gray35")
box(lwd=2)
```

Tendo o modelo ajustado podemos (e devemos) fazer a análise dos resíduos do modelo.

```{r}
# método 1
residuals(mod.temp)

# método 2 para obter os resíduos
mod.temp$residuals
```

Podemos contruir gráficos de diagnóstico do modelo ajustado.

> **Pro tip:** O conjunto de gráficos de diagnóstio ajuda, de uma maneira visual a explicar o quão bem o modelo se ajusta aos dados. Os dois gráficos do lado direito ajudam a análise de resíduos; o objectivo aqui é não encontrar padrões definidos, garantindo assim a homocedasticidade (i.e. variância constante) dos resíduos (com uma linha vermelha horizontal em torno do valor 0). O gráfico do canto superior direito é um gráfico Q-Q ou *quantil-quantil*, este gráfico coloca os resíduos standardizados em função de quantis teóricos calculados considerando uma distribuição Normal; visualmente, se os pontos (ou maioria deles) se ajustarem bem à linha diagonal, pode-se assumir que os dados do modelo têm uma distribuição Normal. Esta ideia pode depois ser corroborada com a aplicação de um teste de normalidade de Shapiro-Wilk. Por fim, o gráfico do canto inferior direito define intervalos para a Distância de Cook. Esta é uma medida utilizada para identificar potenciais observações influentes e/ou candidatos a outliers; visualmente, uma observação seria considerada influente ou outlier caso passe o intervalo 0.5 ou 1.

```{r, fig.width=7, fig.height=7, fig.align='center'}
par(mfrow = c(2, 2))
plot(mod.temp, pch=16, col="coral3")
```

Olhando para os resultados não parece haver nenhum motivo para considerar o modelo como pouco adequado, por isso poderíamos pensar usar o modelo para prever a densidade da espécie em função da temperatura associada. No entanto, é de notar que estamos a trabalhar com um conjunto de dados relativamente pequeno o que dificultaria qualquer que fosse o motivo para considerar o modelo como errado.

<!-- We can look at the diagnostic plots for the fitted model. There is no reason to suspect that the model is not adequate, so we could use this model for predicting `Densidade` as a function of `Temperatura`. Nonetheless, note that with such a small sample size, it would be almost impossible to find any strong reasons to find the model unsuitable, unless these were blatant! -->

Uma análise de correlação das duas variáveis ajudaria a justificar esta relação.

```{r}
cor(ddens$Temperatura, ddens$Densidade)
```

A correlação, $R$, é naturalmente negativa, uma vez que um aumento da temperatura leva à redução da densidade.  
O $R^2$ é aproximadamente `r round(cor(ddens$Temperatura,ddens$Densidade)^2,4)`. Este valor está presente no sumário do modelo enquanto coeficiente de determinação e indica que o modelo explica `r round(with(ddens,cor(Temperatura,Densidade))^2*100,2)`% da variabilidade da densidade da espécie.

```{r}
# mostrar como se calcula o R^2
cor(ddens$Temperatura, ddens$Densidade) ^ 2
```

<br>

***

# Exercício 4.

**Com base nos dados da alínea anterior, efectue uma análise de regressão múltipla, sendo a variável dependente a densidade e as independentes todas as restantes. Comente os resultados e verifique se os pressupostos são cumpridos.**

<!-- ```{r, fig.width=8, fig.height=6, fig.align='center'} -->
<!-- pairs(ddens, pch=16) -->
<!-- ``` -->

Fazemos o modelo com a densidade a depender de todas as variáveis dos dados.

```{r}
mod.dens <- lm(Densidade ~ ., data = ddens)
summary(mod.dens)
```

Quando todas as variáveis são incluídas no modelo apenas a `Turbidez` e a `Salinidade` surgem como variáveis significativas no modelo. A variável `Temperatura`, que no exercício anterior parecia funcionar bem quando isolada aqui pode ser excluída!

```{r}
# modelo só com turbidez e salinidade
mod.TS <- lm(Densidade ~ Turbidez + Salinidade, data = ddens)
summary(mod.TS)
```

Ajustar um modelo só com as duas variáveis faz com elas se tornem ainda mais significativas. Este novo modelo `mod.TS` tem um coeficiente de determinação de 0.9701. Comparando com o estimado no modelo `mod.dens` podemos dizer que retirar as variáveis para a temperatura e oxigénio apenas se perde `r (0.9761-0.9701)*100`% da explicação da variabilidade (o que é uma informação vestigial). E além disso, 0.97 é muito bom!

Vamos então fazer uma análise aos resíduos do modelo. Começando por uma análise visual.

```{r, fig.width=7, fig.height=7, fig.align='center'}
par(mfrow = c(2, 2))
plot(mod.TS, pch=16)
```

Os resíduos parecem bons (com a excepção da observação 14 que pode potencialmente estar a influenciar a tendência do modelo -- experimentar retirar o ponto e voltar a ajustar o modelo pode ser uma boa solução), podemos complementar o gráfico quantil-quantil para testar a sua normalidade (teste de Shapiro-Wilk).

```{r}
shapiro.test(residuals(mod.TS))
```

Com um p-value maior que os niveis de significância usuais não temos qualquer evidência para rejeitar a hipótese de que os residuos sejam provenientes de uma distribuição Gaussiana, pelo que não há evidencias para suspeitar que o modelo não seja adequado.

<br>

***

# Exercício 5.

**Aplique um GLM aos dados `DataTP7anchova.csv`, referentes à abundância de anchova (variável dependente) em função de várias variáveis ambientais (variáveis independentes). Explore os resultados.**

```{r}
danchova <- read.csv2("DataTP7anchova.csv", dec = ".")
```

Dar uma vista de olhos aos dados.

```{r}
# estrutura
str(danchova)

# sumário
summary(danchova)
```

Qual será distribuição adequada para a variável resposta `ANC`? Por vezes a distribuição das observações pode ajudar-nos a ter uma ideia de como se comportam os dados.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
# histograma para tentar ver a dist da veriável resposta
hist(danchova$ANC, ylim = c(0, 20), col = "cyan4", xlab = "Abundância de Anchova", main = "")
box(lwd = 2)
```

O histograma não parece dar grande informação sobre o tipo de distribuição com que estamos a lidar. Podemos ver e construir um modelo linear generalizado (GLM) assumindo que a variável resposta `ANC` tem uma distribuição Gamma com função de ligação logarítmica (`family = Gamma(link = "log")`).

```{r}
# danchova <- glm(ANC ~ ., data = danchova, family = Gamma(link = "log"))
```

No entanto, ajustar o GLM (em cima comentado) devolverá um erro uma vez que existem valores de 0 nas observações da variável resposta, coisa que é incomportável quando se pretende aplicar trabalhar com uma Gamma. Uma maneira de contornar este problema é adicionar uns *pozinhos* à variável resposta (`ANC + 0.001`), assim garantindo que estes zeros passam a observações positivas e o modelo pode ser ajustado sem que os valores sejam significantemente adulterados.

```{r, tidy=T}
fitANC <- glm((ANC + 0.001) ~ ., family = Gamma(link = "log"), data = danchova)
summary(fitANC)
```

Olhando para o sumário dos dados há duas covariáveis que saltam à vista, a salinidade e a turbidez da água, que de um ponto de vista biológico pode fazer sentido. Podemos criar um modelo para a densidade das anchovas incluindo apenas essas duas variáveis.

```{r, tidy=T}
fitANC2 <- glm((ANC + 0.001) ~ Salinity + Turbidity, family = Gamma(link = "log"), data = danchova)
summary(fitANC2)
```

A covariável `Salinity` no modelo comlpeto surgia como *borderline* significativa (ao nível de significância $\alpha$ = 0.05 já não seria considerada significativa). Neste segundo modelo torna-se importante para a descrição da densidade das anchovas. É de notar que as estimativas das covariáveis são baixas (na ordem das milésimas), pelo que apesar de significativas só se sente o seu impacto no modelo devido às suas amplitudes (`Salinity` = (`r range(danchova$Salinity)`) e `Turbidity` = (`r range(danchova$Turbidity)`)).

<!-- ```{r} -->
<!-- # modelo linear -->
<!-- lmANC2 <- lm((ANC + 0.001) ~ Salinity + Turbidity, data = danchova) -->
<!-- summary(lmANC2) -->
<!-- ``` -->

Fazendo a análise dos resíduos do modelo.

```{r, fig.width=7, fig.height=7, fig.align='center'}
par(mfrow = c(2, 2))
plot(fitANC2, pch = 20, col = "cyan4")
```

Temos a observação 58 como um valor algo discrepante das restantes. 

```{r}
# vamos ver o valor dessa observação
danchova$ANC[58]
```

Esta observação é o 0 registado na variável resposta. Podemos construir o mesmo modelo desta vez retirando essa observação nos dados.

```{r, tidy=T}
fitANC2_clean <- glm(ANC ~ Salinity + Turbidity, family = Gamma(link = "log"), data = danchova[-58, ])
summary(fitANC2_clean)
```

Ao retirar esta observação a Salinidade mais uma vez deixa de ser significativa considerando um nivel de significancia de 5% na explicação da densidade de anchovas. Ainda assim podemos fazer a análise de resíduos e comparar as diferenças com o modelo `fitANC2`.

```{r, fig.width=7, fig.height=7, fig.align='center'}
par(mfrow=c(2,2))
plot(fitANC2_clean, pch = 20, col = "cyan4")
```

É possível ver um gráfico quantil-quantil com os pontos mais bem ajustados à linha, ainda que os extremos não o estejam a 100%. O mesmo se pode dizer do gráfico Scale-Location, que já não mostra o ponto 58 a influenciar a linha.  
O modelo pode ser considerado aceitável e o exercício poderia ficar por aqui, interpretando que variáveis podem influenciar a densidade da espécie e em que cenários (com ou sem zeros nas observações). Mas como somos seres curiosos e não custa nada ir um bocado mais além, vamos ajustar um modelo para a densidade das anchovas (sem o 0) apenas com a variável turbidez enquanto variável explicativa.

```{r}
fitANC.T <- glm(ANC ~ Turbidity, family = Gamma(link = "log"), data = danchova[-58, ])
summary(fitANC.T)
```

Modelo bem ajustado e com uma análise de resíduos semelhante às anteriores. Tendo apenas uma covariável podemos estudar visualmente a relação de dependência da variável alearória, incluindo também uma curva para os valores preditivos.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
# criamos sequencia de valores de turbidez onde queremos prever a densidade de anchovas
turb.seq <- seq(40, 3200, 0.2)

# fazemos uso da função predict.glm() para prever a densidade em função da sequancia de valores, tendo em conta o modelo criado
predict.turb <- predict.glm(fitANC.T, newdata = data.frame(Turbidity = turb.seq), type = "response")

# gráfico com os pontos densidade em função da turbidez
plot(danchova$Turbidity, danchova$ANC, xlab = "Turbidez", ylab = "Densidade Anchova", pch = 20, cex = 2, col = "cyan4", xlim = c(0, 3500), ylim = c(0, 1), las = 1)
# linha para valores preditivos em função da sequencia de turbidez
lines(turb.seq, predict.turb, lty = 2, col = "gray35", lwd = 2)
box(lwd = 2)
```

<br>

***

# Exercício 6.

**Aplique um GLM aos dados `DataTP7presenca.csv`, os quais descrevem a abundância de uma espécie em função de algumas variáveis ambientais. Utilize a distribuição Gama como família de distribuição do erro (Gamma) e a função inversa como função de ligação (inverse).**

Ler os dados.

```{r}
dpresenca <- read.csv2("DataTP7presenca.csv")
```

Ver os dados.

```{r}
# estrutura
str(dpresenca)

# sumário
summary(dpresenca)
```

Temos `r ncol(dpresenca)` variáveis, com `r nrow(dpresenca)` observações. Podemos assumir que a variável `presenca` depende de todas as outras (nota: com 10 observações, um modelo com 5 variáveis é demasiado complexo para correr bem, mas vamos experimentar...)

```{r, tidy=T}
fit.pres <- glm(presenca ~ ., family = Gamma(link = "inverse"), data = dpresenca)
summary(fit.pres)
```

Apenas a variável `estradas` parece ter um impacto significativo no modelo. Podemos tentar construir um modelo com essa variável.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
plot(dpresenca$estradas, dpresenca$presenca, xlab = "Estradas", ylab = "Presença", pch = 20, cex = 2, col = "goldenrod", xlim = c(0, 25), ylim = c(0,20), las = 1)
box(lwd = 2)
```

```{r, tidy=T}
fit.estradas <- glm(presenca ~ estradas, family = Gamma(link = "inverse"), data = dpresenca)
summary(fit.estradas)
```

Finalmnete, vamos visualizar esta relação.

```{r, fig.width=8, fig.height=6, fig.align='center', tidy=T}
est.seq <- seq(0, 25, 0.1)
predict.estradas <- predict.glm(fit.estradas, newdata = data.frame(estradas = est.seq), type = "response")

plot(dpresenca$estradas, dpresenca$presenca, xlab = "Estradas", ylab = "Presença", pch = 20, cex = 2, col = "goldenrod", xlim = c(0, 25), ylim = c(0, 20), las = 1)
lines(est.seq, predict.estradas, lty = 2, col = "gray50", lwd = 2)
box(lwd = 2)
```

Na ficha TP7b este conjunto de dados vai voltar a ser explorado, de uma forma mais minuciosa... até lá... bom estudo :)

<br>

***