library(dlnm) ; library(splines) ; library(MASS) ; library(readxl) ;
library(ggplot2) ; library(writexl); library(openxlsx); library(Epi)
library(bbmle); library(openxlsx); library (dplyr); library(car)

library(quantmod)
library(MuMIn)
library(tseries)
library(gridExtra)
library(grid)
library(lattice)
library(qcc)
library(nortest)
library(DT)

library(readxl)
obs <- read_excel("obs.xlsx")
DT::datatable(obs)

obs$fecha=as.Date(obs$fecha)
obs$tipodia=as.factor(obs$tipodia)

argvar1_SO2 <- list(fun="lin")
templin <- list(fun="lin")

MODELO SO2 cardiovascular sin Rezago

maxlag <-0
cb1_SO2 <- crossbasis(obs$SO2,maxlag,argvar1_SO2)
templ <- crossbasis(obs$temperatura,maxlag,templin)

Construcion del MODELO

m1 <- glm(mcardio ~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
          data=obs, family=quasipoisson)

m_pois=glm(mcardio~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
           data=obs,family=poisson)

Cálculo QAIC

dfun <- function(object) {
  with(object,sum((weights * residuals^2)[weights > 0])/df.residual)
}
qAIC(m_pois,dispersion=dfun(m1))

## [1] 10898.37

Cálculo RR e IC 95%

(eff1 <- ci.lin(m1,subset="cb1_SO2",Exp=T))

##             Estimate      StdErr          z         P exp(Est.)      2.5%
## cb1_SO2 -0.001193546 0.002700365 -0.4419944 0.6584933 0.9988072 0.9935348
##            97.5%
## cb1_SO2 1.004107

predicción

varcen <- round(mean(obs$SO2),0)

cp <- crosspred(cb1_SO2,m1,cen=varcen,by=0.1)
cen <- cp$predvar[which.min(cp$allRRfit)]

pred <- crosspred(cb1_SO2,m1,cen=cen, by=0.1)

head(cbind(pred$allRRfit,pred$allRRlow,pred$allRRhigh))

##         [,1]      [,2]     [,3]
## 0.7 1.017457 0.9422953 1.098614
## 0.8 1.017336 0.9426817 1.097902
## 0.9 1.017214 0.9430682 1.097190
## 1   1.017093 0.9434548 1.096478
## 1.1 1.016971 0.9438416 1.095767
## 1.2 1.016850 0.9442286 1.095057

grafico

xlab <- expression(paste("SO2 concentracion mcg/m3"))

par(mar=c(5,4,3,1),mgp=c(3,1,0),las=1,cex.axis=0.9,cex.lab=1)
plot(pred,"overall",col="red",ylim=c(0.9,1.5),axes=T,lab=c(6,5,7),xlab=xlab,
     ylab="RR",main="Mortalidad cardiovascular")

summary(pred)

## PREDICTIONS:
## values: 146 
## centered at: 15.2 
## range: 0.7 , 15.2 
## lag: 0 0 
## exponentiated: yes 
## cumulative: no 
## 
## MODEL:
## parameters: 1 
## class: glm lm 
## link: log

pred$predvar

##   [1]  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.0  2.1
##  [16]  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9  3.0  3.1  3.2  3.3  3.4  3.5  3.6
##  [31]  3.7  3.8  3.9  4.0  4.1  4.2  4.3  4.4  4.5  4.6  4.7  4.8  4.9  5.0  5.1
##  [46]  5.2  5.3  5.4  5.5  5.6  5.7  5.8  5.9  6.0  6.1  6.2  6.3  6.4  6.5  6.6
##  [61]  6.7  6.8  6.9  7.0  7.1  7.2  7.3  7.4  7.5  7.6  7.7  7.8  7.9  8.0  8.1
##  [76]  8.2  8.3  8.4  8.5  8.6  8.7  8.8  8.9  9.0  9.1  9.2  9.3  9.4  9.5  9.6
##  [91]  9.7  9.8  9.9 10.0 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 11.0 11.1
## [106] 11.2 11.3 11.4 11.5 11.6 11.7 11.8 11.9 12.0 12.1 12.2 12.3 12.4 12.5 12.6
## [121] 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 13.5 13.6 13.7 13.8 13.9 14.0 14.1
## [136] 14.2 14.3 14.4 14.5 14.6 14.7 14.8 14.9 15.0 15.1 15.2

