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 cardiorespiratoria sin Rezago

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

Construcion del MODELO

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

m_pois=glm(cardioresp~ 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] 11402.05

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.000908956 0.002297405 -0.3956447 0.6923672 0.9990915 0.9946028
##          97.5%
## cb1_SO2 1.0036

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.013267 0.9492233 1.081632
## 0.8 1.013175 0.9495645 1.081047
## 0.9 1.013083 0.9499058 1.080462
## 1   1.012991 0.9502473 1.079877
## 1.1 1.012899 0.9505889 1.079293
## 1.2 1.012807 0.9509305 1.078709

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 cardiorespiratoria")

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(cardioresp~ 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.9990915 0.9946028 1.003600
## Lag 1 0.9962929 0.9918414 1.000764
## Lag 2 0.9969823 0.9925343 1.001450
## Lag 3 0.9965907 0.9921450 1.001056
## Lag 4 0.9987824 0.9943253 1.003259
## Lag 5 0.9978643 0.9934021 1.002347
## Lag 6 0.9983022 0.9938434 1.002781
## Lag 7 0.9958076 0.9913536 1.000282

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.36493, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 30.925, df = 2, p-value = 1.926e-07

shapiro.test(resi_MODEL)

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

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 2.0584, p-value = 3.093e-05

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.35199, df = 1, p-value = 0.553

MODELO SO2 cardiorespiratoria con 1 Rezago

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

Modelo

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

m_pois=glm(cardioresp~ 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] 11361.29

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.00106174 0.0005209633 -2.038032 0.04154669 0.9989388 0.9979194
##             97.5%
## cb1_SO2 0.9999593

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 cardiorespiratoria")

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(cardioresp~ 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.9990915 0.9946028 1.003600
## Lag 1 0.9962929 0.9918414 1.000764
## Lag 2 0.9969823 0.9925343 1.001450
## Lag 3 0.9965907 0.9921450 1.001056
## Lag 4 0.9987824 0.9943253 1.003259
## Lag 5 0.9978643 0.9934021 1.002347
## Lag 6 0.9983022 0.9938434 1.002781
## Lag 7 0.9958076 0.9913536 1.000282

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.36493, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 30.925, df = 2, p-value = 1.926e-07

shapiro.test(resi_MODEL)

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

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 2.0584, p-value = 3.093e-05

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.35199, df = 1, p-value = 0.553

MODELO SO2 cardiorespiratoria con 2 Rezago

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

MODELO

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

m_pois=glm(cardioresp~ 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] 11391.07

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.003022231 0.002281424 -1.324713 0.1852665 0.9969823 0.9925343
##           97.5%
## cb1_SO2 1.00145

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 cardiorespiratoria")

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(cardioresp~ 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.9990915 0.9946028 1.003600
## Lag 1 0.9962929 0.9918414 1.000764
## Lag 2 0.9969823 0.9925343 1.001450
## Lag 3 0.9965907 0.9921450 1.001056
## Lag 4 0.9987824 0.9943253 1.003259
## Lag 5 0.9978643 0.9934021 1.002347
## Lag 6 0.9983022 0.9938434 1.002781
## Lag 7 0.9958076 0.9913536 1.000282

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.36493, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 30.925, df = 2, p-value = 1.926e-07

shapiro.test(resi_MODEL)

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

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 2.0584, p-value = 3.093e-05

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.35199, df = 1, p-value = 0.553

MODELO SO2 cardiorespiratoria 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(cardioresp ~ cb1_SO2  + diasem + tipodia + ns(fecha,df=round(6*length(fecha)/365.25)) + ns(templ),
          data=obs, family=quasipoisson)

m_pois=glm(cardioresp~ 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] 11398.65

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.003116771 0.002640883 -1.180201 0.2379205 0.9968881 0.9917415
##            97.5%
## cb1_SO2 1.002061

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 cardiorespiratoria")

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(cardioresp~ 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.9990915 0.9946028 1.003600
## Lag 1 0.9962929 0.9918414 1.000764
## Lag 2 0.9969823 0.9925343 1.001450
## Lag 3 0.9965907 0.9921450 1.001056
## Lag 4 0.9987824 0.9943253 1.003259
## Lag 5 0.9978643 0.9934021 1.002347
## Lag 6 0.9983022 0.9938434 1.002781
## Lag 7 0.9958076 0.9913536 1.000282

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.36493, p-value < 2.2e-16
## alternative hypothesis: two-sided

jarque.bera.test(resi_MODEL)

## 
##  Jarque Bera Test
## 
## data:  resi_MODEL
## X-squared = 30.925, df = 2, p-value = 1.926e-07

shapiro.test(resi_MODEL)

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

ad.test(resi_MODEL)

## 
##  Anderson-Darling normality test
## 
## data:  resi_MODEL
## A = 2.0584, p-value = 3.093e-05

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.35199, df = 1, p-value = 0.553

Modelos Cardiorespiratorio

2023-11-02

MODELO SO2 cardiorespiratoria 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 cardiorespiratoria con 1 Rezago

Modelo

CALCULO QAIC

Cálculo RR e IC 95%

predicción

grafico

Tabla con estimaciones

escribir riesgos

Residuales

MODELO SO2 cardiorespiratoria 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 cardiorespiratoria con Promedio rezago 1-2

Modelo

CALCULO QAIC

Cálculo RR e IC 95%

predicción

grafico

Tabla con estimaciones

escribir riesgos

Residuales