Out of Pocket Expenditure by Snakebite victims in Rural Ghana.Rmd

---
date: "`r Sys.Date()`"
author: "Ibrahim Duah Kwaku"
title: "Out of pocket expenditure by snakebite victims in Ghana"
bibliography: snakebites.bib
csl: citestyle.csl
documentclass: article
fontsize: 12pt
output: 
  officedown::rdocx_document:
    tables:
      style: Table
      layout: autofit
      width: 1.0
      caption:
       style: "Table Caption"
       pre: 'Table '
       sep: ': '
    plots:
      style: Normal
      align: left
      caption:
       style: "Image Caption"
       pre: 'Figure '
       sep: ': '
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F, fig.cap = TRUE, warning = F, message = F, error = F)
library(officedown)
library(officer)
library(tidyverse)
library(gtsummary)
library(scales)
library(magrittr)
library(ggpubr)
library(flextable)
# library(gt)
bites <- readRDS("bites 15Jun2022.rds")

theme_gtsummary_compact(font_size = 12)
```

# Table of content

<!---BLOCK_TOC--->

\newpage

# Abstract

\newpage

# Introduction

Worldwide, snakebite victims are often the vulnerable poor rural folks[@Williams2019; @Kasturiratne2017]. Inequitable distribution of healthcare facilities and the lack of appropriate care exacerbate the plight of snakebite victims. An estimated 125,000 to 500,000 cases of snakebites occur annually in Africa even though most of them are thought to be undocumented[@Williams2011]. About 95,000 deaths occur globally annually from snakebites with only about 300,000 survivors[@TheLancet2017]. The survivors are often left with permanent disabilities or disfigurement and are often left stigmatized and destitute[@TheLancet2019]. Life-threatening effects of snakebite envenoming include shock, spontaneous systemic bleeding, paralysis involving respiratory and skeletal muscles, and can also lead to acute renal failure[@WHO2017]. Amputations, disfigurement, mutilations, and tissue necrosis are common complications of snakebite envenoming [@Gutierrez2006; @Habib2013]. The high morbidity and mortality associated with the bites result in high socioeconomic consequences[@WHO2017] for individuals and families.

Victims of snakebites require a range of services, from antivenom administration to supplementary medical interventions such as cardio-respiratory and/or fluid resuscitation; airway intubation; mechanical ventilation; hemodialysis; wound debridement and reconstructive surgery; physiotherapy; and other rehabilitation services[@WHO2017]. Unfortunately, these services are not usually available in primary care health settings (Community-based Health Services and Planning) in rural Ghana where a lot of the cases occur. In situations where they are available, they are not cheap. Snakebite victims therefore often start by going to traditional healers or use ineffective or unproven methods before seeking hospital care, resulting in delays in the administration of antivenoms which results in complications and possible mortality[@Chippaux2017; @Fry2018].

The socioeconomic impact of snakebite is under-appreciated around the world, even though the impacts transcend individuals and families into communities and the health systems[@Habib2018]. It is estimated the burden to some families to be as much as their 12-year salary [@Vaiyapuri2013]. The average cost for an effective treatment, based on recommended doses is USD 124 in Sub-Saharan Africa @Brown2012, 4 times the minimum monthly wage in Ghana then[@Africapay.org]. In Sri Lanka, 79% of snakebite victims suffered an economic loss following a snakebite with a median Out-of-Pocket Payments(OOP) of USD 11 and a median estimated loss of income of USD 28.57 and USD 33.21 for those in employment or self-employment, respectively[@Kasturiratne2017]. The total estimated OOP in the country was USD 1,981,699[@Kasturiratne2017]. Additionally, family members also lost income to help care for patients. In India, 53.5% of snakebite victims spent 1 to 6 months or more at home after the bite, 43.5% of the victims had to sell an asset due to snakebites, with the majority having to sell their farm crops. Four of the victims had to forfeit their education because of the bite, an unfortunate incident that must not happen. The annual estimated total number of DALYs was 11,101 to 15,076 per year for envenoming following snakebite[@Vaiyapuri2013].

The health system cost for Sri Lanka is estimated to be USD 10,260,652 annually[@Kasturiratne2017]. Using the conservative estimate from @Brown2012, then multiplying by the average yearly incidence in Ghana (9600)[@Ardayfio2020], it can be estimated that the government of Ghana spends an average of USD 1,190,400 on antivenoms since they are free in Ghana.

To the best of our knowledge, no studies exist in sub-Saharan Africa which report the OOP experienced by snakebite victims. The study, therefore, reports the direct OOP by snakebite victims in rural Ghana.

# Methods

## Ethics statement

This study has ethical approval from the Ghana Health Service. The reference number for it is GHS-ERC010/03/20. It is part of the Snakebite Incidence Treatment and Effect in Ghana (SnakebITE) project being run by the author[]. Permission was also sought from the administration of the hospital to extract records from their electronic health records (EHR). 

## Data extraction

The OOP by snakebite victims data were extracted from the electronic health records database of the hospital. The hospital has used the system since 2015. Author[] can script in Transactional-Structured Query Language (T-SQL) and extracted the records from the database in the presence of the hospital I.T/Health Information Officer. A total of `r nrow(bites)` were retrieved from the database.

Currency conversion was done with data from [@www.focus-economics.com]. The average  yearly rate was taken from the website and used to compute the OOP for each payment made by the patient while on the admission for the snakebite diagnosis.

## Data analysis

Statistical analysis was done with the R Statistical Software (version 4.2) with the aid of third party packages @gtsummary^,^ @tidyverse^,^,@r42^,^,@officer^,^@ggpubr^,^@flextable. Frequencies and percentages were recorded for the count variables and the median [interquartile range (IQR)] was reported for the continuous variables. A Generalized Linear Model (GLM) was fit on the data to determine the predictors of OOP among the victims of snakebite at the hospital. The payment method (Insured/Cash & Carry: A term in Ghana for patients without an insurance), year, outcome, length of stay was selected through stepwise backwards elimination method until such a time the model AIC failed to reduce significantly. Variable interaction was done for the year and length of stay, and the payment method and the year variables.

# Results and analysis

```{r echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
dem <- bites %>% 
  select(age, gender, year, seas, nhis, occup, act, outcome, los) %>% 
   tbl_summary(type = list(nhis ~ "categorical"))
```

We report cases of snakebites reported in a public hospital from 2016 to 2019 (Table \@ref(tab:tbldem)). A total of `r comma(nrow(bites))` records were retrieved from the EHR. The median (IQR) age of the victims was `r inline_text(dem, variable = age)`. Most, [`r inline_text(dem, variable = gender, level = "Male")`], of the snakebite victims were males. Most of the cases were reported in 2018 `r inline_text(dem, variable = year, level = "2018")` than other years. The primary occupation of most of the victims was farming `r inline_text(dem, variable = occup, level = "Farmer")`, `r inline_text(dem, variable = occup, level = "Student")` were students, and `r inline_text(dem, variable = occup, level = "Other")` engaged in other occupations. The primary occupation of  `r inline_text(dem, variable = occup, level = "Unknown")` of the victims were not indicted in the EHR. A total of `r inline_text(dem, variable = seas, level = "Dry")` were reported in the harmattan season and `r inline_text(dem, variable = act, level = "Irrigation/hunting")` was recorded during periods where there was rather little farming activity and mainly irrigation farming and hunting. Most of the victims, [`r inline_text(dem, variable = outcome, level = "Discharged Alive")`], were treated successfully at the hospital and were discharged. The median (IQR) length of stay of the victims at the hospital was `r inline_text(dem, variable = los)`.


```{r tab.id="tbldem", tab.cap="Demographic characteristics of snakebite victims extracted from EHR"}
dem %>% 
  as_flex_table()
```

## OOP for snakebite victims
```{r}
pay.smry <- bites %>% 
  tbl_continuous(variable = paid, include = c(mode, gender, year, seas, act, outcome)) 

pay.sum <- bites %>% 
  tbl_continuous(variable = paid, include = c(mode, gender, year, seas, act, outcome), statistic = list(everything() ~ "{sum}")) 

pay.ttl <- bites %>% 
  tbl_continuous(variable = paid, include = c(gender, year, seas, act, outcome),  by = mode)

oop.tbl <- tbl_merge(
  list(
    tot = pay.sum,
    smry = pay.smry,
    ttl = pay.ttl
  ),
  tab_spanner = c("Total", "Median(IQR)", "Breakdown")
) 
```

The total OOP from 2016 to 2019 at the hospital was USD `r sum(bites$paid) %>% comma()` (Table \@ref(tab:pay-smry)) of which USD `r inline_text(oop.tbl, variable = mode, level = "NHIS", pattern = "{stat_0_1}")` were payments made by victims that had an insurance cover (presumed to be the Nation Health Insurance Service [NHIS]) at the time of admissions and USD `r inline_text(oop.tbl, variable = mode, level = "Cash & Carry", pattern = "{stat_0_1}")` were paid by victims without insurance cover. The median (IQR) amount paid by NHIS clients was USD `r inline_text(oop.tbl, variable = mode, level = "NHIS", pattern = "{stat_0_2}")` compared to a median (IQR) of USD `r inline_text(oop.tbl, variable = mode, level = "Cash & Carry", pattern = "{stat_0_2}")` by non-insured clients. The total OOP by males was USD `r inline_text(oop.tbl, variable = gender, level = "Male", pattern = "{stat_0_1}")` compared to USD `r inline_text(oop.tbl, variable = gender, level = "Female", pattern = "{stat_0_1}")` among females. The median (IQR) OOP between gender was relatively similar even though it was higher in males [USD `r inline_text(oop.tbl, variable = gender, level = "Male", pattern = "{stat_0_2}")`] compared to females [USD `r inline_text(oop.tbl, variable = gender, level = "Female", pattern = "{stat_0_2}")`]. When the insurance cover of the victims was taken into account, uninsured victims paid more than 4 times as much as those with insurance paid, with very little differences between males and females. From 2016 to 2018, there was a steady increase in the total OOP for snakebite victims at the hospital. However, the total OOP almost tripled, from USD `r inline_text(oop.tbl, variable = year, level = "2018", pattern = "{stat_0_1}")` in 2018 to USD `r inline_text(oop.tbl, variable = year, level = "2019", pattern = "{stat_0_1}")` in 2019. In response, the median (IQR) increased by about a factor of 4 from USD `r inline_text(oop.tbl, variable = year, level = "2018", pattern = "{stat_0_2}")` in 2018 to USD `r inline_text(oop.tbl, variable = year, level = "2019", pattern = "{stat_0_2}")` in 2019. The situation was much more dire for non-insured clients when the median(IQR) increased from USD `r inline_text(oop.tbl, variable = year, level = "2018", pattern = "{stat_1_3}")` in 2018 to USD `r inline_text(oop.tbl, variable = year, level = "2019", pattern = "{stat_1_3}")` in 2019. The insured clients were not spared the surge, paying a median(IQR) of USD `r inline_text(oop.tbl, variable = year, level = "2018", pattern = "{stat_2_3}")` in 2018 to USD `r inline_text(oop.tbl, variable = year, level = "2019", pattern = "{stat_2_3}")`. It can be noted that the absolute difference in the medians between the insured and uninsured victims

For a successful treatment, victims without insurance cover  paid a median(IQR) of USD `r inline_text(oop.tbl, variable = outcome, level = "Discharged Alive", pattern = "{stat_1_3}")` compared to USD `r inline_text(oop.tbl, variable = outcome, level = "Discharged Alive", pattern = "{stat_2_3}")` who had insurance cover.

```{r payment-summary, tab.id="pay-smry", tab.cap = "Summary of OOP experienced by snakebite victims"}
oop.tbl %>% 
  as_flex_table()
```

### OOP by services received

```{r trend-of-oop, eval=FALSE,message=FALSE, warning=FALSE, fig.width=8, fig.height=6, fig.cap="Scatterplot of OOP and age (L) and LoS (R)", fig.id = "oop-trend"}
age.pay <- bites %>% 
  ggscatter(x = "age", y = "paid", color = "mode", alpha = 0.4, ylab = "Amount paid", xlab = "Age", legend.title = "Gender")

los.pay <- bites %>% 
  ggscatter(x = "los", y = "paid", color = "mode", alpha = 0.4, ylab = "Amount paid", xlab = "Length of Stay", legend.title = "Gender")

ggarrange(age.pay, los.pay, common.legend = T) 
```

```{r oop-trend, include=FALSE,warning=FALSE, message=FALSE, fig.cap = "Trend of OOP between victims with and without insurance cover at the time of the bite", fig.id="oop-trend"}
pay.smry <- bites %>%
  group_by(year, mode) %>%
  summarise(
    median = median(paid),
    sec = quantile(paid)[2],
    third = quantile(paid)[4]
  ) %>%
  pivot_longer(c(median, sec,third)) %>%
  mutate(
    name = fct_relevel(name, "third", "median", "sec"),
    name = fct_recode(name, "75%" = "third", "Median" = "median", "25%" = "sec")
  ) 
pay.smry %>%
  ggline(x = "year", y = "value", color = "name", facet.by = "mode", legend.title = "Percentile", xlab = "Year", ylab = "OOP amount (USD)")
```

```{r include=FALSE}
pay.yr.gp <- bites %>%
  ggdotplot(x = "year", y = "paid", fill = "mode", legend.title = "Payment", alpha = .35, ylab = "OOP (USD)", xlab = "Year")
```

```{r cost-summary, include=F}
my_stats <- function(data, full_data, ...) {
  sum_salary <- sum(data$cost, na.rm = TRUE)
  total_sum <- sum(full_data$cost, na.rm = TRUE)
  dplyr::tibble(
    sum = sum_salary,
    total_sum = total_sum,
    p = round(sum_salary / total_sum, 3)
  )
}

mydat <- bites %>%
  select(id, cons:svs) %>%
  sjlabelled::drop_labels() %>%
  pivot_longer(c(cons:svs), names_to = "item", values_to = "cost") %>%
  select(-id) %>%
  mutate(item = fct_recode(
    item,
    "Consultation" = "cons",
    "Drugs" = "drugs",
    "Services" = "svs"
  ))

t1 <- mydat %>%
  tbl_custom_summary(
    include = item,
    stat_fns = ~ my_stats,
    statistic = ~ "{sum} ({p}%)",
    label = list(item ~ "Item")
  )

pay.tbl <-  tbl_merge(list(
  t1,
  mydat %>%
    tbl_continuous(
      variable = cost,
      label = list(item ~ "Item"),
      statistic = everything() ~ "{median} ({p25}, {p75})"
    )
  ),
  tab_spanner = c("Sum(%)", "Median(IQR)")
) %>% 
  modify_footnote(
    update = list(stat_0_2 ~ "Median (IQR")
  )
```

The total (%) OOP on accounts of drugs purchases were USD `r inline_text(pay.tbl, variable = item, level = 'Drugs', column = stat_0_1)` of all OOP with a median (IQR)= USD `r inline_text(pay.tbl, variable = item, level = 'Drugs', column = stat_0_2)` (Table \@ref(tab:svs-smry)). This was followed by the provision of services (e.g., ward admissions, wound dressing, x-rays, labs, etc.)  accounting for USD `r inline_text(pay.tbl, variable = item, level = 'Services', column = stat_0_1)`, with a median (IQR) = `r inline_text(pay.tbl, variable = item, level = 'Drugs', column = stat_0_2)`. Lastly the total OOP [(%), Median (IQR)] for consultations was USD `r inline_text(pay.tbl, variable = item, level = Consultation, pattern = "{stat_0_1}, {stat_0_2}")` of all OOP.

```{r tab.id = "svs-smry", tab.cap = "OOP by services received", cache=FALSE}
  pay.tbl %>% 
  as_flex_table()
```

#### Payments by mode of Insurance status
```{r include=FALSE}
p1 <- bites %>%
  select(paid, cons:exems, mode) %>%
  tbl_summary(by = mode, digits = list(everything() ~ 1))

p2 <- bites %>%
  select(paid, cons:exems, mode) %>%
  tbl_summary(
    by = mode,
    statistic = list(everything() ~ "{sum}")
  ) %>%
  add_overall(last = T,
              col_label = "Total")

ftbl <- tbl_merge(list("Maximum" = p2, "Median" = p1 ), tab_spanner = c("Total", "Summary")) 
```

As indicated earlier, the total OOP at the hospital was USD `r inline_text(ftbl, variable = paid, pattern="{stat_0_1}")`. Of the total OOP,  USD `r inline_text(ftbl, variable = paid, pattern="{stat_1_1}")` were paid by victims without an insurance cover at the time of admission. The median (IQR) for clients with an insurance cover was USD `r inline_text(ftbl, variable = paid, pattern="{stat_1_2}")` and `r inline_text(ftbl, variable = paid, pattern="{stat_2_2}")` those without an insurance (Table \@ref(tab:exptbl)). Victims without a health insurance cover paid a total of USD `r inline_text(ftbl, variable = cons, pattern="{stat_1_1}")` for consultation compared to USD `r inline_text(ftbl, variable = cons, pattern="{stat_2_1}")` among victims with the NHIS, bringing it to a total of `r inline_text(ftbl, variable = cons, pattern="{stat_0_1}")`. The median payment for consultation for victims without insurance cover was USD `r inline_text(ftbl, variable = cons, pattern="{stat_1_2}")` compared to a median(IQR) of USD `r inline_text(ftbl, variable = cons, pattern="{stat_2_2}")` among victims with insurance cover. The hospitals exempted an amount of USD 140 for victims with an insurance cover, reasons for which were not stated in the EHR.

```{r echo=FALSE, tab.cap="Summary of payments by snakebite victims", tab.id="exptbl",  fig.align='center', warning=FALSE, message=FALSE, error=FALSE}
ftbl %>% 
  as_flex_table()
```

## Analysis of data
```{r include=FALSE}
mod <- glm(formula = paid ~ mode*year + outcome + year  + los*mode, data = bites)
out.mod <- tbl_regression(mod, label = list(year ~ "Year"))
```

The fitted regression model shows victims with an insurance cover (NHIS) at the time of the bite experienced a statistically significant average of USD `r inline_text(out.mod, variable = mode, level = 'NHIS')` lesser OOP compared to the victims with an insurance cover. Compared to 2016, the average OOP expenditure in 2017 was lower by an average of USD `r inline_text(out.mod, variable = year, level = "2017")`. The average OOP further declined in 2018 to USD `r inline_text(out.mod, variable = year, level = "2018")` but raised sharply to USD `r inline_text(out.mod, variable = year, level = "2019")` in 2019. The average differences from 2018 to 2019 were statistically significant. Victims that died spent USD `r inline_text(out.mod, variable = outcome, level = 'Discharged Alive')` less compared to those that were discharged alive from the hospital. The average OOP increased by USD `r inline_text(out.mod, variable = los)` for each day on admission.

The interaction of the variables, however, showed a different phenomenon. Victims with an insurance cover spent averagely less than those with an insurance cover in 2017, USD`r inline_text(out.mod, variable = "mode:year", level = "NHIS * 2017")`. Compared to the same reference year (2016), OOP expenditure by victims on the NHIS increased from USD `r inline_text(out.mod, variable = "mode:year", level = "NHIS * 2018")` in 2018 to USD `r inline_text(out.mod, variable = "mode:year", level = "NHIS * 2019")` in 2019, suggesting the NHIS may have exerted a negative influence of the OOP experienced by victims. However, the differences in the averages were not statistical significant from the OOP in 2016. Lastly, for each day on admission, the average OOP for victims with an insurance cover compared to those that did not was USD `r inline_text(out.mod, variable = "mode:los", level = "NHIS * Length of stay")`.

```{r pay-regression, tab.cap="Regression model on predictors of OOP at the hospital", tab.id="reg-mod", message=FALSE, warning=FALSE}
out.mod %>% 
  as_flex_table()
```

## Discussions

The objective of the study was to describe the OOP experienced by victims of snakebites visiting a rural hospital in Ghana using the records available in the Electronic Health Records (EHR) of the hospital from 2016 to 2019. To the best of our knowledge, this is the first of its kind in Ghana and Sub-Saharan Africa. A total of `r nrow(bites)` records retrieved from the EHR and the total OOP of snakebite victims was USD `r sum(bites$paid) %>% comma()`. The victims were mostly males and younger.

```{r}
dif <- bites %>% group_by(mode) %>% 
  summarise(
    pay = median(paid)
  ) %>% pull(pay)

dif.dd <- bites %>% select(mode, paid)

dif.test <- wilcox.test(paid ~ mode, data = dif.dd) %>% broom.mixed::tidy()

sum.paid <- bites %>% 
  select(paid) %>% 
  tbl_summary()
```

This study has demonstrated that snakebite victims in rural Ghana experience a significant OOP. The median (IQR) OOP for the victims was USD `r inline_text(sum.paid, variable = paid)` which is very high given the average monthly income in Ghana was around USD `r 800/7` at the time. Household income in  rural areas is certainly lower than the national average. Given that the major economic activity within the communities is farming, the OOP may have pushe some households into catastrophic health expenditure. 

The victims without an insurance experiencing the worst median OOP expenditure `r median(bites$paid) %>% comma(accuracy = 0.1)` which compares very similarly to OOP in Kenya as reported by @Okumu2019 but harshly with OOP in Sri Lanka according to @Kasturiratne2017. Most of the victims did not have a health insurance cover at the time of the bite and unsurprisingly spent `r (dif[1] - dif[2]) %>% comma()` a median more than their colleagues with an insurance. The median difference in OOP between those with insurance and those without was statistically significant (`r pvalue(dif.test[2] %>% pull())`). 

```{r}
yrly <- bites %>% 
  group_by(year) %>% 
  summarise(pay = median(paid))

yrly.stat <-
  bites %$% 
  glm(paid ~ year*mode) %>% 
  tbl_regression(label = year ~ "Year", add_estimate_to_reference_rows = T) %>% 
  add_global_p()
```

The median cost of care for snakebite victims increased with each year (2016 to 2019). The year of the incidence exhibited a strong linear relationship with OOP `r inline_text(yrly.stat, variable = year, pattern = "{p.value}")`. However, when the insurance status of the victim was taken into account, the statistical difference was no more (`r inline_text(yrly.stat, variable = 3, pattern = "{p.value}")`) demonstrating that the NHIS fails to protect snakebite victims of OOP. It is the opinion of the authors that the surge in OOP in 2019 was due to the hospital reviewing the cost of their drugs. This further suggests the NHIS has played a little role in alleviating OOP among snakebite victims year-on-year.

The biggest driver of OOP among the victims was drugs, accounting for USD `r inline_text(pay.tbl, variable = 1, level="Drugs", pattern="{stat_0_1}")` of the total cost. OOP for services (including laboratory services) followed as the second highest cost. The NHIS did not provide a significant cover for snakebite victims. For consultation fees, NHIS beneficiaries paid USD -2.65 less than their counterparts that did not have an insurance but this difference was not statistically significant. For drugs payments, NHIS clients paid a statistically significant USD -6.28 less than victims with an insurance coverage. It was payment for services, e.g. nursing care, wound dressing, etc where NHIS victims enjoyed the most benefit. NHIS victims paid almost USD 17 less than their counterparts without an insurance.The table below can have 2 plausible interpretations; the NHIS did not provide a significant cover NHIS victims but drugs costs accounted for more than 50% of the total cost of care and it would appear NHIS victims were spared just about USD 7. Another plausible explanation would be that since antivenoms are supposed to be free in Ghana, irrespective of NHIS status, ASV cost was waived for all victims and therefore the NHIS would contribute less to financial protection for NHIS clients.



```{r}
read_rds("nhis item cost-diff.rds")
```





# List of tables

<!---BLOCK_TOC{seq_id: 'tab'}--->
# List of figures

<!---BLOCK_TOC{seq_id: 'fig'}--->

# Reference