README.md
In caldwellst/augury: Provides Streamlined Methods for Data Imputation and Forecasting for WHO DDI Statistics
augury 
The goal of augury is to streamline the process of fitting models to
infill and forecast data, particularly for models used within the WHO’s
Triple Billion framework.
Installation
You can install augury from GitHub.
The package depends on the INLA package
which is not available on CRAN. You will need to separately install that
prior to installing augury, following the code below.
if (!require("INLA")) install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
remotes::install_github("gpw13/augury")
Overview
Most of the functions available through the package are designed to
streamline fitting a model and replacing missing observations within a
single data frame in one pass.
predict_general_mdl() is a function that takes in a data frame,
generic R modelling function (such as stats::lm), and fits a model
and returns newly projected data, error metrics, the fitted model
itself, or a list of all 3.predict_inla() is a function similar to the above, except rather
than taking a generic modelling function it accepts a formula and
additional argument passed to INLA::inla() to perform Bayesian
analysis of structured additive models.predict_forecast() instead uses forecasting methods available from
forecast::forecast() to perform exponential smoothing or other
time series models on a response variable.
Given that models might need to be fit separately to different portions
of a data frame, such as fitting time series models for each country
individually, these functions now accept the argument group_models
that determines whether to fit separate models to each group in the data
frame before joining them back together into the original data frame.
Additional model wrappers
Additional functions are provided as requested to streamlining fitting
of specific models. For general modeling, there are grouped and
individual functions for stats::lm(), stats::glm(), and
lme4::lmer() to fit linear models, generalized linear models, and
linear mixed-effects models respectively. These wrappers are:
predict_lm()predict_glm()predict_lmer()
As well, wrappers are provided around INLA for models currently in place
for the Triple Billion framework. These are a time series model with no
additional covariates (fit individually by country) and a mixed-effects
model using covariates.
predict_inla_me()predict_inla_ts()
And for forecasting, generic functions are provided for simple
exponential smoothing and Holt’s linear trend exponential smoothing.
predict_holt()predict_ses()
Covariates
For convenience, the covariates used within the default INLA modeling
are exported from the package and can be easily joined up to a data
frame for use in modeling. These are available through
augury::covariates_df.
Grouped trend development
While we typically want to fit a model to a data set directly, augury
has a set of predict_..._avg_trend() functions that allows you to
average data by group (e.g. by region), fit the model to the grouped and
averaged data, and extract that trend to apply to the original dataset.
See the section below on Average trend modeling for more details.
Error metrics
When doing model selection, it is always important to use metrics to
evaluate model fit and predictive accuracy. All predict_... functions
in augury use the same evaluation framework, implemented through
model_error(). If a test column is defined, then the evaluation
framework is only applied to this set
Additional functionality
To help streamline other portions of the modeling process, some other
functions are included in the package.
probit_transform() uses a probit transformation on select columns
of a data frame, and inverses it if specified.scale_transform() scales a vector by a single number, very simple,
and can inverse the scaling if specified.predict_simple() performs linear interpolation or flat
extrapolation in the same manner as the other predict_...
functions, but without modeling or confidence bounds.predict_average() performs averaging by groups of columns in the
same manner as the other predict_... functions, but without
modeling or confidence bounds.
Together, they can be used to transform and scale data for better
modeling, and then inverse these to get data back into the original
feature space.
INLA modeling examples
In order to show the INLA modeling wrappers provided in augury, we will
look at two datasets publicly available on the World Health
Organization’s Global Health Observatory. These can be accessed using
the ghost package, which provides
an R interface for the GHO OData API.
Time series modeling
The first indicator will be on safe sanitation. We will also use the
billionaiRe package to
quickly transform the GHO data into the simple format used by augury,
billionaiRe, and other packages.
library(augury)
df <- ghost::gho_data("WSH_SANITATION_SAFELY_MANAGED",
query = "$filter=Dim1 eq 'TOTL'") %>%
billionaiRe::wrangle_gho_data(source = "WHO GHO",
type = "estimated")
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <lgl> <lgl> <lgl> <lgl> <chr> <chr> <lgl>
#> 1 AFG 2000 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 2 AFG 2001 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 3 AFG 2002 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 4 AFG 2003 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 5 AFG 2004 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 6 AFG 2005 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
Now that we have the input data available from the GHO in an easy to use
format, we can now join up with the covariates_df available in augury
and run a time series model to predict sanitation out to 2023. For
simplicity, let’s just look at Albania, with ISO3 code "ALB".
library(dplyr)
df <- left_join(covariates_df,
df,
by = c("iso3", "year")) %>%
filter(iso3 == "ALB")
head(df)
#> # A tibble: 6 × 19
#> iso3 year year_n region sdi sdi_scaled e0 e0_scaled ind value lower
#> <chr> <int> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <lgl>
#> 1 ALB 2000 1 Southe… 0.604 0.650 74.0 0.752 hpop_… 40.2 NA
#> 2 ALB 2001 2 Southe… 0.61 0.657 74.3 0.759 hpop_… 40.5 NA
#> 3 ALB 2002 3 Southe… 0.614 0.661 74.6 0.766 hpop_… 40.8 NA
#> 4 ALB 2003 4 Southe… 0.619 0.667 74.8 0.771 hpop_… 41.1 NA
#> 5 ALB 2004 5 Southe… 0.625 0.673 75.0 0.776 hpop_… 41.4 NA
#> 6 ALB 2005 6 Southe… 0.632 0.681 75.2 0.780 hpop_… 41.9 NA
#> # … with 8 more variables: upper <lgl>, use_dash <lgl>, use_calc <lgl>,
#> # source <chr>, type <chr>, type_detail <lgl>, other_detail <lgl>,
#> # upload_detail <lgl>
Of course, the only “covariate” being used in this time series model is
going to be year_n, but the rest are available if we want to expand to
test other types of modeling. Let’s run the modeling now. We are going
to scale the data and probit transform it before and after the modeling.
We will use the predict_inla_ts() to fit a time series model to the
data.
modeled_df <- df %>%
scale_transform("value") %>%
probit_transform("value") %>%
predict_inla_ts(type_col = "type",
source_col = "source",
source = "augury modeling") %>%
probit_transform(c("value", "pred", "upper", "lower"), inverse = TRUE) %>%
scale_transform(c("value", "pred", "upper", "lower"), divide = FALSE)
# Only look at recent years and relevant columns
modeled_df %>%
filter(year > 2015) %>%
select(iso3, year, value, pred, lower, upper, source, type)
#> # A tibble: 10 × 8
#> iso3 year value pred lower upper source type
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 ALB 2016 46.7 46.7 NA NA WHO GHO estimated
#> 2 ALB 2017 47.1 47.1 NA NA WHO GHO estimated
#> 3 ALB 2018 47.4 47.4 NA NA WHO GHO estimated
#> 4 ALB 2019 47.6 47.6 NA NA WHO GHO estimated
#> 5 ALB 2020 47.7 47.7 NA NA WHO GHO estimated
#> 6 ALB 2021 47.9 47.9 NA NA augury modeling projected
#> 7 ALB 2022 48.1 48.1 NA NA augury modeling projected
#> 8 ALB 2023 48.2 48.2 NA NA augury modeling projected
#> 9 ALB 2024 48.4 48.4 NA NA augury modeling projected
#> 10 ALB 2025 48.6 48.6 NA NA augury modeling projected
And there we go, we have now fit a time series model to our data,
provided new type and source, and merged this into our existing data
frame. However, in this setup, the error calculations returned by
predict_inla_ts() are calculated in the probit space. If we wanted to
scale and probit transform the response variable prior to model fitting,
but still calculate error metrics and automatically return the response
and predicted values back in the original space, we can set
scale = 100 and probit = TRUE within predict_inla_ts().
df %>%
predict_inla_ts(scale = 100,
probit = TRUE,
type_col = "type",
source_col = "source",
source = "augury modeling") %>%
filter(year > 2015) %>%
select(iso3, year, value, pred, lower, upper, source, type)
#> # A tibble: 10 × 8
#> iso3 year value pred lower upper source type
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 ALB 2016 46.7 46.7 NA NA WHO GHO estimated
#> 2 ALB 2017 47.1 47.1 NA NA WHO GHO estimated
#> 3 ALB 2018 47.4 47.4 NA NA WHO GHO estimated
#> 4 ALB 2019 47.6 47.6 NA NA WHO GHO estimated
#> 5 ALB 2020 47.7 47.7 NA NA WHO GHO estimated
#> 6 ALB 2021 47.9 47.9 NA NA augury modeling projected
#> 7 ALB 2022 48.1 48.1 NA NA augury modeling projected
#> 8 ALB 2023 48.2 48.2 NA NA augury modeling projected
#> 9 ALB 2024 48.4 48.4 NA NA augury modeling projected
#> 10 ALB 2025 48.6 48.6 NA NA augury modeling projected
And we can see that the results here are the same as manually scaling
and probit transforming the variables.
Mixed-effects modeling
Now we will look at another indicator, a composite of 13 International
Health Regulations core capacity scores, SPAR version. Since countries
only have two data points at most, we will use mixed-effects modeling to
infill and project the data.
df <- ghost::gho_data("SDGIHR2018") %>%
billionaiRe::wrangle_gho_data(source = "Electronic State Parties Self-Assessment Annual Reporting Tool (e-SPAR)",
type = "reported")
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <lgl> <lgl> <lgl> <lgl> <chr> <chr> <lgl>
#> 1 AFG 2018 espar 35 NA NA TRUE TRUE Elect… repo… NA
#> 2 AFG 2019 espar 43 NA NA TRUE TRUE Elect… repo… NA
#> 3 AFG 2020 espar 47 NA NA TRUE TRUE Elect… repo… NA
#> 4 AGO 2018 espar 59 NA NA TRUE TRUE Elect… repo… NA
#> 5 AGO 2019 espar 63 NA NA TRUE TRUE Elect… repo… NA
#> 6 AGO 2020 espar 65 NA NA TRUE TRUE Elect… repo… NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
With this, let’s go straight into the modeling like last time, except we
will now use predict_inla_me() for mixed-effects modeling using
covariates found in covariates_df. This time, we want to model a first
order auto-regressive process across time rather than a second-order
random walk, so we use the "ar1" model available in INLA.
modeled_df <- df %>%
right_join(covariates_df, by = c("iso3", "year")) %>%
scale_transform("value") %>%
probit_transform("value") %>%
predict_inla_me(model = "ar1",
type_col = "type",
source_col = "source",
source = "WHO DDI Preliminary infilling and projections") %>%
probit_transform(c("value", "pred", "upper", "lower"), inverse = TRUE) %>%
scale_transform(c("value", "pred", "upper", "lower"), divide = FALSE)
# Look at an example for Afghanistan
modeled_df %>%
filter(year > 2017, iso3 == "AFG") %>%
select(iso3, year, value, pred, lower, upper, source, type)
#> # A tibble: 8 × 8
#> iso3 year value pred lower upper source type
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 AFG 2018 35 40.6 NA NA Electronic State Parties Self-Asse… repor…
#> 2 AFG 2019 43 41.5 NA NA Electronic State Parties Self-Asse… repor…
#> 3 AFG 2020 47 42.7 NA NA Electronic State Parties Self-Asse… repor…
#> 4 AFG 2021 43.5 43.5 NA NA WHO DDI Preliminary infilling and … proje…
#> 5 AFG 2022 44.5 44.5 NA NA WHO DDI Preliminary infilling and … proje…
#> 6 AFG 2023 45.4 45.4 NA NA WHO DDI Preliminary infilling and … proje…
#> 7 AFG 2024 46.4 46.4 NA NA WHO DDI Preliminary infilling and … proje…
#> 8 AFG 2025 47.3 47.3 NA NA WHO DDI Preliminary infilling and … proje…
And exactly as we were able to do with the time series modeling, we now
have infilled missing data for this indicator using mixed-effects
modeling in INLA.
Building further
Building further on this work, you can tweak any of the arguments passed
to these INLA models or use the base predict_inla() and other
covariates to test and compare other models. There is much more
functionality to test modeling accuracy and iteratively develop methods
available in this package not shown here, so please continue to explore
and play around.
Forecasting examples
To look at using forecast methods to predict data, we will again be
using the ghost package, which
provides an R interface for the GHO OData API and accessing data on
blood pressure. We will load in data for the USA and Great Britain
initially, which provide full time series from 1975 to 2015.
library(augury)
df <- ghost::gho_data("BP_04", query = "$filter=SpatialDim in ('USA', 'GBR') and Dim1 eq 'MLE' and Dim2 eq 'YEARS18-PLUS'") %>%
billionaiRe::wrangle_gho_data() %>%
dplyr::right_join(tidyr::expand_grid(iso3 = c("USA", "GBR"),
year = 1975:2017))
#> Warning: Some of the rows are missing a source value.
#> Joining, by = c("iso3", "year")
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <chr> <lgl>
#> 1 GBR 1975 bp 37.8 26.7 49.1 TRUE TRUE NA <NA> NA
#> 2 GBR 1976 bp 37.6 27.4 48 TRUE TRUE NA <NA> NA
#> 3 GBR 1977 bp 37.3 27.9 46.8 TRUE TRUE NA <NA> NA
#> 4 GBR 1978 bp 37.1 28.4 45.9 TRUE TRUE NA <NA> NA
#> 5 GBR 1979 bp 36.9 28.8 45.2 TRUE TRUE NA <NA> NA
#> 6 GBR 1980 bp 36.7 29.2 44.4 TRUE TRUE NA <NA> NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
With this data, we can now use the predict_forecast() function like we
would any of the other predict_... functions from augury to forecast
out to 2017. First, we will do this just on USA data and use the
forecast::holt to forecast using exponential smoothing.
usa_df <- dplyr::filter(df, iso3 == "USA")
predict_forecast(usa_df,
forecast::holt,
"value",
sort_col = "year") %>%
dplyr::filter(year >= 2012)
#> Registered S3 method overwritten by 'quantmod':
#> method from
#> as.zoo.data.frame zoo
#> # A tibble: 6 × 16
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <chr> <lgl>
#> 1 USA 2012 bp 15.7 11.7 20.3 TRUE TRUE NA <NA> NA
#> 2 USA 2013 bp 15.5 11.2 20.8 TRUE TRUE NA <NA> NA
#> 3 USA 2014 bp 15.4 10.8 21.3 TRUE TRUE NA <NA> NA
#> 4 USA 2015 bp 15.3 10.4 21.8 TRUE TRUE NA <NA> NA
#> 5 USA 2016 <NA> 15.2 NA NA NA NA NA <NA> NA
#> 6 USA 2017 <NA> 15.1 NA NA NA NA NA <NA> NA
#> # … with 5 more variables: other_detail <lgl>, upload_detail <lgl>, pred <dbl>,
#> # pred_upper <dbl>, pred_lower <dbl>
Of course, we might want to run these models all together for each
country individually. In this case, we can use the group_models = TRUE
function to perform the forecast individually by country. To save a bit
of limited time, let’s use the wrapper predict_holt() to automatically
supply forecast::holt as the forecasting function.
predict_holt(df,
response = "value",
group_col = "iso3",
group_models = TRUE,
sort_col = "year") %>%
dplyr::filter(year >= 2014, year <= 2017)
#> # A tibble: 8 × 16
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <chr> <lgl>
#> 1 GBR 2014 bp 18.5 14 23.3 TRUE TRUE NA <NA> NA
#> 2 GBR 2015 bp 17.9 13 23.2 TRUE TRUE NA <NA> NA
#> 3 GBR 2016 <NA> 17.3 NA NA NA NA NA <NA> NA
#> 4 GBR 2017 <NA> 16.7 NA NA NA NA NA <NA> NA
#> 5 USA 2014 bp 15.4 10.8 21.3 TRUE TRUE NA <NA> NA
#> 6 USA 2015 bp 15.3 10.4 21.8 TRUE TRUE NA <NA> NA
#> 7 USA 2016 <NA> 15.2 NA NA NA NA NA <NA> NA
#> 8 USA 2017 <NA> 15.1 NA NA NA NA NA <NA> NA
#> # … with 5 more variables: other_detail <lgl>, upload_detail <lgl>, pred <dbl>,
#> # pred_upper <dbl>, pred_lower <dbl>
Et voila, we have the same results for the USA and have also ran
forecasting on Great Britain as well. However, you should be careful on
the data that is supplied for forecasting. The forecast package
functions default to using the longest, contiguous non-missing data for
forecasting. augury instead automatically pulls the latest contiguous
observed data to use for forecasting, to ensure that older data is not
prioritized over new data. However, this means any break in a time
series will prevent data before that from being used.
bad_df <- dplyr::tibble(x = c(1:4, NA, 3:2, rep(NA, 4)))
predict_holt(bad_df, "x", group_col = NULL, sort_col = NULL, group_models = FALSE)
#> # A tibble: 11 × 6
#> x pred pred_upper pred_lower upper lower
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 NA NA NA NA NA
#> 2 2 NA NA NA NA NA
#> 3 3 NA NA NA NA NA
#> 4 4 NA NA NA NA NA
#> 5 NA NA NA NA NA NA
#> 6 3 NA NA NA NA NA
#> 7 2 NA NA NA NA NA
#> 8 1.17 1.17 2.55 -0.217 NA NA
#> 9 0.338 0.338 2.33 -1.66 NA NA
#> 10 -0.494 -0.494 2.14 -3.12 NA NA
#> 11 -1.32 -1.32 1.98 -4.63 NA NA
It’s advisable to consider if other data infilling or imputation methods
should be used to generate a full time series prior to the use of
forecasting methods to prevent issues like above from impacting the
predictive accuracy.
Simple prediction methods
The simple methods available in augury are easy to use, but provide the
same functionality of allowing a test column and returning error metrics
as the more complex modeling functions available in the package. Let’s
use data on alcohol from the GHO to demonstrate the functionality.
library(augury)
df <- ghost::gho_data("SA_0000001688",
query = "$filter=Dim1 eq 'BTSX'") %>%
billionaiRe::wrangle_gho_data(source = "WHO GHO",
type = "estimated") %>%
dplyr::arrange(iso3, year)
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <chr> <chr> <lgl>
#> 1 AFG 2000 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 2 AFG 2005 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 3 AFG 2010 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 4 AFG 2015 alcohol 0 0 0 TRUE TRUE WHO GHO esti… NA
#> 5 AFG 2019 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 6 AGO 2000 alcohol 3.3 2.4 4.4 TRUE TRUE WHO GHO esti… NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
Here we can see that data has time series and gaps in years. We can use
linear interpolation and flat extrapolation here to get data out to
2023.
df <- tidyr::expand_grid(iso3 = unique(df$iso3),
year = 2000:2023) %>%
dplyr::left_join(df, by = c("iso3", "year"))
df %>%
dplyr::filter(iso3 == "AFG",
year >= 2010,
year <= 2018) %>%
dplyr::select(iso3,
year,
value)
#> # A tibble: 9 × 3
#> iso3 year value
#> <chr> <int> <dbl>
#> 1 AFG 2010 0
#> 2 AFG 2011 NA
#> 3 AFG 2012 NA
#> 4 AFG 2013 NA
#> 5 AFG 2014 NA
#> 6 AFG 2015 0
#> 7 AFG 2016 NA
#> 8 AFG 2017 NA
#> 9 AFG 2018 NA
Let’s now use our linear interpolation and flat extrapolation on this
data.
pred_df <- predict_simple(df,
group_col = "iso3",
sort_col = "year")
pred_df %>%
dplyr::filter(iso3 == "AFG",
year >= 2010,
year <= 2018) %>%
dplyr::select(iso3, year, value)
#> # A tibble: 9 × 3
#> iso3 year value
#> <chr> <int> <dbl>
#> 1 AFG 2010 0
#> 2 AFG 2011 0
#> 3 AFG 2012 0
#> 4 AFG 2013 0
#> 5 AFG 2014 0
#> 6 AFG 2015 0
#> 7 AFG 2016 0
#> 8 AFG 2017 0
#> 9 AFG 2018 0
And we can see our linear interpolation there. We can also see the flat
extrapolation.
pred_df %>%
dplyr::filter(iso3 == "AFG",
year > 2016) %>%
dplyr::select(iso3, year, value)
#> # A tibble: 7 × 3
#> iso3 year value
#> <chr> <int> <dbl>
#> 1 AFG 2017 0
#> 2 AFG 2018 0
#> 3 AFG 2019 0
#> 4 AFG 2020 0
#> 5 AFG 2021 0
#> 6 AFG 2022 0
#> 7 AFG 2023 0
We can use the predict_average() function in much the same way, except
it is most useful when we have robust series for a set of countries, and
no data for others. We can then use something like the regional average
to infill data for missing countries.
df <- ghost::gho_data("PHE_HHAIR_PROP_POP_CLEAN_FUELS") %>%
billionaiRe::wrangle_gho_data(source = "WHO GHO",
type = "estimated") %>%
dplyr::filter(whoville::is_who_member(iso3))
#> Warning: Some of the rows are missing a ind value.
x <- whoville::who_member_states()
x[!(x %in% df$iso3)]
#> [1] "LBN" "CUB" "BGR" "LBY"
Above, we have 4 missing WHO member states, Lebanon, Cuba, Bulgaria, and
Libya. Let’s use regional averaging to fill in this data. We can use the
most recent World Bank income groups from the whoville
package as our relevant group.
df <- tidyr::expand_grid(iso3 = x,
year = 2000:2018) %>%
dplyr::left_join(df, by = c("iso3", "year")) %>%
dplyr::mutate(region = whoville::iso3_to_regions(iso3, region = "wb_ig"))
predict_average(df,
average_cols = c("region", "year"),
group_col = "iso3",
sort_col = "year",
type_col = "type",
source_col = "source",
source = "WB IG regional averages") %>%
dplyr::filter(iso3 == "LBN")
#> # A tibble: 19 × 15
#> iso3 year ind value lower upper use_dash use_calc source type
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <chr> <chr>
#> 1 LBN 2000 <NA> 67.2 NA NA NA NA WB IG regional … imput…
#> 2 LBN 2001 <NA> 68.2 NA NA NA NA WB IG regional … imput…
#> 3 LBN 2002 <NA> 69.2 NA NA NA NA WB IG regional … imput…
#> 4 LBN 2003 <NA> 70.2 NA NA NA NA WB IG regional … imput…
#> 5 LBN 2004 <NA> 71.2 NA NA NA NA WB IG regional … imput…
#> 6 LBN 2005 <NA> 72.2 NA NA NA NA WB IG regional … imput…
#> 7 LBN 2006 <NA> 73.1 NA NA NA NA WB IG regional … imput…
#> 8 LBN 2007 <NA> 74.0 NA NA NA NA WB IG regional … imput…
#> 9 LBN 2008 <NA> 74.8 NA NA NA NA WB IG regional … imput…
#> 10 LBN 2009 <NA> 75.6 NA NA NA NA WB IG regional … imput…
#> 11 LBN 2010 <NA> 76.3 NA NA NA NA WB IG regional … imput…
#> 12 LBN 2011 <NA> 77.0 NA NA NA NA WB IG regional … imput…
#> 13 LBN 2012 <NA> 77.6 NA NA NA NA WB IG regional … imput…
#> 14 LBN 2013 <NA> 78.2 NA NA NA NA WB IG regional … imput…
#> 15 LBN 2014 <NA> 78.7 NA NA NA NA WB IG regional … imput…
#> 16 LBN 2015 <NA> 79.2 NA NA NA NA WB IG regional … imput…
#> 17 LBN 2016 <NA> 79.7 NA NA NA NA WB IG regional … imput…
#> 18 LBN 2017 <NA> 80.1 NA NA NA NA WB IG regional … imput…
#> 19 LBN 2018 <NA> 80.5 NA NA NA NA WB IG regional … imput…
#> # … with 5 more variables: type_detail <lgl>, other_detail <lgl>,
#> # upload_detail <lgl>, region <chr>, pred <dbl>
Hope these examples have been clear and highlight some of the usefulness
of these simple modelling functions.
Average trend modeling
While we most often want to directly build models on our original
dataset to generate predicted values, we might instead want to generate
average trends across larger groups instead, and then apply this to our
original data. For instance, generating trends by region, and then
applying those regional trends back to the country level. The
predict_...avg_trend() functions in augury allow us to do just that,
applying any of the models we are used to a grouped set of columns.
These work across specific groups, specified by average_cols, and
averaging numeric values specified as the response variable or variables
extracted from a formula. The specified model is then fit to this
averaged data, and the predicted values are joined back up to the
original data frame. Let’s look at an example using blood pressure data,
which has a comprehensive time series.
library(augury)
df <- ghost::gho_data("BP_04", query = "$filter=Dim1 eq 'MLE' and Dim2 eq 'YEARS18-PLUS'") %>%
billionaiRe::wrangle_gho_data() %>%
dplyr::right_join(covariates_df) %>%
dplyr::select(iso3, year, year_n, value) %>%
dplyr::filter(whoville::is_who_member(iso3), # keep WHO member states
year >= 2000, year <= 2023) %>% # get relevant years
dplyr::mutate(who_region = whoville::iso3_to_regions(iso3)) # get WHO regions
#> Warning: Some of the rows are missing a source value.
#> Joining, by = c("iso3", "year")
ur <- unique(df$who_region)
ur
#> [1] "EMR" "AFR" "EUR" "AMR" "WPR" "SEAR"
Alright, so, here we have 6 WHO regions. We will use these regions to
fit a model to and use INLA to predict out to 2023, then apply these
trends to input countries.
pred_df <- df %>%
predict_inla_avg_trend(formula = value ~ f(year_n, model = "rw2"),
average_cols = c("who_region", "year_n"),
group_models = TRUE,
group_col = "iso3",
sort_col = "year_n")
pred_df %>%
dplyr::filter(iso3 == "AFG", year >= 2013)
#> # A tibble: 11 × 10
#> iso3 year year_n value who_region pred pred_upper pred_lower upper lower
#> <chr> <int> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 AFG 2013 14 30.4 EMR 29.4 29.4 29.4 NA NA
#> 2 AFG 2014 15 30.4 EMR 29.4 29.4 29.4 NA NA
#> 3 AFG 2015 16 30.4 EMR 29.4 29.4 29.4 NA NA
#> 4 AFG 2016 17 29.4 EMR 29.4 29.4 29.4 NA NA
#> 5 AFG 2017 18 29.4 EMR 29.4 29.4 29.4 NA NA
#> 6 AFG 2018 19 29.4 EMR 29.4 29.4 29.4 NA NA
#> 7 AFG 2019 20 29.4 EMR 29.4 29.4 29.4 NA NA
#> 8 AFG 2020 21 29.4 EMR 29.4 29.4 29.4 NA NA
#> 9 AFG 2021 22 29.4 EMR 29.4 29.4 29.4 NA NA
#> 10 AFG 2022 23 29.4 EMR 29.4 29.4 29.4 NA NA
#> 11 AFG 2023 24 29.4 EMR 29.4 29.4 29.4 NA NA
Above, we can see we have a generated a model using 2nd order random
walk with INLA, however, the model was generated by averaging data to
WHO regions first, fitting the random walk to each reach (since
group_models = TRUE) and then fitting those trends to the original
data. Note some specifics of what had to be set, as the
predict_..._avg_trend() functions are slightly more complex than
others:
average_cols must contain the sort_col. So, since we use
year_n in the time series rather than year, we will sort by that
this time.average_cols refers to the groupings used for averaging (we take
the average for each WHO region and year in this case). Then, the
model is fit to average_cols that are NOT the sort_col.group_col is the groupings used on the original data frame,
which is still necessary here when applying the trend back to the
original data.- If a variable is in
formula, it must either be in average_cols
or it must be a numeric column that can be averaged. This is because
the formula is applied to the data frame after dplyr::group_by()
and dplyr::summarize() have reduced it.
To highlight this point, in the above example, what’s actually happening
is we are actually fitting the model on the summarized data:
df %>%
dplyr::group_by(who_region, year_n) %>%
dplyr::summarize(value = mean(value, na.rm = T)) %>%
head()
#> `summarise()` has grouped output by 'who_region'. You can override using the `.groups` argument.
#> # A tibble: 6 × 3
#> # Groups: who_region [1]
#> who_region year_n value
#> <chr> <dbl> <dbl>
#> 1 AFR 1 29.7
#> 2 AFR 2 29.6
#> 3 AFR 3 29.5
#> 4 AFR 4 29.4
#> 5 AFR 5 29.3
#> 6 AFR 6 29.2
Since average_cols = c("who_region", "year_n"), we took the mean of
all values in formula not in average_cols, in this case just
value. If for instance, we tried to specify a model using iso3 in
the formula:
predict_inla_avg_trend(df,
formula = value ~ iso3 + f(year_n, model = "rw2"),
average_cols = c("who_region", "year_n"),
group_models = TRUE,
group_col = "iso3",
sort_col = "year_n")
#> Error: iso3 must be numeric columns for use in averaging, or included in `average_cols` for grouping.
We get an error message indicating that iso3 must be numeric or included
in average_cols for grouping. This is because without it being numeric
or in the average_cols, there’s no way
dplyr::group_by() %>% dplyr::summarize() a non-numeric column
automatically (how would we reduce country-level ISO3 codes to the
regional level?).
While slightly complex, ensuring you follow the above means you should
easily and successfully get out meaningful trend predictions for your
data frames using trends generated on grouped data.
caldwellst/augury documentation built on Oct. 10, 2024, 8:20 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
augury
The goal of augury is to streamline the process of fitting models to infill and forecast data, particularly for models used within the WHO’s Triple Billion framework.
Installation
You can install augury from GitHub. The package depends on the INLA package which is not available on CRAN. You will need to separately install that prior to installing augury, following the code below.
if (!require("INLA")) install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
remotes::install_github("gpw13/augury")
Overview
Most of the functions available through the package are designed to streamline fitting a model and replacing missing observations within a single data frame in one pass.
predict_general_mdl()is a function that takes in a data frame, generic R modelling function (such asstats::lm), and fits a model and returns newly projected data, error metrics, the fitted model itself, or a list of all 3.predict_inla()is a function similar to the above, except rather than taking a generic modelling function it accepts a formula and additional argument passed toINLA::inla()to perform Bayesian analysis of structured additive models.predict_forecast()instead uses forecasting methods available fromforecast::forecast()to perform exponential smoothing or other time series models on a response variable.
Given that models might need to be fit separately to different portions
of a data frame, such as fitting time series models for each country
individually, these functions now accept the argument group_models
that determines whether to fit separate models to each group in the data
frame before joining them back together into the original data frame.
Additional model wrappers
Additional functions are provided as requested to streamlining fitting
of specific models. For general modeling, there are grouped and
individual functions for stats::lm(), stats::glm(), and
lme4::lmer() to fit linear models, generalized linear models, and
linear mixed-effects models respectively. These wrappers are:
predict_lm()predict_glm()predict_lmer()
As well, wrappers are provided around INLA for models currently in place for the Triple Billion framework. These are a time series model with no additional covariates (fit individually by country) and a mixed-effects model using covariates.
predict_inla_me()predict_inla_ts()
And for forecasting, generic functions are provided for simple exponential smoothing and Holt’s linear trend exponential smoothing.
predict_holt()predict_ses()
Covariates
For convenience, the covariates used within the default INLA modeling
are exported from the package and can be easily joined up to a data
frame for use in modeling. These are available through
augury::covariates_df.
Grouped trend development
While we typically want to fit a model to a data set directly, augury
has a set of predict_..._avg_trend() functions that allows you to
average data by group (e.g. by region), fit the model to the grouped and
averaged data, and extract that trend to apply to the original dataset.
See the section below on Average trend modeling for more details.
Error metrics
When doing model selection, it is always important to use metrics to
evaluate model fit and predictive accuracy. All predict_... functions
in augury use the same evaluation framework, implemented through
model_error(). If a test column is defined, then the evaluation
framework is only applied to this set
Additional functionality
To help streamline other portions of the modeling process, some other functions are included in the package.
probit_transform()uses a probit transformation on select columns of a data frame, and inverses it if specified.scale_transform()scales a vector by a single number, very simple, and can inverse the scaling if specified.predict_simple()performs linear interpolation or flat extrapolation in the same manner as the otherpredict_...functions, but without modeling or confidence bounds.predict_average()performs averaging by groups of columns in the same manner as the otherpredict_...functions, but without modeling or confidence bounds.
Together, they can be used to transform and scale data for better modeling, and then inverse these to get data back into the original feature space.
INLA modeling examples
In order to show the INLA modeling wrappers provided in augury, we will look at two datasets publicly available on the World Health Organization’s Global Health Observatory. These can be accessed using the ghost package, which provides an R interface for the GHO OData API.
Time series modeling
The first indicator will be on safe sanitation. We will also use the billionaiRe package to quickly transform the GHO data into the simple format used by augury, billionaiRe, and other packages.
library(augury)
df <- ghost::gho_data("WSH_SANITATION_SAFELY_MANAGED",
query = "$filter=Dim1 eq 'TOTL'") %>%
billionaiRe::wrangle_gho_data(source = "WHO GHO",
type = "estimated")
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <lgl> <lgl> <lgl> <lgl> <chr> <chr> <lgl>
#> 1 AFG 2000 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 2 AFG 2001 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 3 AFG 2002 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 4 AFG 2003 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 5 AFG 2004 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> 6 AFG 2005 hpop… NA NA NA TRUE TRUE WHO G… esti… NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
Now that we have the input data available from the GHO in an easy to use
format, we can now join up with the covariates_df available in augury
and run a time series model to predict sanitation out to 2023. For
simplicity, let’s just look at Albania, with ISO3 code "ALB".
library(dplyr)
df <- left_join(covariates_df,
df,
by = c("iso3", "year")) %>%
filter(iso3 == "ALB")
head(df)
#> # A tibble: 6 × 19
#> iso3 year year_n region sdi sdi_scaled e0 e0_scaled ind value lower
#> <chr> <int> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <lgl>
#> 1 ALB 2000 1 Southe… 0.604 0.650 74.0 0.752 hpop_… 40.2 NA
#> 2 ALB 2001 2 Southe… 0.61 0.657 74.3 0.759 hpop_… 40.5 NA
#> 3 ALB 2002 3 Southe… 0.614 0.661 74.6 0.766 hpop_… 40.8 NA
#> 4 ALB 2003 4 Southe… 0.619 0.667 74.8 0.771 hpop_… 41.1 NA
#> 5 ALB 2004 5 Southe… 0.625 0.673 75.0 0.776 hpop_… 41.4 NA
#> 6 ALB 2005 6 Southe… 0.632 0.681 75.2 0.780 hpop_… 41.9 NA
#> # … with 8 more variables: upper <lgl>, use_dash <lgl>, use_calc <lgl>,
#> # source <chr>, type <chr>, type_detail <lgl>, other_detail <lgl>,
#> # upload_detail <lgl>
Of course, the only “covariate” being used in this time series model is
going to be year_n, but the rest are available if we want to expand to
test other types of modeling. Let’s run the modeling now. We are going
to scale the data and probit transform it before and after the modeling.
We will use the predict_inla_ts() to fit a time series model to the
data.
modeled_df <- df %>%
scale_transform("value") %>%
probit_transform("value") %>%
predict_inla_ts(type_col = "type",
source_col = "source",
source = "augury modeling") %>%
probit_transform(c("value", "pred", "upper", "lower"), inverse = TRUE) %>%
scale_transform(c("value", "pred", "upper", "lower"), divide = FALSE)
# Only look at recent years and relevant columns
modeled_df %>%
filter(year > 2015) %>%
select(iso3, year, value, pred, lower, upper, source, type)
#> # A tibble: 10 × 8
#> iso3 year value pred lower upper source type
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 ALB 2016 46.7 46.7 NA NA WHO GHO estimated
#> 2 ALB 2017 47.1 47.1 NA NA WHO GHO estimated
#> 3 ALB 2018 47.4 47.4 NA NA WHO GHO estimated
#> 4 ALB 2019 47.6 47.6 NA NA WHO GHO estimated
#> 5 ALB 2020 47.7 47.7 NA NA WHO GHO estimated
#> 6 ALB 2021 47.9 47.9 NA NA augury modeling projected
#> 7 ALB 2022 48.1 48.1 NA NA augury modeling projected
#> 8 ALB 2023 48.2 48.2 NA NA augury modeling projected
#> 9 ALB 2024 48.4 48.4 NA NA augury modeling projected
#> 10 ALB 2025 48.6 48.6 NA NA augury modeling projected
And there we go, we have now fit a time series model to our data,
provided new type and source, and merged this into our existing data
frame. However, in this setup, the error calculations returned by
predict_inla_ts() are calculated in the probit space. If we wanted to
scale and probit transform the response variable prior to model fitting,
but still calculate error metrics and automatically return the response
and predicted values back in the original space, we can set
scale = 100 and probit = TRUE within predict_inla_ts().
df %>%
predict_inla_ts(scale = 100,
probit = TRUE,
type_col = "type",
source_col = "source",
source = "augury modeling") %>%
filter(year > 2015) %>%
select(iso3, year, value, pred, lower, upper, source, type)
#> # A tibble: 10 × 8
#> iso3 year value pred lower upper source type
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 ALB 2016 46.7 46.7 NA NA WHO GHO estimated
#> 2 ALB 2017 47.1 47.1 NA NA WHO GHO estimated
#> 3 ALB 2018 47.4 47.4 NA NA WHO GHO estimated
#> 4 ALB 2019 47.6 47.6 NA NA WHO GHO estimated
#> 5 ALB 2020 47.7 47.7 NA NA WHO GHO estimated
#> 6 ALB 2021 47.9 47.9 NA NA augury modeling projected
#> 7 ALB 2022 48.1 48.1 NA NA augury modeling projected
#> 8 ALB 2023 48.2 48.2 NA NA augury modeling projected
#> 9 ALB 2024 48.4 48.4 NA NA augury modeling projected
#> 10 ALB 2025 48.6 48.6 NA NA augury modeling projected
And we can see that the results here are the same as manually scaling and probit transforming the variables.
Mixed-effects modeling
Now we will look at another indicator, a composite of 13 International Health Regulations core capacity scores, SPAR version. Since countries only have two data points at most, we will use mixed-effects modeling to infill and project the data.
df <- ghost::gho_data("SDGIHR2018") %>%
billionaiRe::wrangle_gho_data(source = "Electronic State Parties Self-Assessment Annual Reporting Tool (e-SPAR)",
type = "reported")
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <lgl> <lgl> <lgl> <lgl> <chr> <chr> <lgl>
#> 1 AFG 2018 espar 35 NA NA TRUE TRUE Elect… repo… NA
#> 2 AFG 2019 espar 43 NA NA TRUE TRUE Elect… repo… NA
#> 3 AFG 2020 espar 47 NA NA TRUE TRUE Elect… repo… NA
#> 4 AGO 2018 espar 59 NA NA TRUE TRUE Elect… repo… NA
#> 5 AGO 2019 espar 63 NA NA TRUE TRUE Elect… repo… NA
#> 6 AGO 2020 espar 65 NA NA TRUE TRUE Elect… repo… NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
With this, let’s go straight into the modeling like last time, except we
will now use predict_inla_me() for mixed-effects modeling using
covariates found in covariates_df. This time, we want to model a first
order auto-regressive process across time rather than a second-order
random walk, so we use the "ar1" model available in INLA.
modeled_df <- df %>%
right_join(covariates_df, by = c("iso3", "year")) %>%
scale_transform("value") %>%
probit_transform("value") %>%
predict_inla_me(model = "ar1",
type_col = "type",
source_col = "source",
source = "WHO DDI Preliminary infilling and projections") %>%
probit_transform(c("value", "pred", "upper", "lower"), inverse = TRUE) %>%
scale_transform(c("value", "pred", "upper", "lower"), divide = FALSE)
# Look at an example for Afghanistan
modeled_df %>%
filter(year > 2017, iso3 == "AFG") %>%
select(iso3, year, value, pred, lower, upper, source, type)
#> # A tibble: 8 × 8
#> iso3 year value pred lower upper source type
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 AFG 2018 35 40.6 NA NA Electronic State Parties Self-Asse… repor…
#> 2 AFG 2019 43 41.5 NA NA Electronic State Parties Self-Asse… repor…
#> 3 AFG 2020 47 42.7 NA NA Electronic State Parties Self-Asse… repor…
#> 4 AFG 2021 43.5 43.5 NA NA WHO DDI Preliminary infilling and … proje…
#> 5 AFG 2022 44.5 44.5 NA NA WHO DDI Preliminary infilling and … proje…
#> 6 AFG 2023 45.4 45.4 NA NA WHO DDI Preliminary infilling and … proje…
#> 7 AFG 2024 46.4 46.4 NA NA WHO DDI Preliminary infilling and … proje…
#> 8 AFG 2025 47.3 47.3 NA NA WHO DDI Preliminary infilling and … proje…
And exactly as we were able to do with the time series modeling, we now have infilled missing data for this indicator using mixed-effects modeling in INLA.
Building further
Building further on this work, you can tweak any of the arguments passed
to these INLA models or use the base predict_inla() and other
covariates to test and compare other models. There is much more
functionality to test modeling accuracy and iteratively develop methods
available in this package not shown here, so please continue to explore
and play around.
Forecasting examples
To look at using forecast methods to predict data, we will again be using the ghost package, which provides an R interface for the GHO OData API and accessing data on blood pressure. We will load in data for the USA and Great Britain initially, which provide full time series from 1975 to 2015.
library(augury)
df <- ghost::gho_data("BP_04", query = "$filter=SpatialDim in ('USA', 'GBR') and Dim1 eq 'MLE' and Dim2 eq 'YEARS18-PLUS'") %>%
billionaiRe::wrangle_gho_data() %>%
dplyr::right_join(tidyr::expand_grid(iso3 = c("USA", "GBR"),
year = 1975:2017))
#> Warning: Some of the rows are missing a source value.
#> Joining, by = c("iso3", "year")
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <chr> <lgl>
#> 1 GBR 1975 bp 37.8 26.7 49.1 TRUE TRUE NA <NA> NA
#> 2 GBR 1976 bp 37.6 27.4 48 TRUE TRUE NA <NA> NA
#> 3 GBR 1977 bp 37.3 27.9 46.8 TRUE TRUE NA <NA> NA
#> 4 GBR 1978 bp 37.1 28.4 45.9 TRUE TRUE NA <NA> NA
#> 5 GBR 1979 bp 36.9 28.8 45.2 TRUE TRUE NA <NA> NA
#> 6 GBR 1980 bp 36.7 29.2 44.4 TRUE TRUE NA <NA> NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
With this data, we can now use the predict_forecast() function like we
would any of the other predict_... functions from augury to forecast
out to 2017. First, we will do this just on USA data and use the
forecast::holt to forecast using exponential smoothing.
usa_df <- dplyr::filter(df, iso3 == "USA")
predict_forecast(usa_df,
forecast::holt,
"value",
sort_col = "year") %>%
dplyr::filter(year >= 2012)
#> Registered S3 method overwritten by 'quantmod':
#> method from
#> as.zoo.data.frame zoo
#> # A tibble: 6 × 16
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <chr> <lgl>
#> 1 USA 2012 bp 15.7 11.7 20.3 TRUE TRUE NA <NA> NA
#> 2 USA 2013 bp 15.5 11.2 20.8 TRUE TRUE NA <NA> NA
#> 3 USA 2014 bp 15.4 10.8 21.3 TRUE TRUE NA <NA> NA
#> 4 USA 2015 bp 15.3 10.4 21.8 TRUE TRUE NA <NA> NA
#> 5 USA 2016 <NA> 15.2 NA NA NA NA NA <NA> NA
#> 6 USA 2017 <NA> 15.1 NA NA NA NA NA <NA> NA
#> # … with 5 more variables: other_detail <lgl>, upload_detail <lgl>, pred <dbl>,
#> # pred_upper <dbl>, pred_lower <dbl>
Of course, we might want to run these models all together for each
country individually. In this case, we can use the group_models = TRUE
function to perform the forecast individually by country. To save a bit
of limited time, let’s use the wrapper predict_holt() to automatically
supply forecast::holt as the forecasting function.
predict_holt(df,
response = "value",
group_col = "iso3",
group_models = TRUE,
sort_col = "year") %>%
dplyr::filter(year >= 2014, year <= 2017)
#> # A tibble: 8 × 16
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <chr> <lgl>
#> 1 GBR 2014 bp 18.5 14 23.3 TRUE TRUE NA <NA> NA
#> 2 GBR 2015 bp 17.9 13 23.2 TRUE TRUE NA <NA> NA
#> 3 GBR 2016 <NA> 17.3 NA NA NA NA NA <NA> NA
#> 4 GBR 2017 <NA> 16.7 NA NA NA NA NA <NA> NA
#> 5 USA 2014 bp 15.4 10.8 21.3 TRUE TRUE NA <NA> NA
#> 6 USA 2015 bp 15.3 10.4 21.8 TRUE TRUE NA <NA> NA
#> 7 USA 2016 <NA> 15.2 NA NA NA NA NA <NA> NA
#> 8 USA 2017 <NA> 15.1 NA NA NA NA NA <NA> NA
#> # … with 5 more variables: other_detail <lgl>, upload_detail <lgl>, pred <dbl>,
#> # pred_upper <dbl>, pred_lower <dbl>
Et voila, we have the same results for the USA and have also ran
forecasting on Great Britain as well. However, you should be careful on
the data that is supplied for forecasting. The forecast package
functions default to using the longest, contiguous non-missing data for
forecasting. augury instead automatically pulls the latest contiguous
observed data to use for forecasting, to ensure that older data is not
prioritized over new data. However, this means any break in a time
series will prevent data before that from being used.
bad_df <- dplyr::tibble(x = c(1:4, NA, 3:2, rep(NA, 4)))
predict_holt(bad_df, "x", group_col = NULL, sort_col = NULL, group_models = FALSE)
#> # A tibble: 11 × 6
#> x pred pred_upper pred_lower upper lower
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 NA NA NA NA NA
#> 2 2 NA NA NA NA NA
#> 3 3 NA NA NA NA NA
#> 4 4 NA NA NA NA NA
#> 5 NA NA NA NA NA NA
#> 6 3 NA NA NA NA NA
#> 7 2 NA NA NA NA NA
#> 8 1.17 1.17 2.55 -0.217 NA NA
#> 9 0.338 0.338 2.33 -1.66 NA NA
#> 10 -0.494 -0.494 2.14 -3.12 NA NA
#> 11 -1.32 -1.32 1.98 -4.63 NA NA
It’s advisable to consider if other data infilling or imputation methods should be used to generate a full time series prior to the use of forecasting methods to prevent issues like above from impacting the predictive accuracy.
Simple prediction methods
The simple methods available in augury are easy to use, but provide the same functionality of allowing a test column and returning error metrics as the more complex modeling functions available in the package. Let’s use data on alcohol from the GHO to demonstrate the functionality.
library(augury)
df <- ghost::gho_data("SA_0000001688",
query = "$filter=Dim1 eq 'BTSX'") %>%
billionaiRe::wrangle_gho_data(source = "WHO GHO",
type = "estimated") %>%
dplyr::arrange(iso3, year)
head(df)
#> # A tibble: 6 × 13
#> iso3 year ind value lower upper use_dash use_calc source type type_detail
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <chr> <chr> <lgl>
#> 1 AFG 2000 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 2 AFG 2005 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 3 AFG 2010 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 4 AFG 2015 alcohol 0 0 0 TRUE TRUE WHO GHO esti… NA
#> 5 AFG 2019 alcohol 0 0 0.1 TRUE TRUE WHO GHO esti… NA
#> 6 AGO 2000 alcohol 3.3 2.4 4.4 TRUE TRUE WHO GHO esti… NA
#> # … with 2 more variables: other_detail <lgl>, upload_detail <lgl>
Here we can see that data has time series and gaps in years. We can use linear interpolation and flat extrapolation here to get data out to 2023.
df <- tidyr::expand_grid(iso3 = unique(df$iso3),
year = 2000:2023) %>%
dplyr::left_join(df, by = c("iso3", "year"))
df %>%
dplyr::filter(iso3 == "AFG",
year >= 2010,
year <= 2018) %>%
dplyr::select(iso3,
year,
value)
#> # A tibble: 9 × 3
#> iso3 year value
#> <chr> <int> <dbl>
#> 1 AFG 2010 0
#> 2 AFG 2011 NA
#> 3 AFG 2012 NA
#> 4 AFG 2013 NA
#> 5 AFG 2014 NA
#> 6 AFG 2015 0
#> 7 AFG 2016 NA
#> 8 AFG 2017 NA
#> 9 AFG 2018 NA
Let’s now use our linear interpolation and flat extrapolation on this data.
pred_df <- predict_simple(df,
group_col = "iso3",
sort_col = "year")
pred_df %>%
dplyr::filter(iso3 == "AFG",
year >= 2010,
year <= 2018) %>%
dplyr::select(iso3, year, value)
#> # A tibble: 9 × 3
#> iso3 year value
#> <chr> <int> <dbl>
#> 1 AFG 2010 0
#> 2 AFG 2011 0
#> 3 AFG 2012 0
#> 4 AFG 2013 0
#> 5 AFG 2014 0
#> 6 AFG 2015 0
#> 7 AFG 2016 0
#> 8 AFG 2017 0
#> 9 AFG 2018 0
And we can see our linear interpolation there. We can also see the flat extrapolation.
pred_df %>%
dplyr::filter(iso3 == "AFG",
year > 2016) %>%
dplyr::select(iso3, year, value)
#> # A tibble: 7 × 3
#> iso3 year value
#> <chr> <int> <dbl>
#> 1 AFG 2017 0
#> 2 AFG 2018 0
#> 3 AFG 2019 0
#> 4 AFG 2020 0
#> 5 AFG 2021 0
#> 6 AFG 2022 0
#> 7 AFG 2023 0
We can use the predict_average() function in much the same way, except
it is most useful when we have robust series for a set of countries, and
no data for others. We can then use something like the regional average
to infill data for missing countries.
df <- ghost::gho_data("PHE_HHAIR_PROP_POP_CLEAN_FUELS") %>%
billionaiRe::wrangle_gho_data(source = "WHO GHO",
type = "estimated") %>%
dplyr::filter(whoville::is_who_member(iso3))
#> Warning: Some of the rows are missing a ind value.
x <- whoville::who_member_states()
x[!(x %in% df$iso3)]
#> [1] "LBN" "CUB" "BGR" "LBY"
Above, we have 4 missing WHO member states, Lebanon, Cuba, Bulgaria, and Libya. Let’s use regional averaging to fill in this data. We can use the most recent World Bank income groups from the whoville package as our relevant group.
df <- tidyr::expand_grid(iso3 = x,
year = 2000:2018) %>%
dplyr::left_join(df, by = c("iso3", "year")) %>%
dplyr::mutate(region = whoville::iso3_to_regions(iso3, region = "wb_ig"))
predict_average(df,
average_cols = c("region", "year"),
group_col = "iso3",
sort_col = "year",
type_col = "type",
source_col = "source",
source = "WB IG regional averages") %>%
dplyr::filter(iso3 == "LBN")
#> # A tibble: 19 × 15
#> iso3 year ind value lower upper use_dash use_calc source type
#> <chr> <int> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <chr> <chr>
#> 1 LBN 2000 <NA> 67.2 NA NA NA NA WB IG regional … imput…
#> 2 LBN 2001 <NA> 68.2 NA NA NA NA WB IG regional … imput…
#> 3 LBN 2002 <NA> 69.2 NA NA NA NA WB IG regional … imput…
#> 4 LBN 2003 <NA> 70.2 NA NA NA NA WB IG regional … imput…
#> 5 LBN 2004 <NA> 71.2 NA NA NA NA WB IG regional … imput…
#> 6 LBN 2005 <NA> 72.2 NA NA NA NA WB IG regional … imput…
#> 7 LBN 2006 <NA> 73.1 NA NA NA NA WB IG regional … imput…
#> 8 LBN 2007 <NA> 74.0 NA NA NA NA WB IG regional … imput…
#> 9 LBN 2008 <NA> 74.8 NA NA NA NA WB IG regional … imput…
#> 10 LBN 2009 <NA> 75.6 NA NA NA NA WB IG regional … imput…
#> 11 LBN 2010 <NA> 76.3 NA NA NA NA WB IG regional … imput…
#> 12 LBN 2011 <NA> 77.0 NA NA NA NA WB IG regional … imput…
#> 13 LBN 2012 <NA> 77.6 NA NA NA NA WB IG regional … imput…
#> 14 LBN 2013 <NA> 78.2 NA NA NA NA WB IG regional … imput…
#> 15 LBN 2014 <NA> 78.7 NA NA NA NA WB IG regional … imput…
#> 16 LBN 2015 <NA> 79.2 NA NA NA NA WB IG regional … imput…
#> 17 LBN 2016 <NA> 79.7 NA NA NA NA WB IG regional … imput…
#> 18 LBN 2017 <NA> 80.1 NA NA NA NA WB IG regional … imput…
#> 19 LBN 2018 <NA> 80.5 NA NA NA NA WB IG regional … imput…
#> # … with 5 more variables: type_detail <lgl>, other_detail <lgl>,
#> # upload_detail <lgl>, region <chr>, pred <dbl>
Hope these examples have been clear and highlight some of the usefulness of these simple modelling functions.
Average trend modeling
While we most often want to directly build models on our original
dataset to generate predicted values, we might instead want to generate
average trends across larger groups instead, and then apply this to our
original data. For instance, generating trends by region, and then
applying those regional trends back to the country level. The
predict_...avg_trend() functions in augury allow us to do just that,
applying any of the models we are used to a grouped set of columns.
These work across specific groups, specified by average_cols, and
averaging numeric values specified as the response variable or variables
extracted from a formula. The specified model is then fit to this
averaged data, and the predicted values are joined back up to the
original data frame. Let’s look at an example using blood pressure data,
which has a comprehensive time series.
library(augury)
df <- ghost::gho_data("BP_04", query = "$filter=Dim1 eq 'MLE' and Dim2 eq 'YEARS18-PLUS'") %>%
billionaiRe::wrangle_gho_data() %>%
dplyr::right_join(covariates_df) %>%
dplyr::select(iso3, year, year_n, value) %>%
dplyr::filter(whoville::is_who_member(iso3), # keep WHO member states
year >= 2000, year <= 2023) %>% # get relevant years
dplyr::mutate(who_region = whoville::iso3_to_regions(iso3)) # get WHO regions
#> Warning: Some of the rows are missing a source value.
#> Joining, by = c("iso3", "year")
ur <- unique(df$who_region)
ur
#> [1] "EMR" "AFR" "EUR" "AMR" "WPR" "SEAR"
Alright, so, here we have 6 WHO regions. We will use these regions to fit a model to and use INLA to predict out to 2023, then apply these trends to input countries.
pred_df <- df %>%
predict_inla_avg_trend(formula = value ~ f(year_n, model = "rw2"),
average_cols = c("who_region", "year_n"),
group_models = TRUE,
group_col = "iso3",
sort_col = "year_n")
pred_df %>%
dplyr::filter(iso3 == "AFG", year >= 2013)
#> # A tibble: 11 × 10
#> iso3 year year_n value who_region pred pred_upper pred_lower upper lower
#> <chr> <int> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 AFG 2013 14 30.4 EMR 29.4 29.4 29.4 NA NA
#> 2 AFG 2014 15 30.4 EMR 29.4 29.4 29.4 NA NA
#> 3 AFG 2015 16 30.4 EMR 29.4 29.4 29.4 NA NA
#> 4 AFG 2016 17 29.4 EMR 29.4 29.4 29.4 NA NA
#> 5 AFG 2017 18 29.4 EMR 29.4 29.4 29.4 NA NA
#> 6 AFG 2018 19 29.4 EMR 29.4 29.4 29.4 NA NA
#> 7 AFG 2019 20 29.4 EMR 29.4 29.4 29.4 NA NA
#> 8 AFG 2020 21 29.4 EMR 29.4 29.4 29.4 NA NA
#> 9 AFG 2021 22 29.4 EMR 29.4 29.4 29.4 NA NA
#> 10 AFG 2022 23 29.4 EMR 29.4 29.4 29.4 NA NA
#> 11 AFG 2023 24 29.4 EMR 29.4 29.4 29.4 NA NA
Above, we can see we have a generated a model using 2nd order random
walk with INLA, however, the model was generated by averaging data to
WHO regions first, fitting the random walk to each reach (since
group_models = TRUE) and then fitting those trends to the original
data. Note some specifics of what had to be set, as the
predict_..._avg_trend() functions are slightly more complex than
others:
average_colsmust contain thesort_col. So, since we useyear_nin the time series rather thanyear, we will sort by that this time.average_colsrefers to the groupings used for averaging (we take the average for each WHO region and year in this case). Then, the model is fit toaverage_colsthat are NOT thesort_col.group_colis the groupings used on the original data frame, which is still necessary here when applying the trend back to the original data.- If a variable is in
formula, it must either be inaverage_colsor it must be a numeric column that can be averaged. This is because the formula is applied to the data frame afterdplyr::group_by()anddplyr::summarize()have reduced it.
To highlight this point, in the above example, what’s actually happening is we are actually fitting the model on the summarized data:
df %>%
dplyr::group_by(who_region, year_n) %>%
dplyr::summarize(value = mean(value, na.rm = T)) %>%
head()
#> `summarise()` has grouped output by 'who_region'. You can override using the `.groups` argument.
#> # A tibble: 6 × 3
#> # Groups: who_region [1]
#> who_region year_n value
#> <chr> <dbl> <dbl>
#> 1 AFR 1 29.7
#> 2 AFR 2 29.6
#> 3 AFR 3 29.5
#> 4 AFR 4 29.4
#> 5 AFR 5 29.3
#> 6 AFR 6 29.2
Since average_cols = c("who_region", "year_n"), we took the mean of
all values in formula not in average_cols, in this case just
value. If for instance, we tried to specify a model using iso3 in
the formula:
predict_inla_avg_trend(df,
formula = value ~ iso3 + f(year_n, model = "rw2"),
average_cols = c("who_region", "year_n"),
group_models = TRUE,
group_col = "iso3",
sort_col = "year_n")
#> Error: iso3 must be numeric columns for use in averaging, or included in `average_cols` for grouping.
We get an error message indicating that iso3 must be numeric or included
in average_cols for grouping. This is because without it being numeric
or in the average_cols, there’s no way
dplyr::group_by() %>% dplyr::summarize() a non-numeric column
automatically (how would we reduce country-level ISO3 codes to the
regional level?).
While slightly complex, ensuring you follow the above means you should easily and successfully get out meaningful trend predictions for your data frames using trends generated on grouped data.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.