Estrutura do Trabalho

O projeto abaixo é dividido em 2 Etapas:

Etapa 1: Análise de Dados do brasil com base no SUS

Foram coletadas e tratadas diversas bases de dados relacionadas a suícidio através do DATASUS, esta fase foi importante pois as informações disponiveis para download estavam mal organizadas, com linhas de totais e subtotais em todas as colunas e outras informações que comprometeriam as análises feitas a posteriore com os dados.

Etapa 2: Séries Temporais e previsão

Foi feita uma análise temporal com os dados coletados no DATASUS com objetivo de encontrar tendências e padrões ao longo dos anos. Como os dados disponibilizados pelo DATASUS tem um gap de 2 anos foi gerado previsões para os anos faltantes, a medida que no futuro se deseje avaliar se os resultados possuem uma boa aproximação com a realidade.

Carregando os pacotes que seram utilizados para a análise

Bibliotecas

#install.packages("dplyr")
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#install.packages("tidyr")
library(tidyr)
#install.packages("readr")
library(readr)
#install.packages("readr")
#install.packages("visdat")
#install.packages("dplyr")
#install.packages("plotly")
#install.packages("ggplot2")
#install.packages("corrplot")
#install.packages("knitr")
#install.packages("webshot")
library(webshot)
#webshot::install_phantomjs()
library(knitr)
library(readr)
library(visdat)
library(tidyr)
library(dplyr)
library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(corrplot)

## corrplot 0.84 loaded

library(ggplot2)
library(plotly)

setwd("~/Portfolio/TCC/Dados")

#Parte 2: Análise de Dados do brasil com base no SUS

Pré Processamento da da base de dados

Para segunda parte do projeto será necessário a utilização de dados no portal do SUS (DATASUS) do SIM (Sistema de informação de mortalidade) que possuem as informações oficiais do Brasil em relação as mortes por suicídio . Quando analisei as informações no site percebi que para uma pessoa coletar todas as informações precisaria de pelo menos uns 20 arquivos csv’s. Pensando nisso essa parte do projeto foi criada para coletar todos os datasets separados e agrupar suas informações. Essa etapa se faz essêncial para que seja possível analisar as informações do Brasil referêntes a suicídio.

Couteudo = Óbitos p/ocrorrencia Grupo CID10= Lesões autoprovocadas intencionalmente Grande Grupo CID10= Lesões autoprovocadas voluntariamente

Fonte: http://tabnet.datasus.gov.br/cgi/deftohtm.exe?sim/cnv/ext10uf.def.

Para que o leitor tenha uma ideia de como os dados estão organizado, abaixo tem uma amostra de um dos arquivos baixados que trata de informações de suicídio por Estado.

df=read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SucidioPorRegiao.csv", sep = ";")
head(df)

Note que temos colunas como Ano.mês.do.Óbito que está mal formatada e mistura linhas com total no ano com o total para cada mês. Além disso para cada estado foi criado uma coluna o que dificulta as possíveis análises feitas com o tema, embora não seja mostrado no slice acima a base também possuiam textos não relacionados aos dados no início e no final do arquivo que foram retirados.

Para fazer este tratamento foi criado a função abaixo, ela resolverá todos os problemas relatados anteriormente.

trataDados <- function(df,x){
#Reunindo os estados em uma coluna
df=gather(df,UF,Suicidios,2:x)
names(df)[1]= "Ano.Mes"

#Retirando textos que estão no meio da tabela
df=df[!is.na(df$Suicidios),]
df=df[!is.null(df$Suicidios)]
df=df[!(df$Suicidios==""),]

#Retirando valores de totais anuais
df=filter(df,!df$Ano.Mes %in% c(1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,
                               2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019))

#Retirando a linha que contem o total 
df=filter(df, !df$Ano.Mes %in% c("Total"))

#Separando a coluna Ano.Mês.do.Obito em ano e mês

chave=paste(df$Ano.Mes,df$UF,sep="_")
df=separate(df,Ano.Mes,into=c("Mes","Ano"), sep="/")

#Colocando uma coluna de chave para fazer a união das futuras tabelas
df$chave=chave

#Retira os dois pontos que ficou na frente de texto
df$Mes=gsub("\\.\\.","",df$Mes)


#Transforma dados em inteiros, isso vai gerar NAS quando o valor for "-"
df$Suicidios=as.integer(df$Suicidios)

#Substituindo os valores NAs gerados por "0"
df[is.na(df$Suicidios),"Suicidios"]=0

return(df)
}


df= trataDados(df,29)

## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(x)` instead of `x` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

