Modelação estatística na análise em processos ambientais

Carla Sofia Carneiro Gomes da Silva Modelação Estatística na Análise em Processos Ambientais outubro de 2019 UMinho | 2019 Carla Sofia Carneiro Gomes da Silva Modelação Estatística na Análise em Processos Ambientais Universidade do Minho Escola de Ciências Carla Sofia Carneiro Gomes da Silva Modelação Estatística na Análise em Processos Ambientais Dissertação de Mestrado Mestrado em Estatística Trabalho efetuado sob a orientação da Professora Doutora Arminda Manuela Andrade Pereira Gonçalves e da Professora Doutora Susana Margarida Ferreira de Sá Faria Universidade do Minho Escola de Ciências outubro de 2019 ii DIREITOS DE AUTOR E CONDIÇÕES DE UTILIZAÇÃO DO TRABALHO POR TERCEIROS Este é um trabalho académico que pode ser utilizado por terceiros desde que respeitadas as regras e boas práticas internacionalmente aceites, no que concerne aos direitos de autor e direitos conexos. Assim, o presente trabalho pode ser utilizado nos termos previstos na licença abaixo indicada. Caso o utilizador necessite de permissão para poder fazer um uso do trabalho em condições não previstas no licenciamento indicado, deverá contactar o autor, através do RepositóriUM da Universidade do Minho. Licença concedida aos utilizadores deste trabalho Atribuição-NãoComercial-SemDerivações CC BY-NC-ND https://creativecommons.org/licenses/by-nc-nd/4.0/ iii DECLARAÇÃO DE INTEGRIDADE Declaro ter atuado com integridade na elaboração do presente trabalho académico e confirmo que não recorri à prática de plágio nem a qualquer forma de utilização indevida ou falsificação de informações ou resultados em nenhuma das etapas conducente à sua elaboração. Mais declaro que conheço e que respeitei o Código de Conduta Ética da Universidade do Minho. Universidade do Minho, 31 de outubro de 2019, (Carla Sofia Carneiro Gomes da Silva) Agradecimentos “Os meus mestres foram todos os homens e mulheres que me deslumbraram em leitura e n˜ao s´o: em exemplos de vida.” (Nat´alia Correia, 1983) A partilha do conhecimento permite a uni˜ao de pessoas, a sua coopera¸c˜ao faculta a realiza¸c˜ao de um trabalho maior. Agrade¸co a todas as pessoas que permitiram de alguma forma a conce¸c˜ao deste trabalho, momentˆanea ou continuamente. Este trabalho n˜ao seria poss´ıvel sem a orienta¸c˜ao da Professora Doutora A. Manuela Gon¸calves e da Professora Doutora Susana Faria. Toda a disponibilidade, o entusiasmo pela procura constante, a partilha de conhecimento e a tranquilidade determinaram a dire¸c˜ao a seguir. E estou grata pela oportunidade de aprender e de vivenciar ensinamentos das “minhas” Professoras, que acompanharam todo este percurso. Um agradecimento ao Engenheiro Vitorino Jos´e pela colabora¸c˜ao prestada na perce- ¸c˜ao das ferramentas inform´aticas, disponibilizadas pela Agˆencia Portuguesa do Ambiente (APA). N˜ao podendo esquecer, pelo tempo e pela brevidade na resposta, na cr´onica da fronteira da RH3, o Doutor Lu´ıs Margalho. Pelo apoio incessante e pela paciˆencia incans´avel durante todo este percurso, agrade¸co `a minha fam´ılia e a todos os meus amigos. v vi Resumo A degrada¸c˜ao do ambiente ´e atualmente um tema de grande importˆancia, quer pela dificuldade na recupera¸c˜ao e na reabilita¸c˜ao, quer tamb´em pelas gravosas consequˆencias sociais e econ´omicas. Em parte, a crise ambiental ´e o somat´orio de muitos erros cometidos pelo Homem e que ainda hoje ´e poss´ıvel observar. Investiga¸c˜oes realizadas no intuito de minimizar ou estimar problemas ambientais tˆem levado a estudos mais aprofundados de m´etodos, que possibilitem uma melhor perce¸c˜ao dos dados associados a estes problemas. Neste estudo, no contexto de um problema de monitoriza¸c˜ao de Qualidade da ´ Agua de superf´ıcie de uma bacia hidrogr´afica, prop˜oe-se uma abordagem baseada em modelos espaciais e temporais com o objetivo de analisar e avaliar a evolu¸c˜ao de s´eries temporais de vari´aveis ambientais. Os dados dizem respeito `a bacia hidrogr´afica do rio Douro localizada no Norte de Portugal. Para o processo de modela¸c˜ao, consideraram-se as s´eries temporais relativas `a vari´avel de qualidade de Oxig´enio Dissolvido (OD), medido mensalmente no per´ıodo de mar¸co de 2002 a fevereiro de 2013. Com o objetivo de obter estimativas de valores mensais de precipita¸c˜ao, em ´area, nas esta¸c˜oes de amostragem de qualidade (onde n˜ao h´a medi¸c˜oes de precipita¸c˜ao), ´e desenvolvida uma metodologia com recurso a processos estoc´asticos espaciais (Kriging), a ser aplicada aos dados de precipita¸c˜ao existentes nesta bacia. Os valores estimados v˜ao representar o fator hidrometeorol´ogico nas esta¸c˜oes de qualidade, para o processo de modela¸c˜ao do Oxig´enio Dissolvido. Para o processo de modela¸c˜ao do Oxig´enio Dissolvido foram estabelecidos Modelos de Efeitos Mistos (ou Modelos Lineares Generalizados de Efeitos Mistos), pois mostram versatilidade e flexibilidade para a inclus˜ao de efeitos aleat´orios, incorpora¸c˜ao de componentes de tendˆencia e de sazonalidade, de covari´aveis (como o fator hidrometeorol´ogico e outras vari´aveis de Qualidade da ´ Agua de superf´ıcie), bem como da estrutura de correla¸c˜ao temporal pr´opria das s´eries ambientais. Foi efetuado um estudo comparativo dos diversos modelos estabelecidos, considerando crit´erios e m´etricas de qualidade de ajustamento. Palavras-chave: Bacia Hidrogr´afica; rio Douro; Qualidade da ´ Agua; Geoestat´ıstica; Modelos de Efeitos Mistos. vii viii Abstract Environmental degradation is nowadays a critical issue, both due to the difficulty of restoration and rehabilitation and to the serious social and economic consequences. The environmental crisis is partially the result of many man-made mistakes that still remain visible today. Investigations aimed at curbing or estimating environmental problems have led to more in-depth study of methods to better understand the data associated with these problems. This study investigates a problem in the context of surface water quality monitoring in a watershed, and we propose an approach based on spatial and temporal models in order to analyze and evaluate the time series evolution of environmental variables. The data refer to the Douro watershed located in northern Portugal and for the modeling process we considered time series relative to the Dissolved Oxygen (DO) quality variable measured monthly from March 2002 to February 2013. In order to obtain estimates of monthly precipitation values, in area, in the quality sampling stations (where there are no precipitation measurements), we developed a methodology using spatial stochastic processes (Kriging) to be applied to the precipitation data extant in this basin. The estimated values will represent the hydrometeorological factor in the quality sampling stations for the Dissolved Oxygen modeling process. For the Dissolved Oxygen modeling process we established Mixed Effects Models (or Generalized Linear Mixed Effects Models) as they show versatility and flexibility in including random effects, in incorporating trend and seasonality components, covariates (such as the hydrometeorological factor and other surface water quality variables), as well as the temporal correlation structure typical of the environmental series. A comparative study of the various established models was performed considering criteria and quality adjustment metrics. Key-words: Watershed; Douro River; Water quality; Geostatistics; Mixed Effects Models. ix 5.8 Mapa das superf´ıcies dos desvio padr˜ao estimados de precipita¸c˜ao para o mˆes de janeiro, nos anos de 2002 at´e 2013, atrav´es da metodologia Kriging Universal. .................................... 95 5.9 Esquerda: Esta¸c˜ao hidro-qualidade autom´atica, de dezembro de 2003, em Ermida-Corgo (DSRH/INAG). Direita: Sonda de Qualidade da ´ Agua e de n´ıvel danificada (DSRH/INAG). . . . . . . . . . . . . . . . . . . . . . . . . 96 5.10 Representa¸c˜ao da Regi˜ao Hidrogr´afica do Douro e as localiza¸c˜oes das esta- ¸c˜oes de medi¸c˜ao de Qualidade da ´ Agua selecionadas para o estudo. . . . . 96 5.11 Esquerda: Diagramas em caixa de bigodes; Direita: As principais m´etricas da vari´avel OD..................................101 5.12 Representa¸c˜ao gr´afica do OD, em fun¸c˜ao das covari´aveis CBO5, Clorofila, pH,Temperatura,CH eALB.........................102 5.13 Perfil temporal da vari´avel resposta, OD, em fun¸c˜ao do Tempo, no per´ıodo observado.....................................103 5.14 Representa¸c˜ao gr´afica dos intervalos de confian¸ca para os coeficientes relativamente `a constante (esquerda) e `a vari´avel tempo (direita), relativamente ao ajustamento linear para cada de amostragem. . . . . . . . . . . . . . . . 103 5.15 Esquerda: Res´ıduos padronizados do Modelo vs Valores Ajustados; Centro: Valores observados vs Valores Ajustados; Direita: Q-Q plot dos res´ıduos normalizados,Modelo3. ............................111 5.16 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 112 A.1 Diagrama em caixa de bigodes das s´erie de Precipita¸c˜ao, nas 18 esta¸c˜oes de amostragem, no per´ıodo observado. . . . . . . . . . . . . . . . . . . . . 123 A.2 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 124 A.3 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 125 xvi A.4 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 126 A.5 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 127 A.6 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 128 A.7 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 129 A.8 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 130 A.9 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 131 A.10 Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. . 132 A.11 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de fevereiro, nos anos de 2002 at´e 2013. . . . . . . . . . . . . . . . . . . . . . . . . 133 A.12 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de mar¸co, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 A.13 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de abril, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 A.14 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de maio, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 A.15 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de junho, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 A.16 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de julho, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 A.17 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de agosto, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 xvii A.18 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de setembro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . 140 A.19 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de outubro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . . . 141 A.20 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de novembro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . 142 A.21 Representa¸c˜oes das superf´ıcies estimadas da precipita¸c˜ao, no mˆes de dezembro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . 143 A.22 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de fevereiro, nos anos de 2002 at´e 2013. . . . . . . . . . . . . . . . . . . . . 144 A.23 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de mar¸co, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . 145 A.24 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de abril, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . 146 A.25 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de maio, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . 147 A.26 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de junho, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . 148 A.27 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de julho, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . . 149 A.28 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de agosto, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . . 150 A.29 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de setembro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . 151 A.30 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de outubro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . . 152 A.31 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de novembro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . 153 xviii A.32 Representa¸c˜oes das superf´ıcies de erros estimados da precipita¸c˜ao, no mˆes de dezembro, nos anos de 2002 at´e 2012. . . . . . . . . . . . . . . . . . . . 154 B.1 Diagrama em caixa de bigodes das s´erie de Oxig´enio Dissolvido, nas 36 esta¸c˜oes de amostragem, no per´ıodo observado. . . . . . . . . . . . . . . . 155 B.2 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................156 B.3 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................157 B.4 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................158 B.5 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................159 B.6 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................160 B.7 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................161 B.8 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................162 B.9 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................163 xix B.10 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................164 B.11 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................165 B.12 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................166 B.13 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................167 B.14 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................168 B.15 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................169 B.16 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................170 B.17 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................171 B.18 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................172 xx B.19 Representa¸c˜oes gr´aficas das s´eries temporais da Oxig´enio Dissolvido e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. .....................................173 B.20 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 174 B.21 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 175 B.22 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 176 B.23 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 177 B.24 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 178 B.25 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 179 B.26 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 180 B.27 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 181 B.28 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 182 B.29 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 183 B.30 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 184 B.31 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 185 xxi B.32 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 186 B.33 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 187 B.34 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 188 B.35 Representa¸c˜oes da s´erie original e dos os valores estimados, da FAC, da FACP e do Q-Q plot dos res´ıduos do modelo, nas esta¸c˜oes. . . . . . . . . . 189 xxii Lista de Tabelas 5.1 As esta¸c˜oes de amostragem da precipita¸c˜ao, na Bacia Hidrogr´afica do rio Douro, e o respetivo per´ıodo observado. . . . . . . . . . . . . . . . . . . . . 86 5.2 Medidas descritivas da vari´avel da precipita¸c˜ao, no per´ıodo observado. . . . 87 5.3 Valores estimados para cada semivariograma e m´etricas utilizadas para a Valida¸c˜aoCruzada................................ 91 5.4 Localiza¸c˜oes das esta¸c˜oes de amostragem de Qualidade da ´ Agua e respetivo per´ıodo observado, na bacia hidrogr´afica do rio Douro. . . . . . . . . . . . 97 5.5 Vari´aveis em estudo e respetiva descri¸c˜ao. . . . . . . . . . . . . . . . . . . 98 5.6 Medidas descritivas da vari´avel do Oxig´enio Dissolvido, no per´ıodo observado.100 5.7 Medidas Descritivas sobre Oxig´enio Dissolvido, em fun¸c˜ao do mˆes, no per´ıodoobservado..................................101 5.8 Estimativas dos coeficientes da parte fixa, do Modelo Completo 1. . . . . 107 5.9 Estimativas dos coeficientes da parte fixa, do Modelo Completo 2. . . . . 108 5.10 Estimativas dos coeficientes da parte fixa, do Modelo Completo 3. . . . . 108 5.11 Estimativas dos coeficientes da parte fixa, do Modelo Completo 4. . . . . 109 5.12 Teste da Raz˜ao de Verosimilhan¸ca aplicado aos modelos em an´alise. . . . . 109 5.13 Estimativas dos coeficientes dos Modelos 1, 2, 3 e 4. . . . . . . . . . . . . . 110 xxiii xxiv Lista de Siglas/Acr´onimos AIC — Akaike Information Criterion (em portuguˆes, Crit´erio de Informa¸c˜ao de Akaike) APA — Agˆencia Portuguesa do Ambiente AR — Autoregressive (em portuguˆes, Autorregressivo) ARIMA — Autoregressive Integrated Moving Average (em portuguˆes, Autorregressivo Integrado de M´edias M´oveis) ARMA — Autoregressive Moving Average (em portuguˆes, Autorregressivo de M´edias M´oveis) BIC — Bayesian Information Criterion (em portuguˆes, Crit´erio de Informa¸c˜ao Bayesiano) BLE — Best Linear Estimator (em portuguˆes, Melhor Estimador Linear) BLUE — Best Linear Unbiased Estimator (em portuguˆes, Melhor Estimador Linear N˜ao Enviesado) BLUP — Best Linear Unbiased Prediction (em portuguˆes, Melhor Preditor Linear N˜ao Enviesado) CBO5 — Carˆencia Bioqu´ımica de Oxig´enio, a 5 dias CH — Coeficiente Hidrometeorol´ogico DSRH — Dire¸c˜ao dos Servi¸cos de Recursos H´ıdricos EQM — Erro Quadr´atico M´edio EQMN — Erro Quadr´atico M´edio Normalizado ET — Estat´ıstica de Teste xxv Cap´ıtulo 1. Introdu¸c˜ao 4 Cap´ıtulo 2 Geoestat´ıstica 2.1 Interpola¸c˜ao Espacial Os M´etodos de Interpola¸c˜ao Espacial correspondem a procedimentos de estima¸c˜ao do valor de um atributo em locais onde n˜ao est˜ao dispon´ıveis observa¸c˜oes do mesmo, a partir de pontos em que tenham sido registadas observa¸c˜oes do atributo em causa. Assim, a interpola¸c˜ao ´e uma t´ecnica cujo objetivo ´e a estima¸c˜ao de valores desconhecidos de uma fun¸c˜ao, a partir de valores conhecidos da mesma fun¸c˜ao. Quando a informa¸c˜ao dispon´ıvel, proveniente de uma amostra recolhida, n˜ao cobre todo o dom´ınio espacial, a interpola¸c˜ao ´e uma op¸c˜ao para completar os valores em falta. Existem v´arios m´etodos de interpola¸c˜ao, como os M´etodos Determin´ısticos e os M´etodos Estoc´asticos. 2.1.1 M´etodos Determin´ısticos Os M´etodos Determin´ısticos, que continuam a ter uma grande importˆancia e aplica¸c˜ao em ´areas de f´enomenos espaciais, v˜ao ser apresentados de um modo resumido. Os Pol´ıgonos de Thiessen visam a subdivis˜ao do dom´ınio espacial em ´areas de influˆencia (pol´ıgonos de influˆencia) das observa¸c˜oes dispon´ıveis (Thiessen, 1911). Assim, qualquer localiza¸c˜ao no espa¸co tem o valor estimado igual ao valor observado mais pr´oximo, que ´e o do centro do pol´ıgono em que a localiza¸c˜ao est´a contida. Os Pol´ıgonos de Voronoy recorrem a m´etodos que tamb´em consistem na divis˜ao geom´etrica do espa¸co em ´areas de influˆencia (poliedros convexos) e utilizam esta decomposi¸c˜ao para o c´alculo do peso de cada valor observado na interpola¸c˜ao. O M´etodo das M´edias M´oveis estima os valores numa determinada localiza¸c˜ao, pela determina¸c˜ao da m´edia aritm´etica dos valores observados nas localiza¸c˜oes mais pr´oximas. No M´etodo da M´edia Aritm´etica, o valor estimado num local ´e calculado pela m´edia aritm´etica de todas as observa¸c˜oes. 5 Cap´ıtulo 2. Geoestat´ıstica A Interpola¸c˜ao Quadr´atica determina o valor estimado num local a partir da soma ponderada dos valores observados, em que a contribui¸c˜ao de cada valor ´e inversamente proporcional ao quadrado da distˆancia ao ponto a estimar. Na Interpola¸c˜ao Multiquadr´atica, o valor ´e estimado com base na pondera¸c˜ao das distˆancias desse local aos locais de observa¸c˜ao (mais uma constante), em que os pesos s˜ao tais que a superf´ıcie de interpola¸c˜ao obtida passa exatamente pelos valores observados. O M´etodo de Ajustamento de uma Superf´ıcie consiste em ajustar os valores observados a uma superf´ıcie polinomial (splines). As dificuldades principais na aplica¸c˜ao dos M´etodos Determin´ısticos traduzem-se na quantifica¸c˜ao da estrutura espacial da grandeza em estudo e na avalia¸c˜ao da incerteza associada `a caracteriza¸c˜ao do fen´omeno espacial. 2.1.2 M´etodos Estoc´asticos Os M´etodos Estoc´asticos pressup˜oem que os fen´omenos se distribuam no espa¸co de uma forma aleat´oria, com uma determinada estrutura de correla¸c˜ao e, assim, com um grau de incerteza associado aos fen´omenos, resultante da falta de informa¸c˜ao dispon´ıvel. Estes consideram os dados como realiza¸c˜oes de um determinado processo aleat´orio e consistem na modela¸c˜ao da estrutura de varia¸c˜ao do processo e utilizam-na para construir um estimador para os valores n˜ao observados. No contexto de um processo de modela¸c˜ao espacial, dado um conjunto de dados provenientes das amostras experimentais, inicia-se pela conce¸c˜ao de um processo aleat´orio que caracteriza o conjunto de dados, sendo considerada a sele¸c˜ao de um n´umero restrito de parˆametros que, sob determinadas hip´oteses, permitem a inferˆencia espacial. Em 1951, o engenheiro de minas sul-africano, D. G. Krige, desenvolveu um m´etodo para estimar o teor em min´erio de um subsolo a partir de amostras extra´ıdas. Com base nas ideias de D. G. Krige, Matheron (1963) estabelece o termo Geoestat´ıstica em que “a Geoestat´ıstica ´e a aplica¸c˜ao do formalismo das fun¸c˜oes aleat´orias ao reconhecimento e a estima¸c˜ao de fen´omenos naturais”. Ao mesmo tempo que s˜ao desenvolvidas as t´ecnicas geoestat´ısticas na ´area da Engenharia Mineira com G. Matheron, as mesmas ideias s˜ao desenvolvidas na ´area da Meteorologia com L. S. Gandin, na Uni˜ao Sovi´etica, sob o nome de “an´alise objetiva” e “interpola¸c˜ao ´otima” (Lef`evre, 1997). Mercer & Hall (1911) consideram algumas caracter´ısticas da Geoestat´ıstica moderna a dependˆencia espacial, a correla¸c˜ao e o efeito de pepita (a variabilidade `a pequena escala). A representa¸c˜ao da correla¸c˜ao espacial, reconhecida como variograma, ´e desenvolvida por Kolmogorov (1941), assim como o m´etodo de interpola¸c˜ao (Ripley, 1981). Mat´ern (1960) desenvolve algumas fun¸c˜oes que permitem descrever a covariˆancia espacial. Jowett (1955) tamb´em estuda e apresenta algumas fun¸c˜oes, posteriormente denominadas como 6 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais variogramas, que expressam a dependˆencia espacial entre amostras vizinhas. A Geoestat´ıstica permite, tamb´em, fornecer estimativas de erros de estima¸c˜ao, a partir de t´ecnicas num´ericas que caracterizam atributos espaciais (Olea, 2012). A Geoestat´ıstica oferece uma forma de descrever a continuidade espacial dos fen´omenos naturais, adaptando t´ecnicas cl´assicas de regress˜ao, aproveitando a continuidade espacial (Isaacs & Srivastava, 1989). Alguns exemplos s˜ao aplica¸c˜oes relacionadas com a Meteorologia (Cressie & Huang, 1999; Kyriakidis et al., 2001) ou Hidrologia (Goovaerts, 2000). Constata-se um grande desenvolvimento na associa¸c˜ao da dimens˜ao temporal `a dimens˜ao espacial, em diversos autores como, por exemplo, Cressie (1993), Cressie & Wikle (2015), Cressie et al. (2019), Diggle & Giorgi (2019) e Goovaerts (1997). B´ardossy & Pegram (2009) e Gr¨ aler (2014) defendem que a covariˆancia tem um papel preponderante na evolu¸c˜ao da Geoestat´ıstica, atrav´es de campos aleat´orios espa¸cotemporais, o que permite uma maior flexibilidade na modela¸c˜ao dos dados. As implementa¸c˜oes destes m´etodos com ferramentas computacionais s˜ao bastante recentes. Exemplos disso s˜ao a package gstat (Pebesma, 2004) e a package spacetime (Pebesma, 2012), atrav´es da extens˜ao para a Geoestat´ıstica espa¸co-temporal proposta por Gr¨ aler (Pebesma & Heuvelink, 2016), em ambiente R. Cressie & Wikle (2015) explicam e abordam explicitamente estat´ısticas de dados espa¸co-temporais, com exemplos pr´aticos. 2.2 Processos Aleat´orios Considerando os dados como uma s´erie espacial associada a nlocaliza¸c˜oes espaciais {s1, s2, ..., sn}e os valores de uma vari´avel cont´ınua {z(s1), z(s2), ..., z(sn)}, observados nestas localiza¸c˜oes. Cada valor observado z(si), i= 1, . . . , n, ´e considerado como uma realiza¸c˜ao particular de uma determinada vari´avel aleat´oria Z(s), em que svaria numa regi˜ao do espa¸co real de dimens˜ao finita positiva, D⊆Rr, e, usualmente, r= 2,3 (espa¸co real bidimensional ou tridimensional). Este conjunto de vari´aveis aleat´orias (geralmente correlacionadas) ´e denominado por processo aleat´orio, campo aleat´orio ou fun¸c˜ao aleat´oria, sendo definido por {Z(s) : s∈D}(2.1) e tem de satisfazer as condi¸c˜oes de simetria e de consistˆencia (Yaglom, 1962). Um processo aleat´orio {Z(s) : s∈D}´e usualmente caracterizado atrav´es da fun¸c˜ao distribui¸c˜ao cumulativa FZ(s1),...,Z(sn)(s1, . . . , sn) = P(Z(s1)≤z1, . . . , Z(sn)≤zn).(2.2) Para cada s∈D,Z(s) ´e uma vari´avel, assim, define-se como fun¸c˜ao valor m´edio ou 7 Cap´ıtulo 2. Geoestat´ıstica momento de primeira ordem ∀s∈D, E[Z(s)] = µZ(s),(2.3) quando a esperan¸ca existe. Tamb´em se pode especificar a covariˆancia do processo aleat´orio, se existir, e ´e tal que Cov(Z(sj), Z(sk)) = E[(Z(sj)−µZ(sj)(Z(sk)−µZ(sk)] sj, sk∈D, j, k = 1, ..., n, (2.4) em particular, a variˆancia ´e tal que Cov(Z(sj), Z(sj)) = E[(Z(sj)−µZ(sj))2] = V ar[Z(sj)], ∀sj∈D, com j= 1, . . . , n. Na maioria dos casos pr´aticos n˜ao se conhece a lei de probabilidade que define o processo aleat´orio. Deve-se inferir a distribui¸c˜ao ou alguns dos seus momentos, o que requer v´arias realiza¸c˜oes do processo Z(s). Na teoria dos processos aleat´orios, a hip´otese (restri- ¸c˜ao) usual ´e a da estacionaridade (ligada `a no¸c˜ao intuitiva de homogeneidade espacial). Um processo aleat´orio espacial {Z(s) : s∈D}diz-se processo estacion´ario de primeira ordem ou intrinsecamente estacion´ario em Dse para qualquer conjunto de localiza¸c˜oes s1, ..., sn∈D, a distribui¸c˜ao conjunta ´e invariante com respeito a qualquer transla¸c˜ao nas localiza¸c˜oes. Isto ´e, para quaisquer n≥1, h∈Rres1+h, . . . , sn+h∈Das distribui¸c˜oes de (Z(s1+h), ..., Z(sn+h)) e (Z(s1), ..., Z(sn)) s˜ao idˆenticas, ou seja, ∀h∈Rr, Fs1+h,...,sn+h(z1, ..., zn) = Fs1,...,sn(z1, ..., zn).(2.5) Na pr´atica, no entanto, a lei de distribui¸c˜ao n˜ao ´e conhecida, pois os dados s˜ao insuficientes para a inferir. Assim, em algumas situa¸c˜oes ´e desej´avel disponibilizar-se de um conceito de estacionaridade menos restritivo, envolvendo apenas os dois primeiros momentos que s˜ao suficientes para aproximar corretamente a solu¸c˜ao do problema. Um processo aleat´orio {Z(s) : s∈D}diz-se estacion´ario de segunda em D⊆Rr, se ∀s∈D, E[Z(s)] = µZ(s) = µZ(2.6) e, para cada duas vari´aveis aleat´orias Z(u) e Z(v), a fun¸c˜ao covariˆancia existe e apenas depende da diferen¸ca entre uev, i.e., ∀u, v ∈D, Cov(Z(u), Z(v)) = CZ(u−v),(2.7) em que CZdesigna-se por covariograma ou fun¸c˜ao de covariˆancia estacion´aria do processo aleat´orio Z(s). A estacionaridade da covariˆancia implica que a variˆancia V ar[Z(s)] existe e n˜ao de8 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais pende de s(implica a estacionaridade da variˆancia). Isto ´e, V ar[Z(s)] = CZ(0), ∀s∈D. Em particular, se o processo espacial considerado ´e tal que CZ(0) >0 e ´e estacion´ario de segunda ordem, a fun¸c˜ao CZ(.) ´e designada por correlograma (ou fun¸c˜ao de correla¸c˜ao estacion´aria) e denominada por ρZ(.), tal que ∀s, u ∈D, ρ(s−u) = CZ(s−u) CZ(0) ∈[−1,1].(2.8) O correlograma, bem como o covariograma, modela a estrutura de dependˆencia espacial do processo aleat´orio (estacion´ario de segunda ordem). Um processo aleat´orio {Z(s) : s∈D}diz-se intrinsecamente estacion´ario (ou de estacionaridade intr´ınseca) em D⊆Rrse o valor m´edio do processo existe e ´e constante em D, isto ´e, ∀s∈D, E[Z(s)] = µZ(s) = µZ,(2.9) em que a variˆancia V ar[Z(u)−Z(v)] existe para todo u, v ∈De depende apenas da diferen¸ca u−v, isto ´e, ∀u, v ∈D, V ar[Z(u)−Z(v)] = E[(Z(u)−Z(v))2] = 2γZ(u−v) (2.10) em que se designa 2γZa fun¸c˜ao variograma e γZ´e denominada fun¸c˜ao semi-variograma do processo Z(.) (salvaguarda-se a existˆencia de autores que definem outra nomenclatura). A defini¸c˜ao de variograma como a variˆancia dos acr´escimos (incrementos) espaciais de um processo aleat´orio faz com que se verifiquem algumas propriedades. Se γZ(.) ´e o semivariograma de um processo aleat´orio intrinsecamente estacion´ario Z(.), ent˜ao ∀n≥1,∀λ1, . . . , λn∈R,∀s1, ..., sn∈Rr, n X i=1 n X j=1 λiλjCz(si−sj)≥0.(2.11) Al´em disso, se o covariograma CZresulta dum processo aleat´orio estacion´ario de segunda ordem Z(s), tem-se que ∀s∈D, CZ(0) = V ar(Z(s)) ≥0,(2.12) ∀u, v ∈D, CZ(u−v) = CZ(v−u) (2.13) e, pela desigualdade de Cauchy-Schwarz, tem-se que ∀u, v ∈D, |CZ(u−v)| ≤ CZ(0).(2.14) 9 Cap´ıtulo 2. Geoestat´ıstica No caso do processo aleat´orio com apenas estacionaridade intr´ınseca, o semivariograma existe mas o covariograma pode n˜ao existir. Se o semivariograma γZ(.) de um processo intrinsecamente estacion´ario Z(s), ent˜ao apresenta algumas propriedades, tais como γZ(0) = 0,(2.15) ∀u, v ∈D, γZ(u−v) = γZ(v−u),(2.16) ∀u, v ∈D, u 6=v, γZ(u−v)>0 (2.17) e ∀u, v ∈D, lim ||u−v||→∞ γZ(u−v) = c0,(2.18) designado por efeito de pepita. A raz˜ao do crescimento de um semivariograma de um processo aleat´orio Z(s) ´e dada por lim ||u−v||→∞ γZ(u−v) ||u−v||2(2.19) e pode ser indicador se o processo ´e intrinsecamente estacion´ario ou n˜ao. Caso o semivariograma γZ(.) tiver um crescimento mais lento que ||u−v||2, ent˜ao o limite anterior tende para zero e o processo Z(.) ´e intrinsecamente estacion´ario. Se o crescimento de γZ(.) for mais r´apido que ||u−v||2, ent˜ao a hip´otese intr´ınseca n˜ao ´e v´alida. Se a estacionaridade do processo aleat´orio ´e de segunda ordem, ent˜ao o variograma e o covariograma existem e s˜ao estruturalmente equivalentes, cuja rela¸c˜ao estrutural ´e dada por CZ(u−v) = CZ(0) −γZ(u−v)⇔γZ(u−v) = CZ(0) −CZ(u−v) (2.20) e se for verificado que lim||u−v||→∞ CZ(u−v) = 0, ent˜ao lim ||u−v||→∞ γZ(u−v)) = CZ(0) = V ar[Z(.)],(2.21) em que CZ(0) ´e designado por patamar do semivariograma. O patamar parcial ´e determinado pela diferen¸ca entre o patamar e o efeito de pepita, CZ(0) −C0. Um processo aleat´orio {Z(s), s ∈D}pode ser descrito como uma combina¸c˜ao de processos aleat´orios n˜ao correlacionados {Zi(s), s ∈D}, com i= 1, . . . , k, ou seja, formalmente ∀λ1, . . . , λk∈R+ 0, Z(s) = k X i=1 λiZi(s).(2.22) Se Zi(s) ´e um processo estacion´ario de segunda ordem com covariograma CZi(s), para 10 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais cada i= 1, . . . , k, ent˜ao Z(s) tamb´em ´e estacion´ario de segunda ordem com covariograma dado por CZ(u−v) = Pk i=1 λ2 iCZi(u−v), quaisquer que sejam u, v ∈D. Se Zi(s) ´e um processo intrinsecamente estacion´ario, com semivariograma γZi(.), para cada i= 1, . . . , k, ent˜ao Z(s) ´e tamb´em intrinsecamente estacion´ario com semivariograma dado por γZ(u−v) = Pk i=1 λ2 iγZi(u−v), quaisquer que sejam u, v ∈D. 2.3 Continuidade Espacial Habitualmente, o estudo de uma caracter´ıstica revela que se a distˆancia entre dois valores for reduzida (pontos pr´oximos), ent˜ao s˜ao mais semelhantes do que aqueles com distˆancias maiores (pontos afastados). A continuidade espacial, em maior ou menor grau, permite a descri¸c˜ao do quanto os valores se dispersam espacialmente, e de que forma variam com as diferentes dire¸c˜oes do espa¸co (anisotropia). Os conceitos que ser˜ao abordados visam descrever e quantificar a continuidade espacial. 2.3.1 Variograma, Covariograma e Correlograma Na Geoestat´ıstica, as fun¸c˜oes para a modela¸c˜ao da dependˆencia espacial e/ou temporal recorrentemente utilizadas s˜ao o variograma, o covariograma e o correlograma. O variograma pode ser classificado de acordo com a sua natureza: o variograma emp´ırico ´e resultante do conjunto de observa¸c˜oes da amostra em estudo, o variograma te´orico ´e o modelo de variograma de referˆencia e o variograma verdadeiro ´e o variograma real e ´e desconhecido. De forma simplista, o variograma quantifica a dispers˜ao natural das vari´aveis e a variabilidade espacial entre pares de valores separados por uma distˆancia previamente estabelecida ||d||, em que d=si−sj. Existem v´arios m´etodos para o c´alculo do semivariograma. Por exemplo, pelo m´etodo dos momentos, o semivariograma ´e calculado pela m´edia aritm´etica do quadrado das diferen¸cas de todos os pares de pontos que est˜ao separados de um vetor d(Matheron, 1963), tal que ˆγZ(d) = 1 2|N(d)|X (i,j)∈N(d)Z(si)−Z(sj)2,(2.23) em que N(d) = (i, j) : si−sj=h, i, j ∈1, . . . , n, e|N(d)|= #N(d) ´e o n´umero de pares de pontos ||d|| distanciados e alinhados segundo a dire¸c˜ao do vetor d. A representa- ¸c˜ao dos pares de valores ||d||,ˆγz(d)num sistema de eixos representa o semivariograma experimental. A aplica¸c˜ao destes m´etodos tem como pressuposto a malha amostral ser regular. No caso contr´ario, deve-se proceder a uma regulariza¸c˜ao angular e por classes de 11 Cap´ıtulo 2. Geoestat´ıstica distˆancias. Um processo aleat´orio designa-se por isotr´opico se o respetivo semivariograma ou o covariograma depender do vetor dapenas na sua norma, n˜ao podendo depender da dire¸c˜ao angular desse mesmo vetor. Um processo espacial intrinsecamente estacion´ario diz-se isotr´opico quando o seu variograma γZ(d) ´e fun¸c˜ao apenas de ||d||, ou seja ∀d, γZ(d) = γZ(||d||).(2.24) O processo que n˜ao verifique a condi¸c˜ao supracitada ´e denominado como anisotr´opico (depende de ||d|| e da dire¸c˜ao de d,dir), ou seja, ∀d, γZ(d) = γZ(||d||, dir).(2.25) A anisotropia pode ser entendida como a variabilidade espacial dependente das dire¸c˜oes do espa¸co. A modela¸c˜ao de fen´omenos isotr´opicos tem como objetivo reduzir as estruturas de continuidade das diferentes dire¸c˜oes a um s´o modelo. Esta modela¸c˜ao ´e conseguida geralmente atrav´es de um conjunto de transformadas geom´etricas do sistema de coordenadas, de modo a que os diferentes semivariogramas nas diferentes dire¸c˜oes sejam equivalentes a um mesmo modelo ou representando separadamente cada um das variabilidades direcionais consideradas. Os dois modelos mais comuns de anisotropia s˜ao a anisotropia geom´etrica e a anisotropia zonal. Diz-se que um processo espacial intrinsecamente estacion´ario {Z(s) : s∈D} exibe anisotropia geom´etrica se tal anisotropia pode ser reduzida a uma isotropia atrav´es de uma transforma¸c˜ao linear das coordenadas, ou seja, se existir uma matriz invert´ıvel Ab×btal que Z(As) ´e isotr´opico. A anisotropia zonal corresponde ao caso em que os semivariogramas ajustados nas diferentes dire¸c˜oes apresentam diferentes caracter´ısticas de variabilidade (diferentes patamares), podendo ter valores de amplitude tamb´em diferentes (amplitude ´e o valor de ||d|| para o qual o semivariograma se estabiliza, isto ´e, γZ(||d||) = C1, onde C1´e uma constante). Mais pormenores sobre anisotropia podem ser consultados em Soares (2000). O objetivo da corre¸c˜ao de anisotropia ´e obter um ´unico semivariograma isotr´opico que possa modelar a variabilidade espacial do fen´omeno em estudo. A isotropia ´e muito importante porque ´e f´acil de interpretar e esta caracter´ıstica ajuda na compreens˜ao do processo e na interpreta¸c˜ao do modelo e, al´em disso, reduz a carga dos c´alculos computacionais. 2.3.2 Modelos Te´oricos de Semivariogramas Os modelos te´oricos de semivariogramas definidos positivos (modelos de transi¸c˜ao) utilizam um n´umero restrito de fun¸c˜oes ou combina¸c˜oes de fun¸c˜oes, definidas positivas, de 12 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais forma a satisfazerem as condi¸c˜oes de positividade, para interpolar os valores experimentais dos variogramas. A utiliza¸c˜ao de um n´umero restrito de fun¸c˜oes, definidas positivas, para interpolar os valores experimentais dos variogramas ´e uma das poss´ıveis formas de satisfazer as condi¸c˜oes de positividade. Estes modelos s˜ao independentes da dire¸c˜ao, is´otropos, simples e s˜ao fun¸c˜oes do escalar ||d||, sendo classificados em modelos de transi¸c˜ao e em modelos n˜ao estacion´arios. Se CZk(.), ∀k∈N, ´e um covariograma v´alido em Rre limk→∞ CZk(d) = CZ(d), ∀d∈Rr, ent˜ao CZ(.) ´e um covariograma v´alido em Rr. A constru¸c˜ao de semivariogramas/covariogramas v´alidos e o estabelecimento de condi¸c˜oes necess´arias e suficientes para que a fun¸c˜ao seja um semivariograma/covariograma v´alido ´e um tema extensamente estudado, embora neste trabalho se limite a expor os modelos de semivariogramas mais cl´assicos. Num processo aleat´orio intrinsecamente estacion´ario, {Z(s), s ∈D}, com semivariograma γZ(.), se lim||d||→∞ γZ(d) = CZ(0) = σ2 Z<∞, ent˜ao o processo aleat´orio Z(.) designa-se por fen´omeno de transi¸c˜ao e a γZ(.) denomina-se por modelo de transi¸c˜ao. σ2 Z+τ2 Z´e, nesta defini¸c˜ao, o valor do patamar do semivariograma γZ(.) e ´e caracterizado pela altura m´axima, atingida pela curva do semivariograma. O menor valor de ||φ|| para o qual γZ(φ(1 + )) = CZ(0) = σ2 Z´e denominado por amplitude do semivariograma na dire¸c˜ao φ ||φ|| (φ∈R). Quando o semivariograma do processo aleat´orio Z(.) possui um valor do patamar (o semivariograma ´e limitado), este ´e o valor da variˆancia do processo Z(.) e Z(.) ´e tamb´em um processo estacion´ario de segunda ordem. Assim, um fen´omeno de transi¸c˜ao corresponde a um processo estacion´ario de segunda ordem. Neste caso, o covariograma do processo aleat´orio Z(.) possui tamb´em um valor do patamar que ´e igual a zero e no caso de existˆencia de amplitude ||φ|| ´e obviamente a mesma para o semivariograma e semicovariograma. V˜ao ser apresentados os principais modelos de transi¸c˜ao que abragem a generalidade das situa¸c˜oes de dispers˜ao de fen´omenos espaciais nas Ciˆencias do Ambiente, Figura 2.1. Modelo de Efeito de Pepita No caso de vari´aveis cont´ınuas ´e expect´avel que o variograma passe na origem, contudo na maioria dos casos tal n˜ao se verifica. Esta descontinuidade representa as varia¸c˜oes locais ou de pequena escala, como erros de amostragem ou de an´alise (a variˆancia dos erros de an´alise contribui para este valor) e a poss´ıvel existˆencia de micro regionaliza¸c˜oes desenvolvendo-se a uma escala n˜ao detet´avel pela escala de amostragem adotada. 13 Cap´ıtulo 2. Geoestat´ıstica quadr´atico m´edio ´e igual `a sua variˆancia, isto ´e, EQM(λZ(s)) = EλZ(s)−Z(s0)2 =V arλZ(s)−Z(s0) =λTΣλ−2λTC0+V ar[Z(s0)], (2.43) em que Σ´e a matriz de covariˆancia de Z(s) e C0´e o vetor de dimens˜ao n. Assim, estabelecem-se as seguintes condi¸c˜oes n X i=1 λi= 1 (2.44) e λTΣλ−2λTC0´e m´ınimo.(2.45) Para a obten¸c˜ao da minimiza¸c˜ao da express˜ao (2.45), recorre-se ao multiplicador de Lagrange, 5, e tem-se que λTΣλ−2λTC0−2(1Tλ−1)5 dλ= 0 ⇒2Σλ−2C0−215= 0 ⇒Σλ−C0−15= 0 ⇒Σλ−15=C0 ⇒Σλ=C0+15 ⇒λ=Σ−1C0+Σ−115 ⇒1Tλ=1TΣ−1C0+1TΣ−115. (2.46) Pela condi¸c˜ao 1Tλ=Pn i=1 λi= 1, tem-se que 1TΣ−1C0+1TΣ−115= 1 ⇒ 5 =1−1TΣ−1C0 1TΣ−11 .(2.47) O vetor λ´e estimado atrav´es de ˆ λ=Σ−1C0+Σ−111−1TΣ−1C0 1TΣ−11(2.48) 20 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais e o estimador para Z(s0) obtido por Kriging Ordin´ario ´e tal que ˆ ZKO(s0) = CT 0Σ−1Z(s0) + 1−1TΣ−1C0 1TΣ−111TΣ−1Z(s0),(2.49) com a variˆancia do erro de estima¸c˜ao de Kriging dada por σ2 KO =V ar[ˆ ZKO −Z(s0)] =λTΣλ−2λTC0+V ar[Z(s0)] =λTC0+λT15−2λTC0+V ar[Z(s0)]. (2.50) Sabe-se que λT1= 1, ent˜ao o estimador pode ser escrito atrav´es de σ2 KO(s0) = −λTC0+5+V ar[Z(s0)].(2.51) Considerando o processo estacion´ario de segunda ordem (intrinsecamente estacion´ario), pode-se adaptar o m´etodo de estima¸c˜ao de Kriging Ordin´ario e exprimir as equa¸c˜oes de covariˆancias, em termos do semivariograma, Λij =γZ(si, sj),(2.52) em que Λ´e a matriz de semivariogramas de Z(s), dimens˜ao n×n, com i, j = 1, . . . , n, υ0i=γZ(s0, si),(2.53) em que υ0i´e o vetor de semivariogramas entre Z(s0) e Z(si), e 1´e o vetor composto por valores unit´arios, de dimens˜ao n. Atrav´es da rela¸c˜ao γZ(d) = CZ(0) −CZ(d), tem-se que V ar[Z(s0)] = CZ(0), C0=CZ(0)1−υ0eΣ=CZ(0)11T−Λ. O erro quadr´atico m´edio pode ser expresso EQMλTZ(s)=λT[CZ(0)11T−Λ]λ−2λT[CZ(0)1−υ0] + CZ(0) =CZ(0)λT11Tλ−λTΛλ−2λTCZ(0) + 2λTυ0+CZ(0),(2.54) tendo em conta que λ1= 1, tem-se que EQMλTZ(s)= 2λTυ0+λTΛλ.(2.55) As condi¸c˜oes apresentadas anteriormente podem ser reescritas n X i=1 λi= 1 (2.56) 21 Cap´ıtulo 2. Geoestat´ıstica e 2λTυ0+λTΛλ´e m´ınimo.(2.57) Para a minimiza¸c˜ao de (2.57), utilizando o multiplicador de Lagrange, 5, 2λTυ0−λTΛλ−2(1Tλ−1)5 dλ= 0 ⇒2υ0−2Λλ−21T5= 0 ⇒υ0−Λλ−1T5= 0 ⇒Λλ=υ0−1T5 ⇒λ=Λ−1υ0−Λ−11T5 ⇒1Tλ=1TΛ−1υ0−1TΛ−115. (2.58) Pela condi¸c˜ao 1Tλ=Pn i=1 λi= 1 tem-se que 1Tλ=1TΛ−1υ0−1TΛ−115= 1 ⇒ 5 =−1−1TΛ−1υ0 1TΛ−11. (2.59) O vetor λ´e estimado atrav´es de ˆ λ=Λ−1υ0+Λ−111−1TΛ−1υ0 1TΛ−11,(2.60) o estimador para Z(s0) obtido por Kriging Ordin´ario ´e tal que ˆ ZKO =υ0Λ−1Z(s) + 1−1TΛ−1υ0 1TΛ−111TΛ−1Z(s)(2.61) e a variˆancia do erro de estima¸c˜ao ´e dada por ˆσ2 KO =V arˆ ZKO(s0)−Z(s0) =EQMˆ ZKO(s0) = 2λTυ0−λTΣλ = 2λTυ0−λT(υ0−15) = 2λTυ0−λTυ0+λT15 =λTυ0+5 =υ0TΛ−1υ0−(1 −1TΛ−1υ0)2 1TΛ−11. (2.62) 22 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Propriedades do Estimador Nos m´etodos apresentados, a variˆancia do erro de estima¸c˜ao n˜ao depende apenas da esperan¸ca, nem das observa¸c˜oes, mas depende do segundo momento do processo aleat´orio. Ent˜ao, ´e poss´ıvel conhecer a qualidade do estimador antes de se observar o processo. A estimativa do estimador nos pontos observados ´e igual `a fun¸c˜ao nesses pontos, ˆ Z(si) = Z(si), i= 1, . . . , n, ou seja, ´e um interpolador exato. Considerando que Z(u) e Z(w) s˜ao vari´aveis aleat´orias independentes, com u6=w, e V ar[Z(s)] = σ2 Z, n˜ao depende de s, ou seja, {Z(s), s ∈D}´e um processo de ru´ıdo branco. Ent˜ao, a estima¸c˜ao em qualquer ponto s0´e a m´edia aritm´etica dos valores observados, do processo Z(.), ou seja, ˆ Z(s0) = 1 n n X i=1 Z(si) = ¯ Z, ∀s06=si, i = 1, . . . , n, (2.63) para al´em disso, tem-se que ˆ λ=11 ne5=σ2 Z n(2.64) e a variˆancia de estima¸c˜ao ´e tal que σ2 KO(s0) = V ar[ˆ Z(s0)−Z(s0)] = σ2 Z+σ2 Z n.(2.65) Das propriedades supracitadas pode-se concluir que a superf´ıcie dos ˆ Z(.) pode n˜ao ser cont´ınua. Considerando o modelo Z(s) = m(s) + Y(s), s∈D,Y(s) ´e um processo aleat´orio de covariˆancia estacion´aria, m´edia zero e fun¸c˜ao covariˆancia conhecida. No caso de m(s) = mconhecido, ent˜ao ˆ Y(s0) = CT 0Σ−1Y(s) =CT 0Σ−1Z(s)−1m =CT 0Σ−1Z(s)−CT 0Σ−11m (2.66) e, consequentemente, o estimador ´e dado por ˆ Z(s0) = m+ˆ Y(s0) =CT 0Σ−1Z(s) + (1 −CT 0Σ−11)m, (2.67) 23 Cap´ıtulo 2. Geoestat´ıstica cujo erro quadr´atico m´edio determina-se por V ar[ˆ Z(s0)−Z(s0)] = V arCT 0Σ−1Z(s)−Z(s0).(2.68) No caso de m(s) = mdesconhecido, ent˜ao ˆm=1TΣ−11−11TΣ−1Z(s),(2.69) com variˆancia V ar[ ˆm] = 1TΣ−11−1.(2.70) Desta forma, ˆ Z∗(s0) = CT 0Σ−1Z(s) + 1−CT 0Σ−111TΣ−11−11TΣ−1Z(s),(2.71) cuja variˆancia ´e determinada por V arˆ Z(s0)−Z(s0)=V arˆ Z∗(s0)−Z(s0)+1−1TΣ−1C02V ar[ ˆm],(2.72) a parcela da equa¸c˜ao reflete a perda de precis˜ao quando se estima m. Outros M´etodos de Estima¸c˜ao A restri¸c˜ao de estacionariedade na m´edia desconhecida, mas constante para todo o dom´ınio, nem sempre ´e f´acil de cumprir, nomeadamente, nos fen´omenos n˜ao estacion´arios e, neste caso, ´e utilizado o Kriging Universal (ou Kriging com deriva externa). OKriging ´e um m´etodo que permite a estima¸c˜ao de um processo {Z(s), s ∈D}, D⊆Rr, utilizando os valores observados deste processo, nas localiza¸c˜oes s1, s2, . . . , sn. Por´em, ´e poss´ıvel melhorar a qualidade da estima¸c˜ao atrav´es de observa¸c˜oes de outro processo, que tem em conta observa¸c˜oes de processos que est˜ao correlacionados com o que se pretende estudar. Neste sentido, ´e poss´ıvel aplicar um m´etodo de estima¸c˜ao de um processo {Z1(s), s ∈D},D⊆Rr, recorrendo aos valores observados, nas localiza¸c˜oes s1, s2, . . . , sn, e a outros processos, como Z2(s), s ∈D,D⊆Rr,Z3(s), s ∈D,D⊆Rr, entre outros. A t´ecnica ´e uma solu¸c˜ao poss´ıvel quando se tem um n´umero reduzido de observa¸c˜oes, no processo em an´alise, e outros processos tˆem um maior n´umero de observa¸c˜oes, designando-se por Cokriging (que n˜ao vai ser aplicado nesta disserta¸c˜ao). 24 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais 2.5 Estima¸c˜ao Global Nas Sec¸c˜oes anteriores descreveram-se os conceitos b´asicos para o m´etodo de estima- ¸c˜ao de Kriging pontual (de um valor do atributo do processo). Mas, quando se pretende a estima¸c˜ao do valor m´edio Z(.) numa determinada regi˜ao A⊆D, em que A´e um subconjunto com volume r-dimensional estritamente positivo (|A|>0), o m´etodo denomina-se por Kriging Global. Este pode ser obtido pela m´edia dos valores pontuais estimados pelo Kriging que comp˜oem A ou pode ser estimado diretamente. O valor m´edio do processo aleat´orio Z(.), numa determinada regi˜ao A´e dado por Z(A) = 1 |A|ZA Z(v)dv. (2.73) De acordo com a teoria de processos aleat´orios, define-se o integral de um processo aleat´orio de uma vari´avel aleat´oria como o limite de uma soma de Riemann, tendo como suporte a defini¸c˜ao supracitada, deduzindo-se as express˜oes seguintes EZ(A)=1 |A|ZA E[Z(v)]dv, (2.74) V arZ(A)=CovZ(A), Z(A)=1 |A|2ZAZA CovZ(v), Z(w)dw dv, (2.75) Cov(Z(A), Z(v)) = 1 |A|ZA CovZ(v), Z(u)du, ∀v∈A. (2.76) Na presen¸ca de um processo intrinsecamente estacion´ario, com semivariograma γZ(.), podem-se escrever as equa¸c˜oes anteriores como E[Z(A)] = µZ,(2.77) V arZ(A)=1 |A|ZAZA V arZ(v)dv −1 |A|2ZAZA γZ(v−w)dwdv (2.78) e 1 2V arZ(A)−Z(v)=ZA γZ(v−w)dw −1 |A|2ZAZA γZ(w−u)dudv, ∀v∈A. (2.79) Devido `a complexidade das express˜oes e dispˆendio no c´alculo num´erico, usualmente utilizam-se aproxima¸c˜oes. Neste sentido, numa regi˜ao Acom uma amostra finita de valores de atributos, localizados nessa regi˜ao, os integrais das express˜oes anteriores podem ser substitu´ıdos por m´edias e gerando, desta forma, aproxima¸c˜oes 25 Cap´ıtulo 2. Geoestat´ıstica V ar[Z(B)] ≃1 n2 n X i=1 n X j=1 CovZ(si)−Z(sj)(2.80) ou V ar[Z(B)] ≃1 n n X i=1 V arZ(si)−1 n2 n X i=1 n X j=1 γZ(Z(si)−Z(sj)),(2.81) Cov[Z(B), Z(s)] ≃1 n n X i=1 CovZ(si)−Z(s),∀s∈B, (2.82) e 1 2V arZ(B)−Z(v)=1 n n X i=1 γZ(si−s)−1 2n2 n X i=1 n X j=1 γZ(si−sj),∀v∈A. (2.83) Na estima¸c˜ao global, o estimador ˆ Z(A) , A⊆D, obtido a partir das observa¸c˜oes pontuais, ´e semelhante ao descrito na estima¸c˜ao pontual. O estimador de Kriging Ordin´ario ´e dado por ˆ ZKO(A) = n X i=1 λA,iZ(si) = λT AZ(s),(2.84) em que λT A= (λA,1, λA,2, . . . , λA,n), Z(s)=(Z(s1), Z(s2), . . . , Z(sn)) ´e o vetor das observa¸c˜oes do processo. Relativamente `as covariˆancias, um processo estacion´ario de segunda ordem, em D, tem variˆancia finita e a fun¸c˜ao covariˆancia existe e o vetor λ´e dado por ˆ λA=Σ−1CA+Σ−111−1TΣ−1cA 1TΣ−11,(2.85) em que Σ´e a matriz de covariˆancias de Z(s), CA´e o vetor de covariˆancias entre Z(B) e Z(si), com i= 1, . . . , n, tal que CAi=Cov(Z(B),Z(si)), 1´e o vetor de valores unit´arios. O estimador de Kriging Ordin´ario ´e dado por ˆ ZKO(A) = CT AΣ−1Z(s) + 1−1TΣ−1CA 1TΣ−111TΣ−1Z(s) (2.86) e a sua variˆancia do erro de estima¸c˜ao ´e dada por ˆσKO(A) = −λT ACA+5+V arZ(A).(2.87) As equa¸c˜oes de covariˆancia podem ser formuladas, em termos de semivariograma. Considerando um processo estacion´ario de segunda ordem, intrinsecamente estacion´ario, o m´etodo de estima¸c˜ao de Kriging Ordin´ario ´e caraterizado pela matriz de semivariogramas 26 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais de Z(s), tal que Λij =γZ(si, sj), i, j = 1, . . . , n, (2.88) pelo vetor de semivariogramas entre Z(B) e Z(si), i= 1, . . . , n, υA=γZ(A, si), i = 1, . . . , n, (2.89) e pelo vetor composto por nvalores unit´arios. A variˆancia do erro da estima¸c˜ao ´e dada por ˆσ2 KO =V arˆ ZKO(A)−Z(A) =EQM(ˆ ZKO(A)) =V arhn X i=1 λAiZ(si)−Z(B)i = 2 n X i=1 λA,i 1 2V arZ(A)−Z(si)− n X i=1 n X j=1 λAjλAi 1 2V arZ(sj)−Z(si) (2.90) e, sob a forma matricial, ˆσ2 KO = 2λT AυA−λT AΛλA.(2.91) Assim, as condi¸c˜oes impostas anteriormente reformulam-se para as seguintes n X i=1 λA,i = 1 (2.92) e 2λT AυA+λT AΛλA´e m´ınimo.(2.93) Recorrendo ao multiplicador de Lagrange, tem-se que ˆ λA=Λ−1υA+Λ−111−1TΛ−1υA 1TΛ−11,(2.94) em que 5A=−1−1TΛ−1υA 1TΛ−11,(2.95) o estimador de Kriging Ordin´ario de Z(A) ´e ˆ Z(A) = υT AΛ−1Z(s) + 1−1TΛ−1υA 1TΛ−111TΛ−1Z(s) (2.96) 27 Cap´ıtulo 2. Geoestat´ıstica e a variˆancia do erro de estima¸c˜ao ´e ˆσ2 KO = 2λT AυA−λT AΛλA =λT AυA+5A =λT AυA−1−1TΛ−1υA 1TΛ−11. (2.97) Uma abordagem alternativa ao m´etodo exposto para a estima¸c˜ao global consiste na determina¸c˜ao do estimador pontual para todos os pontos e na obten¸c˜ao da m´edia dos valores estimados obtidos. Sob a hip´otese de estacionariedade na regi˜ao A, pode-se provar que estas duas metodologias s˜ao equivalentes. A abordagem utiliza uma grelha de pontos para a aproxima¸c˜ao num´erica do vetor υAe determina o estimador global, com base nessa aproxima¸c˜ao, o que ´e equivalente a calcular todos os estimadores pontuais na grelha e o valor estimado para Z(A) ser considerado como a m´edia de todos os estimadores pontuais. No entanto, as variˆancias dos erros de estima¸c˜ao pontuais e globais n˜ao tˆem uma rela¸c˜ao t˜ao simplista, porque a variˆancia n˜ao se traduz em combina¸c˜oes lineares. 2.6 Valida¸c˜ao Cruzada O m´etodo de estima¸c˜ao requer que sejam validados os modelos de variograma e as hip´oteses de homogeneidade espacial para todo o campo em an´alise, prevenindo dificuldades do ajustamento, como a presen¸ca de valores discrepantes (outliers). Um m´etodo bastante utilizado ´e a Valida¸c˜ao Cruzada (cross-validation), que utiliza a informa¸c˜ao dispon´ıvel e compara os valores observados e os estimados, nas localiza¸c˜oes das observa¸c˜oes. Com os valores reais e os valores estimados (Z(si),ˆ Z(si)), i= 1, . . . , n, pode-se determinar estat´ısticas das distribui¸c˜oes univariadas das estimativas e dos erros, para se determinar a qualidade do modelo do variograma adotado. O procedimento inicia-se com o c´alculo do preditor, atrav´es de Kriging, e a sua variˆancia de precis˜ao, ou seja, determina-se ˆ Z(si), i= 1, . . . , n, em fun¸c˜ao das vari´aveis Z(s1), . . . , Z(si−1), Z(si+1), . . . , Z(sn), i= 1, . . . , n, atrav´es do semivariograma ajustado e do valor da variˆancia da estima¸c˜ao. Posteriormente, para cada i, calcula-se a diferen¸ca normalizada entre o estimado e o observado, ou seja, calcula-se Z(si)−ˆ Z(si) σ(si). Para a avalia¸c˜ao da qualidade do ajustamento, utilizam-se os histogramas das diferen¸cas normalizadas para dete¸c˜ao de valores discordantes, o valor m´edio das diferen¸cas normalizadas, que, no cen´ario ideal, ´e zero ou 28 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais pr´oximo de zero. Assim, a m´edia do erro de predi¸c˜ao ´e dada por MEP =1 n n X i=1 Z(si)−ˆ Z(si) σ(si)(2.98) e o erro quadr´atico m´edio normalizado, que no cen´ario ideal ´e um ou pr´oximo de um, ´e dado por EQMN =1 n n X i=1 Z(si)−ˆ Z(si) σ(si)2.(2.99) O processo de estima¸c˜ao da Valida¸c˜ao Cruzada ´e influenciado por trˆes fatores, que est˜ao bastante relacionados, que tornam dif´ıcil a perce¸c˜ao da sua influˆencia nos valores das estat´ısticas globais dos desvios: as hip´oteses de estacionaridade/homogeneidade espacial, o modelo de variograma (que se pretende validar) e o pr´oprio processo de estima¸c˜ao. Antes de verificar qual o variograma mais apropriado, deve ser primeiramente validada a hip´otese de estacionariedade/ homogeneidade espacial, em segundo lugar h´a que julgar se os desvios n˜ao tˆem a ver com o tipo de estimador usado. O imbricamento destes fatores torna qualquer destas an´alises extremamente dif´ıcil. Se se considerar uma pequena ´area com elevada variabilidade local, ent˜ao v˜ao ser gerados grandes desvios entre os valores reais e os estimados. Esta dificuldade pode ser ultrapassada pelo aumento do Efeito de Pepita e, assim, conseguir uma atenua¸c˜ao dos desvios locais, com base no aumento da m´edia global, relativamente aos valores locais. No entanto, um variograma que tem boas m´etricas na Valida¸c˜ao Cruzada n˜ao ´e condi¸c˜ao suficiente para que o modelo seja o mais adequado. 29 Cap´ıtulo 3. S´eries Temporais e a autocorrela¸c˜ao ρkestimada por ˆρk=ˆγk ˆγ0 =Pn−k t=1 (Yt−¯ Y)(Yt+k−¯ Y) Pn t=1(Yt−¯ Y)2.(3.21) Para al´em da correla¸c˜ao global, tamb´em se pode recorrer `a correla¸c˜ao parcial entre Yt eYt+k, quando s˜ao fixadas as vari´aveis interm´edias Yt+1, . . . , Yt+k−1.´ E poss´ıvel obter a correla¸c˜ao parcial com base na regress˜ao linear, tal que Yt+k=φk1Yt+k−1+···+φkkYt+t+k,(3.22) em que φkj s˜ao os coeficientes do modelo, com j= 1, . . . , k, uma vez que os erros seguem uma distribui¸c˜ao Normal. O valor de φkk ´e o coeficiente de correla¸c˜ao do Modelo de Regress˜ao Linear, onde {t, t ∈Z}, s˜ao independentes e Gaussianos, com m´edia nula e variˆancia σ2. Este coeficiente exprime a varia¸c˜ao entre tet+k, quando os restantes coeficientes s˜ao constantes. A partir da equa¸c˜ao (3.22), multiplicando por Yt+k−j,j= 1, . . . , k, aplicando o valor esperado e dividindo por ρ0, tem-se que ρj=φk1ρj−1+···+φkkYk−jj= 1, . . . , k, (3.23) resolvendo em ordem a φkj, recorrendo `a regra de Cramer, consegue-se obter a fun¸c˜ao de autocorrela¸c˜ao parcial, φkk. A fun¸c˜ao de autocorrela¸c˜ao parcial pode-se definir por φkk =Cor[Yt, Yt+k|Yt+1, Yt+2, . . . , Yt+k−1] = |P∗ k| |Pk|,(3.24) em que Pk´e a matriz k×kdada por Pk=          1ρ1ρ2. . . ρk−1 ρ11ρ1. . . ρk−2 ρ2ρ11. . . ρk−3 . . .. . .. . ..... . . ρk−1ρk−2ρk−3. . . 1          (3.25) eP∗ k´e a matriz k×kde autocorrela¸c˜oes em que a ´ultima coluna ´e substitu´ıda por [ρ1, ρ2, . . . , ρk]T. Sabe-se que φ11 =ρ1, φ22 =ρ2−ρ2 1 1−ρ2 1 eφ33 =ρ3(1 −ρ2 1) + ρ1(ρ2 1+ρ2 2−2ρ2) (1 −ρ2)(1 −ρ2−2ρ2 1).(3.26) Um processo de ru´ıdo branco ´e um processo estoc´astico caracterizado pela sucess˜ao 36 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais de vari´aveis aleat´orias independentes e identicamente distribu´ıdas com m´edia e variˆancia constantes, com as seguintes propriedades ∀t, E[t] = µt,(3.27) ∀t, V ar(t) = σ2 ,(3.28) ∀t, ∀k=±1,±2, . . . , Cov(t, t+k) = γk= 0.(3.29) Se, al´em disso, as vari´aveis aleat´orias seguirem uma distribui¸c˜ao Normal, designa-se o processo por ru´ıdo branco gaussiano. Um ru´ıdo branco ´e um processo estacion´ario cujas fun¸c˜oes de autocorrela¸c˜ao (FAC) e autocorrela¸c˜ao parcial (FACP) s˜ao nulas para todo o k6= 0. 3.2.1 Processo Estoc´astico N˜ao Estacion´ario No contexto ambiental, as s´eries geralmente s˜ao n˜ao estacion´arias. Um processo pode ser n˜ao estacion´ario, na medida em que a m´edia e/ou a variˆancia s˜ao fun¸c˜oes do tempo e n˜ao constantes. Numa primeira an´alise, a s´erie pode ser transformada de forma a obter-se uma s´erie estacion´aria (estabilizar a m´edia e/ou a variˆancia). No caso de uma s´erie n˜ao estacion´aria em m´edia e em variˆancia, deve-se estabilizar a variˆancia e s´o depois a m´edia (Murteira et al., 1993; Caiado, 2016). Mas tamb´em existem m´etodos que extraem a tendˆencia e a sazonalidade na s´erie temporal, fazendo com que haja estacionariedade, com base na decomposi¸c˜ao das suas componentes. Transforma¸c˜oes para a Estacionariedade Em muitos processos, quando se pretende estabilizar a m´edia, utiliza-se a diferencia¸c˜ao, atrav´es da aplica¸c˜ao do operador diferen¸ca ∆. Assim, a s´erie, Yt, n˜ao estacion´aria pode ser sujeita a uma diferencia¸c˜ao de primeira ordem, tal que ∆Yt=Yt−Yt−1, t = 2,3, . . . , n. (3.30) Se, ap´os a aplica¸c˜ao da diferencia¸c˜ao de primeira ordem, n˜ao se atingir a estacionariedade, aplica-se a diferencia¸c˜ao de segunda ordem, isto ´e, ∆2Yt= ∆(∆Yt) = ∆(Yt−Yt−1) = Yt−2Yt−1+Yt−2, t = 3,4, . . . , n. (3.31) 37 Cap´ıtulo 3. S´eries Temporais O operador de diferencia¸c˜ao de ordem ddefine-se como ∆dYt= ∆(∆d−1Yt), d ≥1 e t=d+ 1, . . . , n. (3.32) ´ E de evitar a sobrediferencia¸c˜ao na medida em que a variˆancia aumenta com a ordem de diferencia¸c˜ao. Uma boa pr´atica ´e a diferencia¸c˜ao at´e `a primeira ou `a segunda ordem para a obten¸c˜ao de uma s´erie estacion´aria. Para estabilizar a variˆancia de uma s´erie n˜ao estacion´aria pode recorrer-se a transforma¸c˜oes param´etricas, como ´e exemplo a transforma¸c˜ao de Box-Cox, dada pela seguinte express˜ao Zt=T(Yt) =    Yλ t λ, λ 6= 0, log(Yt), λ = 0 , λ ∈[−1,1].(3.33) Habitualmente, este tipo de transforma¸c˜oes est˜ao definidas para s´eries temporais de valores positivos. Para ultrapassar esta dificuldade, nas s´eries que apresentam valores negativos, utiliza-se a adi¸c˜ao de uma constante cque a torne positiva e s´o depois se recorre `as transforma¸c˜oes, como o logaritmo. De notar que, ap´os a transforma¸c˜ao dos dados, os valores ajustados pelo modelo estar˜ao nas unidades transformadas, isto significa que ´e necess´ario reverter as transforma¸c˜oes de modo a obter as previs˜oes nas unidades originais (Jebb et al., 2015). Passeio Aleat´orio Um passeio aleat´orio ´e caracterizado por movimentos de tendˆencia crescente ou decrescente, em per´ıodos longos, seguidos de mudan¸cas abruptas imprevis´ıveis e define-se Yt=Yt−1+t,(3.34) com t´e um ru´ıdo branco. No caso em que Y0´e conhecido, o modelo de passeio aleat´orio representa-se por Yt=Y0+ t X i=1 i.(3.35) Por´em, o modelo de passeio aleat´orio ´e caracterizado por ser um processo n˜ao estacion´ario, na medida em que a sua variˆancia depende do tempo t(Enders, 2015). 38 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Passeio Aleat´orio com Drift Um passeio aleat´orio com drift ´e uma extens˜ao do modelo de passeio aleat´orio com a adi¸c˜ao de um termo constante, a0, podendo ser formulado por Yt=a0+Yt−1+t,(3.36) em que t´e o ru´ıdo branco. A tendˆencia ´e descrita por parcelas determin´ısticas e parcelas estoc´asticas. Assim, o valor m´edio depende do tempo t, no caso em que Y0´e conhecido, tem-se o modelo de passeio aleat´orio com drift formulado por Yt=Y0+a0t+ t X i=1 i.(3.37) Ao calcular-se o valor esperado de Yt,Y0+a0t, ´e percet´ıvel que depende de te, assim, este processo ´e n˜ao estacion´ario. Estacionariedade A an´alise pr´evia da estacionariedade pode ser realizada pela representa¸c˜ao gr´afica da s´erie, ao longo do tempo. Posteriormente, ´e necess´ario utilizar testes estat´ısticos de forma a realizar um estudo formal. Existem v´arios testes que permitem a avalia¸c˜ao da estacionariedade da s´erie, nomeadamente, o teste de Dickey-Fuller, o teste de Dickey-Fuller Aumentado (Augmented Dickey Fuller), o teste Phillips-Perron e o teste de KwiatkowskiPhillips-Schmidt-Shin. Os trˆes primeiros testes tˆem como hip´otese testar a presen¸ca de uma raiz unit´aria (n˜ao estacionariedade) e na sua n˜ao rejei¸c˜ao os testes fornecem informa¸c˜ao sobre o n´umero de diferencia¸c˜oes necess´arias para atingir a estacionariedade. A hip´otese nula do ´ultimo teste considera que a s´erie temporal ´e estacion´aria. Mais detalhes sobre estes testes podem ser consultados em Dickey & Fuller (1979); Said & Dickey (1984); Phillips & Perron (1988); Kwiatkowski et al. (1992). 3.3 Metodologia Box-Jenkins Nesta Sec¸c˜ao apresenta-se uma breve s´umula sobre as metodologias a adotar no estudo de s´eries temporais, no ˆambito deste estudo. Dar-se-´a ˆenfase aos modelos AR, utilizados, posteriormente, no Cap´ıtulo 4. Em 1970, Box & Jenkins desenvolveram o seu trabalho sobre os modelos SARIMA (Seasonal Autoregressive Integrated Moving Average), com o objetivo da modela¸c˜ao e da previs˜ao de s´eries temporais estacion´arias e n˜ao estacion´arias. Estes modelos descrevem a s´erie Ytcomo uma fun¸c˜ao dos seus valores passados e como combina¸c˜ao linear de uma 39 Cap´ıtulo 3. S´eries Temporais sucess˜ao de choques aleat´orios. Dentro destes modelos, os mais simples s˜ao os Modelos Autorregressivos (AR), os Modelos de M´edias M´oveis (MA) e os Modelos Autorregressivos de M´edias M´oveis (ARMA). O primeiro descreve o comportamento da s´erie `a custa dos seus valores passados, o segundo atrav´es de uma sucess˜ao de choques aleat´orios, ao longo do tempo. O modelo ARMA ´e a combina¸c˜ao dos dois modelos anteriores. Estes modelos s˜ao ´uteis para s´eries estacion´arias. Todavia, os modelos mencionados, quando aplicados a processos n˜ao estacion´arios, n˜ao revelam um bom ajustamento, neste caso ´e aconselh´avel recorrer aos Modelos Autorregressivos Integrados de M´edias M´oveis (ARIMA). Este tipo de modelos s˜ao designados por modelos integrados, uma vez que o modelo estacion´ario, que ´e ajustado aos dados diferenciados, deve ser somado ou integrado para fornecer um modelo para os dados n˜ao estacion´arios. ` A semelhan¸ca dos modelos ARMA, estes modelos podem ser generalizados para incluir termos sazonais dando origem aos Modelos Autorregressivos Integrados de M´edias M´oveis Sazonais (SARIMA). Quando se utiliza este tipo de modelos recorre-se `a metodologia Box-Jenkins para a sua sele¸c˜ao. Esta metodologia implica um processo iterativo constitu´ıdo por trˆes fases: identifica¸c˜ao do modelo, estima¸c˜ao dos parˆametros e an´alise de diagn´ostico. A ideia base da identifica¸c˜ao do modelo ´e que se uma s´erie temporal ´e gerada a partir de um processo SARIMA, ent˜ao deve ter algumas propriedades te´oricas de autocorrela¸c˜ao. Box & Jenkins (1970) propuseram, ent˜ao, usar a fun¸c˜ao de autocorrela¸c˜ao (FAC) e a fun¸c˜ao de autocorrela¸c˜ao parcial (FACP) como ferramentas b´asicas para identificar as ordens do Modelo Autorregressivo Integrado de M´edias M´oveis Sazonal (SARIMA). A estima¸c˜ao permite a obten¸c˜ao dos parˆametros do modelo escolhido (pelo M´etodo de M´axima Verosimilhan¸ca e pelo M´etodo dos M´ınimos Quadrados) e ´e feita a sua avalia¸c˜ao na an´alise diagn´ostico. A fase de diagn´ostico engloba duas etapas: a avalia¸c˜ao da qualidade das estimativas obtidas e a avalia¸c˜ao da qualidade do ajustamento do modelo `as observa¸c˜oes da s´erie em estudo (deve-se proceder `a an´alise dos res´ıduos que devem ter um comportamento semelhante a um ru´ıdo branco). No contexto de dados ambientais existem v´arios estudos que mostraram que o processo para melhor descrever o comportamento dos res´ıduos ´e o processo AR. 3.3.1 Processo Autorregressivo de ordem p, AR(p) Considerando o processo autorregressivo de ordem p, AR(p), sabe-se que este tem como suporte o facto de que a observa¸c˜ao da vari´avel no instante testar relacionada linearmente com as observa¸c˜oes nos instantes anteriores. O processo Ytdiz-se um processo autorregressivo de ordem p, AR(p), quando satisfaz a equa¸c˜ao 40 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Yt=φ1Yt−1+φ2Yt−2+···+φpYt−p+t,(3.38) em que t´e um ru´ıdo branco com m´edia nula e independente de Yt−k,∀k≥1. A vari´avel Ytpode ser vista como uma vari´avel dependente que ´e explicada atrav´es de uma regress˜ao linear m´ultipla, em que as observa¸c˜oes em pinstantes anteriores funcionam como vari´aveis explicativas e φis˜ao os coeficientes de cada Yt−i. Existe outra forma de representa¸c˜ao deste processo atrav´es do operador atraso Φp(B)Yt=t,(3.39) em que Φp(B) = 1 −φ1B− ··· − φpBP´e o polin´omio autorregressivo de ordem p. A fatoriza¸c˜ao deste polin´omio ´e poss´ıvel atrav´es das pra´ızes, G−1 1, . . . , G−1 p, da caracter´ıstica, Φp(B) = 0, tornando-se poss´ıvel fatorizar o polin´omio autorregressivo do seguinte modo Φp(B)Yt= p Y i=1 (1 −GiB).(3.40) Se o m´odulo de cada ra´ız da equa¸c˜ao caracter´ıstica for inferior a um, ou de forma equivalente, |Gi|<1 com i= 1, . . . , p, ent˜ao ´e condi¸c˜ao necess´aria e suficiente para que o processo seja estacion´ario. Garantida esta condi¸c˜ao, sabe-se que o processo ´e invert´ıvel, o que significa que a dependˆencia do passado vai sendo menor `a medida que o passado se torna mais remoto. Graficamente, a FACP de um processo AR(p) revela uma queda brusca para zero a partir do lag p+ 1 e a FAC apresenta um decaimento exponencial ou sinusoidal amortecido para zero. 3.3.2 Sazonalidade Os dados ambientais tipicamente apresentam uma forte sazonalidade e esta ´e poss´ıvel ser integrada no processo de modela¸c˜ao. Indicadores Sazonais Nas s´eries temporais com sazonalidade, a sazonalidade pode ser modelada atrav´es da especifica¸c˜ao de um modelo de regress˜ao que inclua uma vari´avel indicatriz para representar cada um dos sper´ıodos sazonais, isto ´e, Yt=tβ +D1γ1+···+Dsγs+t, t = 1, . . . , n, (3.41) em que tβ representa a tendˆencia, γ1, . . . , γss˜ao os coeficientes que representam os s efeitos sazonais e Dks˜ao as vari´aveis indicatrizes, que representam os diferentes per´ıodos 41 Cap´ıtulo 3. S´eries Temporais sazonais: tomam o valor 1 quando o tempo tpertence ao per´ıodo ke 0 nos restantes casos. Por exemplo, para dados mensais, se D1corresponder `as ocorrˆencias no mˆes de janeiro (ou seja, 1, se tocorre em janeiro e 0, caso contr´ario), ent˜ao γ1s´o ´e tido em considera¸c˜ao para observa¸c˜oes registadas nesse mˆes. Sazonalidade Harm´onica Na integra¸c˜ao da sazonalidade no modelo, pode ser considerada uma vari´avel indicatriz por cada per´ıodo sazonal considerado. Na pr´atica, os efeitos sazonais s˜ao refletidos de forma cont´ınua e suave, o que leva a considerar outros tipos de representa¸c˜ao da sazonalidade. Cowpertwait & Melcalfe (2009) apresentam uma alternativa atrav´es de um modelo sazonal harm´onico, recorrendo a fun¸c˜oes trigonom´etricas (seno e cosseno), para incorporar as oscila¸c˜oes observadas. De forma simples, pode-se descrever uma onda sinusoidal por Asen(2πft +φ) = αccos(2πft) + αssen(2πft),(3.42) em que f´e a frequˆencia dos ciclos, Arepresenta a amplitude, φ´e a constante de fase, αs=Acos(φ) e αc=Asen(φ). O modelo sazonal harm´onico pode ser definido por Yt=tβ + s/2 X k=1 "α1kcos(2πkt s) + α2ksen(2πkt s)#+t,(3.43) em que tβ representa a tendˆencia, α1keα2ks˜ao os parˆametros desconhecidos de interesse, s´e o per´ıodo sazonal (s= 12 para dados mensais), k´e um ´ındice que varia entre 1 e s/2, et´e uma vari´avel codificada que representa o tempo (por exemplo, no caso em estudo, t= 1, . . . , 132 para 132 observa¸c˜oes igualmente espa¸cadas). 42 Cap´ıtulo 4 Modelos Muitos estudos estat´ısticos tˆem como objetivo principal o estudo da rela¸c˜ao entre vari´aveis ou, em particular, a an´alise da influˆencia que uma ou mais vari´aveis (explicativas), medidas em indiv´ıduos ou objetos, tˆem sobre uma vari´avel de interesse, que se denomina por vari´avel resposta (Turkman & Silva, 2000). O modo como o estat´ıstico aborda tal problema ´e atrav´es do modelo de regress˜ao. 4.1 Modelos Lineares Generalizados O Modelo Linear Normal (MLN) ´e o mais usado na modela¸c˜ao estat´ıstica. Este modelo tem v´arias limita¸c˜oes: a rela¸c˜ao ´e descrita atrav´es de uma fun¸c˜ao linear; exige a independˆencia das respostas e a vari´avel dependente condicionada aos valores das vari´aveis explicativas segue a Distribui¸c˜ao Normal, com variˆancia constante (condicionada aos valores das vari´aveis explicativas). O MLN pode ser expresso pela seguinte equa¸c˜ao Y=β0+β1X1+···+βpXp+, (4.1) em que Y´e a vari´avel resposta, X1, . . . , Xps˜ao as vari´aveis explicativas, p´e o n´umero de vari´aveis explicativas, β0, β1, . . . , βps˜ao os parˆametros do modelo e ´e o erro aleat´orio e n˜ao observ´avel e assume-se que E[] = 0 e V ar[] = σ2. Outra poss´ıvel representa¸c˜ao do MLN ´e Y|X∼N(µ, σ2), E[Y|X] = µ=β0+β1X1+···+βpXp. (4.2) Em dados reais, os pressupostos do MLN s˜ao dif´ıceis de verificar e muitas vezes recorrese a transforma¸c˜oes. Box & Cox (1964) prop˜oem uma transforma¸c˜ao que tem como objetivo verificar os pressupostos da normalidade, da variˆancia constante e da linearidade. 43 Cap´ıtulo 4. Modelos Esta transforma¸c˜ao faz com que a vari´avel resposta seja alterada, podendo at´e deixar de ser definida no espa¸co amostral original. Com a evolu¸c˜ao dos modelos, associada ao desenvolvimento computacional, foi poss´ıvel estabelecer uma extens˜ao do MLN a distribui¸c˜oes n˜ao normais, os Modelos Lineares Generalizados (MLG) apresentados por Nelder & Wedderburn (1972). Os MLG tiveram um impacto consider´avel na evolu¸c˜ao da estat´ıstica aplicada, mais concretamente, nas ´ultimas duas d´ecadas, existiu um grande avan¸co que permitiu a acessibilidade e o dinamismo destes modelos. Turkman & Silva (2000) descrevem esta importˆancia: “Do ponto de vista te´orico a sua importˆancia adv´em, essencialmente, do facto de a metodologia destes modelos constituir uma abordagem unificada de muitos procedimentos estat´ısticos correntemente usados nas aplica¸c˜oes e promover o papel central da verosimilhan¸ca na teoria da inferˆencia”. As vantagens principais dos MLG s˜ao a possibilidade de admitir v´arias distribui¸c˜oes para a vari´avel resposta, atrav´es da fam´ılia exponencial de distribui¸c˜oes, e a flexibilidade para a rela¸c˜ao funcional entre o valor esperado da vari´avel resposta (µ) e o preditor linear (η=β0+β1X1+···+βpXp). Algumas das limita¸c˜oes dos MLG s˜ao o facto de manterem uma estrutura de linearidade, das distribui¸c˜oes da vari´avel resposta se restringirem `a fam´ılia exponencial e de exigirem a independˆencia das respostas. 4.1.1 Nota¸c˜ao e Terminologia Os indiv´ıduos s˜ao considerados as unidades do estudo. O valor da vari´avel resposta, do indiv´ıduo i, denomina-se por yie ´e uma realiza¸c˜ao da vari´avel resposta (ou vari´avel dependente) Yi, em que ivaria de 1 a n, em que n´e o n´umero total dos indiv´ıduos em estudo. A vari´avel resposta, para um indiv´ıduo i, representa-se pelo vetor das vari´aveis resposta, resultante das medi¸c˜oes, isto ´e, Y=      Y1 Y2 . . . Yn       ,(4.3) ou, em alternativa, YT= (Y1, Y2, . . . , Yn). Cada vari´avel aleat´oria Yiassociada a cada iest´a o vetor das pcovari´aveis (ou vari´aveis explicativas ou vari´aveis independentes) com dimens˜ao p×1, isto ´e, XT= (X1, ..., Xp). Assim, para Y, a matriz da vari´avel resposta, de ordem n×p´e dada por 44 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais X=      x11 . . . x1p x21 . . . x2p . . ..... . . xn1. . . xnp       ,(4.4) associado ao vetor de parˆametros desconhecidos β= (β1, β2, . . . , βp)T. 4.1.2 Fam´ılia Exponencial Nos Modelos Lineares Generalizados, a vari´avel resposta segue uma distribui¸c˜ao que pertence `a fam´ılia exponencial. Uma vari´avel aleat´oria Ytem distribui¸c˜ao pertencente `a fam´ılia exponencial de dispers˜ao (ou simplesmente fam´ılia exponencial) se a sua fun¸c˜ao densidade de probabilidade (f.d.p.) ou a sua fun¸c˜ao massa de probabilidade (f.m.p.) se puder escrever da seguinte forma f(y|θ, φ) = exp(yθ −b(θ) a(φ)+c(y, φ)),(4.5) em que θ´e o parˆametro de localiza¸c˜ao e φ´e o parˆametro de dispers˜ao. a(.), b(.), c(.) s˜ao fun¸c˜oes reais espec´ıficas para cada distribui¸c˜ao. Uma descri¸c˜ao mais pormenorizada desta fam´ılia pode ser consultada em Cox & Hinkley (1974). Turkman & Silva (2000) apontam que quando φfor conhecido, tem-se uma distribui¸c˜ao da fam´ılia exponencial com parˆametro can´onico θ. No caso de φser desconhecido, a distribui¸c˜ao pode ou n˜ao fazer parte da fam´ılia exponencial, a fun¸c˜ao b(.) ´e diferenci´avel e que o suporte da distribui¸c˜ao n˜ao depende dos parˆametros. Habitualmente, tem-se que a(φ) = φ w, em que w´e uma constante conhecida e revela o peso da observa¸c˜ao. Considerando `(θ, φ, y) = log(f(y|θ, φ)), em que `representa o logaritmo da fun¸c˜ao de verosimilhan¸ca, define-se por fun¸c˜ao Score S(θ) = ∂`(θ, φ, y) ∂θ .(4.6) Sabe-se que no caso das fam´ılias regulares, tem-se que E[S(θ)] = 0,(4.7) E[S2(θ)] = E"∂`(θ, φ, y) ∂θ 2#,(4.8) 45 Cap´ıtulo 4. Modelos 1. Dado ˆ β(k)(com k= 0) determina-se u(k) i, que ´e um vetor com elemento gen´erico u(k) i=η(k) i+ (yi−µ(k) i)∂ηi(k) ∂µ(k) i (4.37) e calcula-se W(k)utilizando a equa¸c˜ao (4.36); 2. A nova itera¸c˜ao β(k+1) ´e calculada usando ˆ β(k+1) = (ZTW(k)Z)−1ZTW(k)u(k).(4.38) O processo iterativo termina quando ´e atingido o seguinte crit´erio ||ˆ β(k+1) −ˆ β(k)|| ||ˆ β(k)|| ≤, (4.39) para algum valor  > 0 previamente definido. Salienta-se que a equa¸c˜ao (4.37) ´e idˆentica `a que se obteria para os estimadores de M´etodo dos M´ınimos Quadrados Ponderados se, em cada itera¸c˜ao, se calculasse a regress˜ao linear de u(k)em vez de Z, em que W(k)´e uma matriz de pesos. Na express˜ao (4.38), a matriz Wcont´em o parˆametro de dispers˜ao φ, mas este ´ultimo n˜ao faz parte dos c´alculos de β(k+1) e torna-se irrelevante para a determina¸c˜ao deste parˆametro. Por este motivo, assumindo φ= 1 n˜ao h´a perda de generalidade, na determina¸c˜ao das estimativas para β. Parˆametro de Dispers˜ao Para a estima¸c˜ao do parˆametro de dispers˜ao ´e poss´ıvel recorrer ao M´etodo de M´axima Verosimilhan¸ca. Mas existem m´etodos mais simples, que tamb´em apresentam bons resultados, como o m´etodo que assenta na distribui¸c˜ao de amostragem para grandes valores de n, da Estat´ıstica de Pearson Generalizada, ˆ φ=1 n−p n X i=1 wi(yi−ˆµi)2 V(µi),(4.40) em que ˆ φ´e um estimador consistente de φ. Quando a amostra ´e muito grande (elevado valor n), esta estat´ıstica pode ser aproximada `a distribui¸c˜ao Qui-quadrado, com n−p graus de liberdade. 52 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Propriedades dos Estimadores de M´axima Verosimilhan¸ca O estimador m´axima verosimilhan¸ca, o ˆ β, ´e um estimador assintoticamente centrado, ou seja, E[ˆ β]≃β.(4.41) A matriz de covariˆancia de ˆ β´e aproximadamente igual ao inverso da matriz de informa¸c˜ao de Fisher, isto ´e, cov(ˆ β)≃E[(ˆ β−β)(ˆ β−β)T] = I−1(β).(4.42) Al´em disso, a distribui¸c˜ao assint´otica de ˆ β´e normal (p+1)-variada, com valor esperado βe matriz de covariˆancia I−1(β), ˆ β∼Np(β, I−1(β)).(4.43) A distribui¸c˜ao assint´otica da estat´ıtica de Wald (ˆ β−β)TI(β)(ˆ β−β) tem uma distribui¸c˜ao assint´otica de um χ2, com p+ 1 graus de liberdade. Estes resultados s˜ao importantes para fazer inferˆencia sobre β, nomeadamente, nos testes de hip´oteses e intervalos de confian¸ca. 4.1.5 Inferˆencia A escolha das vari´aveis apropriadas ´e uma etapa extremamente importante na modela¸c˜ao estat´ıstica. A an´alise da rela¸c˜ao entre a vari´avel resposta e as covari´aveis permite ter uma ideia das vari´aveis que ser˜ao importantes para um modelo final. Testes de Hip´oteses Os testes de hip´oteses s˜ao procedimentos de valida¸c˜ao estat´ıstica que possibilitam avaliar hip´oteses sobre determinadas caracter´ısticas da popula¸c˜ao, sujeitos a um determinado n´ıvel de risco de falharem associado (n´ıvel de confian¸ca). Os testes s˜ao divididos em param´etricos ou n˜ao param´etricos. Os testes param´etricos baseiam-se em parˆametros ou caracter´ısticas quantitativas da vari´avel dependente, sob o pressuposto de que as vari´aveis seguem uma Distribui¸c˜ao Normal. Os testes n˜ao param´etricos s˜ao utilizados quando as vari´aveis s˜ao ordinais ou categ´oricas e os pressupostos param´etricos n˜ao se verifiquem. Estes testes n˜ao s˜ao t˜ao potentes comparativamente com os testes param´etricos. 53 Cap´ıtulo 4. Modelos Teste de Wald O Teste de Wald ´e apropriado para testar uma componente do vetor parˆametro, com as hip´oteses de teste H0:Cβ=ξ H1:Cβ6=ξ, (4.44) em que C´e uma matriz q×pde caracter´ıstica completa q, o vetor ξtem dimens˜ao qe ´e especificado pelo investigador. A estat´ıstica de teste (ET), denominada por Estat´ıstica de Wald, ´e baseada na normalidade assint´otica do estimador de MV, ´e dada por W= (Cˆ β−ξ)T[CI−1(ˆ β)CT]−1(Cˆ β−ξ) (4.45) e tem distribui¸c˜ao assint´otica χ2, com qgraus de liberdade. No caso de ET > χ2 1−α,q, em que α´e o n´ıvel de significˆancia, ent˜ao a hip´otese nula ´e rejeitada. Teste da Raz˜ao de Verosimilhan¸ca O Teste da Raz˜ao de Verosimilhan¸ca ´e apropriado para testar H0:Cβ=ξ H1:Cβ6=ξ, (4.46) em que C´e uma matriz q×pde caracter´ıstica completa q, o vetor ξtem dimens˜ao qe ´e especificado pelo investigador. A estat´ıstica de teste (ET), denominada por estat´ıstica de Wilks ou estat´ıstica da Raz˜ao de Verosimilhan¸cas, ´e baseada na distribui¸c˜ao assint´otica da raz˜ao do m´aximo das verosimilhan¸cas sob as hip´oteses H0eH0∪H1, ´e dada por ET =−2`(˜ β)−`(ˆ β)∼χ2 q.(4.47) No caso de ET > χ2 1−α,q, em que α´e o n´ıvel de significˆancia, ent˜ao a hip´otese nula ´e rejeitada. De notar que tamb´em existem outras estat´ısticas que n˜ao s˜ao mencionadas neste trabalho, como a estat´ıstica de Rao (ou estat´ıstica Score), baseada nas propriedades assint´oticas da fun¸c˜ao Score. Sele¸c˜ao de Vari´aveis Existem trˆes procedimentos habitualmente utilizados na sele¸c˜ao de vari´aveis: backward elimination,forward selection estepwise selection. Estes tˆem como base um algoritmo que 54 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais permite perceber se a vari´avel ´e importante ou n˜ao para o modelo em estudo, recorrendo `a significˆancia estat´ıstica do coeficiente da vari´avel em an´alise. No procedimento backward elimination, ou processo de sele¸c˜ao regressiva, come¸ca-se com um modelo que inclui todas as covari´aveis, modelo completo. Com base no n´ıvel de confian¸ca do teste estat´ıstico, elimina-se a vari´avel menos significativa. Ajusta-se novamente o modelo, excluindo a vari´avel definida, e repete-se o procedimento at´e existirem no modelo covari´aveis todas significativas. O procedimento forward selection, ou processo de sele¸c˜ao progressiva, ´e o inverso do processo anterior. Inicia-se com o modelo sem covari´aveis e adiciona-se a vari´avel com o menor valor de prova. Repete-se o processo at´e n˜ao existirem mais covari´aveis significativas. O procedimento stepwise selection ´e a combina¸c˜ao das duas metodologias anteriores. Uma determinada covari´avel pode ser adicionada ou removida, a cada passo testa-se se as restantes s˜ao significativas para verificar se deve ser removida. 4.1.6 Qualidade de Ajustamento Ap´os a elabora¸c˜ao de um modelo, ´e necess´ario verificar a qualidade do seu ajustamento. Para avaliar essa qualidade, pode-se recorrer `a fun¸c˜ao desvio e `a Estat´ıstica de Pearson Generalizada. Fun¸c˜ao Desvio Como j´a foi mencionado, o logaritmo da fun¸c˜ao de verosimilhan¸ca de um MLG ´e dado por `(β) = logL(β) = n X i=1 wiyiθi−b(θi) φ+c(yi, φ, wi).(4.48) Para se comparar o modelo em estudo com o modelo completo (ou saturado), recorre-se `a estat´ıstica de raz˜ao de verosimilhan¸cas, obtendo-se D∗(y; ˆµ) = −2`M(ˆ βM)−`S(ˆ βS)=D(y; ˆµ) φ,(4.49) em que D∗(y; ˆµ) designa-se por desvio reduzido, ˆ βM´e o vetor de parˆametros, para o modelo em investiga¸c˜ao, ˆ βS´e o vetor de parˆametros, para o modelo saturado (modelo com todas as covari´aveis em estudo), e D(y; ˆµ) ´e o desvio para o modelo em estudo e ´e dado por 55 Cap´ıtulo 4. Modelos D(y; ˆµ) = n X i=1 2winyiq(yi)−q(ˆµi)−q(yi)+bq(ˆµi)o= n X i=1 di,(4.50) em que di´e a diferen¸ca entre os logaritmos das fun¸c˜oes de verosimilhan¸ca da observada e ajustada em cada observa¸c˜ao. A estat´ıstica da raz˜ao de verosimilhan¸ca para comparar dois modelos, M1 e M2, ´e dada por D(y; ˆµ1)−D(y; ˆµ2) φ∼χ2 p1−p2,(4.51) em que p1ep2´e a dimens˜ao do vetor dos parˆametros β, para os Modelos 1 e 2, respetivamente. A compara¸c˜ao de modelos encaixados (modelos em que um ´e submodelo do outro) pode ser realizada atrav´es da diferen¸ca dos desvios. Estat´ıstica de Pearson Generalizada Como anteriormente visto, a Estat´ıstica de Pearson Generalizada ´e dada por ˆ φ=1 n−p n X i=1 wi(yi−ˆµi)2 V(µi),(4.52) em que ˆ φ´e um estimador consistente de φ. Quando a amostra ´e bastante grande (elevado valor n), esta estat´ıstica pode ser aproximada `a distribui¸c˜ao Qui-quadrado, com n−p graus de liberdade. Sele¸c˜ao e Compara¸c˜ao de Modelos A sele¸c˜ao e a compara¸c˜ao de modelos permite que se encontre o modelo mais parcimonioso, ou seja, o modelo que explique o comportamento da vari´avel resposta, com o m´ınimo de parˆametros poss´ıveis. A compara¸c˜ao da qualidade do ajustamento de dois modelos aninhados ´e geralmente realizada a partir de testes de hip´oteses, como o Teste da Raz˜ao de Verosimilhan¸ca. O Teste da Raz˜ao de Verosimilhan¸ca ´e apropriado para testar dois modelos, MpeMq, com p eqn´umero de vari´aveis respetivamente, desde que os mesmos sejam modelos aninhados, com as hip´oteses de teste H0: As q−pvari´aveis no modelo n˜ao s˜ao significativas H1: As q−pvari´aveis no modelo n˜ao s˜ao significativas,(4.53) sob a hip´otese nula, a estat´ıstica de teste e distribui¸c˜ao respetiva ´e 56 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais −2loghLMp(β) LMq(β)i∼χ2 q−p,(4.54) em que LMp(β)´e a fun¸c˜ao verosimilhan¸ca do modelo MpeLMq(β)´e a fun¸c˜ao verosimilhan¸ca do modelo Mq. Na maioria dos casos pr´aticos, existem v´arias vari´aveis que s˜ao potenciais candidatas para explicar a variabilidade da vari´avel resposta. Tal tem como consequˆencia a possibilidade de diversos modelos com combina¸c˜oes diferentes das covari´aveis para explicar o fen´omeno em estudo, o que leva a que o processo da sele¸c˜ao seja dif´ıcil e moroso. Usualmente recorre-se ao m´etodo de sele¸c˜ao stepwise porque facilita a sele¸c˜ao dos modelos. Quando se pretende comparar modelos n˜ao encaixados (isto ´e, no caso em que o Teste da Raz˜ao de Verosimilhan¸ca n˜ao ´e aplic´avel) usam-se os crit´erios de informa¸c˜ao. Os crit´erios mais utilizados para a sele¸c˜ao de modelos s˜ao o Crit´erio de Informa¸c˜ao de Akaike (AIC) e o Crit´erio Bayesiano de Schwarz (BIC). Nestes crit´erios, a compara¸c˜ao dos modelos ´e feita a partir da maximiza¸c˜ao do logaritmo da verosimilhan¸ca, acrescentando uma penalidade ao n´umero de parˆametros de modelo. Neste sentido, o modelo selecionado ´e o modelo de menor valor do crit´erio de informa¸c˜ao. Akaike (1974) prop˜oe o Crit´erio de Informa¸c˜ao de Akaike que tem a seguinte f´ormula AIC =−2`(ˆ β,ˆα)+2npar,(4.55) em que `´e a fun¸c˜ao logaritmo de verosimilhan¸ca, p´e o n´umero de parˆametros, npar ´e o n´umero de parˆametros do modelo em estudo. Schwarz (1978) prop˜oe o Crit´erio de Informa¸c˜ao Bayesiano, estabelecendo que BIC =−2`(ˆ β,ˆα)+2nparln(N),(4.56) em que nnpar ´e o n´umero de parˆametros do modelo em estudo e N´e o n´umero total de observa¸c˜oes. 4.1.7 An´alise Diagn´ostico A an´alise de diagn´ostico tem como objetivo averiguar a existˆencia de desvios isolados do modelo, ou seja, a existˆencia de uma ou mais observa¸c˜oes mal ajustadas, n˜ao tendo o mesmo comportamento que as restantes observa¸c˜oes. Os desvios sistem´aticos podem ser provocados pela sele¸c˜ao inadequada da fun¸c˜ao de variˆancia, da fun¸c˜ao de liga¸c˜ao e da matriz do modelo, ou pela defini¸c˜ao da escala da vari´avel resposta ou das covari´aveis. As discrepˆancias isoladas podem ser resultantes dos pontos estarem nos extremos da amplitude de validade da covari´avel ou dos pontos estarem errados, devido a uma leitura 57 Cap´ıtulo 4. Modelos errada ou a uma transcri¸c˜ao mal realizada ou fatores n˜ao controlados. As metodologias de an´alise de res´ıduos do MLG s˜ao semelhantes `as dos Modelos Lineares, sofrendo algumas modifica¸c˜oes. A variˆancia residual ´e alterada por uma estimativa consistente do parˆametro φ. A an´alise da variˆancia requer um cuidado especial, pois deve ser adequada `a distribui¸c˜ao em estudo. O seu comportamento deve ser constante e unit´ario. O conceito de res´ıduos, ri=Yi−ˆ Yi, do Modelo Linear tem que ser adaptado `a vari´avel dependente ajustada e ao preditor linear. Um res´ıduo tem como objetivo descrever a discrepˆancia entre o valor observado e o valor ajustado pelo modelo e pode ser calculado de v´arias formas. Os res´ıduos ordin´arios s˜ao definidos por ri=yi−ˆµi.(4.57) Os res´ıduos de Pearson determinam-se por rpi =(yi−ˆµi)√wi pV(ˆµi).(4.58) Os res´ıduos de Pearson padronizados calculam-se por rP pi =rpi qˆ φ(1 −hii) ,(4.59) em que hii ´e o i-´esimo elemento da diagonal principal da matriz H, dada por H=W1/2Z(ZTWZ)−1ZTW1/2.(4.60) A matriz Hdepende das vari´aveis explicativas, da fun¸c˜ao de liga¸c˜ao e da fun¸c˜ao de variˆancia, tornando mais dif´ıcil a interpreta¸c˜ao da medida de alavanca. Os res´ıduos deviance (desvio residual) s˜ao determinados por rDi =sinal(yi−ˆµi)pdi,(4.61) em que d(y, ˆµ) = Pn i=1 di. Os res´ıduos deviance padronizados s˜ao dados por rP Di =rDi qˆ φ(1 −hii) .(4.62) Se os res´ıduos de Pearson e deviance revelam a variˆancia aproximadamente constante, ent˜ao ´e ind´ıcio que o modelo ´e adequado. 58 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Turkman & Silva (2000) sugerem a representa¸c˜ao gr´afica dos res´ıduos padronizados vs preditores lineares ou alguma transforma¸c˜ao adequada dos preditores lineares. Algumas transforma¸c˜oes s˜ao apresentadas por McCullagh & Nelder (1989). Se n˜ao existirem anomalias, os res´ıduos devem predispor-se em torno de zero, sem ordem, para os diferentes valores de ˆµ. Se existirem anomalias, estas podem dever-se `a escolha errada da fun¸c˜ao de liga¸c˜ao, `a escolha errada da escala de uma ou mais covari´aveis ou `a omiss˜ao de um termo quadr´atico. Os gr´aficos adequados s˜ao as representa¸c˜oes dos res´ıduos padronizados vs preditores lineares ou fun¸c˜ao dos valores ajustados ˆµou ´ındices. Nesta avalia¸c˜ao deve-se encontrar um padr˜ao nulo, os res´ıduos devem estar dispostos em torno de zero com uma amplitude constante para diferentes valores de ˆµ. Para avaliar a presen¸ca de correla¸c˜ao entre observa¸c˜oes deve-se recorrer a gr´aficos dos res´ıduos do modelo vs a ordem da observa¸c˜ao. 4.2 Modelo de Efeitos Mistos A evolu¸c˜ao da ciˆencia dos dados permitiu o aumento da complexidade na an´alise estat´ısticas e das metodologias estat´ıstica (Kass et al., 2016). Quando cada indiv´ıduo apresenta medidas repetidas ao longo do tempo duma determinada caracter´ıstica, sendo o pr´oprio tempo um fator de interesse, diz-se que se tem dados longitudinais (Molenberghs & Verbeke, 2005). A sua obten¸c˜ao poder´a ser prospetiva (os registos s˜ao obtidos atrav´es dos ind´ıviduos, em estudo, ao longo do tempo) ou retrospetiva, atrav´es de um hist´orico de onde s˜ao extra´ıdas as medi¸c˜oes dos indiv´ıduos em estudo (Diggle et al., 2002). Os estudos longitudinais analisam medi¸c˜oes repetidas ao longo do tempo sobre o mesmo indiv´ıduo, que permite separar o que no contexto de um estudo populacional se chama o efeito de coorte do efeito da idade (Diggle et al., 2002). Outra particularidade dos dados longitudinais ´e o facto de serem agrupados, ou seja, cada grupo ´e o resultado das medi¸c˜oes repetidas de um determinado indiv´ıduo, ao longo do tempo. Estas medi¸c˜oes tˆem uma ordem cronol´ogica e, consequentemente, existe uma correla¸c˜ao entre as observa¸c˜oes de um determinado indiv´ıduo. Cada indiv´ıduo tem um vetor resposta com todas as observa¸c˜oes ao longo do tempo que usualmente est˜ao correlacionadas, o que leva a que a estrutura de autocorrela¸c˜ao tenha um papel fulcral na modela¸c˜ao e estima¸c˜ao de cada parˆametro do modelo (Diggle et al., 2002). Um estudo longitudinal tem como objetivo caracterizar a altera¸c˜ao da vari´avel resposta, ao longo do tempo, assim como a rela¸c˜ao entre as covari´aveis e a vari´avel resposta. Consideram-se alguns fatores que fazem com que esta an´alise estat´ıstica seja complexa, nomeadamente, a estrutura de autocorrela¸c˜ao revela uma especial importˆancia no ajusta59 Cap´ıtulo 4. Modelos mento, a existˆencia de variabilidade entre indiv´ıduos distintos, o n´umero de observa¸c˜oes de cada indiv´ıduo pode ser diferente e as covari´aveis tamb´em se podem modificar ao longo do tempo. A principal vantagem ´e a utiliza¸c˜ao da totalidade dos dados que leva a evidenciar altera¸c˜oes dentro do mesmo indiv´ıduo, o aumento da potˆencia estat´ıstica (´e poss´ıvel distinguir os erros de medi¸c˜ao dos erros aleat´orios) e a redu¸c˜ao do enviesamento. Existem v´arios tipos de modelos para a an´alise de dados longitudinais, nomeadamente o Modelo Marginal e o Modelo de Efeitos Aleat´orios. O Modelo Marginal (Populationaverage) tem como objetivo inferir sobre o valor m´edio populacional. Num Modelo Marginal, o valor esperado marginal ´e modelado como fun¸c˜ao das covari´aveis. Como as medi¸c˜oes s˜ao repetidas, em cada indiv´ıduo, n˜ao tˆem tendˆencia a ser independentes, a an´alise marginal tem que incluir pressupostos em rela¸c˜ao `a correla¸c˜ao. O Modelo Marginal tem a vantagem do valor esperado da vari´avel resposta e a covariˆancia serem modelados separadamente (Diggle et al., 2002) Os Modelos de Efeitos Aleat´orios (Subject-specific) permitem descrever as altera¸c˜oes da resposta m´edia da vari´avel resposta de cada indiv´ıduo e a rela¸c˜ao destas com as covari´aveis possibilitam realizar inferˆencias sobre o indiv´ıduo, a modela¸c˜ao da sobredispers˜ao e da correla¸c˜ao intr´ınseca a cada indiv´ıduo, atrav´es da incorpora¸c˜ao dos efeitos aleat´orios. Estes modelos podem ser aplicados a vari´aveis resposta Distribui¸c˜ao Normal ou n˜ao. Os efeitos aleat´orios incorporam a heterogeneidade entre indiv´ıduos e s˜ao representados por vari´aveis aleat´orias que usualmente seguem a Distribui¸c˜ao Normal (Diggle et al., 1994). 4.2.1 Terminologia Consideram-se algumas nota¸c˜oes e conceitos b´asicos utilizados no presente Cap´ıtulo. Os indiv´ıduos s˜ao considerados unidades do estudo. As observa¸c˜oes no mesmo indiv´ıduo is˜ao medi¸c˜oes, ao longo do tempo, em que ivaria de 1 a n. O valor da vari´avel resposta, do indiv´ıduo ino tempo t, denomina-se por yit e ´e uma realiza¸c˜ao da vari´avel resposta Yit, em que tvaria de 1 a mi. A vari´avel resposta de um indiv´ıduo irepresenta-se pelo vetor das vari´aveis resposta, resultante das medi¸c˜oes, isto ´e, Yi=      Yi1 Yi2 . . . Yimi       (4.63) ou, em alternativa, YT i= [Yi1, Yi2, . . . , Yimi], sendo as medi¸c˜oes obtidas no indiv´ıduo i designadas por yT i= (yi1, yi2, . . . , yimi). Cada vari´avel aleat´oria Yit tem associado o vetor das pcovari´aveis com dimens˜ao p×1, isto ´e, xT it = (x1 it, ..., xp imi). Assim, para cada Yi, 60 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais est´a associada a matriz de ordem mi×pe ´e dada por Xi=      x1 i1. . . xp i1 x1 i2. . . xp i2 . . .. . .. . . x1 imi. . . xp imi       .(4.64) O valor esperado e a variˆancia de Yis˜ao designados por E(Yi) = µieV ar(Yi) = vie o vetor p×1 dos parˆametros desconhecidos, β= (β1, . . . , βp). A vari´avel resposta ´e dada pelo vetor de dimens˜ao M×1, Y=      YT 1 YT 2 . . . YT mi       =                       Y11 . . . Y1m1 Y21 . . . Y2m2 . . . Yn1 . . . Ynmn                       ,(4.65) em que M=Pn i=1 mi. A respetiva matriz de desenho ´e dada por X=                               x1 11 . . . xp 11 x1 12 . . . xp 12 . . .. . .. . . x1 1m1. . . xp 1m1 x1 21 . . . xp 21 x1 22 . . . xp 22 . . .. . .. . . x1 2m2. . . xp 2m2 . . .. . .. . . x1 n1. . . xp n1 x1 n2. . . xp n2 . . .. . .. . . x1 nmn. . . xp nmn                               ,(4.66) associado ao vetor de parˆametros desconhecidos β= (β1, β2, . . . , βp)T. 61 Cap´ıtulo 4. Modelos seguinte equa¸c˜ao ˆ βREML =n X i=1 XT iV−1 i(αREML)Xi−1n X i=1 XT iV−1 i(αREML)Yi =XTV−1 REMLXi−1XTV−1 REMLY. (4.93) Relativamente aos m´etodos apresentados, verifica-se que os estimadores resultantes n˜ao dependem de BeβREML ´e diferente de βML. A partir de W, este m´etodo permite que n˜ao haja perda de informa¸c˜ao sobre α(Patterson & Thompson, 1971). O m´etodo REML depende do estimador βobtido atrav´es do m´etodo ML. O termo −1 2logPn i=1 XT iV−1 i(α)Xrevela que quaisquer modifica¸c˜oes na matriz de desenho originam reparametriza¸c˜oes nos efeitos fixos e, como consequˆencia, n˜ao ´e poss´ıvel comparar diferentes Modelos de Efeitos Mistos com diferentes componentes fixas (Pinheiro & Bates, 2000). Os m´etodos conduzem a resultados semelhantes (Zuur et al., 2009), contudo as estimativas obtidas por ML das componentes da variˆancia s˜ao menores do que as obtidas por REML. 4.2.4 Predi¸c˜ao dos Efeitos Aleat´orios Os efeitos aleat´orios uirepresentam a variabilidade dentro do mesmo indiv´ıduo (evolu¸c˜ao do indiv´ıduo), relativamente `a m´edia da popula¸c˜ao Xiβ. Nos Modelos de Efeitos Mistos ´e importante perceber quais os efeitos a incluir. A distin¸c˜ao dos efeitos aleat´orios dos fixos, a sele¸c˜ao e a inferˆencia n˜ao ´e f´acil. Um cuidado adicional ´e a especifica¸c˜ao dos modelos aninhados, integrando os efeitos aleat´orios. Na pr´atica, os efeitos aleat´orios ser˜ao determinados pelo estudo definido ou pela amostra recolhida (Schielzeth & Nakagawa, 2013). Bates et al. (2015) fornecem um bom guia para determinar, de forma iterativa, a complexidade ideal da estrutura de efeitos aleat´orios. Considerando que uis˜ao vari´aveis aleat´orias, recorre-se `a abordagem Bayesiana para obter a melhor predi¸c˜ao. Anteriormente, viu-se que Yi|ui∼N(Xiβ+Ziui,Σi) (4.94) e que a distribui¸c˜ao a priori ´e dada por ui∼N(0,D).(4.95) A sua distribui¸c˜ao a posteriori pode ser determinada a partir da distribui¸c˜ao de ui, condicionada a Yi=yi. Assim, a distribui¸c˜ao a posteriori ´e dada por 68 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais f(ui|yi) = f(ui|yi) = f(yi|ui)f(ui) Rf(yi|ui)f(ui)dui .(4.96) Neste sentido, a distribui¸c˜ao conjunta ´e Normal Multivariada (Azzalini, 1996), ui Yi!∼N 0 Xiβ!,"D DZT i ZiD ZiDZT i+Σi#! (4.97) e, consequentemente, ui|yi∼NDZT i(ZiDZT i+Σi)−1(yi−Xiβ),V,(4.98) em que V=D−DZT i(ZiDZT i+Σi)−1ZiD. O valor esperado para a distribui¸c˜ao a posteriori pode ser escrito da seguinte forma E(ui|Yi=yi) = Zuif(ui|yi)dui =DZT i(ZiDZT i+Σi)−1(yi−Xiβ) =DZT iV−1 i(yi−Xiβ). (4.99) Considerando que α´e conhecido, o melhor preditor linear centrado para a vari´avel que representa os efeitos aleat´orios ´e determinado com base na equa¸c˜ao (4.99), em que β ´e aproximado por ˆ β(α)=(XV−1X)−1XTV−1y(McCulloch & Searle, 2001), em que se tem uBLUP,i(α) = DZT iV−1 i(yi−Xiˆ β(α)),(4.100) ou, alternativamente, uBLUP (α) = DZTV−1(y−Xˆ β(α)).(4.101) O preditor da combina¸c˜ao linear de v=vT Bβ+vT uui, condicionado a α´e dado por vBLUP (α) = vT βˆ β(α) + vT uˆuBLUP (α),(4.102) em que vβrepresenta o vetor p×1 dos efeitos fixos e vurepresenta o vetor p×1 dos efeitos aleat´orios. Harville (1976) e Searle et al. (1992) provam que vBLUP (α) ´e o melhor preditor linear centrado para v. Equa¸c˜oes A solu¸c˜ao do sistema de equa¸c˜oes lineares, `as quais se d´a o nome de equa¸c˜oes do Modelo de Efeitos Mistos, resultam em estimadores dos efeitos fixos βestimados e preditores dos efeitos aleat´orios u. Sabe-se que Y|u∼N(Xβ+Zu,Σ) e que u∼N(0,F), assim a 69 Cap´ıtulo 4. Modelos densidade conjunta de Yeu´e dada por f(y,u) = f(y|b)f(b) = (2π)−N 2|Σ−1 2| ×exp(−1 2(y−Xβ−Zu)TΣ−1(y−Xβ−Zu)) ×(2π)−qn 2|F−1 2|expn−1 2uTF−1uo = (2π)−N+qn 2|Σ−1 2||F−1 2| ×exp(−1 2h(y−Xβ−Zu)TΣ−1(y−Xβ−Zu) + uTF−1ui). (4.103) Aplicando o logaritmo na equa¸c˜ao anterior tem-se que logf(y,u) = −N+qn 2log(2π)−1 2log|Σ|− 1 2log|F| −1 2h(y−Xβ−Zu)TΣ−1(y−Xβ−Zu) + uTF−1ui. (4.104) Derivando em rela¸c˜ao aos parˆametros de interesse e igualando a zero, tem-se o seguinte sistema de equa¸c˜oes      ∂logf(y,b) ∂β= 0 ∂logf(y,b) ∂u= 0 ⇔         XTΣ−1(y−Xβ−Zu) = 0 ZTΣ−1(y−Xβ−Zu)−Fu = 0 ,(4.105) alternativamente, "XTΣ−1X XTΣ−1Z ZTΣ−1X ZTΣ−1Z+F#"β u#="XTΣ−1y ZTΣ−1y#.(4.106) Considerando FeΣconhecidos, tem-se que ˆ β(α) = (XTV−1X)−1XTV−1y =n X i=1 XT iV−1 iXi−1n X i=1 XT iV−1 iyi (4.107) e uBLUP =FZTV−1y−Xβ(α),(4.108) em alternativa, uBLUP,i =FZT iV−1 iyi−Xiβ(α).(4.109) 70 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais O estimador β, resultante das equa¸c˜oes do Modelo de Efeitos Mistos, ´e igual ao obtido pelo m´etodo ML. M´etodos Num´ericos As equa¸c˜oes de maximiza¸c˜ao do logaritmo da fun¸c˜ao verosimilhan¸ca exigem M´etodos Num´ericos para a sua otimiza¸c˜ao (Pinheiro & Bates, 1995). Existem v´arios M´etodos Num´ericos de otimiza¸c˜ao, como s˜ao exemplo o Algoritmo Expected-Maximization e o M´etodo de Newton-Raphson. Estes dois est˜ao implementados no R. O Algoritmo Expected-Maximization ´e um m´etodo iterativo e cada itera¸c˜ao consiste em dois passos. O primeiro determina o valor esperado do logaritmo da fun¸c˜ao verosimilhan¸ca e o segundo a sua maximiza¸c˜ao. Para garantir a existˆencia de convergˆencia no algoritmo, ´e necess´ario ter aten¸c˜ao aos valores iniciais. A cada itera¸c˜ao existe um aumento da verosimilhan¸ca (Dempster et al., 1977; David e Giltinian, 1993). O M´etodo de Newton-Raphson recorre `a derivada do logaritmo da fun¸c˜ao verosimilhan¸ca, fun¸c˜ao Score, que, igualando a zero, resulta num sistema de equa¸c˜oes, geralmente n˜ao lineares. Em cada itera¸c˜ao ´e necess´ario o c´alculo da fun¸c˜ao Score e da sua derivada. O c´alculo desta ´ultima ´e bastante complexo, existem estudos que tentam reduzir essa complexidade, nomeadamente o M´etodo Quasi-Newton (Thisted, 1988). No Rest´a implementada uma metodologia em que se inicia com as estimativas utilizando o algoritmo EM e, quando perto do ´otimo, utiliza-se o M´etodo de Newton-Raphson. 4.2.5 Matriz Variˆancia-Covariˆancia dos Erros Aleat´orios O Modelo Misto permite flexibilizar a estrutura dos efeitos aleat´orios, sob a condi¸c˜ao da estrutura dos erros aleat´orios Σi=σ2Imi. Contudo, esta condi¸c˜ao nem sempre se verifica, pois, na maioria dos casos, as observa¸c˜oes est˜ao correlacionadas. A heterocedasticidade dos erros aleat´orios pode ser modelada se for considerado outro tipo de estrutura para Σi. Considerando a equa¸c˜ao do Modelo de Efeitos Mistos, a estrutura para os erros aleat´orios ´e dada por i∼N(0, σ2Λi)i= 1, . . . , n, (4.110) em que Λi´e uma matriz definida positiva parametrizada por um n´umero geralmente pequeno de parˆametros que se designa por λ. Dado que admite raiz quadrada invert´ıvel, Λ−1/2 i, tem-se que Λi= (Λ1/2 i)TΛ1/2 iΛ−1 i=Λ−1/2 i(Λ−1/2 i)T.(4.111) 71 Cap´ıtulo 4. Modelos Considerando a reparametriza¸c˜ao Y∗ i= (Λ−1/2 i)TYi,X∗ i= (Λ−1/2 i)TXi,Z∗ i= (Λ−1/2 i)TZie∗ i= (Λ−1/2 i)Ti, o modelo pode ser escrito da seguinte forma Y∗ i=X∗ iβ+Z∗ iu∗ i+∗ ii= 1, . . . , n, (4.112) em que ∗ i∼N(0, σ2I) e u∗ i∼N(0,D). De notar que E[∗ i]=(Λ−1/2 i)TE[i] e var(∗ i) = Λ−1/2 ivar[i](Λ−1/2 i)T=σ2I. A fun¸c˜ao de verosimilhan¸ca para o modelo ´e dado por LML(y;β, θ, σ2, λ) = n Y i=1 f(yi;β, θ, σ2, λ) = n Y i=1 f(y∗ i;β, θ, σ2, λ)Λ−1/2 i =LML(y∗;β, θ, σ2, λ) n Y i=1 Λ−1/2 i, (4.113) e a fun¸c˜ao de verosimilhan¸ca restrita ´e dada por LREML(y;β, θ, σ2, λ) = LREML(y∗;β, θ, σ2, λ) n Y i=1 Λ−1/2 i.(4.114) A matriz de variˆancia-covariˆancia Λipode ser decomposta num produto de matrizes Λi=WiCiWi,(4.115) em que WieCis˜ao uma matriz diagonal (variˆancia) e uma matriz de correla¸c˜ao, respetivamente. Para que a matriz Wiseja ´unica, ´e necess´ario impor que Witenha todos os elementos da diagonal principal positivos, por outro lado, que var(it) = σ2[Wi]2 tt e que corr(it, il)=[Ci]2 tl pelo que Widescreve a variˆancia e Cidescreve a correla¸c˜ao dos erros identro do grupo. Esta decomposi¸c˜ao permite a flexibilidade nestas estruturas e, consequentemente, ´e poss´ıvel definir uma estrutura de correla¸c˜ao e modelar a dependˆencia. Cressie (1993) mostra que a estrutura da modela¸c˜ao da dependˆencia dos erros aleat´orios dentro do grupo ´e isotr´opica. Considerando dois erros, a sua correla¸c˜ao depende da distˆancia entre eles, isto ´e, corr(it, il) = h[d(pit,pil),ρ]i= 1, . . . , n t = 1, . . . , tmi,(4.116) em que ρ´e um vetor de parˆametros de correla¸c˜ao e h(.) ´e a fun¸c˜ao de correla¸c˜ao (ou fun¸c˜ao de autocorrela¸c˜ao) e varia entre -1 e 1 e tal que h(0,ρ) = 1. A fun¸c˜ao de autocorrela¸c˜ao emp´ırica ´e uma estimativa n˜ao param´etrica da fun¸c˜ao h(.). 72 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Supondo que rit = (yit −ˆyit)/ˆσit, em que ˆσit ´e o estimador da variˆancia do erro aleat´orio do indiv´ıduo i, no tempo t. Define-se a fun¸c˜ao de autocorrela¸c˜ao emp´ırica no espa¸camento l, dado por ˆρ(l) = Pn i=1 Ptmi t=1 ritri(t+l)/N(l) Pn i=1 Ptmi t=1 r2 it/N(0) ,(4.117) em que N(l) representa o n´umero de pares de res´ıduos que s˜ao utilizados no somat´orio do numerador da fun¸c˜ao. De notar que quando os valores se aproximam gradualmente de zero, o processo pode ser identificado como autoregressivo e, se a fun¸c˜ao supracitada for consistente no intervalo ±z1−α/2 pN(l)z1−α/2 , ap´os o espa¸camento 2, pode ser um poss´ıvel processo de m´edias m´oveis 1 ou 2. Existem v´arias estruturas de correla¸c˜ao temporal, como a correla¸c˜ao n˜ao estruturada, a simetria composta e a autorregressiva. N˜ao estruturada Na estrutura de correla¸c˜ao n˜ao estruturada ´e assumido um parˆametro para cada correla¸c˜ao entre as observa¸c˜oes, dada por h(k, ρ) = ρk, k = 1,2, . . . (4.118) Esta correla¸c˜ao ´e ´util apenas quando se tem poucos dados. Simetria Composta Na estrutura de correla¸c˜ao sim´etrica composta assume-se uma correla¸c˜ao igual entre todos os erros do mesmo grupo, dada por corr(it, il) = ρ, ∀t6=l, h(k, ρ) = ρ, k = 1,2, . . .,(4.119) em que ρrepresenta o coeficiente de correla¸c˜ao intragrupo e ´e ´unico. Autorregressiva Na estrutura de correla¸c˜ao autorregressiva ´e assumido um parˆametro que depende do espa¸camento, isto ´e, t=φ1t−1+···+φpt−p+at,(4.120) em que at´e um processo de ru´ıdo branco gaussiano. Este processo denota-se por AR(p), em que p´e a ordem do processo, tal que Φ= (φ1, . . . , φp). 73 Cap´ıtulo 4. Modelos Um exemplo ´e o processo autorregressivo de ordem 1, AR(1), em que os erros no tempo ts˜ao modelados em fun¸c˜ao dos erros no tempo t−1, tal que t=φ1t−1+at,|φ|<1.(4.121) A fun¸c˜ao de correla¸c˜ao ´e dada por h(k, φ) = φk, k = 0,1, . . . , (4.122) ou seja, decresce exponencialmente em valor absoluto com o espa¸camento (lag). Este processo ´e descrito com maior detalhe no Cap´ıtulo 2. 4.2.6 Inferˆencia Estat´ıstica O ajustamento do modelo marginal aos dados permite a inferˆencia dos parˆametros do modelo de modo a que os resultados obtidos possam ser generalizados para a popula¸c˜ao a partir da qual a amostra foi obtida. Distribui¸c˜ao Assint´otica Considerando o modelo marginal (4.74), condicionada a αconhecido, sabe-se que o estimador β ˆ β(α) = n X i=1 XT iV−1 iXi!−1n X i=1 XT iV−1 iYi(4.123) tem Distribui¸c˜ao Normal Multivariada com valor esperado E(ˆ β(α)) = n X i=1 XT iV−1 iXi!−1n X i=1 XT iV−1 iE[Yi] = β(4.124) e matriz variˆancia-covariˆancia dada por var(ˆ β(α)) = n X i=1 XT iV−1 iXi!−1n X i=1 XT iV−1 ivar(Yi)V−1 iXi × n X i=1 XT iV−1 iXi!−1 = n X i=1 XT iV−1 iXi!−1 . (4.125) No caso em que α´e desconhecido, ´e necess´ario considerar a sua distribui¸c˜ao assint´otica. 74 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Sob os pressupostos apresentados por Pinheiro & Bates (1995, 2000), os estimadores ML s˜ao consistentes e a sua distribui¸c˜ao assint´otica ´e Normal Multivariada. A inversa da matriz de informa¸c˜ao de Fisher aproxima a matriz de variˆancia-covariˆancia dos estimadores (Cox & Hinkley, 1974; Pinheiro & Bates, 2000). Tendo em conta que E"∂2` ∂β∂θT#= 0 E"∂2` ∂β∂σ2#= 0,(4.126) os estimadores dos efeitos fixos, atrav´es de M´etodo de M´axima Verosimilhan¸ca, n˜ao s˜ao assintoticamente correlacionados com os estimadores ML de θe de σ2. A sua distribui¸c˜ao assint´otica ´e dada por ˆ β∼Nβ, σ2XTM−1(θ)X−1(4.127) e ˆ θ logˆσ!∼N θ logσ!,I−1(θ, σ)!,(4.128) em que `(θ, σ2) ´e o logaritmo da fun¸c˜ao verosimilhan¸ca marginal dos efeitos fixos e I(θ, σ) ´e a matriz emp´ırica de Informa¸c˜ao de Fisher e ´e dada por I(θ, σ) =     ∂2`(θ, σ2) ∂θ∂θT ∂2`(θ, σ2) ∂logσ∂θT ∂2`(θ, σ2) ∂θ∂logσ ∂2`(θ, σ2) ∂2logσ    .(4.129) Recorre-se a log σ em vez de σpara simplifica¸c˜ao de parametriza¸c˜ao. Tal transforma- ¸c˜ao faz com que a aproxima¸c˜ao `a Distribui¸c˜ao Normal seja mais evidente. ` A semelhan¸ca dos estimadores de M´axima Verosimilhan¸ca, os obtidos pelo REML s˜ao tamb´em consistentes e com Distribui¸c˜ao Assint´otica Normal Multivariada. Pinheiro (1994) prova esta conclus˜ao a partir de equa¸c˜oes semelhantes `as anteriores, com a diferen¸ca no logaritmo da fun¸c˜ao verosimilhan¸ca, a qual passa a ser restrita. Na maioria dos casos, os parˆametros s˜ao desconhecidos e procede-se `a sua implementa¸c˜ao por aproxima¸c˜ao (4.128). Esta constata¸c˜ao ´e importante, na medida em que suporta os testes e intervalos de confian¸ca para os efeitos fixos. Efeitos fixos O Teste da Raz˜ao de Verosimilhan¸ca ´e aplic´avel na compara¸c˜ao de modelos encaixados, apenas definindo a estrutura de efeitos fixos. A estat´ıstica de teste ´e dada por 2logL1 L0= 2(logL1−logL0),(4.130) 75 Cap´ıtulo 4. Modelos em que L1´e a verosimilhan¸ca do modelo geral (com mais parˆametros) e L0´e a verosimilhan¸ca do modelo encaixado. Considerando a qualidade do ajustamento, para os dois modelos, a estat´ıstica de teste segue assintoticamente a distribui¸c˜ao do Qui-Quadrado com k1−k0graus de liberdade, em que k1−k0´e a diferen¸ca entre o n´umero de parˆametros de cada modelo. O presente teste s´o ´e poss´ıvel ser efetuado se os estimadores dos modelos forem obtidos a partir do M´etodo de M´axima Verosimilhan¸ca, na medida em que o logaritmo da fun¸c˜ao verosimilhan¸ca restrita se altera se as condi¸c˜oes dos efeitos fixos forem alterados. Dado que este teste conduz a valor de prova inferiores ao verdadeiro, esta imprecis˜ao faz com que autores citem outros testes para avalia¸c˜ao da significˆancia dos efeitos fixos, como o teste-t e o teste-F aproximado. Atrav´es do teste-t aproximado pretende-se testar H0:βj= 0 H1:βj6= 0 , j = 1, . . . , p, (4.131) com base na estat´ıstica de teste ET =ˆ βj ˆσREMLv u u t"XT iM−1 i(ˆ θ)Xi−1#jj ,(4.132) sob hip´otese nula, a estat´ıstica de teste tem distribui¸c˜ao assint´otica `a t de Student com glj. Este teste permite a avalia¸c˜ao do coeficiente do efeito fixo j, quando os restantes est˜ao presentes no modelo. Atrav´es do teste-F aproximado pretende-se testar H0:Lβ=0H1:Lβ6=0,(4.133) com base na estat´ıstica de teste ET = ˆ βTLThLˆσREMLXT iM−1 i(ˆ θ)Xi−1LTiLˆ β c(L),(4.134) sob hip´otese nula, a estat´ıstica de teste tem distribui¸c˜ao assint´otica F de Snedcor com (l, k) graus de liberdade, em que lrepresenta o n´umero de graus de liberdade do numerador dado pela caracter´ıstica da matriz L,c(L). Este teste permite a avalia¸c˜ao dos coeficientes dos efeitos fixos, presentes no modelo. A constru¸c˜ao dos intervalos de confian¸ca aproximados para os efeitos fixos ´e baseada nos testes-t aproximados. Considerando gljo n´umero de graus de liberdade do denominador teste-t aproximado, relativamente ao efeito fixo j, ao n´ıvel de confian¸ca 1 −α, o 76 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais intervalo de confian¸ca para βj´e dado por #ˆ βj−tβj,1−α/2ˆσREMLv u u t"XT iM−1 i(ˆ θ)Xi−1#jj , ˆ βj+tβj,1−α/2ˆσREMLv u u t"XT iM−1 i(ˆ θ)Xi−1#jj", (4.135) em que tβj,1−α/2representa o quantil 1 −α/2 da distribui¸c˜ao t de Student com graus de liberdade glj. Componentes de Variˆancia Para al´em de estudar o comportamento m´edio da popula¸c˜ao, existe a necessidade de clarificar quais os efeitos aleat´orios que devem ser inclu´ıdos no modelo e a estrutura de correla¸c˜ao a adotar. O Teste da Raz˜ao de Verosimilhan¸ca ´e utilizado quando os parˆametros s˜ao estimados, usando o m´etodo REML, j´a que a estrutura fixa ´e a mesma nos dois modelos a comparar. Uma das condi¸c˜oes exigidas ´e o facto da estat´ıstica de teste ter distribui¸c˜ao assint´otica Qui-Quadrado com graus de liberdade igual `a diferen¸ca entre o especificado nas hip´oteses alternativa e nula. No Rest´a implementada uma metodologia sugerida por Pinheiro & Bates (2000). Esta solu¸c˜ao calcula os valores de prova de forma sobrevalorizada, o que conduz a uma an´alise conservativa. O intervalo de confian¸ca aproximado ao n´ıvel αpara o desvio padr˜ao σ´e dado por #ˆσexpn−z1−α/2p[I−1]σσ,ˆσexpnz1−α/2p[I−1]σσo",(4.136) em que z1−α/2representa o quartil 1 −α/2 da distribui¸c˜ao Normal padr˜ao. Para as componentes da matriz de variˆancia-covariˆancia dos efeitos aleat´orios, os intervalos de confian¸ca aproximados s˜ao um pouco mais dif´ıcieis de construir e s˜ao estimados com menor precis˜ao do que os dos efeitos fixos e o do desvio padr˜ao dentro dos grupos. O aumento da precis˜ao na estima¸c˜ao dos primeiros s´o ´e poss´ıvel com o aumento do n´umero de grupos estudados (Pinheiro & Bates, 2000). Efeitos Aleat´orios O valor esperado do BLUP de ui´e dado por E[uBLUP (α)] = E[DZT iV−1 i(Yi−Xiˆ β)] = DZT iV−1 i(Xiβ−Xiˆ β)] = 0 (4.137) 77 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais gelo encontrarem temperaturas positivas, os flocos de neve derretem e transformam-se em chuva. A varia¸c˜ao da precipita¸c˜ao depende da regi˜ao, como a exposi¸c˜ao geogr´afica, a altitude, a proximidade do mar e a distribui¸c˜ao de zonas de press˜ao e de frente polar. A precipita¸c˜ao ´e uma vari´avel que apresenta uma enorme variabilidade espacial e temporal. A precipita¸c˜ao ´e um fen´omeno espacialmente distribu´ıdo de natureza cont´ınua que apenas ´e avaliado em localiza¸c˜oes pontuais atrav´es das esta¸c˜oes climatol´ogicas. A medi¸c˜ao da precipita¸c˜ao tem como intuito a recolha de dados sobre a quantidade de precipita¸c˜ao, durante um determinado intervalo de tempo, quando existe queda de ´agua. A amostra recolhida deve ser representativa da quantidade real de precipita¸c˜ao, que ocorreu na regi˜ao em estudo. A unidade de medida da intensidade de precipita¸c˜ao ´e o mil´ımetro (mm) que corresponde `a altura de ´agua de 1 litro por metro quadrado. Para determinar a intensidade de precipita¸c˜ao utilizam-se aparelhos para a medi¸c˜ao de precipita¸c˜ao, como os ud´ometros (Figura 5.3). Como qualquer sistema de medi¸c˜ao, os ud´ometros est˜ao sujeitos a certos tipos de erros que afetam a fiabilidade da estima¸c˜ao de precipita¸c˜ao. Os ud´ometros s˜ao constitu´ıdos por um cilindro com uma abertura superior e um recipiente de recolha de ´agua, que originam usualmente boas estimativas pontuais. No entanto, sob fatores extremos, como precipita¸c˜ao intensa ou ventos fortes, as medi¸c˜oes podem conter erros de medi¸c˜ao. A sua localiza¸c˜ao ´e bastante importante, normalmente localizam-se em locais descampados e abrigados do vento, parcialmente enterrados no solo ou fixos por estruturas met´alicas. A prote¸c˜ao do aparelho ´e importante para reduzir os efeitos aerodinˆamicos que levam a uma subestima¸c˜ao ou sobreestima¸c˜ao da medi¸c˜ao. Em Portugal, a rede udom´etrica ´e bastante limitada, na medida em que a cobertura de uma regi˜ao por ud´ometros ´e bastante dispendiosa e de manuten¸c˜ao e opera¸c˜ao complexas. A quantifica¸c˜ao espacial e temporal da precipita¸c˜ao sobre uma ´area geogr´afica, em particular, ´e imprescind´ıvel no c´alculo dos balan¸cos h´ıdricos para a estima¸c˜ao indireta de caudais em cursos de ´agua e para o desenvolvimento do estudo de carga dos aqu´ıferos. ´ E ainda essencial `a modela¸c˜ao de diversos fen´omenos ambientais, entre os quais a varia¸c˜ao da Qualidade da ´ Agua numa bacia hidrogr´afica fluvial. A Qualidade da ´ Agua, num determinado local, ´e o reflexo das condi¸c˜oes dominantes da bacia de alimenta¸c˜ao desse local, nomeadamente de fatores hidrometeorol´ogicos. A varia¸c˜ao espa¸co-temporal de uma vari´avel de qualidade est´a associada `a varia¸c˜ao do caudal e esta acompanha geralmente a varia¸c˜ao sazonal da precipita¸c˜ao. Da´ı a necessidade de estimar um fator hidrometeorol´ogico (via precipita¸c˜ao) nas esta¸c˜oes de amostragem de Qualidade da ´ Agua na bacia do rio Douro, onde n˜ao h´a nem valores de caudais nem valores de precipita¸c˜ao. A bacia hidrogr´afica do rio Douro ´e monitorizada por v´arias esta¸c˜oes de amostragem 84 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Figura 5.3: Esquerda: Ud´ometro, utilizado para medir a precipita¸c˜ao total que caiu num determinado per´ıodo de tempo, 28 de setembro de 1992, em Vale dos Camelos (DSRH/INAG). Direita: Ud´ometro entupido e com ´agua acumulada da Esta¸c˜ao Autom´atica sem telemetria, 28 de dezembro de 2012, em Laranjal, Ponte Sˆor (DMSDIH/APA). de precipita¸c˜ao distribu´ıdas pelo rio Douro e pelos seus principais afluentes. No entanto, s´o foram consideradas 18 esta¸c˜oes de amostragem devido `a enorme falta de dados (dados omissos) em grande parte delas. As esta¸c˜oes selecionadas apresentam uma percentagem de dados omissos inferior a cerca de 20 %. Nas Tabelas 5.1 e 5.2 apresentam-se as caracter´ısticas e o per´ıodo de observa¸c˜ao destas esta¸c˜oes de amostragem e a sua localiza¸c˜ao na Figura 5.4. O acesso aos dados mensais foi realizado a partir do Sistema Nacional de Informa¸c˜ao de Recursos H´ıdricos (SNIRH). 5.1.1 An´alise Descritiva A Tabela 5.2 sintetiza algumas estat´ısticas descritivas da precipita¸c˜ao no per´ıodo observado. ´ E poss´ıvel verificar que quanto maior a m´edia, maior ´e o valor do desvio padr˜ao. ´ E tamb´em vis´ıvel que o valor m´ınimo ´e zero em todas as esta¸c˜oes, o que significa que existiu pelo menos um mˆes em que n˜ao ocorreu precipita¸c˜ao em cada esta¸c˜ao de amostragem. O maior valor de precipita¸c˜ao observado (408,50 mm) ´e identificado na esta¸c˜ao de amostragem de Santa Marta da Montanha, em novembro de 2002. A esta¸c˜ao de amostragem de Vilar Formoso ´e a que apresenta menor desvio padr˜ao (31,74 mm) e menor coeficiente de varia¸c˜ao (84 %). A esta¸c˜ao de amostragem de Santa Maria da Montanha ´e a que apresenta maior desvio padr˜ao (93,34 mm). A esta¸c˜ao de amostragem de Fonte Longa ´e a que revela maior coeficiente de varia¸c˜ao (112 %). No Apˆendice A procedeu-se `a representa¸c˜ao gr´afica da evolu¸c˜ao da precipita¸c˜ao ao longo do tempo observado em cada 85 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais Tabela 5.1: As esta¸c˜oes de amostragem da precipita¸c˜ao, na Bacia Hidrogr´afica do rio Douro, e o respetivo per´ıodo observado. C´odigo Nome Latitude (oN) Longitude (oO) Per´ıodo Observado 1 10O/02UG Almeidinha 40,59400 -7,12900 12/2003 - 1/2013 2 04R/01UG Argozelo 41,63800 -6,60100 1/2003 - 2/2013 3 07O/05UG Castelo Melhor 41,01546 -7,06808 3/2002 - 11/2009 4 08J/06G Castro Daire (Lamelas) 40,92400 -7,94000 3/2002 - 9/2011 5 02R/02G Deil˜ao 41,84700 -6,58600 10/2002 - 2/2013 6 08P/02G Escalh˜ao 40,94800 -6,92400 3/2002 - 2/2013 7 06N/03UG Fonte Longa 41,23200 -7,26800 3/2002 - 10/2012 8 08I/01UG Mosteiro de Cabril 40,94700 -8,10000 3/2002 - 2/2012 9 04R/02G Pinelo 41,63500 -6,55200 3/2002 - 2/2013 10 09O/01G Pinhel 40,77100 -7,06100 3/2002 - 2/2013 11 04K/02G Santa Maria da Montanha 41,50075 -7,74599 3/2002 - 9/2011 12 04O/01G Torre de Dona Chama 41,65654 -7,11589 3/2002 - 9/2010 13 08K/01UG Touro 40,89700 -7,74800 3/2002 - 12/2009 14 03N/01G Travancas 41,82797 -7,30561 3/2002 - 6/2012 15 05M/04UG Vales (Valpa¸cos) 41,46547 -7,35143 3/2002 - 5/2012 16 05L/01UG Vila Pouca de Aguiar 41,49806 -7,63600 6/2003 - 5/2012 17 10Q/01UG Vilar Formoso 40,60900 -6,83100 3/2002 - 5/2012 18 02O/02UG Vinhais 41,82798 -6,99384 11/2003 - 2/2013 Figura 5.4: Representa¸c˜ao da bacia hidrogr´afica do rio Douro e localiza¸c˜oes das esta¸c˜oes de medi¸c˜ao de precipita¸c˜ao. esta¸c˜ao de amostragem, bem como os diagramas em caixa de bigodes e os histogramas. Pode observar-se que existem esta¸c˜oes com outliers, valores discrepantes que representam meses com grande intensidade de precipita¸c˜ao. Estas observa¸c˜oes n˜ao s˜ao exclu´ıdas, na medida em que mostram o comportamento extremo de precipita¸c˜ao que foi atingido. 86 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Tabela 5.2: Medidas descritivas da vari´avel da precipita¸c˜ao, no per´ıodo observado. C´odigo N´umero Dados M´ınimo M´aximo 1.oQuartil Mediana 3.oQuartil M´edia Desvio Coeficiente Omissos (mm) (mm) (mm) (mm) (mm) (mml) Padr˜ao de Varia¸c˜ao (mm) (%) 10O/02UG 110 24 0,00 196,00 11,70 32,15 54,88 41,30 39,96 0,97 04R/01UG 98 36 0,00 219,70 11,15 34,25 66,45 46,99 46,17 0,98 07O/05UG 95 39 0,00 157,40 11,40 25,70 44,00 34,87 33,88 0,97 08J/06G 107 27 0,00 366,90 25,20 63,00 139,10 93,25 91,34 0,98 02R/02G 99 35 0,00 227,20 12,35 33,60 69,55 49,49 48,08 0,97 08P/02G 126 8 0,00 239,00 14,50 33,15 59,35 42,84 39,41 0,92 06N/03UG 106 28 0,00 233,00 7,50 26,65 51,45 39,57 44,42 1,12 08I/01UG 102 32 0,00 329,80 23,98 52,30 107,40 78,20 74,76 0,96 04R/02G 124 10 0,00 201,00 7,30 29,30 58,20 43,53 47,44 1,09 09O/01G 130 4 0,00 166,20 7,85 25,50 52,92 35,68 36,40 1,02 04K/02G 105 29 0,00 408,50 23,50 57,50 129,90 92,98 93,57 1,01 04O/01G 104 30 0,00 139,70 12,55 26,60 49,95 36,63 33,30 0,91 08K/01UG 95 39 0,00 223,10 12,30 28,40 70,95 47,73 48,82 1,02 03N/01G 103 31 0,00 269,00 13,90 37,20 74,95 56,08 59,08 1,05 05M/04UG 109 25 0,00 172,90 7,80 23,70 48,10 36,26 39,31 1,08 05L/01UG 96 38 0,00 354,60 16,85 40,65 100,55 67,21 71,17 1,06 10Q/01UG 112 22 0,00 128,00 13,75 30,65 50,70 37,72 31,74 0,84 02O/02UG 109 25 0,00 211,20 11,40 29,40 56,60 41,72 44,40 1,06 5.1.2 An´alise da Continuidade Espacial No presente estudo existe a necessidade de dispor de medi¸c˜oes mensais m´edias em ´area, dependentes da quantidade de precipita¸c˜ao que cai sobre uma certa regi˜ao geogr´afica, que influenciam um determinado local da bacia hidrogr´aficado rio Douro e que funcionam como um fator hidrometeorol´ogico na modela¸c˜ao (espacial e temporal) da Qualidade da ´ Agua. Pretende-se a identifica¸c˜ao de modelos que possibilitem a obten¸c˜ao dessas medi¸c˜oes, construindo a distribui¸c˜ao espacial da precipita¸c˜ao, na bacia do rio Douro, em locais onde n˜ao h´a valores observados (nas esta¸c˜oes de amostragem de Qualidade, onde foram efetuadas medi¸c˜oes das vari´aveis de Qualidade da ´ Agua de superf´ıcie). Foi realizado o estudo dos dados de precipita¸c˜ao espaciais atrav´es de gr´aficos que resumem as observa¸c˜oes de acordo com os seus principais atributos (gr´afico precipita- ¸c˜ao/Latitude e gr´afico Longitude/precipita¸c˜ao). A an´alise da variabilidade espacial da precipita¸c˜ao ´e realizada separadamente para cada mˆes do ano e foi necess´ario transformar as coordenadas devido `a presen¸ca de uma componente de tendˆencia polinomial. Ap´os a transforma¸c˜ao, as observa¸c˜oes transformadas apresentam-se j´a de forma aleat´oria e os gr´aficos tamb´em j´a indicam estacionaridade na m´edia e distribui¸c˜oes gaussianas da precipita¸c˜ao mensal, ao longo dos anos observados. Falta verificar a condi¸c˜ao de isotropia dos processos subjacentes `a precipita¸c˜ao. A isotropia verifica-se quando os valores do processo estoc´astico dependem apenas do m´odulo do vetor da distˆancia entre eles, n˜ao podendo depender da dire¸c˜ao angular desses mesmos vetores. Foram calculados os variogramas direcionais e estes apresentam uma forma semelhante de conduta, dependendo da dire¸c˜ao angular, at´e aproximadamente `a distˆan87 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais cia m´axima considerada. As dire¸c˜oes 0oe 45otˆem valores de semicovariˆancia um pouco superiores a zero e as dire¸c˜oes 90oe 135otˆem valor aproximadamente de um. Os dados exp˜oem-se de forma aleat´oria, na superf´ıcie em estudo; nos variogramas direcionais, a magnitude da semicovariˆancia ´e bastante reduzida, comparativamente ao variograma com tendˆencia na m´edia. Considera-se que o pressuposto de isotropia ´e cumprido. Assumem-se as hip´oteses de homogeneidade do processo (na regi˜ao em estudo, o processo ´e estacion´ario de segunda ordem, logo intrinsecamente estacion´ario) e a isotropia. Os dados s˜ao analisados atrav´es de diferentes medidas de continuidade espacial (o covariograma, o semivariograma e o correlograma emp´ıricos). Modela¸c˜ao da Continuidade Espacial O variograma emp´ırico ´e bastante usado na an´alise explorat´oria de dados, para al´em de poder ser utilizado para descobrir o modelo de correla¸c˜ao espacial. Para a constru¸c˜ao de um variograma definem-se duas regras muito utilizadas na pr´atica: a distˆancia m´axima ´e considerada cerca de 60 % da distˆancia m´axima das localiza¸c˜oes e cada ponto do semivariograma tem que ser calculado atrav´es de, pelo menos, 30 pares de observa¸c˜oes. A no¸c˜ao de continuidade espacial adv´em da hip´otese de que as observa¸c˜oes pr´oximas no espa¸co tendem a ser mais semelhantes do que as mais afastadas. Para a an´alise desta associa¸c˜ao entre as medi¸c˜oes existem diversos m´etodos como o correlagrama, o covariograma e o semivariograma. No presente trabalho recorre-se ao semivariograma emp´ırico determinado pelas equa¸c˜oes mencionadas no Cap´ıtulo 2, supondo que cada medi¸c˜ao ao longo do tempo s˜ao r´eplicas independentes do mesmo processo. De notar que a estima¸c˜ao do variograma pode ser prejudicada pela presen¸ca de outliers, mas opta-se por manter os outliers. O variograma de Matheron estima os valores com base no c´alculo de mais de 30 pares de pontos. Ap´os a determina¸c˜ao do variograma emp´ırico para cada mˆes do ano, o objetivo ´e ajustar o modelo te´orico que mais se adeque. Esta etapa consiste em estimar os parˆametros desconhecidos dos modelos te´oricos. No ajustamento das observa¸c˜oes aos modelos de transi¸c˜ao foram considerados o Modelo Exponencial, o Modelo Esf´erico e o Modelo Gaussiano. Apresentam-se graficamente o semivariograma estimado, bem como os ajustamentos para cada um dos modelos de transi¸c˜ao anteriores (Figuras 5.5 e 5.6). O ajustamento de todos os modelos de transi¸c˜ao considerados foi efetuado pelo M´etodo dos M´ınimos Quadrados. 88 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Figura 5.5: Representa¸c˜oes gr´aficas dos semivariogramas emp´ıricos e dos ajustamentos aos modelos te´oricos, pelo M´etodo dos M´ınimos Quadrados. 89 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais Figura 5.6: Representa¸c˜oes gr´aficas dos semivariogramas emp´ıricos e dos ajustamentos aos modelos te´oricos, pelo M´etodo dos M´ınimos Quadrados. 90 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Tabela 5.3: Valores estimados para cada semivariograma e m´etricas utilizadas para a Valida¸c˜ao Cruzada. Mˆes Modelo ˆ C0ˆ C1ˆ C2ˆ MEP ˆ EQMN janeiro Exponencial 2694,0494 3755,9527 1,4050 -0,7073 0,9935 Esf´erico 3188,9815 808,3802 0,0065 15,9762 0,9976 Gaussiano 2969,6354 31464,2498 3,8027 -2,7856 0,9937 fevereiro Exponencial 1680,5434 1713681,7978 1901,9569 -2,3646 0,9947 Esf´erico 2103,4303 130,2867 0,0100 2,3457 0,9940 Gaussiano 1811,9229 1685,9462 1,2240 -4,1105 0,9955 mar¸co Exponencial 1623,2178 2662,1474 2,0588 0,1131 0,9937 Esf´erico 1624,4941 1924,5520 2,5823 -0,0198 0,9937 Gaussiano 1759,6503 1759,6503 1,0821 -0,7046 0,9938 abril Exponencial 1414,0899 1949226,7488 2771,4596 -0,0128 0,9937 Esf´erico 343,8506 1493,8159 0,0000 7,3619 0,9959 Gaussiano 1500,5562 1364,5816 1,2145 0,2645 0,9937 maio Exponencial 772,1211 350925,3785 2283,3108 -0,3200 0,9940 Esf´erico 801,4512 68,2032 0,0096 -0,0286 0,9940 Gaussiano 795,5256 65286,5162 20,7456 -0,6595 0,9940 junho Exponencial 540,7291 282418,4835 2827,5728 -1,1848 0,9944 Esf´erico 601,1948 0,1673 0,0100 -2,1623 0,9946 Gaussiano 552,1412 22443,4841 14,4906 -1,3728 0,9943 julho Exponencial 156,6149 87427,6303 3284,5601 0,5631 0,9942 Esf´erico 156,8556 17,0238 0,0100 0,6902 0,9944 Gaussiano 156,7661 17,1760 0,0430 0,7238 0,9945 agosto Exponencial 1369,6930 1194,9599 2,5556 1,6822 0,9944 Esf´erico 1424,5478 203,8984 0,0071 1,7706 0,9945 Gaussiano 1440,1172 116374,8005 17,4713 2,1366 0,9944 setembro Exponencial 853,8326 386767,6486 2267,6718 -0,9458 0,9940 Esf´erico 167,0694 789,8636 0,0010 -1,1823 0,9939 Gaussiano 742,3707 214,5623 0,0010 -1,1823 0,9939 outubro Exponencial 3967,3306 2333780,5287 1058,1762 2,6941 0,9941 Esf´erico 4961,8748 356,8896 0,0105 11,1138 0,9971 Gaussiano 4253,6523 384908,2541 12,9947 2,0860 0,9937 novembro Exponencial 4267,3646 2065591,9849 1455,6130 1,8141 0,9937 Esf´erico 4818,1061 324,7365 0,0103 8,7028 0,9963 Gaussiano 4481,9407 157796,2363 10,6039 -0,5441 0,9936 dezembro Exponencial 3280,7657 3133396,5178 1620,6543 -1,5792 0,9935 Esf´erico 3291,6828 1231,9447 0,0007 3,1170 0,9943 Gaussiano 3484,0075 2622,0005 0,9267 -2,0149 0,9935 91 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais Valida¸c˜ao Cruzada Na presen¸ca de v´arios modelos te´oricos para ajustar um variograma te´orico aos dados emp´ıricos, o M´etodo de Valida¸c˜ao Cruzada permite diagnosticar eventuais problemas com o ajustamento dos modelos escolhidos, indicando qual ser´a o mais adequado. Assim sendo, para avaliar a qualidade de ajustamento de cada um dos variogramas obtidos, analisam-se os valores da m´edia dos erros de predi¸c˜ao (MEP) e do erro quadr´atico m´edio normalizado (EQMN). Os dois primeiros valores dever˜ao ser aproximadamente zero e o ´ultimo dever´a ser aproximadamente um. Por exemplo, no mˆes de janeiro, na Tabela 5.3, verifica-se que no modelo exponencial, a medida de avalia¸c˜ao MEP ´e pr´oximo de zero e o EQMN ´e pr´oximo de um, relativamente aos restantes modelos. Com base nestas m´etricas, o modelo te´orico que melhor se ajusta ao variograma emp´ırico ´e o Modelo Exponencial, atrav´es do M´etodo dos M´ınimos Quadrados, para o mˆes de janeiro. Nos restantes meses, o Modelo Exponencial, ´e tamb´em o modelo te´orico que melhor se ajusta ao variograma emp´ırico. Assim, opta-se pelo Modelo Exponencial para a modela¸c˜ao espacial, em cada mˆes. 5.1.3 Predi¸c˜ao Pontual e Global Na predi¸c˜ao global (em ´area), o processo que se segue diz respeito `a previs˜ao da precipita¸c˜ao numa certa ´area geogr´afica. Uma vez que a predi¸c˜ao se refere a inferˆencias sobre valores do processo estoc´astico que n˜ao foram observados, ap´os a escolha do modelo te´orico que melhor se ajusta ao variograma emp´ırico, ser´a aplicado o modelo Kriging sobre uma grelha de 200 por 200 (aproximadamente 40000 pontos), delimitada pela fronteira de RH3. Considerando a Latitude e a Longitude onde se localiza a Regi˜ao Hidrogr´afica 3, sabe-se que 1 grau de Longitude ´e aproximadamente 1,11 km, ou seja, o espa¸camento entre dois pontos da grelha ´e aproximadamente 1,11 km. Na aplica¸c˜ao da t´ecnica de Kriging ´e necess´ario ter em conta o modelo te´orico escolhido (neste estudo, o Exponencial), as estimativas dos valores dos parˆametros e o tipo de Kriging utilizado (neste caso, o Universal). O ajustamento do variograma emp´ırico e a interpola¸c˜ao realizada pelo m´etodo de Kriging Universal possibilitaram a constru¸c˜ao de dois gr´aficos, os quais indicam os valores estimados da precipita¸c˜ao na ´area da bacia hidrogr´afica e os desvios padr˜ao associados a essas estimativas. Os resultados da aplica¸c˜ao do m´etodo de Kriging Universal em ´area s˜ao apresentados para o mˆes de janeiro na Figura 5.7 e para os restantes meses nas Figuras A.11 a A.21 (Apˆendice A) e os desvios padr˜ao das estimativas da precipita¸c˜ao na bacia hidrogr´afica do rio Douro s˜ao apresentados para o janeiro na Figura 5.8 e para os restantes meses nas Figuras A.22 a A.32 (Apˆendice A). 92 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Estes gr´aficos sugerem que os n´ıveis de precipita¸c˜ao estimados mais baixos se encontram na zona Noroeste, assinalada pela cor verde e, `a medida que se avan¸ca para Este, surgem n´ıveis de precipita¸c˜ao mais elevados, estando os valores m´aximos assinalados pela cor cinzenta. As linhas desenhadas a preto representam as localidades“contornadas”pelas estimativas do mesmo n´ıvel, ou seja, as zonas interiores `as linhas de contorno assumem estimativas de concentra¸c˜ao de precipita¸c˜ao na ordem do valor apresentado nessa mesma linha. As superf´ıcies s˜ao bastante semelhantes e a estima¸c˜ao parece bem conseguida. O processo de Kriging utilizado permitiu, tamb´em, o mapeamento dos desvios padr˜ao, resultando na representa¸c˜ao dos pontos verde como sendo as zonas de menor desvio padr˜ao, as quais correspondem `as zonas onde foram amostrados mais pontos. Note-se que a existˆencia de desvios padr˜ao relativamente mais elevados pode dever-se ao facto de existirem observa¸c˜oes mais extremas, que resultaram numa estima¸c˜ao menos precisa e, consequentemente, num erro de estima¸c˜ao maior. O objetivo final da an´alise espacial efetuada foi a estima¸c˜ao pontual da precipita¸c˜ao, por mˆes, nas 36 esta¸c˜oes de amostragem de Qualidade da ´ Agua de superf´ıcie da bacia hidrogr´afica do rio Douro. Estes valores, em n´umero elevado, n˜ao s˜ao apresentados na disserta¸c˜ao, pois s˜ao estimativas mensais para os anos de 2002 at´e 2013. O processo de estima¸c˜ao para a precipita¸c˜ao deu origem a uma nova “vari´avel”, CHt, que representa o fator hidrometeorol´ogico e pretende “traduzir” o valor aproximado da medida de um “volume m´edio” de ´agua que passa, por mˆes, na sec¸c˜ao de determinada esta¸c˜ao de amostragem de Qualidade da ´ Agua de superf´ıcie da bacia hidrogr´afica do rio Douro, e que ser´a utilizado nos instantes t,t−1 e t−2 (isto ´e, nos dois meses antes, no mˆes anterior e no pr´oprio mˆes). Consideram-se os trˆes instantes de tempo, porque existem estudos que comprovam que os meses anteriores em que houve precipita¸c˜ao tˆem maior influˆencia na Qualidade da ´ Agua do que o pr´oprio mˆes. 93 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais Tabela 5.6: Medidas descritivas da vari´avel do Oxig´enio Dissolvido, no per´ıodo observado. C´odigo N´umero Dados M´ınimo M´aximo 1.oQuartil Mediana 3.oQuartil M´edia Desvio Coeficiente Omissos (mg/l) (mg/l) (mg/l) (mg/l) (mg/l) (mg/l) Padr˜ao de Varia¸c˜ao (mg/l) (%) 05K/02 125 7 6,50 11,90 8,30 9,20 10,20 9,15 1,21 13,22 04P/01 127 5 6,00 12,20 8,60 9,20 10,00 9,24 1,06 11,45 05R/01 126 6 6,90 12,20 8,30 8,95 9,70 9,08 1,10 12,11 05P/01 126 6 7,10 12,10 8,40 9,20 10,00 9,18 1,03 11,24 07G/04 130 2 4,70 11,40 7,87 8,70 9,88 8,82 1,33 15,11 06N/02 129 3 6,60 11,10 8,50 9,10 9,80 9,16 1,01 11,02 06N/01 127 5 6,00 12,30 8,30 9,20 10,00 9,16 1,11 12,16 10P/02 116 16 3,30 11,70 7,89 8,75 9,83 8,67 1,50 17,33 07M/01 130 2 6,00 11,90 8,20 8,96 9,70 9,04 1,03 11,42 11O/02 117 15 6,20 12,00 8,20 8,70 9,60 8,77 1,07 12,22 08P/02 116 16 5,37 11,00 7,70 8,80 9,43 8,64 1,24 14,31 02Q/02 126 6 5,80 12,70 8,60 9,22 10,14 9,32 1,13 12,17 06K/04 126 6 7,20 11,30 8,50 9,20 10,09 9,25 0,93 10,11 06H/01 129 3 5,80 12,50 8,20 9,20 10,00 9,10 1,20 13,20 06P/02 126 6 7,00 12,20 8,31 9,20 10,10 9,21 1,12 12,20 09O/03 116 16 5,20 12,00 7,80 8,62 9,60 8,62 1,30 15,13 06M/04 128 4 6,80 11,70 8,20 8,90 9,70 8,99 1,08 12,00 08K/01 115 17 4,90 12,40 8,30 9,20 10,10 9,08 1,35 14,89 03M/03 115 17 4,78 12,00 7,10 8,30 9,55 8,25 1,53 18,54 07H/05 130 2 5,10 11,70 8,12 9,05 10,00 9,10 1,26 13,81 04N/05 120 12 5,30 13,80 8,54 9,50 10,16 9,39 1,32 14,03 07K/01 128 4 4,80 12,80 8,40 9,40 10,20 9,30 1,26 13,55 07H/06 129 3 5,40 12,50 7,80 8,90 9,70 8,73 1,34 15,34 08H/02 129 3 4,10 12,30 8,60 9,50 10,30 9,43 1,22 12,98 07G/05 129 3 6,30 11,20 8,00 8,90 9,90 8,92 1,22 13,66 06G/07 118 14 6,80 11,40 8,43 9,30 10,08 9,25 0,99 10,73 07L/01 125 7 5,60 12,70 8,20 9,00 9,90 9,02 1,27 14,08 07K/04 128 4 6,00 12,40 7,97 8,75 10,00 8,88 1,34 15,11 04L/01 112 20 3,10 12,50 7,95 9,25 9,94 8,81 1,76 20,01 07J/02 128 4 5,40 11,90 8,10 9,15 10,01 9,07 1,30 14,38 07H/04 129 3 3,60 12,20 8,80 9,60 10,40 9,58 1,17 12,21 04N/01 127 5 6,40 13,00 8,45 9,60 10,30 9,39 1,24 13,19 06I/04 130 2 6,70 12,80 8,40 9,40 10,30 9,40 1,20 12,76 04N/06 128 4 5,30 12,50 8,40 9,35 10,40 9,33 1,42 15,22 06G/06 127 5 6,70 11,80 8,60 9,30 10,00 9,31 0,96 10,30 04J/09 131 1 7,40 13,50 8,90 9,80 10,40 9,73 1,11 11,37 `a contamina¸c˜ao das massas de ´agua por polui¸c˜ao de origem urbana, industrial e agr´ıcola e `a contamina¸c˜ao de ´aguas subterrˆaneas. No PGRH3 constata-se que nas zonas junto `as localidades de Bragan¸ca e R´egua a falta de ´agua no Ver˜ao ´e mais acentuada e que as sub-bacias, onde as necessidades de ´agua para a ind´ustria s˜ao mais elevadas, s˜ao as do Douro e Costeiras entre o Douro e o Vouga. As massas de ´agua em mau estado/incumprimento, na sua maioria, situam-se nas zonas m´edias e inferiores das principais bacias hidrogr´aficas da RH3, designadamente perto do litoral, nas sub-bacias do Douro, do Tˆamega e do Cˆoa. Nestas localiza¸c˜oes verificam-se as maiores densidades populacionais e ´areas de ocupa¸c˜ao urbana, o que se reflete no valor m´edio do OD observado nas esta¸c˜oes de amostragem perto do litoral e na regi˜ao Nordeste. 100 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Tabela 5.7: Medidas Descritivas sobre Oxig´enio Dissolvido, em fun¸c˜ao do mˆes, no per´ıodo observado. C´odigo N´umero Dados M´ınimo M´aximo 1.ˇz Quartil Mediana 3.ˇz Quartil M´edia Desvio Coeficiente Omissos (mg/l) (mg/l) (mg/l) (mg/l) (mg/l) (mg/l) Padr˜ao de Varia¸c˜ao (mg/l) (%) janeiro 372 24 8,10 12,70 10.00 10.40 10.72 10.38 0.63 9.67 janeiro 372 24 8,10 12,70 10,00 10,40 10,72 10,38 0,63 9,67 fevereiro 381 15 8,30 13,50 10,00 10,40 10,80 10,47 0,72 0,04 mar¸co 384 12 8,10 12,80 9,60 10,00 10,40 9,99 0,70 0,01 abril 369 27 6,30 11,70 9,10 9,60 10,00 9,54 0,76 8,35 maio 388 8 6,50 11,53 8,50 8,90 9,30 8,90 0,78 8,45 junho 382 14 4,78 11,30 7,60 8,30 8,90 8,25 0,92 10,18 julho 375 21 3,10 10,80 7,40 8,00 8,50 7,90 0,94 10,21 agosto 353 43 4,00 10,60 7,50 8,00 8,50 7,97 0,86 9,79 setembro 376 20 3,30 11,20 7,60 8,20 8,60 8,02 1,01 11,01 outubro 374 22 3,60 13,80 8,00 8,50 9,00 8,46 1,05 11,42 novembro 389 7 6,40 11,32 8,90 9,30 9,90 9,30 0,80 9,27 dezembro 355 41 7,31 12,80 9,50 10,00 10,40 10,00 0,83 9,13 Para o estudo da sazonalidade da s´erie recorre-se `a representa¸c˜ao das m´edias mensais da vari´avel OD, por mˆes (Figura 5.11). A evolu¸c˜ao do Oxig´enio Dissolvido em fun¸c˜ao dos meses ´e not´oria: os valores mais elevados s˜ao observados em dezembro, janeiro, fevereiro e mar¸co e os valores mais baixos nos meses de julho, agosto e setembro. Existem trabalhos que apontam para a existˆencia de sazonalidade para este tipo de vari´aveis ambientais (Cabecinha et al., 2009; Costa & Gon¸calves, 2011; Costa & Gon¸calves, 2012; Gon¸calves e Costa, 2013). No per´ıodo em estudo, as s´eries temporais apresentaram valores em falta, mas n˜ao foi realizada qualquer a¸c˜ao para imputar dados. Na Figura 5.13 s˜ao apresentadas as s´eries temporais de cada esta¸c˜ao numa ´unica representa¸c˜ao para a concentra¸c˜ao de Oxig´enio Dissolvido, o perfil temporal da vari´avel OD, onde se verifica um padr˜ao ao longo do tempo. Figura 5.11: Esquerda: Diagramas em caixa de bigodes; Direita: As principais m´etricas da vari´avel OD. 101 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais Figura 5.12: Representa¸c˜ao gr´afica do OD, em fun¸c˜ao das covari´aveis CBO5, Clorofila, pH,Temperatura,CH eALB. 102 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Figura 5.13: Perfil temporal da vari´avel resposta, OD, em fun¸c˜ao do Tempo, no per´ıodo observado. Figura 5.14: Representa¸c˜ao gr´afica dos intervalos de confian¸ca para os coeficientes relativamente `a constante (esquerda) e `a vari´avel tempo (direita), relativamente ao ajustamento linear para cada de amostragem. Na Figura 5.12 representa-se graficamente a rela¸c˜ao entre a concentra¸c˜ao de OD e as 103 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais outras covari´aveis. A inspe¸c˜ao gr´afica sugere que o aumento da Temperatura diminui a concentra¸c˜ao da vari´avel OD. N˜ao ´e clara uma existˆencia de uma associa¸c˜ao entre as vari´aveis ALB,CBO5, Clorofila,pH,Temperatura,CH e o OD. Com o objetivo de identificar os efeitos aleat´orios a considerar, faz-se um ajustamento linear individual e analisam-se quais os parˆametros que mais variam de esta¸c˜ao para esta¸c˜ao. Esta an´alise ´e feita com base no gr´afico da estimativa dos intervalos de confian¸ca para os parˆametros do modelo ajustado a cada esta¸c˜ao (Figura 5.14). Como se pode verificar, h´a variabilidade nas estimativas quer da intersec¸c˜ao quer do declive. Na constru¸c˜ao dos modelos ser˜ao considerados os efeitos aleat´orios no termo constante e no termo da covari´avel Tempo. 5.2.3 Formula¸c˜ao dos Modelos Alpuim & El-Shaarawi (2008) e Gon¸calves & Alpuim (2011) sugerem que uma boa op¸c˜ao para modelar a estrutura de correla¸c˜ao temporal dos erros aleat´orios ´e um processo autorregressivo de ordem 1 (AR(1)), no contexto de dados ambientais. Este pressuposto prende-se com o facto da Qualidade da ´ Agua ser influenciada por condi¸c˜oes que dependem do mˆes anterior. Os modelos em estudo consideram esta estrutura e as representa¸c˜oes do comportamento da FAC e da FACP dos res´ıduos mostrarem que existe uma correla¸c˜ao temporal fraca na s´erie, ou seja, os res´ıduos comportam-se como um processo autorregressivo AR(1). Na parte fixa consideram-se dois casos: sem intera¸c˜ao entre as covari´aveis em estudo e a vari´avel Tempo e com intera¸c˜ao entre as covari´aveis em estudo e a vari´avel Tempo. No segundo caso, se a intera¸c˜ao for significativa, ent˜ao o efeito da covari´avel no valor esperado da concentra¸c˜ao de OD depende do tempo em estudo. Os efeitos sazonais usualmente variam de forma suave e cont´ınua, pelo que faz sentido a considera¸c˜ao de fun¸c˜oes de suaviza¸c˜ao, recorrendo a fun¸c˜oes seno e cosseno para descrever as oscila¸c˜oes observadas ao longo do tempo, ajustando um modelo sazonal harm´onico. Os Modelos Sazonais Harm´onicos consideram um somat´orio, que descreve a varia¸c˜ao sazonal dos dados, ao longo de cada per´ıodo (s = 12, para dados mensais): •cosx´e o cosseno de per´ıodo x, onde xassume os valores 12; 6; 4; 3; 2,4 e 2. Por exemplo, cos6´e a codifica¸c˜ao para o cosseno de per´ıodo 6; •senx´e o seno de per´ıodo x, onde xassume os valores 12; 6; 4; 3 e 2,4. De notar que nos Modelos Sazonais Harm´onicos (com erros correlacionados), quando uma curva cosseno com um certo per´ıodo ´e significativa, esta deve ser inclu´ıda juntamente com o correspondente seno do mesmo per´ıodo e vice-versa. 104 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Considerando que irepresenta a esta¸c˜ao de amostragem em estudo; Tempoit designa o tempo desde mar¸co de 2002 at´e fevereiro de 2013, em meses; CBO5it refere-se `a vari´avel de Carˆencia Bioqu´ımica de Oxig´enio, em mg/l, na esta¸c˜ao de amostragem ino instante t;Clorofilait retrata o valor dos n´ıveis de Clorofila, em ug/l, na esta¸c˜ao de amostragem ino instante t;ALBit designa se a esta¸c˜ao de amostragem i´e ou n˜ao uma albufeira; pHit representa o valor de pH, na esta¸c˜ao de amostragem ino instante t;Temperaturait designa o valor de temperatura da amostra, em ◦C, na esta¸c˜ao de amostragem ino instante t;CHit,CHi(t−1) eCHi(t−2) referem-se ao fator hidrometeorol´ogico, na esta¸c˜ao de amostragem i, nos instantes t,t−1 e t−2, respetivamente; s´e o per´ıodo de sazonalidade (s= 12). Os quatro modelos analisados s˜ao: (i) O Modelo 1 consiste no modelo com um efeito aleat´orio no termo constante (associado ao parˆametro β0), ODit =β0+β1Tempoit +β2CBO5it +β3Clorofilait +β4pHit +β5ALBit +β6Temperaturait +β7CHit +β8CHi(t−1) +β9CHi(t−2) + s/2 X k=1 "α1icos2πkTempoit s+α2isen2πkTempoit s#+u1it +it i=1, . . . , 36, t = 1, . . . , 132 (5.1) em que it ´e o erro aleat´orio tal que it =φ1it−1+ait, com ait ∼N(0, σ2); u1i representa o efeito aleat´orio e u1it ∼N(0, d2 11); (ii) O Modelo 2 consiste no modelo com efeito aleat´orio no termo constante e na vari´avel tempo, isto ´e, ODit =β0+β1Tempoit +β2CBO5it +β3Clorofilait +β4pHit +β5ALBit +β6Temperaturait +β7CHit +β8CHi(t−1) +β9CHi(t−2) + s/2 X k=1 "α1icos2πkTempoit s+α2isen2πkTempoit s# +u1it +u2itTempoit +it i=1, . . . , 36, t = 1, . . . , 132 (5.2) em que it ´e o erro aleat´orio tal que it =φ1it−1+ait, com ait ∼N(0, σ2); u1it,u2it representam os efeitos aleat´orios e 105 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais u1i u2i!∼N(0,D),D="d2 11 d12 d21 d2 22#; (iii) O Modelo 3 consiste no modelo com um efeito aleat´orio no termo constante e um termo de intera¸c˜ao ODit =β0+β1Tempoit +β2CBO5it +β3Clorofilait +β4pHit +β5ALBit +β6Temperaturait +β7CHit +β8CHi(t−1) +β9CHi(t−2) +β10CBO5it :Tempoit +β11Clorofilait :Tempoit +β12pHit :Tempoit +β10ALBit :Tempoit +β13Temperaturait :Tempoit +β14CHit :Tempoit +β15CHi(t−1) :Tempoit +β16CHi(t−2) :Tempoit + s/2 X k=1 "α1icos2πkTempoit s+α2isen2πkTempoit s# +u1it +it i=1, . . . , 36, t = 1, . . . , 132 (5.3) em que it ´e o erro aleat´orio tal que it =φ1it−1+ait, com ait ∼N(0, σ2); u1it representa o efeito aleat´orio e u1it ∼N(0, d2 11); (iv) O Modelo 4 consiste no modelo com um efeito aleat´orio no termo constante e na vari´avel tempo e um termo de intera¸c˜ao ODit =β0+β1Tempoit +β2CBO5it +β3Clorofilait +β4pHit +β5ALBit +β6Temperaturait +β7CHit +β8CHi(t−1) +β9CHi(t−2) +β10CBO5it :Tempoit +β11Clorofilait :Tempoit +β12pHit :Tempoit +β10ALBit :Tempoit +β13Temperaturait :Tempoit +β14CHit :Tempoit +β15CHi(t−1) :Tempoit +β16CHi(t−2) :Tempoit + s/2 X k=1 "α1icos2πkTempoit s+α2isen2πkTempoit s# +u1it +u2itTempoit +it i=1, . . . , 36, t = 1, . . . , 132 (5.4) 106 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais em que it ´e o erro aleat´orio tal que it =φ1it−1+ait, com ait ∼N(0, σ2); u1it,u2it representam os efeitos aleat´orios e u1i u2i!∼N(0,D),D="d2 11 d12 d21 d2 22#. Com base na formula¸c˜ao dos modelos estimaram-se os modelos completos. As estimativas dos efeitos fixos, os erros padr˜ao, os valores da estat´ıstica de teste e os valores de prova s˜ao apresentados nas Tabelas 5.8 a 5.11. O M´etodo de M´axima Verosimilhan¸ca Restrito (REML) foi utilizado para o ajustamento dos diferentes modelos. Os modelos finais s˜ao apresentados na Tabela 5.13. Com base no Teste da Raz˜ao de Verosimilhan¸ca apresentado na Tabela 5.12, onde se comparam os modelos encaixados (Modelo 1 vs Modelo 3 e Modelo 2 vs Modelo 4), pode-se verificar que os valores de prova s˜ao inferiores a 5 %, rejeitando-se a hip´otese nula. A compara¸c˜ao entre o Modelo 3 e o Modelo 4 ´e feita com base nos crit´erios de informa¸c˜ao AIC eBIC, uma vez que estes n˜ao s˜ao aninhados. O Modelo 4 ´e o modelo que apresenta maior valor de AIC eBIC (Tabela 5.12) e, por essa raz˜ao, o Modelo 3 ´e o modelo selecionado. Tabela 5.8: Estimativas dos coeficientes da parte fixa, do Modelo Completo 1. Vari´avel Estimativa Erro Padr˜ao Estat´ıstica t Valor de Prova Constante 7,7822 0,3145 24,7482 <0,0001 Tempoit 0,0018 0,0005 3,2851 0,0010 CBO5 0,0607 0,0181 3,3605 0,0008 Clorofila 0,0007 0,0018 0,3586 0,7200 pH 0,2732 0,0408 6,6895 <0,0001 Temperatura -0,0591 0,0061 -9,6186 <0,0001 ALB -0,1393 0,1164 -1,1968 0,2399 CHt-0,0002 0,0005 -0,3102 0,7565 CHt−10,0018 0,0005 3,3020 0,0010 CHt−20,0007 0,0005 1,2282 0,2195 cos12 0,8606 0,0553 15,5681 <0,0001 cos60,0997 0,0262 3,8072 <0,0001 cos4-0,0373 0,0209 -1,7866 0,0741 cos30,0414 0,0189 2,1899 0,0286 cos2,40,0053 0,0174 0,3028 0,7621 cos20,0203 0,0132 1,5305 0,1260 sen12 -0,0434 0,0352 -1,2356 0,2167 sen6-0,0273 0,0264 -1,0348 0,3009 sen4-0,0415 0,0204 -2,0374 0,0417 sen30,0119 0,0191 0,6229 0,5334 sen2,4-0,0175 0,0183 -0,9566 0,3388 107 Cap´ıtulo 5. Aplica¸c˜ao aos Dados Ambientais Tabela 5.9: Estimativas dos coeficientes da parte fixa, do Modelo Completo 2. Vari´avel Estimativa Erro Padr˜ao Estat´ıstica t Valor de Prova Constante 7,7822 0,3145 24,7482 <0,0001 Tempoit 0,0018 0,0005 3,2851 0,0010 CBO5 0,0607 0,0181 3,3605 0,0008 Clorofila 0,0007 0,0018 0,3586 0,7200 pH 0,2732 0,0408 6,6895 <0,0001 Temperatura -0,0591 0,0061 -9,6186 <0,0001 Albufeira -0,1393 0,1164 -1,1968 0,2399 CHt-0,0002 0,0005 -0,3102 0,7565 CHt−10,0018 0,0005 3,3020 0,0010 CHt−20,0007 0,0005 1,2282 0,2195 cos12 0,8606 0,0553 15,5681 <0,0001 cos60,0997 0,0262 3,8072 0,0001 cos4-0,0373 0,0209 -1,7866 0,0741 cos30,0414 0,0189 2,1899 0,0286 cos2,40,0053 0,0174 0,3028 0,7621 cos20,0203 0,0132 1,5305 0,1260 sen12 -0,0434 0,0352 -1,2356 0,2167 sen6-0,0273 0,0264 -1,0348 0,3009 sen4-0,0415 0,0204 -2,0374 0,0417 sen30,0119 0,0191 0,6229 0,5334 sen2,4-0,0175 0,0183 -0,9566 0,3388 Tabela 5.10: Estimativas dos coeficientes da parte fixa, do Modelo Completo 3. Vari´avel Estimativa Erro Padr˜ao Estat´ıstica t Valor de Prova Constante 7,9377 0,4665 17,0165 <0,0001 Tempoit 0,0007 0,0065 0,1112 0,9115 CBO5 0,0222 0,0334 0,6669 0,5049 Clorofila 0,0048 0,0033 1,4753 0,1402 pH 0,3082 0,0631 4,8802 <0,0001 Temperatura -0,0781 0,0085 -9,1897 <0,0001 ALB -0,2068 0,1305 -1,5840 0,1227 CHt-0,0002 0,0008 -0,1981 0,8430 CHt−10,0015 0,0009 1,6558 0,0979 CHt−20,0002 0,0009 0,2239 0,8229 cos12 0,8652 0,0557 15,5285 <0,0001 cos60,0950 0,0265 3,5863 0,0003 cos4-0,0346 0,0209 -1,6535 0,0983 cos30,0412 0,0190 2,1677 0,0303 cos2,40,0059 0,0177 0,3338 0,7386 cos20,0201 0,0135 1,4951 0,1350 sen12 -0,0441 0,0361 -1,2210 0,2222 sen6-0,0382 0,0266 -1,4342 0,1516 sen4-0,0451 0,0204 -2,2051 0,0275 sen30,0142 0,0192 0,7382 0,4605 sen2,4 -0,0186 0,0184 -1,0123 0,3115 CBO5 : Tempoit 0,0007 0,0006 1,3007 0,1935 Clorofila :Tempoit -0,0001 0,0001 -1,4800 0,1390 pH -0,0008 0,0009 -0,8707 0,3840 Temperatura :Tempoit 0,0003 0,0001 3,1540 0,0016 ALB :Tempoit 0,0011 0,0011 1,0459 0,2957 CHt:Tempoit 0,0000 0,0000 0,0790 0,9371 CHt−1:Tempoit 0,0000 0,0000 0,0642 0,9488 CHt−2:Tempoit 0,0000 0,0000 0,3525 0,7245 108 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Tabela 5.11: Estimativas dos coeficientes da parte fixa, do Modelo Completo 4. Vari´avel Estimativa Erro Padr˜ao Estat´ıstica t Valor de Prova Constante 7,9281 0,4733 16,7521 <0,0001 Tempoit 0,0008 0,0067 0,1194 0,9049 CBO5 0,0218 0,0334 0,6521 0,5144 Clorofila 0,0048 0,0033 1,4708 0,1415 pH 0,3112 0,0641 4,8562 <0,0001 Temperatura -0,0787 0,0085 -9,2484 <0,0001 ALB -0,2080 0,1337 -1,5559 0,1293 CHt-0,0002 0,0008 -0,2041 0,8383 CHt−10,0015 0,0009 1,6470 0,0997 CHt−20,0002 0,0009 0,2132 0,8312 cos12 0,8635 0,0557 15,4974 <0,0001 cos60,0952 0,0265 3,5980 0,0003 cos4-0,0345 0,0209 -1,6497 0,0991 cos30,0410 0,0190 2,1569 0,0311 cos2,40,0059 0,0177 0,3345 0,7380 cos20,0202 0,0135 1,5023 0,1331 sen12 -0,0434 0,0361 -1,2021 0,2294 sen6-0,0384 0,0266 -1,4412 0,1496 sen4-0,0452 0,0204 -2,2136 0,0269 sen30,0143 0,0192 0,7420 0,4582 sen2,4-0,0187 0,0184 -1,0162 0,3096 CBO5 : Tempoit 0,0007 0,0006 1,3183 0,1875 Clorofila :Tempoit -0,0001 0,0001 -1,4727 0,1410 pH :Tempoit -0,0008 0,0009 -0,8859 0,3757 Temperatura :Tempoit 0,0003 0,0001 3,2109 0,0013 ALB :Tempoit 0,0011 0,0011 1,0168 0,3093 CHt:Tempoit 0,0000 0,0000 0,0889 0,9292 CHt−1:Tempoit 0,0000 0,0000 0,0768 0,9388 CHt−2:Tempoit 0,0000 0,0000 0,3598 0,7190 Tabela 5.12: Teste da Raz˜ao de Verosimilhan¸ca aplicado aos modelos em an´alise. Modelo AIC BIC loglik Teste Valor de Prova Modelo 1 6492,439 6592,938 -3229,220 Modelo 3 6496,439 6608,761 -3229,220 1 vs 3 0,0178 Modelo 2 6497,727 6604,130 -3230,863 Modelo 4 6501,467 6619,694 -3230,734 2 vs 4<0,0001 109 Bibliografia [12] Burnham, K. P., & Anderson, D. R. (2004). Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods & Research, 33(2), 261 –304. [13] Burnham, K. P., & Anderson, D.R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. New York: Springer-Verlag, 2a ed. [14] Cabecinha, E., Cortes, R., Pardal, M. ˆ A., & Cabral, J. A. (2009). A Stochastic Dynamic Methodology (StDM) for reservoirs water quality management: Validation of a multi-scale approach in a south european basin (Douro, Portugal). Ecological Indicators, 9(2), 329 –345. [15] Caiado, J. (2016). M´etodos de Previs˜ao em Gest˜ao - Com Aplica¸c˜oes em Excel. Lisboa: Edi¸c˜oes S´ılabo, 2a ed. [16] Chatfield, C. (2000). Time-Series Forecasting. Chapman & Hall/CRC, 1a ed. [17] Chatfield, C. (2004). The Analysis of Time Series: An Introduction. Chapman & Hall/CRC, 5a ed. [18] Chow, V. T., Maidment, D. R., & Mays, L. W. (1988). Applied Hidrology. McGrawHill Series in Water Resources and Environmental Engineering, McGraw-Hill, 1a ed. [19] Costa, M., & Gon¸calves, A. M. (2011). Clustering and forecasting of dissolved oxygen concentration on a river basin. Stochastic Environmental Research and Risk Assessment, 25(2), 151 –163. [20] Costa, M., & Gon¸calves, A. M. (2012). Combining statistical methodologies in water quality monitoring in a hydrological Basin-Space and time approaches. Water Quality Monitoring and Assessment, 121 –142. [21] Cowpertwait, P., & Metcalfe, A. (2009). Introductory Time Series with R. New York: Springer, 1a ed. [22] Cox, D. R., & Hinkley, D. V. (1974). Theoretical statistics. New York: Chapman&Hall/CRC, 1a ed. [23] Cressie, N. A. C. (1993). Statistics for Spatial Data. New York: Wiley, 2a ed. [24] Cressie, N., & Huang, H. C. (1999). Classes of nonseparable, spatio-temporal stationary covariance functions. Journal of the American Statistical Association, 94(448), 1330 –1339. 116 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais [25] Cressie, N., & Wikle, C. K. (2015). Statistics for spatio-temporal data. John Wiley & Sons. [26] Cressie, N., Zammit-Mangion, A., & Wikle, C. K. (2019). Spatio-Temporal Statistics with R. Chapman & Hall/CRC. [27] D0Odorico, P., Carr, J., Dalin, C., DellAngelo, J., Konar, M., Laio, F., & Tuninetti, M. (2019). Global virtual water trade and the hydrological cycle: patterns, drivers, and socio-environmental impacts. Environmental Research Letters, 14(5), 053001. [28] Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1), 1 –22. [29] Dickey, D., & Fuller, W. (1979). Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74, 427 –431. [30] Diggle, P. J., & Giorgi, E. (2019). Model-based Geostatistics for Global Public Health: Methods and Applications. Chapman & Hall/CRC Interdisciplinary Statistics, 1a ed. [31] Diggle, P. J., Heagerty, P., Liang, K. Y., Heagerty, P. J., & Zeger, S. (2002). Analysis of longitudinal data. Oxford University Press. [32] Diggle, P. J., Liang, K. Y., & Zeger, S. L. (1994). Analysis of Longitudinal Data, Oxford: Oxford University Press. [33] Enders, W. (2015). Applied Econometric Time Series. Alabama: John Wiley and Sons, 4a ed. [34] Fausto, M., Carneiro, M., Antunes, C., Pinto, J., & Colosimo, E. (2008). O modelo de regress˜ao linear misto para dados longitudinais: uma aplica¸c˜ao na an´alise de dados antropom´etricos desbalanceados. Cadernos de Sa´ude P´ublica, 24. [35] Gon¸calves, A. M., & Alpuim, T. (2011). Water quality monitoring using cluster analysis and linear models. Environmetrics, 22(8), 933 –945. [36] Gon¸calves, A. M., & Costa, M. (2013). Predicting seasonal and hydro-meteorological impact in environmental variables modelling via Kalman filtering. Stochastic Environmental Research and Risk Assessment, 27(5), 1021 –1038. [37] Goovaerts, P. (1997). Geostatistics for natural resources evaluation. Oxford University Press on Demand. 117 Bibliografia [38] Gr¨ aler, B. (2014). Modelling skewed spatial random fields through the spatial vine copula. Spatial Statistics, 10, 87 –102. [39] Harville, D. (1976). Extension of the Gauss-Markov theorem to include the estimation of random effects. The Annals of Statistics, 4(2), 384 –395. [40] Harville, D. A. (1974). Bayesian inference for variance components using only error contrasts. Biometrika, 61(2), 383 –385. [41] Harville, D. A. (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American statistical association, 72(358), 320 –338. [42] Ho, L., Pham, D., Van Echelpoel, W., Muchene, L., Shkedy, Z., Alvarado, A., Espinoza-Palacios, J., Arevalo-Durazno, M., Thas, O. & Goethals, P. (2018). A closer look on spatiotemporal variations of dissolved oxygen in waste stabilization ponds using mixed models. Water, 10(2), 201. [43] Isaacs, E. H., & Srivastava, R. M. (1989), An Introduction to Applied Geostatistics, Oxford University Press. [44] Jebb, A., Tay, L., Wang, W., & Huang, Q. (2015). Time series analysis for psychological research: examining and forecasting change. Frontiers in Psychology, 6, 727. [45] Jowett, G. H. (1955). Sampling properties of local statistics in stationary stochastic series. Biometrika, 42(1/2), 160 –169. [46] Kass, R. E., Caffo, B. S., Davidian, M., Meng, X. L., Yu, B., Reid, N. (2016). Ten simple rules for effective statistical practice. PLOS Computational Biology, 12(6). [47] Kirchg¨ assner, G., & Wolters, J. (2008). Introduction to modern time series analysis. Springer-Verlag. [48] Kolmogorov, A. N. (1941). Interpolation and Extrapolation of Stationary Sequences. Izvestiya the Academy of Sciences of the USSR Series Mathematic, 5, 3-14. [49] Krige, D. G. (1951). A statistical approach to some basic mine valuation problems on the Witwatersrand. Journal of the Southern African Institute of Mining and Metallurgy, 52(6), 119 –139. [50] Kwiatkowski, D., Phillips, P., Schmidt, P., & Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root. Journal of Econometrics, 54, 159 –178. 118 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais [51] Kyriakidis, P. C., Kim, J., & Miller, N. L. (2001). Geostatistical mapping of precipitation from rain gauge data using atmospheric and terrain characteristics. Journal of Applied Meteorology, 40(11), 1855 –1877. [52] Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963 –974. [53] Lefˆevre, C. (1997). ”Kriging”e ”Cokringing”e aplica¸c˜oes ao Radar Meteorol´ogico. Master Thesis, University of Lisbon, Lisbon. [54] Makridakis, S., Wheelwright, S., & Hyndman, R. (1998). Forecasting: Methods and Applications. New York: John Wiley and Sons, 3a ed. [55] Matern, B. (1960). Spatial Variation. Meddelanden Statens fr˚an Skogsforskningsinstitut Stockholm, 49 (5), 1 –144. [2nd edition (1986). Spatial variation. Lecture Notes in Statistics, No. 36, Springer, New York, 2a ed.] [56] Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58(8), 1246 –1266. [57] Moshogianis, A. (2015). A Statistical Model for the Prediction of Dissolved Oxygen Dynamics and the Potential for Hypoxia in the Mississippi Sound and Bight. Master Thesis. The University of Southern Mississippi, Mississippi. [58] McCulloch, C. E., & Searle, S. R. (2001). Generalized, Linear, and Mixed Models. Wiley Series in Probability and Statistics. [59] McCulloch, P., Nelder, J.A. (1989) Generalized Linear Models. London: Chapman and Hall, 2a Ed. [60] Mercer, W. B., & Hall, A. D. (1911). The experimental error of field trials. The Journal of Agricultural Science, 4(2), 107 –132. [61] Min, S. K., Zhang, X., Zwiers, F. W., & Hegerl, G. C. (2011). Human contribution to more-intense precipitation extremes. Nature, 470(7334), 378. [62] Molenberghs, G., & Verbeke, G. (2005). Model for Discrete Longitudinal Data. New York: Springer. [63] Murteira, B., Muller, D., & Turkman, K. (1993). An´alise de Sucess˜oes Cronol´ogicas. Lisboa: McGraw-Hill. [64] Nelder, J. A., & Wedderburn, R. W. (1972). Generalized linear models. Journal of the Royal Statistical Society: Series A (General), 135(3), 370 –384. 119 Bibliografia [65] Olea, R. A. (2012). Geostatistics for engineers and earth scientists. Springer Science & Business Media. [66] Patterson, H. D., & Thompson, R. (1971). Recovery of inter-block information when block sizes are unequal. Biometrika, 58(3), 545 –554. [67] Pebesma, E. J. (2004). Multivariable geostatistics in S: the gstat package. Computers & Geosciences, 30(7), 683 –691. [68] Pebesma, E. J. (2012). spacetime: Spatio-temporal data in R. Journal of Statistical Software, 51(7), 1 –30. [69] Pebesma, E., & Heuvelink, G. (2016). Spatio-temporal interpolation using gstat.RFID Journal, 8(1), 204 –218. [70] Phillips, P., & Perron, P. (1988). Testing for unit roots in time series regression. Biometrika, 75, 335 –346. [71] Pinheiro, J. C. (1994). Topics in mixed effects models. Ph. D. Thesis, University of Wisconsin, Madison. [72] Pinheiro, J. C., & Bates, D. M. (1995). Approximations to the log-likelihood function in the nonlinear mixed-effects model. Journal of computational and Graphical Statistics, 4(1), 12 –35. [73] Pinheiro, J. C., & Bates, D. M. (2000). Linear mixed-effects models: basic concepts and examples. Mixed-effects models in S and S-Plus, 3 –56. [74] R Core Team (2017). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.Rproject.org/ [75] Ripley, B. D. (1981). Spatial Statistics. New York: Wiley & Sons, 252 Rouhani, S., & Wackernagel, H. (1990). Multivariate geostatistical approach to spacetime data analysis. Water Resources Research, 26(4), 585 –591. [76] Said, S., & Dickey, D. (1984). Testing for unit roots in autoregressive moving-average models with unknown order.Biometrika, 71, 599 –607. [77] Schielzeth, H., Nakagawa, S. (2013). Nested by design: model fitting and interpretation in a mixed model era. Methods in Ecology Evolution, 4(1), 14 –24. [78] Schwarz, G. (1978). Estimating the dimension of a model. Annals of Statistics, 6, 461 –464. 120 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais [79] Searle, S. R., Casella, G., & McCulloch, C. E. (1992). Variance components. Hoboken: John Wiley & Sons. [80] Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3/4), 591 –611. [81] Silva-Santos, P., Pardal, M. ˆ A., Lopes, R. J., M´urias, T., & Cabral, J. A. (2008). Testing the Stochastic Dynamic Methodology (StDM) as a management tool in a shallow temperate estuary of south Europe (Mondego, Portugal). Ecological Modelling, 210(4), 377 –402. [82] Thiessen, A. (1911). Precipitation Averages for Large Areas. Monthly Weather Review, 39(7), 1082 –1084. [83] Thisted, R. A. (1988). Elements of Statistical Computing. London: Chapman & Hall. [84] Trigo, R. M., PozoV´azquez, D., Osborn, T. J., CastroD´ıez, Y., G´amizFortis, S., & EstebanParra, M. J. (2004). North Atlantic Oscillation influence on precipitation, river flow and water resources in the Iberian Peninsula. International Journal of Climatology: A Journal of the Royal Meteorological Society, 24(8), 925 –944. [85] Turkman, M. A. A., & Silva, G. L. (2000). Modelos Lineares Generalizados da teoria a pr´atica. Lisboa: VIII Congresso Anual da Sociedade Portuguesa de Estat´ıstica. [86] van Dijk, G. M., van Liere, L., Admiraal, W., Bannink, B. A., & Cappon, J. J. (1994). Present state of the water quality of european rivers and implications for management. Science of the Total Environment, 145(1-2), 187 –195. [87] Verbeke, G., & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data, New York: Springer. [88] Voronoi, G. (1908). Nouvelles applications des parametres continus a la theorie des formes quadratiques. Deuxieme Memoire: Recherche sur les paralleloedres primitifs. Journal Reine Angew. Math., 134, 198 –287. [89] Wu, L. (2009). Mixed effects models for complex data. Chapman & Hall/CRC. [90] Yaglom, A. M. (1962), An introduction to the theory of stationary random functions: Englewood Cliffs, Pretice-Hall, Inc. [91] Zuur A. F., Ieno E. N., Walker N. J., Saveliev A. A., Smith G.M. (2009). Mixed Effects Models and Extensions in Ecology with R. New York: Springer. 121 Bibliografia 122 Apˆendice A Geoestat´ıstica A.1 Representa¸c˜oes Gr´aficas da Precipita¸c˜ao em cada Esta¸c˜ao de Amostragem Figura A.1: Diagrama em caixa de bigodes das s´erie de Precipita¸c˜ao, nas 18 esta¸c˜oes de amostragem, no per´ıodo observado. 123 Apˆendice A. Geoestat´ıstica Figura A.2: Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. 124 Modela¸c˜ao Estat´ıstica na An´alise em Processos Ambientais Figura A.3: Representa¸c˜oes gr´aficas das s´eries temporais da precipita¸c˜ao e dos respetivos diagramas em caixas de bigodes e histogramas, no per´ıodo observado. 125