#----------------------------------------------------------------------------------
#              PRÁTICA 5 -  ANÁLISE DE SÉRIES TEMPORAIS - MINAS GERAIS
#----------------------------------------------------------------------------------

library(RcmdrMisc)
library(epiDisplay)
library(tidyverse)

# Importe o banco de dados "Complica.xls" e dê o nome de dados

dados <- readXL("C:/A DATA/PSP5103/2021/Dados/Complica com extras.xls", 
                rownames=FALSE, header=TRUE, na="", sheet="Data", stringsAsFactors=FALSE)

names(dados)
dados

plot(Mg~Ano, xlab="Tempo (anos)", ylab="CM/100.000 h", col="blue",pch=16, data=dados)

dados$Ano
dados$Ano2 <- dados$Ano^2
dados$Ano3 <- dados$Ano^3

cor(dados[,c("Ano","Ano2","Ano3")], use="complete")

scatterplotMatrix(~Ano+Ano2+Ano3, regLine=TRUE, smooth=FALSE, 
                  diagonal=list(method="none"), data=dados)

# Problema: multicolinearidade quase perfeita
# Não é o caso de duas variáveis diferentes que tem colinearidade - menos grave.
# porque estou o ajustando um polinomio. 

dados$Ano0 <- dados$Ano - 1998              # Ano começando no zero
dados$Ano0
dados$Ano02 <- dados$Ano0^2
dados$Ano03 <- dados$Ano0^3

cor(dados[,c("Ano0","Ano02","Ano03")], use="complete")

scatterplotMatrix(~Ano0+Ano02+Ano03, regLine=TRUE, smooth=FALSE, 
                  diagonal=list(method="none"), data=dados)


mean(dados$Ano)

dados$Anoc <- dados$Ano - mean(dados$Ano)   # Ano centralizado
dados$Anoc

mean(dados$Anoc)

dados$Anoc2 <- dados$Anoc^2
dados$Anoc3 <- dados$Anoc^3

cor(dados[,c("Anoc","Anoc2","Anoc3")], use="complete")

scatterplotMatrix(~Anoc+Anoc2+Anoc3, regLine=TRUE, smooth=FALSE, 
                  diagonal=list(method="none"), data=dados)


# Minas 
#---------

# Diagramas de dispersão

par(mfrow=c(1,2))
plot(Mg~Ano, xlab="Tempo (anos)", ylab="CM/100.000 h", col="blue",pch=16, data=dados)
plot(Mg.mm~Ano, xlab="Tempo (anos)", ylab="Média móvel de CM/100.000 h",pch=16 ,col="blue",data=dados)

# Verificando a suposição de Normalidade

normalityTest(~Mg, test="shapiro.test", data=dados)
normalityTest(~Mg, test="lillie.test", data=dados)

normalityTest(~Mg.mm, test="shapiro.test", data=dados)
normalityTest(~Mg.mm, test="lillie.test", data=dados)


par(mfrow=c(1,2))
qqPlot(dados$Mg, dist="norm", pch=16, id=list(method="y", n=2, labels=rownames(dados)))
qqPlot(dados$Mg.mm, dist="norm", pch=16, id=list(method="y", n=2, labels=rownames(dados)))

# Modelos com diferentes parametrizações para o Ano

mod.Mg.Ano <- lm(Mg~Ano, data=dados)
summary(mod.Mg.Ano)

mod.Mg.Ano0 <- lm(Mg~Ano0, data=dados)     # Utilizando o ano começando no zero
summary(mod.Mg.Ano0)                       # Agora, bo é o CM médio no ano zero (1988)

mod.Mg.Anoc <- lm(Mg~Anoc, data=dados)     # Utilizando o ano centralizado
summary(mod.Mg.Anco)                       # Agora, bo é o CM médio do período
mean(dados$Mg)

par(mfrow=c(1,3))

plot(Mg~Ano, xlab="Tempo (anos)", ylab="CM/100.000 h", col="blue",data=dados)
abline(mod.Mg.Ano)
plot(Mg~Ano0, xlab="Tempo (anos)", ylab="CM/100.000 h", col="blue",data=dados)
abline(mod.Mg.Ano0)
plot(Mg~Anoc, xlab="Tempo (anos)", ylab="CM/100.000 h", col="blue",data=dados)
abline(mod.Mg.Anoc)

