Load

## Load Harmonized Data ----
data <- read.csv("harmonized_data.csv")
head(data)
##   X year month       date    province dengue_cases population   tasmax
## 1 1 2023     1 2023-01-01  Xaisomboun            1     114000 23.89677
## 2 2 2023     1 2023-01-01      Attapu           10     166000 30.78710
## 3 3 2023     1 2023-01-01   Champasak            5     772000 30.78710
## 4 4 2023     1 2023-01-01      Xekong            3     134000 30.63548
## 5 5 2023     1 2023-01-01     Salavan            9     457000 29.26774
## 6 6 2023     1 2023-01-01 Savannakhet           15    1102000 28.32581
##      tasmin      tas prlr
## 1  7.970968 15.93387  8.8
## 2 20.403226 25.59516  0.8
## 3 18.535484 24.66129 32.8
## 4 16.403226 23.51935  5.0
## 5 17.719355 23.49355 19.2
## 6 15.461290 21.89355  0.7
## Load shapefile
laos<- st_read("gadm41_LAO1_cleaned.shp")
## Reading layer `gadm41_LAO1_cleaned' from data source 
##   `/home/akawieck/Documents/projects/dhis2.personal/ghrsuite_laos/gadm41_LAO1_cleaned.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 18 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 100.0868 ymin: 13.90968 xmax: 107.635 ymax: 22.5004
## Geodetic CRS:  WGS 84
ggplot() +
  geom_sf(data = laos)

# Calculate Pseudo-anomaly ----

data <- data %>%
  group_by(month, province) %>%
  mutate(tas.avg = mean(tas, na.rm = TRUE), 
         tasmax.avg = mean(tasmax, na.rm = TRUE),
         tasmin.avg = mean(tasmin, na.rm = TRUE),
         prlr.avg = mean(prlr, na.rm = TRUE)) %>%
  mutate(tas.diff = tas - tas.avg, 
         tasmax.diff = tasmax - tasmax.avg,
         tasmin.diff = tasmin - tasmin.avg,
         prlr.diff = prlr - prlr.avg) %>%
  ungroup() %>%
  # person time = 100000 person-month
  mutate(dengue_incidence=(dengue_cases / population) * 100000 ) # calculate incidence

GHRexplore description

GHRexplore plot functions:

plot_correlation

plot_bivariate

plot_timeseries

plot_heatmap

plot_seasonality

plot_map

GHRexplore plot management functions:

plot_multiple

Outputs a list with multiple plots, each representing one variable. Variable names and palette can be customized for each plot, otherwise parameters will be the same for each variable plotted (options depend on the plot type chosen). Returns a list of plots, each list element named after the variable plotted

plot_function indicates which of the ghr plot types to use Options are: ‘plot_timeseries’, ‘plot_heatmap’, ‘plot_seasonality’, ‘plot_bivariate’, ‘plot_map’.

data an object of class “data.frame” containing covariates or cases

var the name of the variables identifying the covariates to be plotted. Must be a character vector except if plot_function = plot_bivariate, in which case it must be a list of 2 character vectors (a list of variable pairs).

plot_combine

Combines plots, each representing one variable, into a single plot. This function takes any input from the cowplot::plot_grid function to customize the organization of the plots. To customize organization of a common plot legend to the combined plots, use parameters with the suffix “_l”.

plot_list requires as input a list of plots

combine_legend logical. If TRUE, assumes the legend of all plots is the same as the legend of the first plot in plot_list and final plot shows only one instance of the common legend. Default is FALSE.

combine_xaxis logical. If TRUE, removes x axis labels from all but the last plot. Default is FALSE.

ncol (from cowplot) number of colums in the plot grid. Default 1.

align (from cowplot) specifies how plots should be aligned. Options are “none”, “hv” (align in both directions), “h”, and “v” (default).

ncol_l when combine_legend = TRUE. Number of colums in which to align plots and the common legend. Default 2

