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.
Feito este preliminar, passemos à resolução.
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.
## 'data.frame': 10 obs. of 5 variables:
## $ mortalidade : int 2 5 6 7 12 15 36 45 58 124
## $ densidade.aero : int 5 25 9 36 11 9 8 7 41 25
## $ altitude : int 150 236 364 489 684 789 1254 1564 1789 1987
## $ proximidade.aero: int 3 5 6 8 9 12 25 69 154 351
## $ estradas : int 21 0 17 1 18 21 23 14 1 2## mortalidade densidade.aero altitude proximidade.aero
## Min. : 2.00 Min. : 5.00 Min. : 150.0 Min. : 3.0
## 1st Qu.: 6.25 1st Qu.: 8.25 1st Qu.: 395.2 1st Qu.: 6.5
## Median : 13.50 Median :10.00 Median : 736.5 Median : 10.5
## Mean : 31.00 Mean :17.60 Mean : 930.6 Mean : 64.2
## 3rd Qu.: 42.75 3rd Qu.:25.00 3rd Qu.:1486.5 3rd Qu.: 58.0
## Max. :124.00 Max. :41.00 Max. :1987.0 Max. :351.0
## estradas
## Min. : 0.00
## 1st Qu.: 1.25
## Median :15.50
## Mean :11.80
## 3rd Qu.:20.25
## Max. :23.00Todas as colunas apresentam valores numéricos.
Podemos agora ver como se comporta a variável mortalidade em função da 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)
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.
Podemos ver o nosso modelo graficamente adicionando uma linha (função abline()) ao gráfico anterior.
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()).
##
## Call:
## lm(formula = mortalidade ~ altitude, data = dpeol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.975 -10.707 -4.129 7.699 39.671
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -15.978726 10.139500 -1.576 0.153701
## altitude 0.050482 0.008995 5.612 0.000503 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.1 on 8 degrees of freedom
## Multiple R-squared: 0.7975, Adjusted R-squared: 0.7721
## F-statistic: 31.5 on 1 and 8 DF, p-value: 0.0005029Pro 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.
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?
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 (-15.9787257) e da variável altitude (0.0504822) 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.
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.
## [1] 59.74457Fazendo esta previsão visualmente.
# 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.
# valores observados e previstos de mortalidade para cada registo de altitude
cbind(altitude = dpeol$altitude,
mort.obs = dpeol$mortalidade,
mort.prev = predict(mod.alt))## altitude mort.obs mort.prev
## 1 150 2 -8.406397
## 2 236 5 -4.064929
## 3 364 6 2.396791
## 4 489 7 8.707065
## 5 684 12 18.551092
## 6 789 15 23.851722
## 7 1254 36 47.325940
## 8 1564 45 62.975419
## 9 1789 58 74.333912
## 10 1987 124 84.329385# 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))## altitude mort.getMort mort.prev
## 1 150 -8.406396 -8.406397
## 2 236 -4.064927 -4.064929
## 3 364 2.396795 2.396791
## 4 489 8.707070 8.707065
## 5 684 18.551099 18.551092
## 6 789 23.851730 23.851722
## 7 1254 47.325953 47.325940
## 8 1564 62.975435 62.975419
## 9 1789 74.333930 74.333912
## 10 1987 84.329406 84.329385# 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.
## 1
## 34.50346## [1] 34.50347Exercí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ávelproximidade.aerotambém poderá ter um efeito positivo e linear na mortalidade.
# 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)##
## Call:
## lm(formula = mortalidade ~ ., data = dpeol)
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10
## -2.146 1.872 -2.155 3.056 -1.753 -1.536 7.734 -2.491 -4.593 2.010
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.418208 10.009923 0.741 0.491957
## densidade.aero -0.394911 0.303776 -1.300 0.250299
## altitude 0.018310 0.004458 4.108 0.009286 **
## proximidade.aero 0.252193 0.029437 8.567 0.000357 ***
## estradas -0.228630 0.446504 -0.512 0.630427
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.872 on 5 degrees of freedom
## Multiple R-squared: 0.9908, Adjusted R-squared: 0.9835
## F-statistic: 135 on 4 and 5 DF, p-value: 2.804e-05Vendo o sumário do modelo podemos criar um modelo mais simples apenas com as covariáveis altitude e proximidade.aero.
##
## Call:
## lm(formula = mortalidade ~ altitude + proximidade.aero, data = dpeol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3555 -0.9453 -0.3324 1.3840 8.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.350733 3.477927 -0.388 0.70928
## altitude 0.017811 0.004552 3.912 0.00581 **
## proximidade.aero 0.245732 0.027434 8.957 4.4e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.48 on 7 degrees of freedom
## Multiple R-squared: 0.9837, Adjusted R-squared: 0.9791
## F-statistic: 211.8 on 2 and 7 DF, p-value: 5.474e-07Exercí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.
# 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.
## 'data.frame': 14 obs. of 5 variables:
## $ Densidade : num 51.4 72.1 53.2 83.2 57.4 66.5 98.3 74.8 92.2 97.9 ...
## $ Turbidez : num 0.2 1.9 0.2 10.7 6.8 10.6 9.6 6.3 10.8 9.6 ...
## $ Salinidade : num 17.8 29.4 17.2 30.2 15.3 17.6 35.6 28.2 34.7 35.8 ...
## $ Temperatura: num 24.6 20.7 18.5 10.6 8.9 11.1 10.6 8.8 11.9 10.8 ...
## $ Oxigenio : num 18.9 8 22.6 7.1 27.3 20.8 5.6 13.1 5.9 5.5 ...## Densidade Turbidez Salinidade Temperatura
## Min. :51.40 Min. : 0.200 Min. :15.30 Min. : 6.70
## 1st Qu.:63.73 1st Qu.: 2.000 1st Qu.:18.93 1st Qu.:10.60
## Median :78.20 Median : 8.200 Median :28.80 Median :11.40
## Mean :76.74 Mean : 7.171 Mean :26.99 Mean :14.39
## 3rd Qu.:91.17 3rd Qu.:10.575 3rd Qu.:33.58 3rd Qu.:19.70
## Max. :98.30 Max. :20.500 Max. :37.90 Max. :26.50
## Oxigenio
## Min. : 0.50
## 1st Qu.: 6.20
## Median : 9.00
## Mean :11.96
## 3rd Qu.:17.75
## Max. :27.30Temos 14 observações e 5 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.
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.
##
## Call:
## lm(formula = Densidade ~ Temperatura, data = ddens)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.317 -10.394 4.093 10.744 16.053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 97.6494 9.7515 10.014 3.53e-07 ***
## Temperatura -1.4531 0.6233 -2.331 0.038 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.3 on 12 degrees of freedom
## Multiple R-squared: 0.3117, Adjusted R-squared: 0.2544
## F-statistic: 5.435 on 1 and 12 DF, p-value: 0.03798Pelo 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.
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.
## 1 2 3 4 5 6
## -10.5040934 4.5289612 -17.5677773 0.9530256 -27.3171813 -15.0204429
## 7 8 9 10 11 12
## 16.0530256 -10.0624876 11.8420074 15.9436382 7.4513948 6.8860802
## 13 14
## 3.6567261 13.1571234## 1 2 3 4 5 6
## -10.5040934 4.5289612 -17.5677773 0.9530256 -27.3171813 -15.0204429
## 7 8 9 10 11 12
## 16.0530256 -10.0624876 11.8420074 15.9436382 7.4513948 6.8860802
## 13 14
## 3.6567261 13.1571234Podemos 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.
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.
Uma análise de correlação das duas variáveis ajudaria a justificar esta relação.
## [1] -0.5583191A correlação, \(R\), é naturalmente negativa, uma vez que um aumento da temperatura leva à redução da densidade.
O \(R^2\) é aproximadamente 0.3117. Este valor está presente no sumário do modelo enquanto coeficiente de determinação e indica que o modelo explica 31.17% da variabilidade da densidade da espécie.
## [1] 0.3117202Exercí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.
Fazemos o modelo com a densidade a depender de todas as variáveis dos dados.
##
## Call:
## lm(formula = Densidade ~ ., data = ddens)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1977 -1.3324 -0.7577 1.9245 4.0316
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -30.8800 37.3544 -0.827 0.42979
## Turbidez 2.0845 0.4562 4.570 0.00135 **
## Salinidade 2.5951 0.7353 3.529 0.00642 **
## Temperatura 0.6453 0.4584 1.408 0.19281
## Oxigenio 1.1147 0.7595 1.468 0.17622
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.079 on 9 degrees of freedom
## Multiple R-squared: 0.9761, Adjusted R-squared: 0.9655
## F-statistic: 91.83 on 4 and 9 DF, p-value: 2.728e-07Quando 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!
# modelo só com turbidez e salinidade
mod.TS <- lm(Densidade ~ Turbidez + Salinidade, data = ddens)
summary(mod.TS)##
## Call:
## lm(formula = Densidade ~ Turbidez + Salinidade, data = ddens)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6288 -1.5121 -0.7341 2.2797 4.8000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.8189 3.1332 7.921 7.17e-06 ***
## Turbidez 1.4783 0.1548 9.548 1.17e-06 ***
## Salinidade 1.5306 0.1152 13.282 4.07e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.114 on 11 degrees of freedom
## Multiple R-squared: 0.9701, Adjusted R-squared: 0.9647
## F-statistic: 178.5 on 2 and 11 DF, p-value: 4.13e-09Ajustar 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 0.6% 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.
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).
##
## Shapiro-Wilk normality test
##
## data: residuals(mod.TS)
## W = 0.94903, p-value = 0.5458Com 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.
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.
Dar uma vista de olhos aos dados.
## 'data.frame': 58 obs. of 6 variables:
## $ ANC : num 0.974 0.96 0.956 0.938 0.873 ...
## $ Salinity : num 11.09 8.74 8.35 5.53 6.14 ...
## $ Turbidity : num 745 300 930 1077 1508 ...
## $ Flow : num 1055 560 560 831 295 ...
## $ Temperature: num 10.4 21.6 21.6 15.9 25.1 10.4 10.8 9.5 21.6 11.3 ...
## $ NAO : num 0.23 -0.2 -0.2 0.92 -0.07 0.23 0.69 -0.83 -0.2 0.23 ...## ANC Salinity Turbidity Flow
## Min. :0.0000 Min. : 0.170 Min. : 41.93 Min. : 260.0
## 1st Qu.:0.2379 1st Qu.: 2.712 1st Qu.: 237.69 1st Qu.: 378.0
## Median :0.4383 Median : 6.030 Median : 364.23 Median : 691.9
## Mean :0.4822 Mean : 5.831 Mean : 629.78 Mean : 816.7
## 3rd Qu.:0.7014 3rd Qu.: 8.477 3rd Qu.: 762.75 3rd Qu.:1055.1
## Max. :0.9737 Max. :12.720 Max. :3198.20 Max. :2557.0
## Temperature NAO
## Min. : 8.2 Min. :-1.030000
## 1st Qu.:10.8 1st Qu.:-0.250000
## Median :11.7 Median :-0.030000
## Mean :15.1 Mean : 0.007241
## 3rd Qu.:21.6 3rd Qu.: 0.630000
## Max. :25.2 Max. : 0.920000Qual 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.
# 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")).
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.
##
## Call:
## glm(formula = (ANC + 0.001) ~ ., family = Gamma(link = "log"),
## data = danchova)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0642 -0.4833 -0.1056 0.3404 0.7205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.329e+00 3.361e-01 -3.954 0.000233 ***
## Salinity 3.825e-02 1.985e-02 1.927 0.059454 .
## Turbidity 3.118e-04 9.714e-05 3.210 0.002279 **
## Flow -9.458e-05 1.606e-04 -0.589 0.558554
## Temperature 1.470e-02 1.407e-02 1.045 0.300901
## NAO 1.402e-01 1.258e-01 1.115 0.270087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2537669)
##
## Null deviance: 26.544 on 57 degrees of freedom
## Residual deviance: 22.318 on 52 degrees of freedom
## AIC: 15.988
##
## Number of Fisher Scoring iterations: 6Olhando 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.
fitANC2 <- glm((ANC + 0.001) ~ Salinity + Turbidity, family = Gamma(link = "log"), data = danchova)
summary(fitANC2)##
## Call:
## glm(formula = (ANC + 0.001) ~ Salinity + Turbidity, family = Gamma(link = "log"),
## data = danchova)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1168 -0.5019 -0.1258 0.3641 0.7860
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.189e+00 1.402e-01 -8.478 1.46e-11 ***
## Salinity 4.157e-02 1.903e-02 2.184 0.03324 *
## Turbidity 2.992e-04 9.477e-05 3.157 0.00259 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2521114)
##
## Null deviance: 26.544 on 57 degrees of freedom
## Residual deviance: 23.011 on 55 degrees of freedom
## AIC: 11.874
##
## Number of Fisher Scoring iterations: 6A 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 = (0.17, 12.72) e Turbidity = (41.93, 3198.2)).
Fazendo a análise dos resíduos do modelo.
Temos a observação 58 como um valor algo discrepante das restantes.
## [1] 0Esta observação é o 0 registado na variável resposta. Podemos construir o mesmo modelo desta vez retirando essa observação nos dados.
fitANC2_clean <- glm(ANC ~ Salinity + Turbidity, family = Gamma(link = "log"), data = danchova[-58, ])
summary(fitANC2_clean)##
## Call:
## glm(formula = ANC ~ Salinity + Turbidity, family = Gamma(link = "log"),
## data = danchova[-58, ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9952 -0.5121 -0.1601 0.3376 0.7912
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.120e+00 1.385e-01 -8.083 7.24e-11 ***
## Salinity 3.341e-02 1.865e-02 1.791 0.07888 .
## Turbidity 2.900e-04 9.105e-05 3.185 0.00241 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2323509)
##
## Null deviance: 16.239 on 56 degrees of freedom
## Residual deviance: 13.299 on 54 degrees of freedom
## AIC: -8.3999
##
## Number of Fisher Scoring iterations: 6Ao 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.
É 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.
fitANC.T <- glm(ANC ~ Turbidity, family = Gamma(link = "log"), data = danchova[-58, ])
summary(fitANC.T)##
## Call:
## glm(formula = ANC ~ Turbidity, family = Gamma(link = "log"),
## data = danchova[-58, ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8574 -0.5423 -0.2000 0.3388 0.9060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.163e-01 8.992e-02 -10.190 2.83e-14 ***
## Turbidity 2.913e-04 9.508e-05 3.063 0.00339 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2536607)
##
## Null deviance: 16.239 on 56 degrees of freedom
## Residual deviance: 14.010 on 55 degrees of freedom
## AIC: -7.3104
##
## Number of Fisher Scoring iterations: 6Modelo 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.
# 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)
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.
Ver os dados.
## 'data.frame': 10 obs. of 5 variables:
## $ presenca : int 2 3 4 3 2 1 12 15 17 16
## $ altitude : int 300 256 246 178 231 125 569 785 569 654
## $ floresta : int 5 9 7 2 8 9 25 36 41 25
## $ temperatura: int 25 22 21 24 22 28 10 8 12 14
## $ estradas : int 21 17 14 18 23 21 0 1 1 2## presenca altitude floresta temperatura
## Min. : 1.00 Min. :125.0 Min. : 2.00 Min. : 8.0
## 1st Qu.: 2.25 1st Qu.:234.8 1st Qu.: 7.25 1st Qu.:12.5
## Median : 3.50 Median :278.0 Median : 9.00 Median :21.5
## Mean : 7.50 Mean :391.3 Mean :16.70 Mean :18.6
## 3rd Qu.:14.25 3rd Qu.:569.0 3rd Qu.:25.00 3rd Qu.:23.5
## Max. :17.00 Max. :785.0 Max. :41.00 Max. :28.0
## estradas
## Min. : 0.00
## 1st Qu.: 1.25
## Median :15.50
## Mean :11.80
## 3rd Qu.:20.25
## Max. :23.00Temos 5 variáveis, com 10 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…)
##
## Call:
## glm(formula = presenca ~ ., family = Gamma(link = "inverse"),
## data = dpresenca)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7
## -0.13372 0.14525 0.30615 0.20438 0.02682 -0.73211 -0.12128
## 8 9 10
## -0.00656 0.04631 0.03274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2824582 0.1578231 1.790 0.1335
## altitude -0.0001675 0.0001378 -1.216 0.2783
## floresta -0.0009144 0.0015712 -0.582 0.5858
## temperatura -0.0090641 0.0066797 -1.357 0.2328
## estradas 0.0207172 0.0053274 3.889 0.0115 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1078066)
##
## Null deviance: 8.76092 on 9 degrees of freedom
## Residual deviance: 0.72915 on 5 degrees of freedom
## AIC: 45.853
##
## Number of Fisher Scoring iterations: 5Apenas a variável estradas parece ter um impacto significativo no modelo. Podemos tentar construir um modelo com essa variável.
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)
fit.estradas <- glm(presenca ~ estradas, family = Gamma(link = "inverse"), data = dpresenca)
summary(fit.estradas)##
## Call:
## glm(formula = presenca ~ estradas, family = Gamma(link = "inverse"),
## data = dpresenca)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.72360 -0.10694 0.06731 0.16441 0.36233
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.050970 0.010765 4.735 0.001473 **
## estradas 0.018495 0.002807 6.588 0.000171 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.09651818)
##
## Null deviance: 8.76092 on 9 degrees of freedom
## Residual deviance: 0.98988 on 8 degrees of freedom
## AIC: 42.953
##
## Number of Fisher Scoring iterations: 5Finalmnete, vamos visualizar esta relação.
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 :)