plot(pred, "slices", type="p", pch=19, cex=1.5, lag=c(0,0,0),
     var=c(1.1,11.1,15.2), ci="bars",ci.arg=list(col=grey(0.8)),
     ylab="RR",main="Lag-specific effects")

Tabla con estimaciones

tablag <- matrix(NA,7+1,3,dimnames=list(paste("Lag",0:7),
                                        c("RR","ci.low","ci.hi")))

# Correr el ciclo
for(i in 0:7) {
  # LAG de las variables i=0, 1 ,2, 3...7
  SO2lag=Lag(obs$SO2,i)
  temlag=Lag(obs$temperatura,i)  # el i puede fijarse 
  m1 <- glm(mcardio~ SO2lag+ diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(temlag),
            data=obs, family=quasipoisson)
  tablag[i+1,] <- ci.lin(m1,subset="SO2lag",Exp=T)[5:7]
}
tablag

##              RR    ci.low     ci.hi
## Lag 0 0.9988072 0.9935348 1.0041075
## Lag 1 0.9934055 0.9881862 0.9986524
## Lag 2 0.9942683 0.9890516 0.9995125
## Lag 3 0.9941522 0.9889385 0.9993933
## Lag 4 0.9984665 0.9932279 1.0037326
## Lag 5 0.9957767 0.9905404 1.0010407
## Lag 6 0.9977365 0.9924994 1.0030014
## Lag 7 0.9943632 0.9891330 0.9996210

escribir riesgos

cardiores.p8 <- data.frame(pred$predvar, pred$allRRfit, pred$allRRlow, pred$allRRhigh)
#write_xlsx(cardiores.p8, "cardiores.p8.xlsx")
DT::datatable(cardiores.p8)

Residuales

resi_MODEL <- residuals(m1,type="response")
############gr?fico residuos
plot(resi_MODEL,ylim=c(-50,130),pch=19,cex=0.4,col=grey(0.6),
     main="Residuales en el tiempo",ylab="Residuales (observado-ajustado)",xlab="Días")
abline(h=1,lty=2,lwd=2)

#Pruebas de normalidad y graficos de normalidad

hist(resi_MODEL)

g2 = resi_MODEL |> hist(freq=FALSE,
                        main="Histograma de los residuales",
                        xlab="Residuales",
                        ylab="Densidad",
                        col ="blue")
resi_MODEL |> density() |> lines()

g1 = data.frame(resi_MODEL) |> ggplot(aes(sample=resi_MODEL))+
  stat_qq(size = 2) +
  stat_qq_line(distribution = stats::qnorm)+
  labs(x = "Cuantil teórico",
       y = "Residuales")+
  theme_minimal()
g1

#test para contrastar hipotesis de normalidad
# test de kolmogorov smirnov muestras mayores 50 datos

ks.test(resi_MODEL,pnorm)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resi_MODEL
## D = 0.34236, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 32.69, df = 2, p-value = 7.969e-08

shapiro.test(resi_MODEL)

## 
##  Shapiro-Wilk normality test
## 
## data:  resi_MODEL
## W = 0.99432, p-value = 1.938e-06

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 1.7514, p-value = 0.0001748

graficos de autocorrelacion y autocorrelacion parcial

layout(matrix(c(1,2),ncol=2,byrow=T))

resi_MODEL_serie=ts(resi_MODEL)
plot(acf(resi_MODEL_serie,plot = FALSE, lag.max = 30))
plot(pacf(resi_MODEL_serie,plot = FALSE, lag.max = 30))

# test-verificarrrrrrrrr 
#Prueba de Ljung-Box para verificar autocorrelación 
Box.test(resi_MODEL)

