Séries temporais: Coronavirus Colômbia
Danna Lesley Cruz Reyes
Departamento de Estatística, Universidade Federal de Minas Gerais, Minas Gerais, Brazil
library(rjags)
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
library(coda)
library(reshape2)
library(ggplot2)data = read.csv("casos.csv", sep=";")
data$Fecha <- as.Date(data$Fecha,format = "%d/%m/%y") Grafica Séries temporais
plot1<-ggplot(data,aes(Fecha, acumulada, group=1)) +
geom_point() +
geom_line(color="red") +
labs(x = "Fecha", y = "Casos", title = "Séries temporais Colômbia") +
scale_x_date(date_breaks="1 day", date_labels="%b %d") + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5))
plot1
set.seed(50)
k=30
T <- length(data$t)
dat <- list("T" = T, "y" = data$acumulada, "k"=k) # names list of numbers
##### Initial values
iniciais <- function(){
inits = list("a"=rgamma(1,0.01,0.01),"b"=rgamma(1,0.2,0.2),"c"=rgamma(1,0.01,0.01))
}
jags.m <- jags.model(file = "serieTiempoModel.bug", data=dat, inits=iniciais, n.chains=3)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 22
## Unobserved stochastic nodes: 33
## Total graph size: 431
##
## Initializing model
update(jags.m,5000)
params <- c("a", "b", "c", "yfuturo", "mufuturo")
samps <- coda.samples(jags.m, params, n.iter=20000)
out <- do.call(rbind.data.frame, samps) #para poder leer las muestras
a<-out$a
b<-out$b
c<-out$c
lim_mu<-4+k-1
ini_fu<-(4+k)
fin_fu<-dim(out)[2]
mufuturo<-out[,4:lim_mu]
futuro<-out[,ini_fu: fin_fu]
# Burn in of 5000. Start at 5001.
#estimacion de a b c
res<-cbind(a,b,c)
esti =round(t(apply(res,2, quantile,prob=c(0.025,0.5,0.975))),4)
plot.ts(cbind(a,b,c))
hist(c)
Especificação da assíntota (total de casos)
assint = (a/b)
est.Assin=quantile(assint,prob=c(0.025,0.5,0.975))#a/b
est.Assin
## 2.5% 50% 97.5%
## 806.4550 879.8265 973.4716As curvas usam valores de (a,b,c) associados a esse quantil da assintota.
Esses valores são obtidos com os seguintes comandos no R
rr = sort(assint)
q1 = round(0.025*length(rr))
q2 = round(0.5*length(rr))
q3 = round(0.975*length(rr))
limites = matrix(c(
a[which(assint==rr[q1])[[1]]],
b[which(assint==rr[q1])[[1]]],
c[which(assint==rr[q1])[[1]]],
a[which(assint==rr[q2])[[1]]],
b[which(assint==rr[q2])[[1]]],
c[which(assint==rr[q2])[[1]]],
a[which(assint==rr[q3])[[1]]],
b[which(assint==rr[q3])[[1]]],
c[which(assint==rr[q3])[[1]]]),3,3)
colnames(limites)=c("lim Inf", "Mediana","lim sup")
rownames(limites)=c("a", "b","c");limites
## lim Inf Mediana lim sup
## a 3.634273054 4.509661887 4.916957799
## b 0.004506539 0.005125629 0.005050956
## c 0.334065606 0.312115365 0.298724305
mu = matrix(0,k,3)
# Estimativas para mi considerando a assintota
for(l in 1:3)
{
for(j in 1:k){
mu[j,l] = limites[1,l]*exp(limites[3,l]*(j))/(1+limites[2,l]*exp(limites[3,l]*(j))) }
}
datas = seq(as.Date("2020-03-06"), by="days",length.out = nrow(mu))
data_res<-as.data.frame(mu, datas)
my_ggplot<- ggplot(data_res, aes(datas)) +
geom_line(aes(y = V1, colour = "Quantile 2.5%"), linetype="dashed", size = 1)+
geom_point(aes(y = V2), size=0.5) +
geom_line(aes(y = V2, colour = "Mediana"), size=0.2)+
geom_line(aes(y = V3, colour = "Quantile 97.5%"), linetype="dashed", size = 1) +
labs(x = "Fecha", y = "Numero acumulado de casos", title = "Séries temporais Colômbia") + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5))+ scale_x_date(date_breaks="2 day", date_labels="%b %d")
my_ggplot + scale_color_discrete(name = " ")
#GRAFICA NMNC
dif = apply(mu,2,diff) # mi_t-mi_t-1
data_nmnc<-as.data.frame(dif, datas[-1])
my_ggplot<- ggplot(data_nmnc, aes(datas[-1])) +
geom_line(aes(y = V1, colour = "Quantile 2.5%"), linetype="dashed", size = 1)+
geom_point(aes(y = V2), size=0.5) +
geom_line(aes(y = V2, colour = "Mediana"), size=0.2)+
geom_line(aes(y = V3, colour = "Quantile 97.5%"), linetype="dashed", size = 1) +
labs(x = "Fecha", y = "Numero medio de novos casos", title = "Grafo NMNC ") + theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5))+ scale_x_date(date_breaks="2 day", date_labels="%b %d")
my_ggplot + scale_color_discrete(name = " ")
#MATRIZ
set.seed(50)
k=214
T <- length(data$t)
dat <- list("T" = T, "y" = data$acumulada, "k"=k) # names list of number
##### Initial values
iniciais <- function(){
inits = list("a"=rgamma(1,0.01,0.01),"b"=rgamma(1,0.2,0.2),"c"=rgamma(1,0.01,0.01))
}
jags.m <- jags.model(file = "serieTiempoModel.bug", data=dat, inits=iniciais, n.chains=3)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 22
## Unobserved stochastic nodes: 217
## Total graph size: 2087
##
## Initializing model
update(jags.m,5000)
params <- c("a", "b", "c", "yfuturo", "mufuturo")
samps <- coda.samples(jags.m, params, n.iter=20000)
out <- do.call(rbind.data.frame, samps)
a<-out$a
b<-out$b
c<-out$c
assint = (a/b)
est.Assin=quantile(assint,prob=c(0.025,0.5,0.975))#a/b
# as curvas usam valores de (a,b,c) associados a esse quantil da assint
# esses valores são obtidos com os seguintes comandos no R
rr = sort(assint)
q1 = round(0.025*length(rr))
q2 = round(0.5*length(rr))
q3 = round(0.975*length(rr))
limites = matrix(c(a[which(assint==rr[q1])[[1]]],
b[which(assint==rr[q1])[[1]]],
c[which(assint==rr[q1])[[1]]],
a[which(assint==rr[q2])[[1]]],
b[which(assint==rr[q2])[[1]]],
c[which(assint==rr[q2])[[1]]],
a[which(assint==rr[q3])[[1]]],
b[which(assint==rr[q3])[[1]]],
c[which(assint==rr[q3])[[1]]]),3,3)
colnames(limites)=c("lim Inf", "Mediana","lim sup")
rownames(limites)=c("a", "b","c");limites
## lim Inf Mediana lim sup
## a 4.754139599 4.118622677 4.788867683
## b 0.005906009 0.004685883 0.004911434
## c 0.315484548 0.320166546 0.302332244
mu = matrix(0,k,3)
for(l in 1:3)
{
for(j in 1:k){
mu[j,l] = limites[1,l]*exp(limites[3,l]*(j))/(1+limites[2,l]*exp(limites[3,l]*(j))) }
}
MATRIX<-matrix(ncol=216, nrow=3)
MATRIX[1, ]<-c(data[dim(data)[1], 3], est.Assin[1], mu[,1 ])
MATRIX[2, ]<-c(data[dim(data)[1], 3], est.Assin[2], mu[,2])
MATRIX[3, ]<-c(data[dim(data)[1], 3], est.Assin[3], mu[,3])
colnames(MATRIX)<-c("ultimo dado", "assintota", seq( from = 1, to = 14), seq( from = 1, to = 200))
rownames(MATRIX)<-c("lim Inf", "Mediana","lim sup")
print(MATRIX)
## ultimo dado assintota 1 2 3 4 5
## lim Inf 798 804.9697 6.465215 8.836988 12.06570 16.44970 22.38162
## Mediana 798 878.9451 5.636436 7.744644 10.63179 14.57723 19.95318
## lim sup 798 975.0453 6.436616 8.688563 11.71886 15.78875 21.24093
## 6 7 8 9 10 11 12 13
## lim Inf 30.37027 41.06070 55.24579 73.85757 97.92053 128.4459 166.2497 211.6981
## Mediana 27.24922 37.09791 50.29599 67.81125 90.75803 120.3180 157.5808 203.2915
## lim sup 28.51988 38.19333 50.97117 67.71494 89.42666 117.2007 152.1191 195.0751
## 14 1 2 3 4 5 6 7
## lim Inf 264.4272 323.1365 385.5824 448.8543 509.8856 566.0253 615.4541 657.3249
## Mediana 257.5282 319.3945 386.8704 456.9598 526.1693 591.1760 649.4289 699.4696
## lim sup 246.5275 306.2226 372.9723 444.5999 518.1449 590.3168 658.0627 719.0525
## 8 9 10 11 12 13 14 15
## lim Inf 691.6481 719.0350 740.4208 756.8405 769.2845 778.6228 785.5788 790.7317
## Mediana 740.9187 774.2283 800.3519 820.4507 835.6872 847.1088 855.5988 861.8701
## lim sup 771.9295 816.2956 852.5093 881.4095 904.0611 921.5654 934.9447 945.0856
## 16 17 18 19 20 21 22 23
## lim Inf 794.5332 797.3293 799.3813 800.8848 801.9851 802.7896 803.3774 803.8068
## Mediana 866.4812 869.8601 872.3297 874.1316 875.4445 876.4001 877.0953 877.6007
## lim sup 952.7232 958.4479 962.7234 965.9080 968.2752 970.0323 971.3351 972.3002
## 24 25 26 27 28 29 30 31
## lim Inf 804.1202 804.3490 804.5160 804.6379 804.7268 804.7916 804.8389 804.8735
## Mediana 877.9679 878.2348 878.4286 878.5694 878.6717 878.7459 878.7998 878.8390
## lim sup 973.0147 973.5435 973.9347 974.2241 974.4380 974.5962 974.7132 974.7996
## 32 33 34 35 36 37 38 39
## lim Inf 804.8986 804.9170 804.9304 804.9402 804.9473 804.9525 804.9563 804.9591
## Mediana 878.8674 878.8881 878.9030 878.9139 878.9218 878.9276 878.9317 878.9347
## lim sup 974.8635 974.9108 974.9457 974.9715 974.9906 975.0047 975.0151 975.0228
## 40 41 42 43 44 45 46 47
## lim Inf 804.9611 804.9626 804.9636 804.9644 804.9650 804.9654 804.9657 804.9659
## Mediana 878.9369 878.9385 878.9397 878.9405 878.9411 878.9416 878.9419 878.9421
## lim sup 975.0285 975.0327 975.0358 975.0381 975.0398 975.0410 975.0420 975.0427
## 48 49 50 51 52 53 54 55
## lim Inf 804.9661 804.9662 804.9663 804.9664 804.9664 804.9664 804.9665 804.9665
## Mediana 878.9423 878.9424 878.9425 878.9426 878.9426 878.9427 878.9427 878.9427
## lim sup 975.0432 975.0435 975.0438 975.0440 975.0442 975.0443 975.0444 975.0444
## 56 57 58 59 60 61 62 63
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0445 975.0445 975.0445 975.0445 975.0446 975.0446 975.0446 975.0446
## 64 65 66 67 68 69 70 71
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 72 73 74 75 76 77 78 79
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 80 81 82 83 84 85 86 87
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 88 89 90 91 92 93 94 95
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 96 97 98 99 100 101 102 103
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 104 105 106 107 108 109 110 111
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 112 113 114 115 116 117 118 119
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 120 121 122 123 124 125 126 127
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 128 129 130 131 132 133 134 135
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 136 137 138 139 140 141 142 143
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 144 145 146 147 148 149 150 151
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 152 153 154 155 156 157 158 159
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 160 161 162 163 164 165 166 167
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 168 169 170 171 172 173 174 175
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 176 177 178 179 180 181 182 183
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 184 185 186 187 188 189 190 191
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 192 193 194 195 196 197 198 199
## lim Inf 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665 804.9665
## Mediana 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427 878.9427
## lim sup 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446 975.0446
## 200
## lim Inf 804.9665
## Mediana 878.9427
## lim sup 975.0446usethis::use_git_config(user.name = “Danna Cruz”, # Seu nome user.email = “dcruzreyes@gmail.com”)
install.packages(“usethis”) library(usethis)
usethis::edit_r_environ()