nrow_l when combine_legend = TRUE. Number of rows in which to align plots and the common legend. Default NULL.

rel_widths_l when combine_legend = TRUE. Vector of widths in which to align plots and the common legend. Default = c(3, 1)

rel_heights_l when combine_legend = TRUE. Vector of heights in which to align plots and the common legend. Default = c(1, 1)

ncol_legend when combine_legend = TRUE number of columns the legend should be distributed in. Default 1 column

plot_compare

Combines multiple plots based on specified variables and plot type. This function takes any input arguments from ghrexplore::plot_combine() and from ghrexplore::plot_multiple() to customize the plots and their organization in a grid.

plot_function indicates which of the ghr plot types to use Options are: ‘plot_timeseries’, ‘plot_heatmap’, ‘plot_seasonality’, ‘plot_bivariate’, ‘plot_map’.

data an object of class “data.frame” containing covariates or cases

var the name of the variables identifying the covariates to be plotted. Must be a character vector except if plot_function = plot_bivariate, in which case it must be a list of 2 character vectors (a list of variable pairs).

Explore Laos data

Correlation

p_correlation <- plot_correlation(
  data= data,
  var = c("dengue_incidence", "tas", "tasmax", "tasmin", "prlr", 
          "tas.diff", "tasmax.diff", "tasmin.diff", "prlr.diff"),
  custom_label = c("Dengue Incidence",
                   "Average Temp.", "Max Temp.",
                   "Min Temp.", "Precipitation", 
                   "Average Temp. Anomaly", "Max Temp. Anomaly",
                   "Min Temp. Anomaly", "Precipitation Anomaly"),
  method = "pearson",
  palette = "BlueRed")
##                  Dengue_incidence         Tas      Tasmax      Tasmin
## Dengue_incidence     1.0000000000  0.17596357  0.09712722  0.22035502
## Tas                  0.1759635746  1.00000000  0.91632658  0.94122193
## Tasmax               0.0971272221  0.91632658  1.00000000  0.72721721
## Tasmin               0.2203550157  0.94122193  0.72721721  1.00000000
## Prlr                 0.2214521989  0.26724158  0.02647315  0.43559499
## Tas.diff            -0.0192117969  0.25138631  0.33410614  0.14866831
## Tasmax.diff         -0.0117324572  0.21859051  0.37464078  0.05822584
## Tasmin.diff         -0.0003059163  0.20513451  0.13927784  0.23406056
## Prlr.diff            0.0157685627 -0.07027154 -0.11537276 -0.02308505
##                         Prlr   Tas.diff Tasmax.diff   Tasmin.diff   Prlr.diff
## Dengue_incidence  0.22145220 -0.0192118 -0.01173246 -0.0003059163  0.01576856
## Tas               0.26724158  0.2513863  0.21859051  0.2051345150 -0.07027154
## Tasmax            0.02647315  0.3341061  0.37464078  0.1392778392 -0.11537276
## Tasmin            0.43559499  0.1486683  0.05822584  0.2340605605 -0.02308505
## Prlr              1.00000000 -0.1260153 -0.16124929 -0.0101481744  0.52818645
## Tas.diff         -0.12601529  1.0000000  0.88449686  0.7106647648 -0.22901375
## Tasmax.diff      -0.16124929  0.8844969  1.00000000  0.3082084670 -0.29212521
## Tasmin.diff      -0.01014817  0.7106648  0.30820847  1.0000000000 -0.03074776
## Prlr.diff         0.52818645 -0.2290138 -0.29212521 -0.0307477594  1.00000000
p_correlation

ggsave("p_correlation.jpg", plot = p_correlation, width = 15, height = 10, units = "cm", dpi = 300)

The correlation analysis suggests that the most precipitation, minimum temperature and average temperature might be the most relevant variables affecting dengue incidence.

Bivariate

# plot bivariate plots for all variables of interest into a list of plots