## 
##  Box-Pierce test
## 
## data:  resi_MODEL
## X-squared = 0.12588, df = 1, p-value = 0.7227

MODELO SO2 cardiovascular con 1 Rezago

maxlag <-c(1,7)
cb1_SO2 <- crossbasis(obs$SO2,maxlag,argvar1_SO2)
templ <- crossbasis(obs$temperatura,maxlag,templin)

Modelo

m1 <- glm(mcardio ~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
          data=obs, family=quasipoisson)

m_pois=glm(mcardio~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
           data=obs,family=poisson)

CALCULO QAIC

dfun <- function(object) {
  with(object,sum((weights * residuals^2)[weights > 0])/df.residual)
}
qAIC(m_pois,dispersion=dfun(m1))

## [1] 10867.37

Cálculo RR e IC 95%

(eff1 <- ci.lin(m1,subset="cb1_SO2",Exp=T))

##             Estimate       StdErr         z           P exp(Est.)      2.5%
## cb1_SO2 -0.001585641 0.0006131706 -2.585971 0.009710518 0.9984156 0.9972164
##             97.5%
## cb1_SO2 0.9996162

predicción

varcen <- round(mean(obs$SO2),0)

cp <- crosspred(cb1_SO2,m1,cen=varcen,by=0.1)
cen <- cp$predvar[which.min(cp$allRRfit)]

pred <- crosspred(cb1_SO2,m1,cen=cen, by=0.1)

DT::datatable(cbind(pred$allRRfit, pred$allRRlow,pred$allRRhigh))

grafico

xlab <- expression(paste("SO2 concentracion mcg/m3"))

par(mar=c(5,4,3,1),mgp=c(3,1,0),las=1,cex.axis=0.9,cex.lab=1)
plot(pred,"overall",col="red",ylim=c(0.9,1.5),axes=T,lab=c(6,5,7),xlab=xlab,
     ylab="RR",main="Mortalidad cardiovascular")

pred$predvar

##   [1]  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.0  2.1
##  [16]  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9  3.0  3.1  3.2  3.3  3.4  3.5  3.6
##  [31]  3.7  3.8  3.9  4.0  4.1  4.2  4.3  4.4  4.5  4.6  4.7  4.8  4.9  5.0  5.1
##  [46]  5.2  5.3  5.4  5.5  5.6  5.7  5.8  5.9  6.0  6.1  6.2  6.3  6.4  6.5  6.6
##  [61]  6.7  6.8  6.9  7.0  7.1  7.2  7.3  7.4  7.5  7.6  7.7  7.8  7.9  8.0  8.1
##  [76]  8.2  8.3  8.4  8.5  8.6  8.7  8.8  8.9  9.0  9.1  9.2  9.3  9.4  9.5  9.6
##  [91]  9.7  9.8  9.9 10.0 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 11.0 11.1
## [106] 11.2 11.3 11.4 11.5 11.6 11.7 11.8 11.9 12.0 12.1 12.2 12.3 12.4 12.5 12.6
## [121] 12.7 12.8 12.9 13.0 13.1 13.2 13.3 13.4 13.5 13.6 13.7 13.8 13.9 14.0 14.1
## [136] 14.2 14.3 14.4 14.5 14.6 14.7 14.8 14.9 15.0 15.1 15.2

plot(pred, "slices", type="p", pch=19, cex=1.5, lag=c(1,1,1),
     var=c(1.1,5, 15.2), ci="bars",
     ci.arg=list(col=grey(0.8)), 
     ylab="RR",main="Lag-specific effects")

plot(pred, "slices", var=c(1.1,5, 15.2), lag=c(1,1,1), col=4,
     ci.arg=list(density=40,col=grey(0.7)))

Tabla con estimaciones

tablag <- matrix(NA,7+1,3,dimnames=list(paste("Lag",0:7),
                                        c("RR","ci.low","ci.hi")))


