london dataset
As an example for this vignette, we will use the built-in dataset
london which contains daily temperature and mortality data
from 2000 to 2012.
head(london)
#> date year dow tmean mort_00_74 mort_75plus mort
#> 1 2000-01-01 2000 Sat 7.642016 101 226 327
#> 2 2000-01-02 2000 Sun 8.351041 100 191 291
#> 3 2000-01-03 2000 Mon 8.188446 107 205 312
#> 4 2000-01-04 2000 Tue 6.320066 97 191 288
#> 5 2000-01-05 2000 Wed 6.550981 92 213 305
#> 6 2000-01-06 2000 Thu 8.728774 82 213 295
Using Bayesian Distributed Non-Linear Models (B-DLNM) we will
estimate the lagged effect of the temperature exposure to mortality in
the population with more than 75 years old, taking 21 days of lags, in
the city of London (2000-2012).
Set DLNM parameters
We will start by building the crossbasis with the dlnm
package.
We will use a natural spline with three knots for the
exposure–response function, placed at the 10th, 75th, and 90th
percentiles of daily mean temperature, and a natural spline for the
lag–response function, with knots equally spaced on the log scale up to
a maximum lag of 21 days:
# Exposure-response and lag-response spline parameters
dlnm_var <- list(
var_prc = c(10, 75, 90),
var_fun = "ns",
lag_fun = "ns",
max_lag = 21,
lagnk = 3
)
# Cross-basis parameters
argvar <- list(fun = dlnm_var$var_fun,
knots = quantile(london$tmean,
dlnm_var$var_prc/100, na.rm = TRUE),
Bound = range(london$tmean, na.rm = TRUE))
arglag <- list(fun = dlnm_var$lag_fun,
knots = logknots(dlnm_var$max_lag, nk = dlnm_var$lagnk))
# Create crossbasis
cb <- crossbasis(london$tmean, lag = dlnm_var$max_lag, argvar, arglag)
Let’s define the seasonality of the mortality time series using a
natural spline with 8 degrees of freedom per year, which flexibly
controls for long-term and seasonal trends in mortality:
# Seasonality of mortality time series
seas <- ns(london$date, df = round(8 * length(london$date) / 365.25))
Finally, we also define the temperature values for which predictions
will be generated. These correspond to an evenly spaced grid with a
0.1ºC step covering the full range of observed temperatures in the
data:
# Prediction values (equidistant points)
temp <- round(seq(min(london$tmean), max(london$tmean), by = 0.1), 1)
Fit the model and posterior samples
Fit the Bayesian DLNM using the function bdlnm(). In
this example, we specify a Poisson model including the cross-basis term,
a seasonal component, and an indicator for the day of the week. We draw
1000 samples from the posterior distribution and set a seed for
reproducibility.
mod <- bdlnm(mort_75plus ~ cb + factor(dow) + seas, data = london, family = "poisson", sample.arg = list(n = 1000, seed = 5243))
#> Warning in inla.posterior.sample(n, rfake, intern = intern, use.improved.mean =
#> use.improved.mean, : Since 'seed!=0', parallel model is disabled and serial
#> model is selected, num.threads='1:1'
The returned object is of class bdlnm and contains the fitted INLA
model, the crossbasis objects included in the formula, and posterior
samples of the estimated coefficients together with their summaries
(mean, standard deviation, and quantiles).
str(mod, max.level = 1)
#> List of 4
#> $ model :List of 48
#> ..- attr(*, "class")= chr "inla"
#> $ basis :List of 1
#> $ coefficients : num [1:123, 1:1000] 5.31646 -0.1515 -0.00617 -0.05908 0.05625 ...
#> ..- attr(*, "dimnames")=List of 2
#> $ coefficients.summary: num [1:123, 1:6] 5.20642 -0.13531 -0.00608 -0.05008 0.0505 ...
#> ..- attr(*, "dimnames")=List of 2
#> - attr(*, "n_sim")= num 1000
#> - attr(*, "class")= chr "bdlnm"
Minimum mortality temperature
We estimate the minimum mortality temperature (MMT), defined as the
temperature at which the overall mortality risk is minimized. This
optimal value will later be used to center the estimated relative
risks.
mmt <- optimal_exposure(mod, exp_at = temp)
#> Registered S3 method overwritten by 'crs':
#> method from
#> predict.gsl.bs np
str(mmt)
#> List of 2
#> $ est : Named num [1:1000] 18.8 18.9 18.8 18.6 19 19 18.8 18.9 19.1 19 ...
#> ..- attr(*, "names")= chr [1:1000] "sample1" "sample2" "sample3" "sample4" ...
#> $ summary: Named num [1:6] 18.896 0.139 18.6 18.9 19.2 ...
#> ..- attr(*, "names")= chr [1:6] "mean" "sd" "0.025quant" "0.5quant" ...
#> - attr(*, "exp_at")= num [1:326] -3.4 -3.3 -3.2 -3.1 -3 -2.9 -2.8 -2.7 -2.6 -2.5 ...
#> - attr(*, "lag_at")= num [1:22] 0 1 2 3 4 5 6 7 8 9 ...
#> - attr(*, "which")= chr "min"
#> - attr(*, "class")= chr "optimal_exposure"
We can visualize the posterior distribution of the MMT estimates. It
is useful to inspect this distribution to assess whether multiple
candidate optimal exposure values exist and to verify that the median
provides a reasonable centering value:
plot(mmt, xlab = "Temperature (ºC)", main = paste0("MMT (Median = ", round(mmt$summary[["0.5quant"]], 1), "ºC)"))
To make the predictions, we will center the risk at the median of
these values:
cen <- mmt$summary[["0.5quant"]]
cen
#> [1] 18.9
Predict exposure-lag-response effects
Predict the exposure–lag–response association between temperature and
mortality from the fitted model at the supplied temperature grid,
centering the predictions at the MMT value:
cpred <- bcrosspred(mod, exp_at = temp, cen = cen)
A bcrosspred object is returned containing posterior
samples of the estimated exposure–lag–response association for the
supplied temperature values and lags, together with their summaries.
str(cpred)
#> List of 16
#> $ exp_at : num [1:326] -3.4 -3.3 -3.2 -3.1 -3 -2.9 -2.8 -2.7 -2.6 -2.5 ...
#> $ lag_at : num [1:22] 0 1 2 3 4 5 6 7 8 9 ...
#> $ cen : num 18.9
#> $ coefficients : num [1:20, 1:1000] -0.1515 -0.00617 -0.05908 0.05625 -0.04859 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:20] "cbv1.l1" "cbv1.l2" "cbv1.l3" "cbv1.l4" ...
#> .. ..$ : chr [1:1000] "sample1" "sample2" "sample3" "sample4" ...
#> $ matfit : num [1:326, 1:22, 1:1000] -0.176 -0.175 -0.174 -0.173 -0.172 ...
#> ..- attr(*, "dimnames")=List of 3
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:22] "lag0" "lag1" "lag2" "lag3" ...
#> .. ..$ : chr [1:1000] "sample1" "sample2" "sample3" "sample4" ...
#> $ allfit : num [1:326, 1:1000] 0.659 0.654 0.649 0.644 0.639 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:1000] "sample1" "sample2" "sample3" "sample4" ...
#> $ matRRfit : num [1:326, 1:22, 1:1000] 0.839 0.84 0.841 0.841 0.842 ...
#> ..- attr(*, "dimnames")=List of 3
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:22] "lag0" "lag1" "lag2" "lag3" ...
#> .. ..$ : chr [1:1000] "sample1" "sample2" "sample3" "sample4" ...
#> $ allRRfit : num [1:326, 1:1000] 1.93 1.92 1.91 1.9 1.9 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:1000] "sample1" "sample2" "sample3" "sample4" ...
#> $ coefficients.summary: num [1:20, 1:6] -0.13531 -0.00608 -0.05008 0.0505 -0.04918 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:20] "cbv1.l1" "cbv1.l2" "cbv1.l3" "cbv1.l4" ...
#> .. ..$ : chr [1:6] "mean" "sd" "0.025quant" "0.5quant" ...
#> $ matfit.summary : num [1:326, 1:22, 1:6] -0.161 -0.16 -0.159 -0.158 -0.157 ...
#> ..- attr(*, "dimnames")=List of 3
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:22] "lag0" "lag1" "lag2" "lag3" ...
#> .. ..$ : chr [1:6] "mean" "sd" "0.025quant" "0.5quant" ...
#> $ allfit.summary : num [1:326, 1:6] 0.606 0.602 0.598 0.594 0.59 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:6] "mean" "sd" "0.025quant" "0.5quant" ...
#> $ matRRfit.summary : num [1:326, 1:22, 1:6] 0.852 0.852 0.853 0.854 0.854 ...
#> ..- attr(*, "dimnames")=List of 3
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:22] "lag0" "lag1" "lag2" "lag3" ...
#> .. ..$ : chr [1:6] "mean" "sd" "0.025quant" "0.5quant" ...
#> $ allRRfit.summary : num [1:326, 1:6] 1.83 1.83 1.82 1.81 1.8 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:326] "-3.4" "-3.3" "-3.2" "-3.1" ...
#> .. ..$ : chr [1:6] "mean" "sd" "0.025quant" "0.5quant" ...
#> $ ci.level : num 0.95
#> $ model.class : chr "inla"
#> $ model.link : chr "log"
#> - attr(*, "class")= chr "bcrosspred"
For instance, the estimated crossbasis coefficients are
stored as:
cpred$coefficients |>
head(c(5, 5))
#> sample1 sample2 sample3 sample4 sample5
#> cbv1.l1 -0.151502754 -0.140536633 -0.118854832 -0.123795275 -0.13058717
#> cbv1.l2 -0.006173712 -0.003457928 0.003135457 -0.007956683 -0.01566591
#> cbv1.l3 -0.059082316 -0.050139201 -0.049991893 -0.045513488 -0.03690102
#> cbv1.l4 0.056250819 0.054350516 0.039827408 0.047954664 0.04874636
#> cbv1.l5 -0.048591640 -0.042215543 -0.047502789 -0.049240365 -0.05852420
The relative risks (RR) for each temperature–lag combination are also
stored in an array for each posterior sample. For example, for the first
sample:
cpred$matRRfit[, , "sample1"] |>
head()
#> lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7
#> -3.4 0.8389132 1.080854 1.149670 1.101739 1.062913 1.041356 1.031633 1.029007
#> -3.3 0.8397120 1.080037 1.148395 1.100927 1.062452 1.041088 1.031456 1.028858
#> -3.2 0.8405115 1.079221 1.147122 1.100116 1.061991 1.040821 1.031278 1.028709
#> -3.1 0.8413116 1.078406 1.145851 1.099305 1.061531 1.040554 1.031101 1.028559
#> -3 0.8421122 1.077592 1.144581 1.098495 1.061071 1.040287 1.030923 1.028410
#> -2.9 0.8429132 1.076779 1.143314 1.097687 1.060611 1.040020 1.030746 1.028261
#> lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15
#> -3.4 1.029047 1.029096 1.028974 1.028697 1.028278 1.027733 1.027074 1.026317
#> -3.3 1.028901 1.028951 1.028830 1.028552 1.028132 1.027585 1.026924 1.026164
#> -3.2 1.028756 1.028807 1.028686 1.028408 1.027986 1.027437 1.026773 1.026010
#> -3.1 1.028610 1.028662 1.028542 1.028263 1.027840 1.027289 1.026623 1.025856
#> -3 1.028465 1.028518 1.028398 1.028118 1.027694 1.027141 1.026472 1.025703
#> -2.9 1.028319 1.028374 1.028254 1.027974 1.027549 1.026993 1.026322 1.025549
#> lag16 lag17 lag18 lag19 lag20 lag21
#> -3.4 1.025476 1.024566 1.023599 1.022592 1.021557 1.020509
#> -3.3 1.025319 1.024404 1.023434 1.022422 1.021382 1.020330
#> -3.2 1.025162 1.024243 1.023268 1.022252 1.021208 1.020151
#> -3.1 1.025005 1.024082 1.023103 1.022082 1.021034 1.019972
#> -3 1.024847 1.023921 1.022938 1.021912 1.020859 1.019793
#> -2.9 1.024690 1.023760 1.022773 1.021743 1.020685 1.019614
The overall cumulative effects for each temperature (summing across
all lags) are stored in:
cpred$allRRfit |>
head(c(5, 5))
#> sample1 sample2 sample3 sample4 sample5
#> -3.4 1.932018 1.743056 1.709815 1.764002 1.920655
#> -3.3 1.922786 1.737823 1.704222 1.756976 1.911513
#> -3.2 1.913600 1.732606 1.698648 1.749979 1.902417
#> -3.1 1.904461 1.727406 1.693093 1.743012 1.893366
#> -3 1.895368 1.722223 1.687558 1.736076 1.884363
Summaries of these effects across posterior samples are also
available:
cpred$matRRfit.summary |>
head(5)
#> , , mean
#>
#> lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7
#> -3.4 0.8515053 1.071324 1.135081 1.094447 1.060246 1.040318 1.030335 1.026475
#> -3.3 0.8522192 1.070717 1.134071 1.093762 1.059822 1.040047 1.030143 1.026317
#> -3.2 0.8529336 1.070110 1.133062 1.093078 1.059400 1.039777 1.029951 1.026159
#> -3.1 0.8536486 1.069504 1.132054 1.092394 1.058977 1.039506 1.029759 1.026001
#> -3 0.8543639 1.068899 1.131047 1.091712 1.058555 1.039236 1.029568 1.025843
#> lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15
#> -3.4 1.025160 1.024241 1.023546 1.023055 1.022750 1.022612 1.022622 1.022762
#> -3.3 1.025016 1.024108 1.023419 1.022933 1.022629 1.022490 1.022496 1.022631
#> -3.2 1.024873 1.023975 1.023293 1.022810 1.022508 1.022367 1.022371 1.022500
#> -3.1 1.024730 1.023842 1.023167 1.022688 1.022387 1.022245 1.022245 1.022370
#> -3 1.024587 1.023709 1.023041 1.022566 1.022266 1.022123 1.022120 1.022239
#> lag16 lag17 lag18 lag19 lag20 lag21
#> -3.4 1.023013 1.023357 1.023776 1.024250 1.024762 1.025292
#> -3.3 1.022876 1.023212 1.023621 1.024085 1.024586 1.025106
#> -3.2 1.022738 1.023066 1.023466 1.023920 1.024411 1.024920
#> -3.1 1.022601 1.022920 1.023311 1.023756 1.024236 1.024734
#> -3 1.022463 1.022775 1.023157 1.023591 1.024061 1.024548
#>
#> , , sd
#>
#> lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7
#> -3.4 1.023770 1.011110 1.012904 1.006822 1.006101 1.006757 1.006194 1.004886
#> -3.3 1.023481 1.010968 1.012745 1.006741 1.006027 1.006674 1.006117 1.004826
#> -3.2 1.023193 1.010827 1.012587 1.006660 1.005954 1.006592 1.006041 1.004765
#> -3.1 1.022907 1.010687 1.012430 1.006580 1.005882 1.006510 1.005966 1.004705
#> -3 1.022624 1.010548 1.012275 1.006500 1.005810 1.006429 1.005891 1.004646
#> lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15
#> -3.4 1.003892 1.003676 1.003872 1.004111 1.004224 1.004174 1.003998 1.003787
#> -3.3 1.003843 1.003629 1.003822 1.004057 1.004169 1.004120 1.003946 1.003738
#> -3.2 1.003794 1.003582 1.003772 1.004004 1.004114 1.004065 1.003894 1.003688
#> -3.1 1.003745 1.003535 1.003722 1.003951 1.004059 1.004012 1.003842 1.003639
#> -3 1.003697 1.003489 1.003673 1.003898 1.004005 1.003958 1.003791 1.003591
#> lag16 lag17 lag18 lag19 lag20 lag21
#> -3.4 1.003684 1.003854 1.004390 1.005268 1.006395 1.007678
#> -3.3 1.003637 1.003805 1.004336 1.005203 1.006316 1.007584
#> -3.2 1.003589 1.003756 1.004281 1.005138 1.006238 1.007490
#> -3.1 1.003542 1.003707 1.004227 1.005074 1.006161 1.007397
#> -3 1.003495 1.003659 1.004173 1.005010 1.006084 1.007305
#>
#> , , 0.025quant
#>
#> lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7
#> -3.4 0.8125568 1.048865 1.106934 1.080210 1.047614 1.026364 1.018071 1.017190
#> -3.3 0.8136784 1.048547 1.106319 1.079719 1.047326 1.026303 1.018106 1.017123
#> -3.2 0.8148014 1.048200 1.105633 1.079262 1.047032 1.026264 1.018110 1.017036
#> -3.1 0.8159837 1.047860 1.105148 1.078775 1.046741 1.026224 1.018052 1.016991
#> -3 0.8170436 1.047518 1.104522 1.078292 1.046450 1.026105 1.017994 1.016948
#> lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15
#> -3.4 1.017417 1.017364 1.016103 1.015198 1.014824 1.014525 1.014577 1.015107
#> -3.3 1.017413 1.017320 1.016122 1.015190 1.014827 1.014487 1.014496 1.015053
#> -3.2 1.017359 1.017236 1.016142 1.015154 1.014830 1.014449 1.014465 1.014998
#> -3.1 1.017277 1.017152 1.016162 1.015118 1.014798 1.014459 1.014461 1.014944
#> -3 1.017237 1.017102 1.016183 1.015083 1.014738 1.014402 1.014459 1.014882
#> lag16 lag17 lag18 lag19 lag20 lag21
#> -3.4 1.015333 1.015670 1.015251 1.014467 1.012532 1.010840
#> -3.3 1.015264 1.015610 1.015211 1.014373 1.012391 1.010867
#> -3.2 1.015195 1.015547 1.015170 1.014305 1.012389 1.010893
#> -3.1 1.015125 1.015482 1.015137 1.014197 1.012391 1.010876
#> -3 1.015057 1.015438 1.015109 1.014088 1.012393 1.010823
#>
#> , , 0.5quant
#>
#> lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7
#> -3.4 0.8520207 1.071217 1.135277 1.094285 1.060433 1.040508 1.030243 1.026470
#> -3.3 0.8527836 1.070636 1.134250 1.093583 1.060016 1.040202 1.030043 1.026297
#> -3.2 0.8536028 1.070016 1.133184 1.092943 1.059589 1.039950 1.029871 1.026126
#> -3.1 0.8543231 1.069410 1.132108 1.092279 1.059166 1.039667 1.029659 1.025972
#> -3 0.8550176 1.068765 1.131050 1.091647 1.058772 1.039345 1.029461 1.025826
#> lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15
#> -3.4 1.025090 1.024037 1.023273 1.022783 1.022499 1.022398 1.022438 1.022718
#> -3.3 1.024961 1.023928 1.023138 1.022672 1.022392 1.022262 1.022325 1.022608
#> -3.2 1.024808 1.023775 1.023006 1.022562 1.022265 1.022154 1.022215 1.022488
#> -3.1 1.024676 1.023655 1.022879 1.022450 1.022164 1.022042 1.022106 1.022364
#> -3 1.024555 1.023542 1.022766 1.022352 1.022054 1.021902 1.021997 1.022237
#> lag16 lag17 lag18 lag19 lag20 lag21
#> -3.4 1.022959 1.023210 1.023706 1.024230 1.024781 1.025265
#> -3.3 1.022817 1.023044 1.023552 1.024094 1.024656 1.025108
#> -3.2 1.022722 1.022899 1.023390 1.023945 1.024443 1.024913
#> -3.1 1.022623 1.022755 1.023224 1.023754 1.024287 1.024670
#> -3 1.022510 1.022627 1.023077 1.023589 1.024119 1.024431
#>
#> , , 0.975quant
#>
#> lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7
#> -3.4 0.8913229 1.095721 1.163639 1.109040 1.072839 1.053700 1.042297 1.036149
#> -3.3 0.8915669 1.094826 1.162214 1.108176 1.072251 1.053347 1.041952 1.035846
#> -3.2 0.8918072 1.093926 1.160791 1.107293 1.071617 1.052917 1.041607 1.035544
#> -3.1 0.8920616 1.093028 1.159367 1.106449 1.071040 1.052426 1.041263 1.035242
#> -3 0.8923106 1.092130 1.157938 1.105605 1.070464 1.051936 1.040919 1.034946
#> lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15
#> -3.4 1.033459 1.032345 1.031874 1.031738 1.031698 1.031445 1.030900 1.030684
#> -3.3 1.033195 1.032127 1.031645 1.031499 1.031491 1.031163 1.030673 1.030433
#> -3.2 1.032930 1.031885 1.031417 1.031261 1.031292 1.030939 1.030447 1.030182
#> -3.1 1.032665 1.031642 1.031189 1.031022 1.031093 1.030749 1.030232 1.029954
#> -3 1.032400 1.031399 1.030961 1.030782 1.030810 1.030559 1.030018 1.029736
#> lag16 lag17 lag18 lag19 lag20 lag21
#> -3.4 1.030421 1.031382 1.032987 1.035557 1.038203 1.041057
#> -3.3 1.030197 1.031142 1.032759 1.035270 1.037851 1.040589
#> -3.2 1.029966 1.030928 1.032530 1.034963 1.037499 1.040210
#> -3.1 1.029750 1.030720 1.032294 1.034632 1.037146 1.039785
#> -3 1.029528 1.030438 1.032002 1.034302 1.036793 1.039466
#>
#> , , mode
#>
#> lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7
#> -3.4 0.8543850 1.066869 1.135865 1.093911 1.061136 1.040744 1.029982 1.027207
#> -3.3 0.8550927 1.066353 1.134937 1.093268 1.060730 1.040435 1.029704 1.027061
#> -3.2 0.8558087 1.065664 1.133798 1.092625 1.060321 1.040224 1.029519 1.026840
#> -3.1 0.8565227 1.065152 1.132666 1.091983 1.059911 1.039914 1.029335 1.026695
#> -3 0.8572315 1.064633 1.131540 1.091445 1.059498 1.039605 1.029151 1.026549
#> lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15
#> -3.4 1.025023 1.022489 1.022406 1.022133 1.021831 1.021776 1.022224 1.022924
#> -3.3 1.024867 1.022409 1.022374 1.022147 1.021783 1.021725 1.022096 1.022845
#> -3.2 1.024769 1.022271 1.022341 1.022091 1.021799 1.021602 1.022031 1.022762
#> -3.1 1.024613 1.022190 1.022306 1.022095 1.021742 1.021545 1.021900 1.022619
#> -3 1.024456 1.022057 1.022269 1.022031 1.021682 1.021484 1.021828 1.022530
#> lag16 lag17 lag18 lag19 lag20 lag21
#> -3.4 1.023118 1.022770 1.023754 1.024518 1.025406 1.027223
#> -3.3 1.022939 1.022630 1.023597 1.024349 1.025161 1.027054
#> -3.2 1.022814 1.022489 1.023440 1.024241 1.024915 1.026882
#> -3.1 1.022689 1.022349 1.023283 1.024050 1.024669 1.026588
#> -3 1.022513 1.022209 1.023125 1.023860 1.024424 1.026411
cpred$allRRfit.summary |>
head(5)
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> -3.4 1.833955 1.065134 1.617029 1.831841 2.077265 1.835723
#> -3.3 1.826291 1.064254 1.612914 1.824431 2.064625 1.828870
#> -3.2 1.818660 1.063380 1.608809 1.817566 2.052884 1.822000
#> -3.1 1.811063 1.062512 1.604714 1.810012 2.041214 1.815474
#> -3 1.803501 1.061652 1.600631 1.802157 2.029614 1.807719
Plot
Several visualizations can be produced from these predictions using
the plot() method.
For example, we can plot the overall cumulative effect for each
temperature value (suming across all lags):
plot(cpred, "overall", xlab = "Temperature (ºC)", ylab = "Relative Risk", col = 4, main="Overall", log = "y")
Alternatively, we can display the posterior sampling variability by
plotting the curves from all posterior samples instead of the credible
interval:
plot(cpred, "overall", xlab = "Temperature (ºC)", ylab = "Relative Risk", col = 4, main="Overall", log = "y", ci = "sampling")
We can also visualize the full exposure–lag–response association
using a 3-D surface:
plot(cpred, "3d", zlab = "Relative risk", col = 4, lphi = 60, cex.axis = 0.6, xlab = "Temperature (ºC)", main = "3D graph of temperature effect")
A contour plot provides a 2-D projection of the same
relationship:
plot(cpred, "contour", xlab = "Temperature (ºC)", ylab = "Lag", main = "Contour plot")
We can also examine the lag–response association for a high
temperature (e.g., the 99th percentile):
htemp <- round(quantile(london$tmean, 0.99), 1)
plot(cpred , "slices", ci = "bars", type = "p", pch = 19, exp_at = htemp, ylab="RR", main=paste0("Association for a high temperature (", htemp, "ºC)"))
Or the exposure–response association at lag 0:
plot(cpred , "slices", lag_at = 0, col=4, ylab="RR", main=paste0("Association at lag 0"))
