2016_gasparrini_AJE_Rcodedata/00.prepdata.R at master · ehm-lab/2016_gasparrini_AJE_Rcodedata · GitHub

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
################################################################################
# Updated version of the R code for the analysis in:
#
#   "Changes in susceptibility to heat within the summer:
#     a multi-country analysis"
#   Antonio Gasparrini and collaborators
#   American Journal of Epidemiology - 2016
#   http://www.ag-myresearch.com/2016_gasparrini_aje.html
#
# Update: 05 December 2017
# * an updated version of this code, compatible with future versions of the
#   software, is available at:
#   https://github.com/gasparrini/2016_gasparrini_AJE_Rcodedata
################################################################################

# This code reproduces the analysis with the subset of data only including UK
# For this reason, the results are slightly different from those published

################################################################################
# PREPARE THE DATA
################################################################################

# LOAD THE PACKAGES
library(dlnm) ; library(mvmeta) ; library(splines)

# CHECK VERSION OF THE PACKAGE
if(packageVersion("dlnm")<"2.2.0")
  stop("update dlnm package to version >= 2.2.0")

# LOAD THE DATASET (INCLUDING THE 10 UK REGIONS ONLY)
orig <- read.csv("regEngWales19902012.csv")
orig$date <- as.Date(orig$date)

# ARRANGE THE DATA AS A LIST OF DATA SETS
regions <- as.character(unique(orig$regnames))
dlist <- lapply(regions,function(x) orig[orig$regnames==x,])
names(dlist) <- regions

# METADATA FOR LOCATIONS
cities <- data.frame(
  city = regions,
  cityname = c("North East","North West","Yorkshire & Humber","East Midlands",
    "West Midlands","East","London","South East","South West","Wales"),
  lat = c(55.00000,54.01670,53.56670,52.98000,52.47000,52.24000,
    51.50720,51.31670,50.96000,52.32044),
  long = c(-1.866700,-2.633300,-1.200000,-0.750000,-2.290000,0.410000,-0.127500,
    -0.500000,-3.220000,-3.732616)
)

# IDENTIFY THE 4 AVERAGE HOTTEST MONTHS IN EACH LOCATION
t(sapply(dlist,function(x)
  sort((1:12)[order(tapply(x$tmean,x$month,mean,na.rm=T))][9:12])))

# RESTRICT TO SUMMER AND CREATE DAY OF SEASON VARIABLE
dlist <- lapply(seq(dlist),function(i) {
  sub <- subset(dlist[[i]],month %in% 6:9)
  sub$dos <- sequence(tapply(sub$date,sub$year,length))
  sub
})

################################################################################

# PARAMETERS FOR THE EXPOSURE-RESPONSE FUNCTION
vararg <- list(fun="bs",degree=2)
varper <- c(75)

# PARAMETERS FOR THE LAG-RESPONSE FUNCTION
lag <- 10
lagnk <- 2
lagdf <- 4

# DEGREES OF FREEDOM FOR SEASONALITY
dfseas <- 4

# DEGREES OF FREEDOM/DECADE FOR TIME TREND
dftrend <- 1

# MODEL FORMULAE
formula <- death ~ cb + dow + ns(dos,df=dfseas):factor(year) +
  ns(date,df=round(length(unique(year))/dftrend/10))
formula1 <- death ~ cb + dow + ns(yday,df=dfseas):factor(year) +
  ns(date,df=round(length(unique(year))/dftrend/10)) + int1
formula2 <- death ~ cb + dow + ns(yday,df=dfseas):factor(year) +
  ns(date,df=round(length(unique(year))/dftrend/10)) + int2

#