## 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
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).
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
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).
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.
# 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.
# 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.
# 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.
# 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.
# 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)
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