Shaping to RR (fitted line to RR)
Nasif Hossain
2025-02-11
Load Library
library(foreign)
library(ggplot2)
library(splines)
library(dlnm)
library(zoo)Read Dataset
data <- readRDS("Bangladeshu.rds")By using geom smooth
ggplot(data, aes(x=Rainfall_MS, y=AllCause))+
geom_smooth(method = "glm", formula = y~ns(x, df=3), se=TRUE)
same model with natural cubic spline function
model_1<-glm(AllCause ~ ns(Rainfall_MS, df=3),
family=quasipoisson(), data=data, na.action="na.exclude")summary(model_1)##
## Call:
## glm(formula = AllCause ~ ns(Rainfall_MS, df = 3), family = quasipoisson(),
## data = data, na.action = "na.exclude")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0825 0.0328 63.495 < 2e-16 ***
## ns(Rainfall_MS, df = 3)1 -0.7908 0.1252 -6.318 3.85e-10 ***
## ns(Rainfall_MS, df = 3)2 -0.5516 0.1229 -4.490 7.88e-06 ***
## ns(Rainfall_MS, df = 3)3 -0.3266 0.1531 -2.133 0.0331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 2.834402)
##
## Null deviance: 3610.0 on 1103 degrees of freedom
## Residual deviance: 3281.2 on 1100 degrees of freedom
## (356 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Predict model 1
dPred <- data.frame(
Rainfall_MS=seq(min(data$Rainfall_MS,na.rm=T),max(data$Rainfall_MS,na.rm=T),
length.out=nrow(data))
)
fitted<-predict(model_1,newdata=dPred,type="response")Visualization of predictive model
plot(data$Rainfall_MS,data$AllCause,pch=19,cex=0.5,col=grey(0.6),
main="Natural cubic spline model \n Precipitation only",
ylab="Daily number of Diarrhoea",
xlab="Precipitation")
lines(dPred$Rainfall_MS,fitted,lwd=2, col="blue")
Cross-basis function, nonlinear but not distruted lag
cb <- crossbasis(data$Rainfall_MS, lag=0,
argvar=list(fun="ns",df=3))Without adjusting seasonality model
To make relative risk plot we will use cross-basis with DLNM package
model with nonlinear function without adjusting seasonality
model_2<-glm(AllCause ~ cb,
family=quasipoisson(), data=data, na.action="na.exclude")summary(model_2)##
## Call:
## glm(formula = AllCause ~ cb, family = quasipoisson(), data = data,
## na.action = "na.exclude")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0825 0.0328 63.495 < 2e-16 ***
## cbv1.l1 -0.7908 0.1252 -6.318 3.85e-10 ***
## cbv2.l1 -0.5516 0.1229 -4.490 7.88e-06 ***
## cbv3.l1 -0.3266 0.1531 -2.133 0.0331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 2.834402)
##
## Null deviance: 3610.0 on 1103 degrees of freedom
## Residual deviance: 3281.2 on 1100 degrees of freedom
## (356 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Model prediction without reference (By default ref. at medium value of predictor)
pred <- crosspred(cb, model_2, by=1,
cumul=TRUE)## centering value unspecified. Automatically set to 1000Model prediction without reference with
pred <- crosspred(cb, model_2, by=1,
cen = as.numeric(names(which.min(pred$allRRfit))),
cumul=TRUE)Visualization of fitted model
## Without specify centering
plot(crosspred(cb, model_2, by=1),
"overall", ci="area", ci.arg=list(col=grey(0.8)),
xlab="Precipitation", ylab="RR", las=1, lwd=2
,col=4,main="Overall effect (Seasonality Unadjusted)")## centering value unspecified. Automatically set to 1000
## centering at MPR (Minimum Precipitation risk)
plot(crosspred(cb, cen=as.numeric(names(which.min(pred$allRRfit))), model_2, by=1),
"overall", ci="area", ci.arg=list(col=grey(0.8)),
xlab="Precipitation", ylab="RR", las=1, lwd=2
,col=4,main="Overall effect (Seasonality Unadjusted)")
Seasonality adjusted model
# same model with natural cubic spline function
model_3<-glm(AllCause ~ ns(Rainfall_MS, df=3)+ns(Time, df=3*7),
family=quasipoisson(), data=data, na.action="na.exclude")
summary(model_3)##
## Call:
## glm(formula = AllCause ~ ns(Rainfall_MS, df = 3) + ns(Time, df = 3 *
## 7), family = quasipoisson(), data = data, na.action = "na.exclude")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.801032 0.144733 12.444 < 2e-16 ***
## ns(Rainfall_MS, df = 3)1 0.193823 0.176523 1.098 0.27245
## ns(Rainfall_MS, df = 3)2 -0.008736 0.200559 -0.044 0.96527
## ns(Rainfall_MS, df = 3)3 -0.030919 0.210006 -0.147 0.88298
## ns(Time, df = 3 * 7)1 0.652058 0.155363 4.197 2.93e-05 ***
## ns(Time, df = 3 * 7)2 0.677218 0.219615 3.084 0.00210 **
## ns(Time, df = 3 * 7)3 0.205496 0.236533 0.869 0.38516
## ns(Time, df = 3 * 7)4 -0.421541 0.235065 -1.793 0.07321 .
## ns(Time, df = 3 * 7)5 -0.211278 0.245201 -0.862 0.38907
## ns(Time, df = 3 * 7)6 0.207646 0.216549 0.959 0.33783
## ns(Time, df = 3 * 7)7 0.096396 0.207597 0.464 0.64250
## ns(Time, df = 3 * 7)8 1.128703 0.197253 5.722 1.36e-08 ***
## ns(Time, df = 3 * 7)9 -0.095038 0.215580 -0.441 0.65941
## ns(Time, df = 3 * 7)10 -0.238636 0.233912 -1.020 0.30786
## ns(Time, df = 3 * 7)11 -0.198851 0.269421 -0.738 0.46063
## ns(Time, df = 3 * 7)12 -0.650349 0.247627 -2.626 0.00875 **
## ns(Time, df = 3 * 7)13 0.432069 0.232518 1.858 0.06341 .
## ns(Time, df = 3 * 7)14 -0.241806 0.230962 -1.047 0.29536
## ns(Time, df = 3 * 7)15 -0.118104 0.238181 -0.496 0.62009
## ns(Time, df = 3 * 7)16 -0.091789 0.243513 -0.377 0.70629
## ns(Time, df = 3 * 7)17 -0.313122 0.320517 -0.977 0.32882
## ns(Time, df = 3 * 7)18 -0.888079 0.295883 -3.001 0.00275 **
## ns(Time, df = 3 * 7)19 -1.316730 0.293456 -4.487 8.00e-06 ***
## ns(Time, df = 3 * 7)20 -0.192950 0.357272 -0.540 0.58926
## ns(Time, df = 3 * 7)21 -1.258974 0.243219 -5.176 2.70e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.872559)
##
## Null deviance: 3610.0 on 1103 degrees of freedom
## Residual deviance: 2275.5 on 1079 degrees of freedom
## (356 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5dPred <- data.frame(
Rainfall_MS=seq(min(data$Rainfall_MS,na.rm=T),max(data$Rainfall_MS,na.rm=T),
length.out=nrow(data)),
Time=seq(1,1104,length.out=nrow(data))
)
fitted<-predict(model_3,newdata=dPred,type="response")Visualization of fitted model
plot(data$Rainfall_MS,data$AllCause,pch=19,cex=0.5,col=grey(0.6),
main="Natural cubic spline model \n Seasonality adjusted",
ylab="Daily number of Diarrhoea",
xlab="Precipitation")
lines(dPred$Rainfall_MS,fitted,lwd=2, col="blue")
# model with nonlinear function with adjusting seasonality
model_4<-glm(AllCause ~ cb+ns(Time, df=3*7),
family=quasipoisson(), data=data, na.action="na.exclude")
summary(model_4)##
## Call:
## glm(formula = AllCause ~ cb + ns(Time, df = 3 * 7), family = quasipoisson(),
## data = data, na.action = "na.exclude")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.801032 0.144733 12.444 < 2e-16 ***
## cbv1.l1 0.193823 0.176523 1.098 0.27245
## cbv2.l1 -0.008736 0.200559 -0.044 0.96527
## cbv3.l1 -0.030919 0.210006 -0.147 0.88298
## ns(Time, df = 3 * 7)1 0.652058 0.155363 4.197 2.93e-05 ***
## ns(Time, df = 3 * 7)2 0.677218 0.219615 3.084 0.00210 **
## ns(Time, df = 3 * 7)3 0.205496 0.236533 0.869 0.38516
## ns(Time, df = 3 * 7)4 -0.421541 0.235065 -1.793 0.07321 .
## ns(Time, df = 3 * 7)5 -0.211278 0.245201 -0.862 0.38907
## ns(Time, df = 3 * 7)6 0.207646 0.216549 0.959 0.33783
## ns(Time, df = 3 * 7)7 0.096396 0.207597 0.464 0.64250
## ns(Time, df = 3 * 7)8 1.128703 0.197253 5.722 1.36e-08 ***
## ns(Time, df = 3 * 7)9 -0.095038 0.215580 -0.441 0.65941
## ns(Time, df = 3 * 7)10 -0.238636 0.233912 -1.020 0.30786
## ns(Time, df = 3 * 7)11 -0.198851 0.269421 -0.738 0.46063
## ns(Time, df = 3 * 7)12 -0.650349 0.247627 -2.626 0.00875 **
## ns(Time, df = 3 * 7)13 0.432069 0.232518 1.858 0.06341 .
## ns(Time, df = 3 * 7)14 -0.241806 0.230962 -1.047 0.29536
## ns(Time, df = 3 * 7)15 -0.118104 0.238181 -0.496 0.62009
## ns(Time, df = 3 * 7)16 -0.091789 0.243513 -0.377 0.70629
## ns(Time, df = 3 * 7)17 -0.313122 0.320517 -0.977 0.32882
## ns(Time, df = 3 * 7)18 -0.888079 0.295883 -3.001 0.00275 **
## ns(Time, df = 3 * 7)19 -1.316730 0.293456 -4.487 8.00e-06 ***
## ns(Time, df = 3 * 7)20 -0.192950 0.357272 -0.540 0.58926
## ns(Time, df = 3 * 7)21 -1.258974 0.243219 -5.176 2.70e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.872559)
##
## Null deviance: 3610.0 on 1103 degrees of freedom
## Residual deviance: 2275.5 on 1079 degrees of freedom
## (356 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Predicting fitted model
pred <- crosspred(cb, model_4, by=1,
cumul=TRUE)## centering value unspecified. Automatically set to 1000pred <- crosspred(cb, model_4, by=1,
cen = as.numeric(names(which.min(pred$allRRfit))),
cumul=TRUE)Visualization of fitted model
# Without specify centering
plot(crosspred(cb, model_4, by=1),
"overall", ci="area", ci.arg=list(col=grey(0.8)),
xlab="Precipitation", ylab="RR", las=1, lwd=2
,col=4,main="Overall effect (Seasonality adjusted)")## centering value unspecified. Automatically set to 1000
# centering at MPR (Minimum Precipitation risk)
plot(crosspred(cb, cen=as.numeric(names(which.min(pred$allRRfit))), model_4, by=1),
"overall", ci="area", ci.arg=list(col=grey(0.8)),
xlab="Precipitation", ylab="RR", las=1, lwd=2
,col=4,main="Overall effect (Seasonality adjusted)")