forked from gasparrini/2013_gasparrini_BMCmrm_Rcodedata
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy path05.mainresults.R
More file actions
110 lines (88 loc) · 3.87 KB
/
Copy path05.mainresults.R
File metadata and controls
110 lines (88 loc) · 3.87 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
###############################################################################
# Updated version of the code for the analysis in:
#
# "Reducing and meta-analyzing estimates from distributed lag non-linear models"
# Gasparrini and Armstrong
# BMC Medical Research Methodology - 2013
# http://www.ag-myresearch.com/2013_gasparrini_bmcmrm.html
#
# Update: 05 December 2017
# * an updated version of this code, compatible with future versions of the
# software, is available at:
# http://www.ag-myresearch.com/2013_gasparrini_bmcmrm.html
###############################################################################
####################################################################
# OVERALL CUMULATIVE SUMMARY ASSOCIATION
# PLOT
pdf("figure2.pdf",height=5,width=13)
par(mar=c(5,4,1,1)+0.1,cex.axis=0.9,mgp=c(2.5,1,0))
layout(matrix(1:2,ncol=2))
plot(cpall,type="n",ylab="RR",ylim=c(.8,2),xlab="Temperature (C)")
for(i in seq(regall)) lines(regall[[i]],ptype="overall",col=grey(0.5),lty=2)
abline(h=1)
lines(cpall,col=2,lwd=2)
mtext("Main model: first-stage and pooled estimates",cex=1)
legend ("top",c("Pooled (with 95%CI)","First-stage region-specific"),
lty=c(1,2),lwd=1.5,col=c(2,grey(0.7)),bty="n",inset=0.1,cex=0.8)
plot(cpall,ylab="RR",col=2,lwd=2,ylim=c(.8,2),xlab="Temperature (C)")
lines(cpall2,col=3,lty=2,lwd=2)
lines(cpall3,col=4,lty=4,lwd=2)
mtext("Comparison of alternative models",cex=1)
legend ("top",c("B-spline of lag 0-21 (with 95%CI)","Constant of lag 0-3",
"Constant of lag 0-21"),lty=c(1,2,4),lwd=1.5,col=2:4,bty="n",inset=0.05,
cex=0.8,title="Function for the lag space:")
dev.off()
# POINT OF MINIMUM MORTALITY
cpall$predvar[which.min(cpall$allRRfit)]
round(sum(regEngWales$tmean<17.1)/nrow(regEngWales)*100,1)
# Q TEST AND I-SQUARE
(qall <- qtest(mvall))
round(((qall$Q-qall$df)/qall$Q)[1]*100,1)
(qall2 <- qtest(mvall2))
round(((qall2$Q-qall2$df)/qall2$Q)[1]*100,1)
(qall3 <- qtest(mvall3))
round(((qall3$Q-qall3$df)/qall3$Q)[1]*100,1)
####################################################################
# PREDICTOR-SPECIFIC SUMMARIES
# PLOT
pdf("figure3.pdf",height=5,width=13)
par(mar=c(5,4,1,1)+0.1,cex.axis=0.9,mgp=c(2.5,1,0))
layout(matrix(1:2,ncol=2))
plot(cphot,type="n",ylab="RR",ylim=c(.95,1.12),xlab="Lag")
for(i in seq(reghot)) lines(reghot[[i]],ptype="overall",col=grey(0.5),lty=2)
abline(h=1)
lines(cphot,col=2,lwd=2)
legend ("top",c("Pooled (with 95%CI)","First-stage region-specific"),
lty=c(1,2),lwd=1.5,col=c(2,grey(0.7)),bty="n",inset=0.1,cex=0.8)
mtext(text=paste("Predictor-specific summary for temperature = ",22,
"C",sep=""),cex=1)
plot(cpcold,type="n",ylab="RR",ylim=c(.95,1.12),xlab="Lag")
for(i in seq(regcold)) lines(regcold[[i]],ptype="overall",col=grey(0.5),lty=2)
abline(h=1)
lines(cpcold,col=2,lwd=2)
legend ("top",c("Pooled (with 95%CI)","First-stage region-specific"),
lty=c(1,2),lwd=1.5,col=c(2,grey(0.7)),bty="n",inset=0.1,cex=0.8)
mtext(text=paste("Predictor-specific summary for temperature = ",0,
"C",sep=""),cex=1)
dev.off()
# OVERALL EFFECTS AT THESE TWO PREDICTOR LEVELS
round(with(cpall,cbind(allRRfit,allRRlow,allRRhigh)["22",]),3)
round(with(cpall,cbind(allRRfit,allRRlow,allRRhigh)["0",]),3)
# TESTS AND STATISTICS
(qhot <- qtest(mvhot))
(qcold <- qtest(mvcold))
round(((qhot$Q-qhot$df)/qhot$Q)[1]*100,1)
round(((qcold$Q-qcold$df)/qcold$Q)[1]*100,1)
####################################################################
# COMPARISON OF RANDOM VS. FIXED EFFECT MODEL FOR SUMMARY AT 22C
mvhot2 <- update(mvhot,method="fixed")
coef(mvhot); coef(mvhot2)
cphot2 <- crosspred(blag,coef=coef(mvhot2),vcov=vcov(mvhot2),
model.link="log",at=0:210/10)
plot(cphot,ci="lines",ylab="RR",ylim=c(.95,1.12),xlab="Lag",col=2)
lines(cphot2,ci="lines",col=4)
mtext(text=paste("Predictor-specific summary for temperature = ",22,
"C",sep=""),cex=1)
legend ("top",c("Random-effects model","Fixed-effects model"),
lty=1,col=c(2,4),bty="n",inset=0.1,cex=0.8)
#