# Correr el ciclo
for(i in 0:7) {
  # LAG de las variables i=0, 1 ,2, 3...7
  SO2lag=Lag(obs$SO2,i)
  temlag=Lag(obs$temperatura,i)  # el i puede fijarse 
  m1 <- glm(mcardio~ SO2lag+ diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(temlag),
            data=obs, family=quasipoisson)
  tablag[i+1,] <- ci.lin(m1,subset="SO2lag",Exp=T)[5:7]
}
tablag

##              RR    ci.low     ci.hi
## Lag 0 0.9988072 0.9935348 1.0041075
## Lag 1 0.9934055 0.9881862 0.9986524
## Lag 2 0.9942683 0.9890516 0.9995125
## Lag 3 0.9941522 0.9889385 0.9993933
## Lag 4 0.9984665 0.9932279 1.0037326
## Lag 5 0.9957767 0.9905404 1.0010407
## Lag 6 0.9977365 0.9924994 1.0030014
## Lag 7 0.9943632 0.9891330 0.9996210

escribir riesgos

cardiores.p8 <- data.frame(pred$predvar, pred$allRRfit, pred$allRRlow, pred$allRRhigh)
#write_xlsx(cardiores.p8, "cardiores.p8.xlsx")
DT::datatable(cardiores.p8)

Residuales

resi_MODEL <- residuals(m1,type="response")
############gr?fico residuos
plot(resi_MODEL,ylim=c(-50,130),pch=19,cex=0.4,col=grey(0.6),
     main="Residuales en el tiempo",ylab="Residuales (observado-ajustado)",xlab="Días")
abline(h=1,lty=2,lwd=2)

#Pruebas de normalidad###
#graficos de normalidad###

##histohgrama# 

hist(resi_MODEL)

g2 = resi_MODEL |> hist(freq=FALSE,
                        main="Histograma de los residuales",
                        xlab="Residuales",
                        ylab="Densidad",
                        col ="blue")
resi_MODEL |> density() |> lines()

g1 = data.frame(resi_MODEL) |> ggplot(aes(sample=resi_MODEL))+
  stat_qq(size = 2) +
  stat_qq_line(distribution = stats::qnorm)+
  labs(x = "Cuantil teórico",
       y = "Residuales")+
  theme_minimal()
g1

#test para contrastar hipotesis de normalidad
# test de kolmogorov smirnov muestras mayores 50 datos

ks.test(resi_MODEL,pnorm)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resi_MODEL
## D = 0.34236, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 32.69, df = 2, p-value = 7.969e-08

shapiro.test(resi_MODEL)

## 
##  Shapiro-Wilk normality test
## 
## data:  resi_MODEL
## W = 0.99432, p-value = 1.938e-06

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 1.7514, p-value = 0.0001748

layout(matrix(c(1,2),ncol=2,byrow=T))
#GRAFICO AUTOCORELACI?N MODELO
resi_MODEL_serie=ts(resi_MODEL)
plot(acf(resi_MODEL_serie,plot = FALSE, lag.max = 30))
plot(pacf(resi_MODEL_serie,plot = FALSE, lag.max = 30))

# test-verificarrrrrrrrr 
#Prueba de Ljung-Box para verificar autocorrelación 
Box.test(resi_MODEL)

## 
##  Box-Pierce test
## 
## data:  resi_MODEL
## X-squared = 0.12588, df = 1, p-value = 0.7227

MODELO SO2 cardiovascular con 2 Rezago

maxlag <-c(2,2)
cb1_SO2 <- crossbasis(obs$SO2,maxlag,argvar1_SO2)
templ <- crossbasis(obs$temperatura,maxlag,templin)

MODELO

m1 <- glm(mcardio ~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
          data=obs, family=quasipoisson)

m_pois=glm(mcardio~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
           data=obs,family=poisson)

Cálculo QAIC

dfun <- function(object) {
  with(object,sum((weights * residuals^2)[weights > 0])/df.residual)
}
qAIC(m_pois,dispersion=dfun(m1))

## [1] 10897.94