p_bivariate_list <- ghrexplore::plot_multiple(
 plot_function = plot_bivariate,
 data = data,
 var = list( c("tas", "dengue_incidence"),
             c("tasmin", "dengue_incidence"),
             c("tasmax", "dengue_incidence"),
             c("prlr", "dengue_incidence")),
 custom_labels = list(c("Average Temp.", "DENV incidence"),
                      c("Minimum Temp.", "DENV incidence"),
                      c("Maximum Temp.", "DENV incidence"),
                      c("Precipitation", "DENV incidence")),
 legend = "Province",
 group= "province",
 palette = c("Qualitative")
)
## [1] "Using palette: Qualitative"
## Warning in check_na(v, data): Missing values found in the ' tas ' column
## [1] "Using palette: Qualitative"
## Warning in check_na(v, data): Missing values found in the ' tasmin ' column
## [1] "Using palette: Qualitative"
## Warning in check_na(v, data): Missing values found in the ' tasmax ' column
## [1] "Using palette: Qualitative"
## Warning in check_na(v, data): Missing values found in the ' prlr ' column
# plot the difference from the mean plots 
p_bivariate_cases.vs.var<- ghrexplore::plot_combine(plot_list = p_bivariate_list[c("prlr_dengue_incidence",
                                                                                   "tas_dengue_incidence",
                                                                                   "tasmin_dengue_incidence",
                                                                                   "tasmax_dengue_incidence"
                                                                                   
                                                                                      )], 
                                                      ncol=2,
                                                      align = "hv", 
                                                      combine_legend = TRUE) 
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 146 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
p_bivariate_cases.vs.var

ggsave("p_bivariate_cases.vs.var.jpg", plot = p_bivariate_cases.vs.var, width = 30, height = 20, units = "cm", dpi = 300)

The bivariate analysis suggests that there might be non-linear effects between precipitation, temperature and dengue incidence, and apparent temperatures at which transmission is optimal.

Time series

# plot timeseries for all variables of interest into a list of plots
p_timeseries_list <- ghrexplore::plot_multiple(
  plot_function = plot_timeseries,
  data = data,
  var = c("tas", "tasmax", "tasmin", "prlr", 
          "tas.diff", "tasmax.diff", "tasmin.diff", "prlr.diff"),
  cases = "dengue_cases",
  pop = "population",
  type = "incidence",
  custom_label = c("Dengue Incidence",
                   "Average Temp.", "Max Temp.",
                   "Min Temp.", "Precipitation", 
                   "Average Temp. Anomaly", "Max Temp. Anomaly",
                   "Min Temp. Anomaly", "Precipitation Anomaly"),
  palette = c("Red",
                   "Purp", "Purp",
                   "Purp","Blue", 
                   "Purp", "Purp",
                   "Purp","Blue"),
  time = "date",
  area = "province", 
  panel =TRUE)
## Warning in plot_function(data = data, cases = cases, var = NULL, pop = pop, :
## More than 10 spatial units detected, evaluate 'highlight' or 'aggregate_space'
## option
## Warning in check_na(var, data): Missing values found in the ' tas ' column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
## Warning in check_na(var, data): Missing values found in the ' tasmax ' column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
## Warning in check_na(var, data): Missing values found in the ' tasmin ' column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
## Warning in check_na(var, data): Missing values found in the ' prlr ' column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
## Warning in check_na(var, data): Missing values found in the ' tas.diff ' column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
## Warning in check_na(var, data): Missing values found in the ' tasmax.diff '
## column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
## Warning in check_na(var, data): Missing values found in the ' tasmin.diff '
## column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
## Warning in check_na(var, data): Missing values found in the ' prlr.diff '
## column
## Warning in plot_function(data = data, var = v, custom_label = var_label, : More
## than 10 time series detected. Try 'highlight' or 'aggregate_space'?
# plot the difference from the mean plots 
p_timeseries_cases.vs.var<- ghrexplore::plot_combine(plot_list = p_timeseries_list[c("dengue_cases","prlr", "dengue_cases","tas"
                                                                                      )], 
                                                      ncol=2,
                                                      align = "hv") 

