-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy path03.firststage.R
152 lines (121 loc) · 5.04 KB
/
03.firststage.R
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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
###############################################################################
# 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
###############################################################################
####################################################################
# FIRST STAGE
# - DEFINE THE CROSS-BASIS MATRICES FOR THE 3 MODELS
# - BUILD OBJECTS TO STORE THE RESULTS
# - RUN THE POISSON TIME SERIES MODELS
# - REDUCE THE FITTED MAIN MODEL TO SUMMARIES
# - STORE THE RESULTS
# COMPUTING TIME IS ~40SEC (IN A 2.66GHz-4GBRAM PC WITH WINDOWS)
####################################################################
####################################################################
# DEFINE THE CROSS-BASIS MATRICES
# NB: THE USER CAN MODIFY THE CHOICES BELOW TO RUN ALTERNATIVE MODELS
# MAIN MODEL
# - PREDICTOR SPACE: QUADRATIC SPLINE WITH SPECIFIC KNOT SELECTION
# - LAG SPACE: NATURAL CUBIC SPLINE WITH DF AT EQUALLY-SPACED LOG-VALUES
lag <- c(0,21)
bound <- colMeans(ranges)
varknots <- equalknots(bound,fun="bs",degree=2,df=4)
lagknots <- logknots(21,df=5,int=T)
argvar <- list(fun="bs",degree=2,knots=varknots,bound=bound)
arglag <- list(fun="ns",knots=lagknots)
# ALTERNATIVE MODELS
# - IDENTICAL BASIS FOR PREDICTOR SPACE
# LAG SPACE: CONSTANT FOR LAG 0-3 AND LAG 0-21
lag2 <- c(0,3)
lag3 <- c(0,21)
arglag2 <- arglag3 <- list(fun="strata",df=1)
####################################################################
# BUILT OBJECTS WHERE RESULTS WILL BE STORED
# y- IS THE MATRIX FOR THE OUTCOME PARAMETERS
# S- IS THE LISTS OF (CO)VARIANCE MATRICES
# OVERALL CUMULATIVE SUMMARIES
yall <- matrix(NA,length(data),4,dimnames=list(regions,paste("b",seq(4),sep="")))
yall2 <- yall3 <- yall
# PREDICTOR-SPECIFIC SUMMARIES FOR MAIN MODEL
yhot <- matrix(NA,length(data),5,dimnames=list(regions,paste("b",seq(5),sep="")))
ycold <- matrix(NA,length(data),5,dimnames=list(regions,paste("b",seq(5),sep="")))
# (CO)VARIANCE MATRICES
Sall <- vector("list",length(data))
names(Sall) <- regions
Shot <- Scold <- Sall2 <- Sall3 <- Sall
# Q-AIC
qaic <- qaic2 <- qaic3 <- 0
####################################################################
# RUN THE MODEL FOR EACH CITY
# WARNING FOR PREDICTION BEYOND BOUNDARIES SUPPRESSED
options(warn=-1)
# LOOP FOR CITIES
for(i in seq(data)) {
# PRINT
cat(i,"")
# LOAD
sub <- data[[i]]
# DEFINE THE CROSS-BASES
cb <- crossbasis(sub$tmean,lag=lag,argvar=argvar,arglag=arglag)
cb2 <- crossbasis(sub$tmean,lag=lag2,argvar=argvar,arglag=arglag2)
cb3 <- crossbasis(sub$tmean,lag=lag3,argvar=argvar,arglag=arglag3)
# SET THE FIRST 21 RECORDS FOR cb2 AS MISSING
# THIS MAKES THE 3 MODELS COMPARABLE THROUGH AIC (SAME OBS)
cb2[0:21,] <- NA
# RUN THE FIRST-STAGE MODELS
mfirst <- glm(death ~ cb+dow+ns(time,df=10*14),family=quasipoisson(),sub)
mfirst2 <- glm(death ~ cb2+dow+ns(time,df=10*14),family=quasipoisson(),sub)
mfirst3 <- glm(death ~ cb3+dow+ns(time,df=10*14),family=quasipoisson(),sub)
####################################################################
# REDUCTION TO SUMMARY ASSOCIATIONS
# TO OVERALL CUMULATIVE SUMMARY
# NB: CENTERING NOT REALLY NEEDED HERE, AS COEF-VCOV (EXTRACTED BELOW) IN THE
# VAR SPACE DO NOT DEPEND ON CENTERING VALUE
crall <- crossreduce(cb,mfirst,cen=17)
crall2 <- crossreduce(cb2,mfirst2,cen=17)
crall3 <- crossreduce(cb3,mfirst3,cen=17)
# TO PREDICTOR-SPECIFIC SUMMARY FOR 22C AND 0C
# NB: CENTERING NEEDED HERE, AS COEF-VCOV (EXTRACTED BELOW) IN THE LAG SPACE
# DO DEPEND ON CENTERING VALUE
crhot <- crossreduce(cb,mfirst,type="var",value=22,cen=17)
crcold <- crossreduce(cb,mfirst,type="var",value=0,cen=17)
####################################################################
# STORE THE RESULTS
# OVERALL CUMULATIVE SUMMARY FOR THE MAIN MODEL
yall[i,] <- coef(crall)
Sall[[i]] <- vcov(crall)
# OVERALL CUMULATIVE SUMMARY FOR THE ALTERNATIVE MODELS
yall2[i,] <- coef(crall2)
yall3[i,] <- coef(crall3)
Sall2[[i]] <- vcov(crall2)
Sall3[[i]] <- vcov(crall3)
# PREDICTOR-SPECIFIC SUMMARY FOR 22C (MAIN MODEL)
yhot[i,] <- coef(crhot)
Shot[[i]] <- vcov(crhot)
# PREDICTOR-SPECIFIC SUMMARY FOR 0C (MAIN MODEL)
ycold[i,] <- coef(crcold)
Scold[[i]] <- vcov(crcold)
# Q-AIC
qaic[i] <- fqaic(mfirst)
qaic2[i] <- fqaic(mfirst2)
qaic3[i] <- fqaic(mfirst3)
}
####################################################################
# TEST: REDUCTION OF ALTERNATIVE MODELS TO THE SPACE OF THE PREDICTOR RETURNS
# THE SAME PARAMETERS APART FROM SCALING (SUMMED UPON 22 LAGS)
coef(crosspred(cb3,mfirst3,cen=17))
coef(crossreduce(cb3,mfirst3,cen=17))/22
# GRAND Q-AIC
sum(qaic) ; sum(qaic2) ; sum(qaic3)
# RESET WARNING
options(warn=0)
#