Cálculo RR e IC 95%

(eff1 <- ci.lin(m1,subset="cb1_SO2",Exp=T))

##             Estimate     StdErr        z          P exp(Est.)      2.5%
## cb1_SO2 -0.005748232 0.00268402 -2.14165 0.03222162 0.9942683 0.9890516
##             97.5%
## cb1_SO2 0.9995125

predicción

varcen <- round(mean(obs$SO2),0)

cp <- crosspred(cb1_SO2,m1,cen=varcen,by=0.1)
cen <- cp$predvar[which.min(cp$allRRfit)]

pred <- crosspred(cb1_SO2,m1,cen=cen, by=0.1)

a=cbind(pred$allRRfit, pred$allRRlow,pred$allRRhigh)
colnames(a)<-c("allRRfit", "allRRlow","allRRhigh")
DT::datatable(a)

grafico

xlab <- expression(paste("SO2 concentracion mcg/m3"))

par(mar=c(5,4,3,1),mgp=c(3,1,0),las=1,cex.axis=0.9,cex.lab=1)
plot(pred,"overall",col="red",ylim=c(0.9,1.5),axes=T,lab=c(6,5,7),xlab=xlab,
     ylab="RR",main="Mortalidad cardiovascular")

plot(pred, "slices", type="p", pch=19, cex=1.5, lag=c(2,2,2),
     var=c(1.1,5, 15.2), ci="bars",
     ci.arg=list(col=grey(0.9)), 
     ylab="RR",main="Lag-specific effects")