p_timeseries_cases.vs.var

ggsave("p_timeseries_cases.vs.var.jpg", plot = p_timeseries_cases.vs.var, width = 30, height = 20, units = "cm", dpi = 300)
# plot the difference from the mean plots 
p_timeseries_cases.vs.diff<- ghrexplore::plot_combine(plot_list = p_timeseries_list[c("dengue_cases","prlr.diff", "dengue_cases","tas.diff"
                                                                                      )], 
                                                      ncol=2,
                                                      align = "hv") 

p_timeseries_cases.vs.diff

ggsave("p_timeseries_cases.vs.diff.jpg", plot = p_timeseries_cases.vs.diff, width = 30, height = 20, units = "cm", dpi = 300)

It is apparent from the seasonal anomaly of average temperature that something happened in the province of Khammouan during 2019 and 2020 (there might have been a data collection error), as those 2 years present monthly temperatures 3 degrees lower than the average during those months.

Heat map

# plot heat map for all variables of interest into a list of plots
p_heatmap_list <- ghrexplore::plot_multiple(
  plot_function = plot_heatmap,
  data = data,
  var = c("tas", "tasmax", "tasmin", "prlr", 
          "tas.diff", "tasmax.diff", "tasmin.diff", "prlr.diff"),
  cases = "dengue_cases",
  pop = "population",
  type = "incidence",
  custom_label = c("Dengue Incidence",
                   "Average Temp.", "Max Temp.",
                   "Min Temp.", "Precipitation", 
                   "Average Temp. Anomaly", "Max Temp. Anomaly",
                   "Min Temp. Anomaly", "Precipitation Anomaly"),
  palette= c("Red",
                   "BlYlRd", "BlYlRd",
                   "BlYlRd", "Water", 
                   "BlueRed", "BlueRed",
                   "BlueRed", "BlueRed"),
  time = "date",
  area = "province")
## Warning in check_na(var, data): Missing values found in the ' tas ' column
## Warning in check_na(var, data): Missing values found in the ' tasmax ' column
## Warning in check_na(var, data): Missing values found in the ' tasmin ' column
## Warning in check_na(var, data): Missing values found in the ' prlr ' column
## Warning in check_na(var, data): Missing values found in the ' tas.diff ' column
## Warning in check_na(var, data): Missing values found in the ' tasmax.diff '
## column
## Warning in check_na(var, data): Missing values found in the ' tasmin.diff '
## column
## Warning in check_na(var, data): Missing values found in the ' prlr.diff '
## column
# plot the difference from the mean plots 
p_heatmap_cases.vs.var<- ghrexplore::plot_combine(plot_list = p_heatmap_list[c("dengue_cases","prlr",
                                                                               "dengue_cases", "tas"
                                                                                      )], 
                                                      ncol=2,
                                                      align = "hv") 

p_heatmap_cases.vs.var

ggsave("p_heatmap_cases.vs.var.jpg", plot = p_heatmap_cases.vs.var, width = 30, height = 20, units = "cm", dpi = 300)
# plot the difference from the mean plots 
p_heatmap_cases.vs.diff<- ghrexplore::plot_combine(plot_list = p_heatmap_list[c("dengue_cases","prlr.diff",
                                                                                "dengue_cases","tas.diff"
                                                                                      )], 
                                                      ncol=2,
                                                      align = "hv") 

p_heatmap_cases.vs.diff

ggsave("p_heatmap_cases.vs.diff.jpg", plot = p_heatmap_cases.vs.diff, width = 30, height = 20, units = "cm", dpi = 300)

These plots alos highlight the temperature anomaly in Khammouan province during 2019 and 2020.

Seasonality

