Análisis de evolución de coeficientes
Claudio Bustos
2020-06-20
Modelo deceso basado en casos
Modelo con lags de 1 a 22
casos.nuevos<-ts(unido$casos.nuevos)
decesos.nuevos<-ts(unido$decesos.nuevos)
dyn.1<-dynlm(decesos.nuevos~L(casos.nuevos,c(1:22)))
summary(dyn.1)
#>
#> Time series regression with "ts" data:
#> Start = 23, End = 109
#>
#> Call:
#> dynlm(formula = decesos.nuevos ~ L(casos.nuevos, c(1:22)))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -95.519 -9.425 -1.157 7.058 177.532
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.916299 6.830481 -0.134 0.89371
#> L(casos.nuevos, c(1:22))1 -0.002099 0.002521 -0.833 0.40815
#> L(casos.nuevos, c(1:22))2 0.006255 0.002583 2.422 0.01829 *
#> L(casos.nuevos, c(1:22))3 -0.002904 0.002469 -1.176 0.24391
#> L(casos.nuevos, c(1:22))4 0.002782 0.014229 0.196 0.84560
#> L(casos.nuevos, c(1:22))5 -0.041102 0.015592 -2.636 0.01051 *
#> L(casos.nuevos, c(1:22))6 0.029277 0.017061 1.716 0.09100 .
#> L(casos.nuevos, c(1:22))7 0.019394 0.017290 1.122 0.26620
#> L(casos.nuevos, c(1:22))8 0.001135 0.016897 0.067 0.94666
#> L(casos.nuevos, c(1:22))9 -0.061245 0.018328 -3.342 0.00139 **
#> L(casos.nuevos, c(1:22))10 0.015036 0.020374 0.738 0.46321
#> L(casos.nuevos, c(1:22))11 0.009427 0.022842 0.413 0.68118
#> L(casos.nuevos, c(1:22))12 -0.025034 0.025414 -0.985 0.32830
#> L(casos.nuevos, c(1:22))13 0.050103 0.025950 1.931 0.05795 .
#> L(casos.nuevos, c(1:22))14 -0.020451 0.025830 -0.792 0.43143
#> L(casos.nuevos, c(1:22))15 -0.012260 0.026425 -0.464 0.64425
#> L(casos.nuevos, c(1:22))16 0.053532 0.027548 1.943 0.05639 .
#> L(casos.nuevos, c(1:22))17 0.050733 0.027033 1.877 0.06512 .
#> L(casos.nuevos, c(1:22))18 0.021920 0.027680 0.792 0.43133
#> L(casos.nuevos, c(1:22))19 0.013397 0.028213 0.475 0.63651
#> L(casos.nuevos, c(1:22))20 -0.087336 0.029065 -3.005 0.00379 **
#> L(casos.nuevos, c(1:22))21 0.064198 0.029987 2.141 0.03610 *
#> L(casos.nuevos, c(1:22))22 -0.045527 0.025196 -1.807 0.07548 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 42.2 on 64 degrees of freedom
#> Multiple R-squared: 0.8413, Adjusted R-squared: 0.7868
#> F-statistic: 15.42 on 22 and 64 DF, p-value: < 2.2e-16dd<-22
res1<-lapply(46:nrow(unido),function(i) {
dl<-unido[1:i,]
dym.l<-dynlm(ts(decesos.nuevos)~L(ts(casos.nuevos),c(1:dd)), data=dl)
#print(dym.l)
s1<-summary(dym.l)
res<-s1$coefficients[,c(1,4)]
colnames(res)<-c("est","p")
res<-res[-1,]
data.frame(coef.n=paste0("c",sprintf("%02d",1:dd)),res,
sign=res[,2]<0.05,
fecha=tail(dl$fecha,1),
observado=tail(dl$decesos.nuevos,1),
predicho=tail(predict(dym.l),1),
r2=s1$adj.r.squared
)
})
res2<-do.call(rbind,res1)obs.pred<-do.call(rbind,lapply(res1,function(x) {tail(x,1)[,c("fecha","observado","predicho"),drop=FALSE]}))
ggplot(melt(obs.pred, id.vars="fecha"),aes(x=fecha,y=value,color=variable))+geom_point()+geom_line()+scale_y_continuous(trans="log10")
Podemos mapear los coeficientes desde el día 30 (1 de Abril) en adelante. Resulta interesante ver que hay lags que resultan significativos en más ocasiones que otros.
t1<-table(res2$coef.n, res2$sign)[,2]
barplot(t1)
Podemos ver también como evolucionan los R² ajustados de la relación entre casos y decesos.
r1<-do.call(rbind,lapply(res1,function(x) {tail(x,1)[,c("fecha","r2"),drop=FALSE]}))
ggplot(r1,aes(x=fecha, y=r2))+geom_point()+geom_line()
paleta<-scale_fill_manual(values = colorRampPalette(brewer.pal(12, "Set3"))(22))
ggplot(res2,aes(x=fecha,y=est,color=coef.n,group=coef.n))+geom_point(aes(alpha=ifelse(sign,1,0.2)),size=2)+geom_line(alpha=0.3)+theme_bw()+paleta
Modelo de camas UCI basado en casos
unido.camas1<-unido[!is.na(unido$camas.uci),]
casos.nuevos.uci<-ts(unido.camas1$casos.nuevos)
camas.nuevas.uci<-ts(unido.camas1$camas.uci)
dyn.uci<-dynlm(camas.nuevas.uci~L(casos.nuevos.uci,c(1:14)))
summary(dyn.uci)
#>
#> Time series regression with "ts" data:
#> Start = 15, End = 79
#>
#> Call:
#> dynlm(formula = camas.nuevas.uci ~ L(casos.nuevos.uci, c(1:14)))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -31.533 -8.538 -2.367 6.336 43.678
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 6.223e+00 3.407e+00 1.826 0.07375 .
#> L(casos.nuevos.uci, c(1:14))1 8.605e-04 7.335e-04 1.173 0.24633
#> L(casos.nuevos.uci, c(1:14))2 6.989e-05 7.372e-04 0.095 0.92486
#> L(casos.nuevos.uci, c(1:14))3 1.311e-03 7.167e-04 1.829 0.07330 .
#> L(casos.nuevos.uci, c(1:14))4 -5.601e-04 4.411e-03 -0.127 0.89947
#> L(casos.nuevos.uci, c(1:14))5 1.303e-02 4.750e-03 2.743 0.00842 **
#> L(casos.nuevos.uci, c(1:14))6 2.486e-03 5.187e-03 0.479 0.63387
#> L(casos.nuevos.uci, c(1:14))7 5.777e-03 4.675e-03 1.236 0.22239
#> L(casos.nuevos.uci, c(1:14))8 1.105e-02 4.560e-03 2.424 0.01901 *
#> L(casos.nuevos.uci, c(1:14))9 -2.280e-03 4.827e-03 -0.472 0.63867
#> L(casos.nuevos.uci, c(1:14))10 -6.934e-03 5.278e-03 -1.314 0.19495
#> L(casos.nuevos.uci, c(1:14))11 -6.518e-03 5.743e-03 -1.135 0.26183
#> L(casos.nuevos.uci, c(1:14))12 -4.156e-03 6.167e-03 -0.674 0.50343
#> L(casos.nuevos.uci, c(1:14))13 -1.110e-02 6.157e-03 -1.802 0.07753 .
#> L(casos.nuevos.uci, c(1:14))14 -2.731e-03 5.790e-03 -0.472 0.63927
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 15.6 on 50 degrees of freedom
#> Multiple R-squared: 0.6194, Adjusted R-squared: 0.5128
#> F-statistic: 5.812 on 14 and 50 DF, p-value: 1.566e-06dd<-14
res1<-lapply(31:nrow(unido.camas1),function(i) {
dl<-unido.camas1[1:i,]
dym.l<-dynlm(ts(camas.uci)~L(ts(casos.nuevos),c(1:dd)), data=dl)
#print(dym.l)
s1<-summary(dym.l)
res<-s1$coefficients[,c(1,4)]
colnames(res)<-c("est","p")
res<-res[-1,]
data.frame(coef.n=paste0("c",sprintf("%02d",1:dd)),res,
sign=res[,2]<0.05,
fecha=tail(dl$fecha,1),
r2=s1$adj.r.squared
)
})
res2<-do.call(rbind,res1)t1<-table(res2$coef.n, res2$sign)[,2]
barplot(t1)
r1<-do.call(rbind,lapply(res1,function(x) {tail(x,1)[,c("fecha","r2"),drop=FALSE]}))
ggplot(r1,aes(x=fecha, y=r2))+geom_point()+geom_line()
paleta<-scale_fill_manual(values = colorRampPalette(brewer.pal(12, "Set3"))(14))
ggplot(res2,aes(x=fecha,y=est,color=coef.n,group=coef.n))+geom_point(aes(alpha=ifelse(sign,1,0.2)),size=2)+geom_line(alpha=0.3)+theme_bw()+paleta
Modelo de decesos basado en camas UCI
unido.camas<-unido[!is.na(unido$camas.uci),]
casos.nuevos.uci<-ts(unido.camas$camas.uci)
decesos.nuevos.uci<-ts(unido.camas$decesos.nuevos)
dyn.uci<-dynlm(decesos.nuevos.uci~L(casos.nuevos.uci,c(1:15)))
summary(dyn.uci)
#>
#> Time series regression with "ts" data:
#> Start = 16, End = 79
#>
#> Call:
#> dynlm(formula = decesos.nuevos.uci ~ L(casos.nuevos.uci, c(1:15)))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -140.05 -48.83 -6.98 35.56 335.09
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.4649 21.2039 0.022 0.9826
#> L(casos.nuevos.uci, c(1:15))1 0.3779 0.8454 0.447 0.6569
#> L(casos.nuevos.uci, c(1:15))2 -0.3822 0.8416 -0.454 0.6518
#> L(casos.nuevos.uci, c(1:15))3 0.2424 0.8457 0.287 0.7756
#> L(casos.nuevos.uci, c(1:15))4 -0.6859 0.9668 -0.709 0.4815
#> L(casos.nuevos.uci, c(1:15))5 1.8750 1.0083 1.860 0.0691 .
#> L(casos.nuevos.uci, c(1:15))6 -0.1752 1.0178 -0.172 0.8640
#> L(casos.nuevos.uci, c(1:15))7 -1.9674 1.0364 -1.898 0.0637 .
#> L(casos.nuevos.uci, c(1:15))8 -0.3513 1.0059 -0.349 0.7284
#> L(casos.nuevos.uci, c(1:15))9 -0.5131 1.0573 -0.485 0.6297
#> L(casos.nuevos.uci, c(1:15))10 1.0965 1.0551 1.039 0.3039
#> L(casos.nuevos.uci, c(1:15))11 1.2459 1.0611 1.174 0.2461
#> L(casos.nuevos.uci, c(1:15))12 -0.2851 1.0308 -0.277 0.7833
#> L(casos.nuevos.uci, c(1:15))13 0.2392 0.9261 0.258 0.7973
#> L(casos.nuevos.uci, c(1:15))14 1.2611 0.9017 1.399 0.1683
#> L(casos.nuevos.uci, c(1:15))15 1.1381 0.9076 1.254 0.2159
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 86.65 on 48 degrees of freedom
#> Multiple R-squared: 0.4508, Adjusted R-squared: 0.2792
#> F-statistic: 2.627 on 15 and 48 DF, p-value: 0.005732
res1.uci<-lapply(31:nrow(unido.camas),function(i) {
dl<-unido.camas[1:i,]
casos.nuevos.uci.l<-ts(dl$camas.uci)
decesos.nuevos.uci.l<-ts(dl$decesos.nuevos)
dym.l<-dynlm(decesos.nuevos.uci.l~L(casos.nuevos.uci.l,c(1:15)))
s1<-summary(dym.l)
res<-s1$coefficients[,c(1,4)]
colnames(res)<-c("est","p")
res<-res[-1,]
data.frame(coef.n=paste0("c",sprintf("%02d",1:15)),
res,
sign=res[,2]<0.05,
fecha=openxlsx::convertToDate(tail(dl$fecha,1)),
r2=s1$adj.r.squared,
observado=tail(dl$decesos.nuevos,1),
predicho=tail(predict(dym.l),1)
)
})
res2.uci<-do.call(rbind,res1.uci)obs.pred<-do.call(rbind,lapply(res1.uci,function(x) {tail(x,1)[,c("fecha","observado","predicho"),drop=FALSE]}))
ggplot(melt(obs.pred, id.vars="fecha"),aes(x=fecha,y=value,color=variable))+geom_point()+geom_line()+scale_y_continuous(trans="log10")
#> Warning in self$trans$transform(x): Se han producido NaNs
#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning in self$trans$transform(x): Se han producido NaNs
#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 1 rows containing missing values (geom_point).
t1.uci<-table(res2.uci$coef.n, res2.uci$sign)[,2]
barplot(t1.uci)
r1<-do.call(rbind,lapply(res1.uci,function(x) {tail(x,1)[,c("fecha","r2"),drop=FALSE]}))
ggplot(r1,aes(x=fecha, y=r2))+geom_point()+geom_line()+ylim(0,1)
#> Warning: Removed 2 rows containing missing values (geom_point).
#> Warning: Removed 1 row(s) containing missing values (geom_path).
paleta<-scale_fill_manual(values = colorRampPalette(brewer.pal(12, "Set3"))(14))
ggplot(res2.uci, aes(x=fecha,y=est,color=coef.n,group=coef.n))+geom_point(aes(alpha=ifelse(sign,1,0.2)),size=2)+geom_line(alpha=0.3)+theme_bw()+paleta
Modelo basado en ventiladores
unido.vent<-unido[!is.na(unido$vent.usados),]
casos.nuevos.vent<-ts(unido.vent$vent.usados)
decesos.nuevos.vent<-ts(unido.vent$decesos.nuevos)
dyn.vent<-dynlm(decesos.nuevos.vent~L(casos.nuevos.vent,c(1:14)))
summary(dyn.vent)
#>
#> Time series regression with "ts" data:
#> Start = 15, End = 66
#>
#> Call:
#> dynlm(formula = decesos.nuevos.vent ~ L(casos.nuevos.vent, c(1:14)))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -153.86 -49.22 -4.23 30.89 347.73
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 32.524663 29.766133 1.093 0.2816
#> L(casos.nuevos.vent, c(1:14))1 0.003204 0.621938 0.005 0.9959
#> L(casos.nuevos.vent, c(1:14))2 -1.031408 0.619519 -1.665 0.1044
#> L(casos.nuevos.vent, c(1:14))3 -0.483016 0.616843 -0.783 0.4386
#> L(casos.nuevos.vent, c(1:14))4 1.205396 0.589695 2.044 0.0481 *
#> L(casos.nuevos.vent, c(1:14))5 -0.362048 0.601397 -0.602 0.5508
#> L(casos.nuevos.vent, c(1:14))6 0.077995 0.607572 0.128 0.8985
#> L(casos.nuevos.vent, c(1:14))7 -1.111373 0.625580 -1.777 0.0839 .
#> L(casos.nuevos.vent, c(1:14))8 -0.454882 0.623691 -0.729 0.4704
#> L(casos.nuevos.vent, c(1:14))9 0.986137 0.610964 1.614 0.1150
#> L(casos.nuevos.vent, c(1:14))10 0.296516 0.632881 0.469 0.6422
#> L(casos.nuevos.vent, c(1:14))11 0.453663 0.640918 0.708 0.4835
#> L(casos.nuevos.vent, c(1:14))12 0.233339 0.645678 0.361 0.7199
#> L(casos.nuevos.vent, c(1:14))13 0.600804 0.637116 0.943 0.3518
#> L(casos.nuevos.vent, c(1:14))14 1.213782 0.632685 1.918 0.0628 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 96.02 on 37 degrees of freedom
#> Multiple R-squared: 0.4392, Adjusted R-squared: 0.227
#> F-statistic: 2.07 on 14 and 37 DF, p-value: 0.0388
res1.vent<-lapply(30:nrow(unido.vent),function(i) {
dl<-unido.vent[1:i,]
casos.nuevos.vent.l<-ts(dl$vent.usados)
decesos.nuevos.vent.l<-ts(dl$decesos.nuevos)
dym.l<-dynlm(decesos.nuevos.vent.l~L(casos.nuevos.vent.l,c(1:14)))
s1<-summary(dym.l)
res<-s1$coefficients[,c(1,4)]
colnames(res)<-c("est","p")
res<-res[-1,]
data.frame(coef.n=paste0("c",sprintf("%02d",1:14)),res,
sign=res[,2]<0.05,fecha=openxlsx::convertToDate(tail(dl$fecha,1)),
r2=s1$adj.r.squared)
})
res2.vent<-do.call(rbind,res1.vent)#t1.vent<-table(res2.vent$coef.n, res2.vent$sign)[,2]
#barplot(t1.uci)r1<-do.call(rbind,lapply(res1.vent,function(x) {tail(x,1)[,c("fecha","r2"),drop=FALSE]}))
ggplot(r1,aes(x=fecha, y=r2))+geom_point()+geom_line()+ylim(0,1)
paleta<-scale_fill_manual(values = colorRampPalette(brewer.pal(12, "Set3"))(14))
ggplot(res2.vent, aes(x=fecha,y=est,color=coef.n,group=coef.n))+geom_point(aes(alpha=ifelse(sign,1,0.2)),size=2)+geom_line(alpha=0.5)+theme_bw()