` ### Tabla con estimaciones

tablag <- matrix(NA,7+1,3,dimnames=list(paste("Lag",0:7),
                                        c("RR","ci.low","ci.hi")))


# Correr el ciclo
for(i in 0:7) {
  # LAG de las variables i=0, 1 ,2, 3...7
  SO2lag=Lag(obs$SO2,i)
  temlag=Lag(obs$temperatura,i)  # el i puede fijarse 
  m1 <- glm(mcardio~ SO2lag+ diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(temlag),
            data=obs, family=quasipoisson)
  tablag[i+1,] <- ci.lin(m1,subset="SO2lag",Exp=T)[5:7]
}
tablag

##              RR    ci.low     ci.hi
## Lag 0 0.9988072 0.9935348 1.0041075
## Lag 1 0.9934055 0.9881862 0.9986524
## Lag 2 0.9942683 0.9890516 0.9995125
## Lag 3 0.9941522 0.9889385 0.9993933
## Lag 4 0.9984665 0.9932279 1.0037326
## Lag 5 0.9957767 0.9905404 1.0010407
## Lag 6 0.9977365 0.9924994 1.0030014
## Lag 7 0.9943632 0.9891330 0.9996210

escribir riesgos

cardiores.p8 <- data.frame(pred$predvar, pred$allRRfit, pred$allRRlow, pred$allRRhigh)
#write_xlsx(cardiores.p8, "cardiores.p8.xlsx")
DT::datatable(cardiores.p8)

Residuales

resi_MODEL <- residuals(m1,type="response")

plot(resi_MODEL,ylim=c(-50,130),pch=19,cex=0.4,col=grey(0.6),
     main="Residuales en el tiempo",ylab="Residuales (observado-ajustado)",xlab="Días")
abline(h=1,lty=2,lwd=2)

Pruebas de normalidad y graficos de normalidad

hist(resi_MODEL)

g2 = resi_MODEL |> hist(freq=FALSE,
                        main="Histograma de los residuales",
                        xlab="Residuales",
                        ylab="Densidad",
                        col ="blue")
resi_MODEL |> density() |> lines()

g1 = data.frame(resi_MODEL) |> ggplot(aes(sample=resi_MODEL))+
  stat_qq(size = 2) +
  stat_qq_line(distribution = stats::qnorm)+
  labs(x = "Cuantil teórico",
       y = "Residuales")+
  theme_minimal()
g1

#test para contrastar hipotesis de normalidad
# test de kolmogorov smirnov muestras mayores 50 datos

ks.test(resi_MODEL,pnorm)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resi_MODEL
## D = 0.34236, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 32.69, df = 2, p-value = 7.969e-08

shapiro.test(resi_MODEL)

## 
##  Shapiro-Wilk normality test
## 
## data:  resi_MODEL
## W = 0.99432, p-value = 1.938e-06

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 1.7514, p-value = 0.0001748

Gráficos de autocorrelacion y autocorrelacion parcial

layout(matrix(c(1,2),ncol=2,byrow=T))
resi_MODEL_serie=ts(resi_MODEL)
plot(acf(resi_MODEL_serie,plot = FALSE, lag.max = 30))
plot(pacf(resi_MODEL_serie,plot = FALSE, lag.max = 30))

# test-verificarrrrrrrrr 
#Prueba de Ljung-Box para verificar autocorrelación 
Box.test(resi_MODEL)

## 
##  Box-Pierce test
## 
## data:  resi_MODEL
## X-squared = 0.12588, df = 1, p-value = 0.7227

MODELO SO2 cardiovascular con Promedio rezago 1-2

maxlag <-c(1,1)
cb1_SO2 <- crossbasis(obs$PSO2,maxlag,argvar1_SO2)
templ <- crossbasis(obs$temperatura,maxlag,templin)

Modelo

m1 <- glm(mcardio ~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
          data=obs, family=quasipoisson)

m_pois=glm(mcardio~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
           data=obs,family=poisson)

CALCULO QAIC

dfun <- function(object) {
  with(object,sum((weights * residuals^2)[weights > 0])/df.residual)
}
qAIC(m_pois,dispersion=dfun(m1))

## [1] 10893.76

Cálculo RR e IC 95%

(eff1 <- ci.lin(m1,subset="cb1_SO2",Exp=T))

##             Estimate      StdErr         z          P exp(Est.)      2.5%
## cb1_SO2 -0.005213656 0.003107266 -1.677891 0.09336829 0.9947999 0.9887599
##            97.5%
## cb1_SO2 1.000877

predicción

varcen <- round(mean(obs$SO2),0)

cp <- crosspred(cb1_SO2,m1,cen=varcen,by=0.1)
cen <- cp$predvar[which.min(cp$allRRfit)]

pred <- crosspred(cb1_SO2,m1,cen=cen, by=0.1)

DT::datatable(cbind(pred$allRRfit, pred$allRRlow,pred$allRRhigh))

grafico

xlab <- expression(paste("SO2 concentracion mcg/m3"))

par(mar=c(5,4,3,1),mgp=c(3,1,0),las=1,cex.axis=0.9,cex.lab=1)
plot(pred,"overall",col="red",ylim=c(0.9,1.5),axes=T,lab=c(6,5,7),xlab=xlab,
     ylab="RR",main="Mortalidad cardiovascular")

summary(pred)

## PREDICTIONS:
## values: 122 
## centered at: 12.9 
## range: 0.8 , 12.9 
## lag: 1 1 
## exponentiated: yes 
## cumulative: no 
## 
## MODEL:
## parameters: 1 
## class: glm lm 
## link: log

plot(pred, "slices", type="p", pch=19, cex=1.5, lag=c(1),
     var=c(1.1), ci="bars",
     ci.arg=list(col=grey(0.8)), 
     ylab="RR",main="Lag-specific effects")

Tabla con estimaciones

tablag <- matrix(NA,7+1,3,dimnames=list(paste("Lag",0:7),
                                        c("RR","ci.low","ci.hi")))


# Correr el ciclo
for(i in 0:7) {
  # LAG de las variables i=0, 1 ,2, 3...7
  SO2lag=Lag(obs$SO2,i)
  temlag=Lag(obs$temperatura,i)  # el i puede fijarse 
  m1 <- glm(mcardio~ SO2lag+ diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(temlag),
            data=obs, family=quasipoisson)
  tablag[i+1,] <- ci.lin(m1,subset="SO2lag",Exp=T)[5:7]
}
tablag

##              RR    ci.low     ci.hi
## Lag 0 0.9988072 0.9935348 1.0041075
## Lag 1 0.9934055 0.9881862 0.9986524
## Lag 2 0.9942683 0.9890516 0.9995125
## Lag 3 0.9941522 0.9889385 0.9993933
## Lag 4 0.9984665 0.9932279 1.0037326
## Lag 5 0.9957767 0.9905404 1.0010407
## Lag 6 0.9977365 0.9924994 1.0030014
## Lag 7 0.9943632 0.9891330 0.9996210

escribir riesgos

cardiores.p8 <- data.frame(pred$predvar, pred$allRRfit, pred$allRRlow, pred$allRRhigh)
#write_xlsx(cardiores.p8, "cardiores.p8.xlsx")
DT::datatable(cardiores.p8)

Residuales

resi_MODEL <- residuals(m1,type="response")
############gr?fico residuos
plot(resi_MODEL,ylim=c(-50,130),pch=19,cex=0.4,col=grey(0.6),
     main="Residuales en el tiempo",ylab="Residuales (observado-ajustado)",xlab="Días")
abline(h=1,lty=2,lwd=2)

#Pruebas de normalidad###
#graficos de normalidad###

##histohgrama# 

hist(resi_MODEL)

g2 = resi_MODEL |> hist(freq=FALSE,
                        main="Histograma de los residuales",
                        xlab="Residuales",
                        ylab="Densidad",
                        col ="blue")
resi_MODEL |> density() |> lines()

g1 = data.frame(resi_MODEL) |> ggplot(aes(sample=resi_MODEL))+
  stat_qq(size = 2) +
  stat_qq_line(distribution = stats::qnorm)+
  labs(x = "Cuantil teórico",
       y = "Residuales")+
  theme_minimal()
g1

#test para contrastar hipotesis de normalidad
# test de kolmogorov smirnov muestras mayores 50 datos

ks.test(resi_MODEL,pnorm)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resi_MODEL
## D = 0.34236, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 32.69, df = 2, p-value = 7.969e-08

shapiro.test(resi_MODEL)

## 
##  Shapiro-Wilk normality test
## 
## data:  resi_MODEL
## W = 0.99432, p-value = 1.938e-06

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 1.7514, p-value = 0.0001748

layout(matrix(c(1,2),ncol=2,byrow=T))
#GRAFICO AUTOCORELACI?N MODELO
resi_MODEL_serie=ts(resi_MODEL)
plot(acf(resi_MODEL_serie,plot = FALSE, lag.max = 30))
plot(pacf(resi_MODEL_serie,plot = FALSE, lag.max = 30))

# test-verificarrrrrrrrr 
#Prueba de Ljung-Box para verificar autocorrelación 
Box.test(resi_MODEL)

## 
##  Box-Pierce test
## 
## data:  resi_MODEL
## X-squared = 0.12588, df = 1, p-value = 0.7227

Modelos Cardiovascular

2023-11-03

MODELO SO2 cardiovascular sin Rezago

Construcion del MODELO

Cálculo QAIC

Cálculo RR e IC 95%

predicción

grafico

Tabla con estimaciones

escribir riesgos

Residuales

graficos de autocorrelacion y autocorrelacion parcial

MODELO SO2 cardiovascular con 1 Rezago

Modelo

CALCULO QAIC

Cálculo RR e IC 95%

predicción

grafico

Tabla con estimaciones

escribir riesgos

Residuales

MODELO SO2 cardiovascular con 2 Rezago

MODELO

Cálculo QAIC

Cálculo RR e IC 95%

predicción

grafico

escribir riesgos

Residuales

Pruebas de normalidad y graficos de normalidad

Gráficos de autocorrelacion y autocorrelacion parcial

MODELO SO2 cardiovascular con Promedio rezago 1-2

Modelo

CALCULO QAIC

Cálculo RR e IC 95%

predicción

grafico

Tabla con estimaciones

escribir riesgos

Residuales