# plot heat map for all variables of interest into a list of plots
p_seasonality_list <- ghrexplore::plot_multiple(
  plot_function = plot_seasonality,
  data = data,
  var = c("tas", "tasmax", "tasmin", "prlr", 
          "tas.diff", "tasmax.diff", "tasmin.diff", "prlr.diff"),
  cases = "dengue_cases",
  pop = "population",
  type = "incidence",
  custom_label = c("Dengue Incidence",
                   "Average Temp.", "Max Temp.",
                   "Min Temp.", "Precipitation", 
                   "Average Temp. Anomaly", "Max Temp. Anomaly",
                   "Min Temp. Anomaly", "Precipitation Anomaly"),
  palette= c("BlYlRd"),
  time = "date",
  area = "province")
## Warning in check_na(var, data): Missing values found in the ' tas ' column
## Warning in check_na(var, data): Missing values found in the ' tasmax ' column
## Warning in check_na(var, data): Missing values found in the ' tasmin ' column
## Warning in check_na(var, data): Missing values found in the ' prlr ' column
## Warning in check_na(var, data): Missing values found in the ' tas.diff ' column
## Warning in check_na(var, data): Missing values found in the ' tasmax.diff '
## column
## Warning in check_na(var, data): Missing values found in the ' tasmin.diff '
## column
## Warning in check_na(var, data): Missing values found in the ' prlr.diff '
## column
# plot the difference from the mean plots 
p_seasonality_cases.vs.var<- ghrexplore::plot_combine(plot_list = p_seasonality_list[c("dengue_cases","prlr",
                                                                                       "dengue_cases","tas"
                                                                                      )], 
                                                      ncol=2,
                                                      align = "hv", 
                                                      combine_legend = TRUE) 
## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_line()`).
p_seasonality_cases.vs.var

ggsave("p_seasonality_cases.vs.var.jpg", plot = p_seasonality_cases.vs.var, width = 30, height = 20, units = "cm", dpi = 300)

Seasonality plots suggest that rainfall and temperature is highly seasonal, following similar patterns in the different provinces of Laos. Dengue incidence also presents seasonal patterns, with increased number of cases in the more recent years (recent years represented in colors toward the red spectrum.

Map

# plot heat map for all variables of interest into a list of plots
p_map_list <- ghrexplore::plot_multiple(
  plot_function = plot_map,
  data = data,
  var = c("tas", "tasmax", "tasmin", "prlr"),
  cases = "dengue_cases",
  pop = "population",
  time = "date",
  area = "province", 
  custom_label = c("Dengue Incidence",
                   "Average Temp.", "Max Temp.",
                   "Min Temp.", "Precipitation"),
  palette= c("BlYlRd"),
  map = laos, 
  map_area = "NAME_1")
## Warning in check_na(var, data): Missing values found in the ' tas ' column
## Warning in check_na(var, data): Missing values found in the ' tasmax ' column
## Warning in check_na(var, data): Missing values found in the ' tasmin ' column
## Warning in check_na(var, data): Missing values found in the ' prlr ' column
# plot the difference from the mean plots 
p_map_cases.vs.var<- ghrexplore::plot_combine(plot_list = p_map_list[c("dengue_cases",
                                                                       "prlr",
                                                                       "dengue_cases",
                                                                       "tas"
                                                                                      )], 
                                                      ncol=2,
                                                      align = "hv") 

p_map_cases.vs.var

ggsave("p_map_cases.vs.var.jpg", plot = p_map_cases.vs.var, width = 30, height = 20, units = "cm", dpi = 300)

Compare

plot_compare(
   plot_function = plot_heatmap,
   data = data,
   var = c("prlr"),
   cases = "dengue_cases",
   type= "incidence",
   pop = "population",
   time = "date",
   area = "province",
   custom_label = c("Dengue Cases", "Precipitation"), 
   palette = c("Red", "BlYlRd"),
   ncol = 2,
   ncol_legend = 1)
## Warning in check_na(var, data): Missing values found in the ' prlr ' column