## Warning in trataDados(df, 29): NAs introduzidos por coerção

df=filter(df, !df$UF %in% c("Total"))

Tratando os os dados relacionado aos obitos por sexo

SuicidioFem= read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorSexoFem.csv",sep=";")
SuicidioMasc= read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorSexoMasc.csv",sep=";")
SuicidioIgn= read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorSexoIgn.csv",sep=";")

SuicidioFem=trataDados(SuicidioFem,29)

## Warning in trataDados(SuicidioFem, 29): NAs introduzidos por coerção

SuicidioMasc=trataDados(SuicidioMasc,29)

## Warning in trataDados(SuicidioMasc, 29): NAs introduzidos por coerção

SuicidioIgn= trataDados(SuicidioIgn,16)

## Warning in trataDados(SuicidioIgn, 16): NAs introduzidos por coerção

SuicidioFem$Sexo="Feminino"
SuicidioMasc$Sexo="Masculino"
SuicidioIgn$Sexo="Ign"
sexo=rbind(SuicidioFem,SuicidioMasc)
sexo=rbind(sexo,SuicidioIgn)

Tratando dados Relacionados a Cor

SuicidioPorCor<- read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorCor.csv", header = T, sep=";")
SuicidioPorCor=trataDados(SuicidioPorCor,8)

## Warning in trataDados(SuicidioPorCor, 8): NAs introduzidos por coerção

names(SuicidioPorCor)[3]="Cor"
SuicidioPorCor=filter(SuicidioPorCor, !SuicidioPorCor$Cor %in% c("Total"))
SuicidioPorCor$chave=NULL

Tratando dados Relacionados a Cor padrão IBGE

# Colocando os dados em Padrão IBGE
PretoPardo=group_by(filter(SuicidioPorCor, SuicidioPorCor$Cor %in% c("Preta","Parda")),Mes,Ano)%>%summarise(Suicidios=sum(Suicidios))

## `summarise()` regrouping output by 'Mes' (override with `.groups` argument)

PretoPardo$Cor="Pretos e Pardos"

SuicidioPorCorIBGE=rbind(filter(SuicidioPorCor, !SuicidioPorCor$Cor %in% c("Preta","Parda")), PretoPardo)

Tratando dados Relacionados a suicidio por Local

SuicidioPorLocal<-read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorLocal.csv", header = T, sep=";")
SuicidioPorLocal=trataDados(SuicidioPorLocal,8)

## Warning in trataDados(SuicidioPorLocal, 8): NAs introduzidos por coerção

names(SuicidioPorLocal)[3]=c("Local")
SuicidioPorLocal=filter(SuicidioPorLocal,!SuicidioPorLocal$Local %in% c("Total"))

Tratamento Dados Relacionados a suicidio por Faixa Etaria

SuicidioPorFaixaEtaria<-read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorFaixaEtaria.csv", header = T, sep=";")
SuicidioPorFaixaEtaria=trataDados(SuicidioPorFaixaEtaria,15)

## Warning in trataDados(SuicidioPorFaixaEtaria, 15): NAs introduzidos por coerção

names(SuicidioPorFaixaEtaria)[3]="FaixaEtaria"
SuicidioPorFaixaEtaria=filter(SuicidioPorFaixaEtaria,!SuicidioPorFaixaEtaria$FaixaEtaria %in% c("Total"))

Tratamento Dados Relacionados a suicidio por Estado Civil