# Vamos trabalhar com o ano centralizado

# Comparando os modelos com os dados originais e com a média móvel

mod.Mgmm.Anoc <- lm(Mg.mm~Anoc, data=dados)  
summary(mod.Mgmm.Anoc)   

par(mfrow=c(1,2))

plot(Mg~Anoc, xlab="Tempo (anos)", ylab="CM/100.000 h", pch=16, col="blue",data=dados)
abline(mod.Mg.Anoc)
plot(Mg.mm~Anoc, xlab="Tempo (anos)", ylab="Média móvel de CM/100.000 h", pch=16, col="blue",data=dados)
abline(mod.Mgmm.Anoc)

# Vamos trabalhar com as médias móveis


# Modelos polinomiais
#----------------------------------------------------------------------------------------

# Nos modelos com médias móveis, para evitar problemas nos gráficos, 
# vamos trabalhar com os dados sem missings, eliminando a primeira e última linha.

dados1 <- na.omit(dados)
View(dados1)

Mg.linear <- lm(Mg.mm ~ Anoc, data=dados1)  
summary(Mg.linear) 
anova(Mg.linear)

Mg.quadratico <- lm(Mg.mm ~ Anoc + Anoc2, data=dados1)  
summary(Mg.quadratico) 
anova(Mg.quadratico)
anova(Mg.linear, Mg.quadratico, test="F")

Mg.cubico <- lm(Mg.mm~Anoc + Anoc2 + Anoc3, data=dados1)  
summary(Mg.cubico) 
anova(Mg.cubico)
anova(Mg.quadratico, Mg.cubico, test="F")

# Quadrático b1-velocidade, b2-aceleração ?
# Cúbico - aceleração não é constante ?


# Utilizando a função poly

Mg.quadratico <- lm(Mg.mm ~ poly (Anoc, 2, raw = T), data=dados1)  # Modelo polinomial (quadrático)
summary(Mg.quadratico) 

Mg.cubico <- lm(Mg.mm~ poly (Anoc, 3, raw = T), data=dados1)  # Modelo polinomial (quadrático)
summary(Mg.cubico) 


# Modelo linear é melhor!

par(mfrow=c(1,3))

plot(Mg.mm~Anoc, xlab="Tempo (anos)", ylab="Média móvel de CM/100.000 h", pch=16, col="blue",data=dados1)
lines(fitted(Mg.linear)~dados1$Anoc, col="red", lwd=1.5)
plot(Mg.mm~Anoc, xlab="Tempo (anos)", ylab="Média móvel de CM/100.000 h", pch=16, col="blue",data=dados1)
lines(fitted(Mg.quadratico)~dados1$Anoc, col="red", lwd=1.5)
plot(Mg.mm~Anoc, xlab="Tempo (anos)", ylab="Média móvel de CM/100.000 h", pch=16, col="blue",data=dados1)
lines(fitted(Mg.cubico)~dados1$Anoc, col="red", lwd=1.5)


#Análise de resíduos

par(mfrow=c(1,4))

plot(Mg.mm~Anoc, xlab="Tempo (anos)", ylab="Média móvel de CM/100.000 h", pch=16, col="blue",data=dados1)
lines(fitted(Mg.linear)~dados1$Anoc, col="red", lwd=1.5)

plot(rstandard(Mg.linear) ~ dados1$Anoc, ylab="Resíduo Padronizado", xlab="Tempo (anos)")
abline(h=0, col="blue"); abline(h=-1.96, col="blue", lty=2); abline(h= 1.96, col="blue", lty=2)

plot(rstandard(Mg.linear) ~ fitted(Mg.linear), ylab="Resíduo Padronizado", xlab="Valores ajustados")
abline(h=0, col="blue"); abline(h=-1.96, col="blue", lty=2); abline(h= 1.96, col="blue", lty=2)

qqPlot(rstandard(Mg.linear), dist="norm", pch=16, id=list(method="y", n=2, labels=rownames(dados1)))