2012_gasparrini_StatMed_Rcodedata/01.mkdata.R at master · ehm-lab/2012_gasparrini_StatMed_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
################################################################################
# Updated version of the code for the analysis in:
#
#   "Multivariate meta-analysis for non-linear and other
#     multi-parameter associations"
#   Gasparrini, Armstrong and Kenward
#   Statistics in Medicine 2012
#   http://www.ag-myresearch.com/2012_gasparrini_statmed.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/2012_gasparrini_StatMed_Rcodedata
################################################################################

################################################################################
# NB: THE EXAMPLE IS DIFFERENT IF COMPARED TO THE PUBLICATION, AS THE ORIGINAL
#     DATA ARE NOT AVAILABLE ANY MORE THROUGH THE PACKAGE NMMAPSlite, NOW
#     ARCHIVED
################################################################################

# LOAD THE PACKAGE
library(dlnm) ; library(splines) ; library(xtable)

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

# LOAD THE DATA
data <- read.csv("regEngWales.csv",row.names=1)

####################################################################
# REGIONS

regions <- as.character(unique(data$regnames))

####################################################################
# LIST OF DATAFRAMES FOR 10 REGIONS

datalist <- lapply(regions, function(region) data[data$regnames==region,])
names(datalist) <- regions

####################################################################
# CITY-LEVEL META-PREDICTORS

lat <- c(54.84815,53.58832,53.72352,52.85539,52.53304,52.03734,51.50583,
  51.24213,51.05361,52.02615)
perclat <- round(quantile(lat,c(1,3)*0.25),1)

####################################################################
# ADDITIONAL INFO

m <- length(datalist)

# MOVING AVERAGE OF TMEAN OVER LAG 0-6
for(i in seq(datalist)) datalist[[i]]$tmean05 <-
  filter(datalist[[i]]$tmean,rep(1,6)/6,side=1)

# TEMPERATURE RANGES (FOR LAG 0-5)
ranges <- t(sapply(datalist,function(x) range(x$tmean05,na.rm=T)))

# COMPUTE 25TH-75TH PERCENTILES OF META-VARIABLES

# DEFINE THE AVERAGE RANGE, CENTERING POINT, DEGREE AND TYPE OF THE SPLINE
# (THESE PARAMETERS CAN BE CHANGED BY THE USER FOR ADDITIONAL ANALYSES)
cen <- 17
bound <- colMeans(ranges)
degree <- 2
type <- "bs"
df <- 6

# DEFINE THE KNOTS AT TEMPERATURE CORRESPONDING TO AVERAGE PERCENTILES
knotperc <- c(5,35,65,95)
knots <- rowMeans(sapply(datalist,function(x)
  quantile(x$tmean05,knotperc/100,na.rm=T)))

# SAVE
#save.image("data.RData")

#