SuicidioPorEstadoCivil<- read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorEstadoCivil.csv", header = T, sep=";")
SuicidioPorEstadoCivil=trataDados(SuicidioPorEstadoCivil,29)

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 168 rows [1, 2,
## 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].

SuicidioPorEstadoCivil$Ano=NULL
names(SuicidioPorEstadoCivil)[1]="EstadoCivil"

Tratamento de Dados Relacionado a suicidio por Tipo de Causa

SuicidioPorCausa<- read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioPorCausa.csv", header = T, sep=";")
SuicidioPorCausa=trataDados(SuicidioPorCausa,5)

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 72 rows [1, 2, 3,
## 4, 5, 7, 9, 10, 11, 13, 14, 15, 17, 18, 19, 21, 22, 23, 26, 27, ...].

## Warning in trataDados(SuicidioPorCausa, 5): NAs introduzidos por coerção

SuicidioPorCausa$Ano=NULL
names(SuicidioPorCausa)[1]="Causa"
names(SuicidioPorCausa)[2]="Sexo"
SuicidioPorCausa=filter(SuicidioPorCausa,!SuicidioPorCausa$Sexo %in% c("Total"))

Tratamento de dados relacionado a suicidio por cor e idade

SuicidioCorIdade<- read.csv("C:/Users/italo/OneDrive/Documents/Portfolio/TCC/Dados/SuicidioCorIdade.csv", header = T, sep=";")
SuicidioCorIdade= trataDados(SuicidioCorIdade,15)

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 84 rows [1, 2, 3,
## 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].

## Warning in trataDados(SuicidioCorIdade, 15): NAs introduzidos por coerção

SuicidioCorIdade$Ano=NULL
names(SuicidioCorIdade)[1]="Cor"
names(SuicidioCorIdade)[2]="Idade"
SuicidioCorIdade=filter(SuicidioCorIdade,!SuicidioCorIdade$Idade %in% c("Total"))
SuicidioCorIdade$Cor=as.factor(SuicidioCorIdade$Cor)
SuicidioCorIdade$Idade=as.factor(SuicidioCorIdade$Idade)
SuicidioCorIdade$chave=NULL

## OBS= A base por cor esta diferente da base por cor e idade

Tratando dados relacionados a suicidio por cor e idade pelo padrão IBGE

# Colocando os dados em Padrão IBGE
PretoPardo=group_by(filter(SuicidioCorIdade, SuicidioCorIdade$Cor %in% c("Preta","Parda")),Idade)%>%summarise(Suicidios=sum(Suicidios))

## `summarise()` ungrouping output (override with `.groups` argument)

PretoPardo$Cor="Pretos e Pardos"

SuicidioCorIdadeIBGE=rbind(filter(SuicidioCorIdade, !SuicidioCorIdade$Cor %in% c("Preta","Parda")), PretoPardo)

Análise estatisitca

hist(df$Suicidios, main = "Histograma taxa de suicídios entre 1996 até 2019")

boxplot(df$Suicidios, main="Boxplot taxa de suicídios entre 1996 até 2019")

summary(df$Suicidios)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    7.00   15.00   28.58   35.00  229.00

Analise Descritiva dos dados

Com os dados devidamente tratados foi feito uma análise exploratória com objetivo de entender melhor o tema no brasil. As divisões abaixo está relacionada a cada base que foi tratada.

Suicidio Por estado

x=aggregate(df$Suicidios, list(local=df$UF), sum)
x=x[order(x[,2],decreasing = T),]
x$indice=c(1:27)
ggplot(x, aes(x=reorder(local,indice), y=x)) + geom_col() + theme_classic() + labs(title ="Suicidios por Estado") + ylab("Número de Suicídios") + xlab("")

Aux=df %>%
  group_by(UF,Ano) %>% summarise(Suicidios=sum(Suicidios))

## `summarise()` regrouping output by 'UF' (override with `.groups` argument)

ggplot(Aux, aes(x=UF, y=Suicidios)) + geom_col() + labs(title ="Taxa de Suicídios antes 2005?") + ylab("Taxa de Suicídios")+ xlab("") + theme_classic()+ theme(legend.direction = "vertical", axis.text.x = element_text(angle=-90)) + facet_wrap(Aux$Ano==c(1996:2005), ncol = 1)

## Warning in Aux$Ano == c(1996:2005): comprimento do objeto maior não é múltiplo
## do comprimento do objeto menor

-Fazer análises

Suicidios Por Causas no geral e para cada sexo

#Ordenando Suicidios pelas principais causas
x=aggregate(SuicidioPorCausa$Suicidios, list(local=SuicidioPorCausa$Causa), sum)
x=x[order(x[,2],decreasing = T),]
x$indice=c(1:25)
x=x[x$indice<11,]
ggplot(x, aes(x=reorder(local,-indice), y=x)) + geom_col() + coord_flip() + theme_classic() + labs(title ="Principais formas de Suicídio no Brasil") + ylab("Número de Suicídios") + xlab("")

Fazer Análises

SuicidioPorCausaMasc=subset(SuicidioPorCausa,SuicidioPorCausa$Sexo=="Masc")
x=aggregate(SuicidioPorCausaMasc$Suicidios, list(local=SuicidioPorCausaMasc$Causa), sum)
x=x[order(x[,2],decreasing = T),]
x$indice=c(1:25)
x=x[x$indice<11,]
ggplot(x, aes(x=reorder(local,-indice), y=x)) + geom_col() + coord_flip() + theme_classic() + labs(title ="Homens cometerem Suicídio no Brasil") + ylab("Número de Suicídios") + xlab("")

Fazer Análises

SuicidioPorCausaFem=subset(SuicidioPorCausa,SuicidioPorCausa$Sexo=="Fem")
x=aggregate(SuicidioPorCausaFem$Suicidios, list(local=SuicidioPorCausaFem$Causa), sum)
x=x[order(x[,2],decreasing = T),]
x$indice=c(1:25)
x=x[x$indice<11,]
ggplot(x, aes(x=reorder(local,-indice), y=x)) + geom_col() + coord_flip() + theme_classic() + labs(title ="Mulheres cometerem Suicídio no Brasil") + ylab("Número de Suicídios") + xlab("")

Fazer Análises

SuicidioPorCausaIgn=subset(SuicidioPorCausa,SuicidioPorCausa$Sexo=="Ign")
x=aggregate(SuicidioPorCausaIgn$Suicidios, list(local=SuicidioPorCausaIgn$Causa), sum)
x=x[order(x[,2],decreasing = T),]
x$indice=c(1:25)
x=x[x$indice<11,]
ggplot(x, aes(x=reorder(local,-indice), y=x)) + geom_col() + coord_flip() + theme_classic() + labs(title ="IGN's cometerem Suicídio no Brasil") + ylab("Número de Suicídios") + xlab("")

Fazer Análises

SUicidio por cor

### Padrão sem ser IBGE
#x= aggregate(SuicidioPorCor$Suicidios, list(local=SuicidioPorCor$Cor), sum)
#x=x[order(x[,2], decreasing = T),]
#x$indice= c(1:6)
#ggplot(x, aes(x=reorder(local,indice), y=x)) + geom_col() + theme_classic() + labs(title #="Nº de Suicidios por Cor") + ylab("Número de Suicídios") + xlab("")
############


x= aggregate(SuicidioPorCorIBGE$Suicidios, list(local=SuicidioPorCorIBGE$Cor), sum)
x=x[order(x[,2], decreasing = T),]
x$indice= c(1:5)
ggplot(x, aes(x=reorder(local,indice), y=x)) + geom_col() + theme_classic() + labs(title ="Nº de Suicidios por Cor (Padrão IBGE)") + ylab("Número de Suicídios") + xlab("")

### TENTAR GRAFICO DE LINHA/ APREDER A TRABALHAR COM DATAS

#install.packages("lubridate")
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

SuicidioPorCor$Ano= as.factor(SuicidioPorCor$Ano)
SuicidioPorCor$Cor=as.factor(SuicidioPorCor$Cor)
ggplot(SuicidioPorCor, aes(x=Cor, y=Suicidios, fill=Cor)) + geom_col() + labs(title ="Nº de Suicidios por Cor e por Ano") + ylab("Número de Suicídios") + xlab("") + theme_classic()+ theme(legend.direction = "vertical", axis.text.x = element_text(angle=-90)) + facet_wrap(SuicidioPorCor$Ano)

#Gráfico padrao IBGE
SuicidioPorCorIBGE$Ano= as.factor(SuicidioPorCorIBGE$Ano)
SuicidioPorCorIBGE$Cor=as.factor(SuicidioPorCorIBGE$Cor)
ggplot(SuicidioPorCorIBGE, aes(x=Cor, y=Suicidios, fill=Cor)) + geom_col() + labs(title ="Nº de Suicidios por Cor e por Ano(Padrão IBGE)") + ylab("Número de Suicídios") + xlab("") + theme_classic()+ theme(legend.direction = "vertical", axis.text.x = element_text(angle=-90)) + facet_wrap(SuicidioPorCorIBGE$Ano)

### Dados Fora do Padrão IBGE
#ggplot(SuicidioCorIdade, aes(x=Cor, y=Suicidios, fill=Cor))+ geom_col() + labs(title ="Nº de #Suicidios por Cor e por Idade") + ylab("Número de Suicídios") + xlab("") + #facet_wrap(SuicidioCorIdade$Idade) + theme_classic()+ theme(legend.direction = "vertical", #axis.text.x = element_text(angle=-90))
########


##Padrão IBGE e retirando idades menores de 10 anos
SuicidioCorIdadeIBGE=filter(SuicidioCorIdadeIBGE,!SuicidioCorIdadeIBGE$Idade %in% c("X1.a.4.anos", "Menor.1.ano"))
ggplot(SuicidioCorIdadeIBGE, aes(x=Cor, y=Suicidios, fill=Cor))+ geom_col() + labs(title ="Nº de Suicidios por Cor e por Idade(Padrão IBGE)") + ylab("Número de Suicídios") + xlab("") + facet_wrap(SuicidioCorIdadeIBGE$Idade) + theme_classic()+ theme(legend.direction = "vertical", axis.text.x = element_text(angle=-90))

Suicidio por Sexo

sexo$Ano= as.factor(sexo$Ano)
ggplot(sexo, aes(x=Sexo, y=Suicidios, fill=Sexo))+ geom_col() + labs(title ="Suicidio por Sexo") + ylab("Número de Suicídios") + xlab("") + theme_classic() +scale_y_continuous(labels = function(x) format(x, scientific = FALSE))

sexo$Ano= as.factor(sexo$Ano)
ggplot(sexo, aes(x=Sexo, y=Suicidios, fill=Sexo))+ geom_col() + labs(title ="Suicidio por Sexo por Ano") + ylab("Número de Suicídios") + xlab("") + theme_classic() +scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) + facet_wrap(sexo$Ano) + coord_flip()+theme(legend.direction = "vertical", axis.text.x = element_text(angle=-90))

Graficos relacionado ao estado civil

x= aggregate(SuicidioPorEstadoCivil$Suicidios, list(local=SuicidioPorEstadoCivil$EstadoCivil), sum)
x=x[order(x[,2], decreasing = T),]
x$indice= c(1:6)
ggplot(x, aes(x=reorder(local,indice), y=x)) + geom_col() + theme_classic() + labs(title ="Suícidio por Estado Civil") + ylab("Número de Suicídios") + xlab("")

x=SuicidioPorEstadoCivil[SuicidioPorEstadoCivil$UF=="PE",]
x= aggregate(x$Suicidios, list(local=x$EstadoCivil), sum)
x=x[order(x[,2], decreasing = T),]
x$indice= c(1:6)
ggplot(x, aes(x=reorder(local,indice), y=x)) + geom_col() + theme_classic() + labs(title ="Suícidio por Estado Civil em PE") + ylab("Número de Suicídios") + xlab("")

Gráfico de local

x= aggregate(SuicidioPorLocal$Suicidios, list(local=SuicidioPorLocal$Local), sum)
x=x[order(x[,2], decreasing = T),]
x$indice= c(1:6)
ggplot(x, aes(x=reorder(local,indice), y=x)) + geom_col() + theme_classic() + labs(title ="Suicidio por local de Óbito") + ylab("Número de Suicídios") + xlab("")

Parte 3- Serie Temporal (Prevendo valores de Suicidio)

#install.packages("forecast")
library(forecast)

## Warning: package 'forecast' was built under R version 4.0.5

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

#install.packages("fpp3")
library(fpp3)

## Warning: package 'fpp3' was built under R version 4.0.5

## -- Attaching packages -------------------------------------------- fpp3 0.4.0 --

## v tibble      3.0.1     v feasts      0.2.1
## v tsibble     1.0.0     v fable       0.3.0
## v tsibbledata 0.3.0

## Warning: package 'tsibble' was built under R version 4.0.4

## Warning: package 'tsibbledata' was built under R version 4.0.4

## Warning: package 'feasts' was built under R version 4.0.4

## Warning: package 'fabletools' was built under R version 4.0.4

## Warning: package 'fable' was built under R version 4.0.4

## -- Conflicts ------------------------------------------------- fpp3_conflicts --
## x lubridate::date()      masks base::date()
## x plotly::filter()       masks dplyr::filter(), stats::filter()
## x fabletools::forecast() masks forecast::forecast()
## x tsibble::intersect()   masks base::intersect()
## x tsibble::interval()    masks lubridate::interval()
## x dplyr::lag()           masks stats::lag()
## x tsibble::setdiff()     masks base::setdiff()
## x tsibble::union()       masks base::union()

###Criando Serie

Decidir utilizar o objeto Tsibble para trabalhar com series temporais, atualmente é o método mais estruturado para trabalhar com esse tipo de problema. Embora não seja o caso desta análise, seria possível trabalhar com diversas series temporais ao mesmo tempo nesse objeto, em uma análise mais completa seria possível criar diversas séries temporais para cada estado do Brasil.

dados=group_by(df,Ano,Mes)%>%summarise(Suicidios=sum(Suicidios))

## `summarise()` regrouping output by 'Ano' (override with `.groups` argument)

dados$Mes=yearmonth(ymd(paste(dados$Ano,dados$Mes,"1",sep="/")))
serie=dados[,c("Mes","Suicidios")] %>%
  as_tsibble(index = Mes)

serie

Grafico de tempo

autoplot(serie,Suicidios) +
  labs(title="Taxa de suicídios no Brasil",
       subtitle = "Periodo: 1996-2019") + theme_bw()

Avaliando a serie temporal sobre suas principais características vemos que ela demonstra uma forte tendência de crescimento ao longo dos anos, parece existir um comportamento sazonal não muito forte e nenhum ciclo aparente no horizonte temporal de dados coletados.

Parcelas Sazonais

Para entender melhor o comportamento da serie, decidir avaliar seu comportamento por mês utilizando um gráfico sazonal para encontrar possíveis padrões nos dados.

serie %>%
  gg_season(Suicidios, labels = "both") + 
  labs(title = "Gráfico Sazonal: Suícidos no Brasil", 
       y="Suicidios/População")+
  theme_bw()

Analisando os dados Sozonais é perceptivel o crescimento das series todos os anos, embora não se pode afirmar que isso é algo puramente devido ao número de suicidios ou se deve também a melhora da forma de coleta dessas informações. Além disso, parece que os dados seguem um padrão de queda e aumento entre março, abril e maio. Também da para se notar um crescimento padrão no mês de setembro,o que justifica setembro ser o mês de prevenção ao suicidio.

Gráfico de Subséries Temporais

#dados=group_by(df,Ano,Mes)%>%summarise(Suicidios=sum(Suicidios))
ts(dados$Suicidios, start = 1996, frequency = 12) %>%
  ggsubseriesplot(ts) + 
  labs(title= "Média de Suicídios por Mês",
       y="Suicídios/População") +theme_bw()

Com esse gráfico fica mais evidente o padrão nos meses de fevereiro,março e abril. Com ele foi possivel perceber que existe uma flutuação em que os mês em média com menor quantidade de suicidios é julho.

Decomposição da série temporal

Para finalizar essa análise exploratória farei a decomposição da série nos seguintes componentes: -Ciclo de tendência -Sazonalidade -Componentes remanescente

Isso ajudará a entender o comportamente da taxa de suicídios no Brasil, mas antes de decompor irei fazer uma transformação e ajustar seus dados para facilitar a análise e posteriores manipulações que serão feitas.

Transformaçao e ajuste

Utilizarei a transformação Box-Cox, onde todos os dados serão colocados em escala logaritmica com base em um parametro lambda. O melhor valor para esse parametro é aquele que busca deixar o padrão da sazonalidade da serie praticamente a mesma, isso simplificará a série e fará com que as futuras previsões feitas sejam menos sensiveis a pontos flutuantes. Para estimar o valor de lambda utilizei o recurso guerrero, ele escolhe o melhor valor para lambda automaticamente e logo em seguida plotei a série ajustada.

serie %>%
  autoplot(Suicidios) + 
           labs(title="Serie ajustada",
           y="Suicídios/população") + theme_bw()

#Colocando uma coluna dos valores em log

Decidir optar pela decomposição aditiva, onde y(t)= S(t) + T(t) + R(t) pois o comportamento da serie não varia ao longo do tempo em relação a sua sazonalidade e tendência. Foi escolhido o método STL para fazer a decomposição opr ele ser o mais utilizado por diversos motivos como ser mais robusto a outline e poder lidar com qualquer tipo de sazonalidade. Definir o ciclo de tendência variando a cada 24 Meses e uma sazonalidade padrão ao longo dos anos.

serie %>%
  model(STL(Suicidios~season(window="periodic")+trend(window=24),robust=TRUE)) %>%
  components() %>%
  autoplot() + theme_bw()

A serie histórica mostra uma forte tendência de crescimento, mas uma baixa sazonalidade sendo seu componente remanescente mostrando que existem variáveis externas que possam estar influênciando os número de suícidios.

Para finalizar decidir encontrar a força da tendência e da sazonalidade.

serie %>%
  features(Suicidios, feat_stl) %>%
  select(c("trend_strength", "seasonal_strength_year"))

A série histórica possui uma alta força de tendência e uma força média de sazonalidade

Outros Recursos

Calculei a entropia espectral (Shannon) que é uma medidade que mostra quão fácil é prever os proximos valores dessa série. Quanto mais próxima de 0 mais fácil será a previsão, o valor para a série estudada deu 0.10, mostrando que é uma série fácil de prever.

serie %>%
  features(Suicidios,feat_spectral)

Escolhendo Modelo de Series temporais

Método de Suavização Exponêncial

Como primeira opção de modelo escolhi o método de suavização exponencial (ETS), para usarmos esse tipo de abordagem é necessário estimar os parametros para o erro (Aditivo ou Multiplicativo), tendência(Nenhuma,Aditiva,Aditiva amortecida) e a sazonalidade(Nenhuma, Atiditivo ou Multiplicativo). Para encontrar esses parametros utilizei a função do R ets() que utiliza o critério de AICa para encontrar um modelo apropriado de forma automatica

fit<-serie %>%
  model(ETS(Suicidios))
report(fit)

## Series: Suicidios 
## Model: ETS(M,Ad,M) 
##   Smoothing parameters:
##     alpha = 0.29069 
##     beta  = 0.005560979 
##     gamma = 0.0001049667 
##     phi   = 0.9714889 
## 
##   Initial states:
##         l        b      s1        s2       s3        s4       s5        s6
##  550.9013 1.779263 1.05667 0.9973694 1.033154 0.9941494 1.000316 0.9363775
##         s7        s8        s9      s10       s11      s12
##  0.9242413 0.9773794 0.9895966 1.056725 0.9703554 1.063666
## 
##   sigma^2:  0.0021
## 
##      AIC     AICc      BIC 
## 3692.164 3694.706 3758.097

O método escolhido foi o ETS(M,Ad,M).

fit<-serie %>%
  model(ETS(Suicidios ~ error("M") + trend("Ad") + season("A")))
report(fit)

## Series: Suicidios 
## Model: ETS(M,Ad,A) 
##   Smoothing parameters:
##     alpha = 0.2629962 
##     beta  = 0.005519299 
##     gamma = 0.0003182844 
##     phi   = 0.9797508 
## 
##   Initial states:
##         l        b       s1        s2       s3        s4        s5        s6
##  553.6935 3.069131 45.62786 -4.783519 28.82009 -1.327922 -2.149453 -44.02996
##         s7        s8        s9      s10       s11      s12
##  -52.39599 -15.00039 -11.42982 42.59666 -24.17899 38.25143
## 
##   sigma^2:  0.0022
## 
##      AIC     AICc      BIC 
## 3701.003 3703.546 3766.936

Inicialmente dividir os dados de treino e teste

# Separando dados em treino e teste
treino=dados[dados$Ano<=2017,c("Mes","Suicidios")] %>%
  as_tsibble(index = Mes)

teste= dados[dados$Ano>2017,c("Mes","Suicidios")] %>%
  as_tsibble(index = Mes)


#Criando um fit do metodo ETS
treino_fit <- treino %>%
  model(Previsao=ETS(Suicidios))

#Prevendo os proximos 24 meses
previsao = treino_fit %>%
  forecast(h=24)

Analise de residuos

treino_fit %>%
  gg_tsresiduals()

residuos=treino_fit %>% residuals()

Visualizando os resultados

autoplot(treino, .vars = Suicidios, color="black") +
  autolayer(teste, .vars=Suicidios, color="black") + theme_bw()+
  autolayer(previsao, level=NULL, color="red") + 
  labs(title = "Previsão da Serie histórica para os próximos 2 Anos", y="Suicidios/100k", subtitle = "2018 Jan - 2019 Dez") + guides(colour= guide_legend(title = "Forecast", label = TRUE))

Avaliando apenas o intervalo de previsão

prev= cbind(previsao$.mean, teste)  
names(prev)<- c("Previsto","Periodo","Real")

ggplot(prev, aes(x=Periodo,y=Real, color="black")) + geom_line() + geom_line(aes(x=Periodo, y=Previsto, color="red")) + theme_bw() + labs(title = "Valores reais e valores previstos ",y="Suicidios/100k", subtitle = "Periodo:2018 Jan - 2019 Dez") + guides(colour=guide_legend(title="Series")) + 
  scale_color_manual(labels = c("Real", "Previsto"), values = c("black", "red"))

autoplot(previsao) + 
  autolayer(teste,.vars=Suicidios) + theme_bw()

Metodo Arima

Teste de normalidade

par(mfrow=c(1,2))
hist(serie$Suicidios, main="Histograma normal")
hist(log(serie$Suicidios), main="Histograma dados transformados")

tseries::jarque.bera.test(serie$Suicidios)

## 
##  Jarque Bera Test
## 
## data:  serie$Suicidios
## X-squared = 16.117, df = 2, p-value = 0.0003164

tseries::jarque.bera.test(log(serie$Suicidios))

## 
##  Jarque Bera Test
## 
## data:  log(serie$Suicidios)
## X-squared = 9.5829, df = 2, p-value = 0.0083

Lags

#par(mfrow=c(2,1))
acf(serie$Suicidios, main="Função de Autocorrelação")

pacf(serie$Suicidios, main="Função de Autocorrelação parcial")

Estacionariedade da Serie

Olhando a serie ao longo do tempo é possível perceber que ela é não estacionária pois seus valores dependem do tempo (Tendência e Sazonalidade). Nessa secção tornarei a serie estacionária, para finalmente utilizar o método Arima.

Será utilizado o teste de raiz unitária, são testes de hipótese que podem ser utilizados para avaliar objetivamente se uma série temporal é estacionária ou não. Existem diversos testes diferentes, para essa análise escolhi o teste Kwiatkowski-Phillips-Schmidt-Shin.

Para esse teste, temos como a hipótese nula que a série é estacionária e a hipotese alternativa de que a série não é estacionária. Será definido o valor de p=0.01

# Colocando uma coluna de log
treino %>%
  features(Suicidios, unitroot_kpss)

A estatistica de testes (4.28) é maior que o valor crítico de 0.01, logo a hipótese nula é rejeitada.

Sabendo disso, é necessário transformar a série em uma série estaciónaria.

Diferenciação ordinária e diferenciação ordinal

Existem duas formas de tratar uma serie não-estacionária:

Diferenciação ordinária: É calculado as diferenças entre as observações consecutivas

-Diferenciação sazonal: É calculado as diferenças entre as observações de um ano para outro.

Segundo o livro FPP3, deve-se escolher inicialmente a diferenciação sazonal se a força sazonal da série for menor que 0.64 nenhuma diferenciação sazonal é sugerida. Anteriormente mostrei que a força sazonal dessa série é de 0.61, devido a isso utilizarei a diferencição ordinária como primeira opção para transformar a série.

Existe uma função no R chamana unitroot_ndiffs() que utiliza os testes KPSS para determinar o número apropriado de diferenças ordinarias a serem utilizados

treino %>%
  features(Suicidios, unitroot_ndiffs)

Como mostrado acima, apenas uma diferença é necessária para transformar os dados em estácionários

treino_arima=treino %>%
  mutate(log_suicidios=log(Suicidios)) %>%
  mutate(dif_suicidios= difference(log(Suicidios)))
?difference

## starting httpd help server ... done

Utilizando o teste de Kwiatkowski-Phillips-Schmidt-Shin para um valor p=0.01 para a série transformada:

treino_arima %>%
  features(dif_suicidios, unitroot_kpss)

Com um valor para o teste de (0.0227), não podemos recusar a hipotese nula de que a série é estacionária. Abaixo segue um plot mostrando as diferentes transformações e as novas séries que foram geradas

treino_arima %>%
  pivot_longer(-Mes,names_to="Type",values_to="Series")%>%
  mutate(
    Type= factor(Type, levels = c("Suicidios",
                                  "log_suicidios",
                                  "dif_suicidios"))
  ) %>%
  ggplot(aes(x=Mes, y=Series)) + geom_line()+
  facet_grid(vars(Type), scales="free_y") + labs(title = "Série História taxa de suícidio no Brasil ", subtitle = "Periodo: 1996 jan - 2017 Dez") + theme_bw()

Método ARIMA

fit_arima=treino_arima %>%
  model(ARIMA(Suicidios))
report(fit_arima)

## Series: Suicidios 
## Model: ARIMA(0,1,3)(0,0,2)[12] w/ drift 
## 
## Coefficients:
##           ma1     ma2      ma3    sma1    sma2  constant
##       -0.7357  0.1067  -0.2743  0.3833  0.2082    1.9275
## s.e.   0.0610  0.0695   0.0665  0.0665  0.0558    0.4161
## 
## sigma^2 estimated as 1797:  log likelihood=-1357.36
## AIC=2728.72   AICc=2729.16   BIC=2753.72

fit_arima %>%
  forecast(h=24) %>%
  autoplot(teste)

rbind(accuracy(fit_arima), accuracy(fit))

arima(serie$Suicidios,c(1,1,3),c(0,0,2,12))

## 
## Call:
## arima(x = serie$Suicidios, order = c(1, 1, 3), seasonal = c(0, 0, 2, 12))
## 
## Coefficients:
##           ar1      ma1     ma2      ma3    sma1    sma2
##       -0.0435  -0.6746  0.0663  -0.1841  0.3942  0.2335
## s.e.   0.3023   0.2973  0.2102   0.0649  0.0596  0.0513
## 
## sigma^2 estimated as 1861:  log likelihood = -1489.34,  aic = 2992.67

#Igual ao prof

Identificação

Estacionariedade Será utilizado o teste de raiz unitária, são testes de hipótese que podem ser utilizados para avaliar objetivamente se uma série temporal é estacionária ou não. Existem diversos testes diferentes, para essa análise escolhi o teste Kwiatkowski-Phillips-Schmidt-Shin.