03-data-manipulation.Rmd

# Data manipulation

## Collapse multiple "no" and "yes" options

Common to have to do this in GlobalSurg projects:

```{r message=FALSE}
library(dplyr)
mydata = tibble(
  ssi.factor = c("No", "Yes, no treatment/wound opened only (CD 1)",    
                 "Yes, antibiotics only (CD 2)", "Yes, return to operating theatre (CD 3)", 
                 "Yes, requiring critical care admission (CD 4)", 
                 "Yes, resulting in death (CD 5)",
                 "Unknown") %>%
    factor(),
  mri.factor = c("No, not available", "No, not indicated", 
                 "No, indicated and facilities available, but patient not able to pay",
                 "Yes", "Unknown", "Unknown", "Unknown") %>% 
    factor()
)

# Two functions make this work
fct_collapse_yn = function(.f){
  .f %>% 
    forcats::fct_relabel(~ gsub("^No.*", "No", .)) %>% 
    forcats::fct_relabel(~ gsub("^Yes.*", "Yes", .))
}

is.yn = function(.data){
  .f = is.factor(.data)
  .yn = .data %>% 
    levels() %>% 
    grepl("(No.+)|(Yes.+)", .) %>% 
    any()
  all(.f, .yn)
}

# Raw variable
mydata %>% 
  count(ssi.factor)

# Collapse to _yn version
mydata %>% 
  mutate(across(where(is.yn), fct_collapse_yn, .names = "{col}_yn")) %>% 
  count(ssi.factor_yn)

```

## Filtering best practice
### From Jamie Farrell 

Particularly useful in OpenSAFELY, but should be adopted by us all. 

This creates inclusion flags, and then summarises those inclusion flags in a table. This allows us to see exactly how many rows are removed at each filter step. 

The inclusion flags are then used to run the final filter. 
```{r message=FALSE, warning=FALSE}
library(finalfit)
library(tidyverse)

# Utility function ----
## Ignore, just needed for the example. 
fct_case_when <- function(...) {
  # uses dplyr::case_when but converts the output to a factor,
  # with factors ordered as they appear in the case_when's  ... argument
  args <- as.list(match.call())
  levels <- sapply(args[-1], function(f) f[[3]])  # extract RHS of formula
  levels <- levels[!is.na(levels)]
  factor(dplyr::case_when(...), levels=levels)
}

# Load data
colon_s = colon_s

# Create inclusion flags:
# Age over 50, male, obstruction data not missing
colon_s = colon_s %>% 
  mutate(
    is_age_abover_50 = age > 50,
    is_male = sex.factor == "Male",
    obstruction_not_missing = !is.na(obstruct.factor)
  )

# Create exclusion table
flowchart = colon_s %>% 
  transmute(
  c0 = TRUE,
  c1 = c0 & is_age_abover_50,
  c2 = c1 & is_male,
  c3 = c2 & obstruction_not_missing,
) %>%
  summarise(
    across(.fns=sum)
  ) %>%
  mutate(pivot_col = NA) %>% 
  pivot_longer(
    cols=-pivot_col,
    names_to="criteria",
    values_to="n"
  ) %>% 
  select(-pivot_col) %>% 
  mutate(
    n_exclude = lag(n) - n,
    pct_all = (n/first(n)) %>% scales::percent(0.1),
    pct_exclude_step = (n_exclude/lag(n)) %>% scales::percent(0.1),
    crit = str_extract(criteria, "^c\\d+"),
    criteria = fct_case_when(
      crit == "c0" ~ "Original dataset",
      crit == "c1" ~ "  with age over 50 years",
      crit == "c2" ~ "  is male",
      crit == "c3" ~ "  obstruction data not missing",
      TRUE ~ NA_character_
    )
  ) %>% 
  select(criteria, n, n_exclude, pct_all, pct_exclude_step)
flowchart
```

Now filter the dataset in one step. 

```{r message=FALSE, warning=FALSE}
# Filtered dataset
colon_filtered = colon_s %>% 
  filter(
    is_age_abover_50,
    is_male,
    obstruction_not_missing
  )
```


## Filter NA: Dropping rows where all specified variables are NA


I want to keep rows that have a value for `important_a` and/or `important_b` (so rows 1, 3, 4).
I don't care whether `whatever_c` is empty or not, but I do want to keep it.

```{r, message=FALSE}
library(tidyverse)

mydata  = tibble(important_a = c("Value", NA, "Value", NA, NA),
                 important_b = c(NA, NA, "Value", "Value", NA),
                 whatever_c  = c(NA, "Value", NA, NA, NA))

mydata
```



Functions for missing values that are very useful, but don't do what I want are:

(1) This keeps complete cases based on all columns:
```{r}
mydata %>% 
  drop_na()
```

(Returns 0 as we don't have rows where all 3 columns have a value).

(2) This keeps complete cases based on specified columns:

```{r}
mydata %>% 
  drop_na(important_a, important_b)
```
This only keeps the row where both a and b have a value.


(3) This keeps rows that have a value in any column:

```{r}
mydata %>% 
  filter_all(any_vars(! is.na(.)))
```


The third example is better achieved using the janitor package:

```{r, message = FALSE}
mydata %>% 
  janitor::remove_empty()
```

Now, (3) is pretty close, but still, I'm not interested in row 2 - where both a and b are empty but c has a value (which is why it's kept).

(4) Simple solution

A quick solution is to use `! is.na()` for each variable inside a `filter()`:

```{r}
mydata %>% 
  filter(! is.na(important_a) | ! is.na(important_b))
```

And this is definitely what I do when I only have a couple of these variables. But if you have tens, then the filtering logic becomes horrendously long and it's easy to miss one out/make a mistake.

(5) Powerful solution:

A scalable solution is to use `filter_at()` with `vars()` with a select helper (e.g., `starts with()`), and then the `any_vars(! is.na(.))` that was introduced in (3).

```{r}
mydata %>% 
  filter_at(vars(starts_with("important_")), any_vars(! is.na(.)))
```

## Vectorising rowwise procedures

We frequently have to aggregate across columns. `dplyr::rowwise()` is the `group_by()` equivalent for rows. But this is painfully slow for large datasets. 

The following works beautifully and allows tidyselect of variable names. 

```{r}
library(tidyverse)
mydata = tibble(a = c(1,2,NA),
       b = c(2,3,NA),
       c = c(3,4,NA),
       d = c(NA,NA,NA))

mydata %>% 
  mutate(s = pmap_dbl(select(., a:d), pmax, na.rm=TRUE))
```

## Multiple imputation and IPW for missing data

Two approaches are commonly used for dealing with missing data. 

Multiple imputation (MI) "fills in" missing values. Inverse probability weighting (IPW) "weights" observed values to reflect characteristics of missing data. 

Missing data is discussed in detail here:https://finalfit.org/articles/missing.html

The IPW methods are taken from here: https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/

### Data

```{r message=FALSE, warning=FALSE}
library(finalfit)
library(tidyverse)
library(mice)

## Inspect dataset
# colon_s %>% ff_glimpse()

## Below we will use these variables to model, 
## so remove missings from them and from outcome for this demonstration
explanatory = c("age", 
                "sex.factor", 
                "obstruct.factor", 
                "perfor.factor",
                "adhere.factor",
                "nodes",
                "differ.factor",
                "rx")


## Remove existing missings
colon_s = colon_s %>% 
  select(mort_5yr, explanatory) %>% 
  drop_na()
```

Here is the truth for mortality at 5 years in these data. 

```{r message=FALSE, warning=FALSE}
## Truth
colon_s %>% 
  summary_factorlist(explanatory = "mort_5yr")
```

Let's conditionally delete some outcomes. 

```{r}
set.seed(1234)
colon_rs = colon_s %>% 
  mutate(mort_5yr = ifelse(mort_5yr == "Alive", 
                           sample(c("Alive", NA_character_), 
                                  size = sum(colon_s$mort_5yr == "Alive"), 
                                  replace = TRUE, prob = c(0.9, 0.1)),
                           sample(c("Died", NA_character_), 
                                  size = sum(colon_s$mort_5yr == "Died"), 
                                  replace = TRUE, prob = c(0.85, 0.15))) %>% 
           factor()
  )
```

The proportion that died, is now lower due to our artifically introduced reporting bias:

```{r message=FALSE, warning=FALSE}
colon_rs %>% 
  summary_factorlist(explanatory = "mort_5yr")
```

### Multiple imputation
First, let's use multiple imputation to "fill in" missing outcome using the existing information we have about the patient. 

```{r}
## Choose which variables in dataset to include in imputation. 
predM = colon_rs %>% 
  select(mort_5yr, explanatory) %>% 
  missing_predictorMatrix(
    drop_from_imputed = c(""), 
    drop_from_imputer = c(""))

# This is needed to drop from imputed
# Include for completeness, but not used here. 
m0 = mice(colon_rs %>% 
            select(mort_5yr, explanatory), maxit=0)
# m0$method[c("study_id")] = ""
# m0$method[c("redcap_data_access_group")] = ""

colon_rs_imputed = colon_rs %>% 
  select(mort_5yr, explanatory) %>% 
  mice(m = 10, maxit = 10, predictorMatrix = predM, method = m0$method,
       print = FALSE) # print false just for example

## Check output here:
# summary(colon_rs_imputed)
# plot(colon_rs_imputed )
```

#### Reduce imputed sets (mean)
As can be seen, the imputed proportion approaches that of the "true" value. 

```{r}
colon_rs_imputed %>% 
  complete("all") %>% 
  map(~ (.) %>% count(mort_5yr) %>% 
        mutate(nn = sum(n),
               prop = n/nn) %>% 
        select(-mort_5yr)) %>% 
  reduce(`+`) / 10
```

### IPW
Rather than filling in the missing values, we can weight the observed values to reflect the characteristics of the missing data. 

```{r}
## Define observed vs non-observed groups
colon_rs = colon_rs %>% 
  mutate(
    group_fu = if_else(is.na(mort_5yr), 0, 1)
  )

colon_rs %>% 
  count(group_fu)
```

This method uses "stabilised weights", meaning weights should have a mean of 1. 

```{r}
# Numerator model
fit_num = colon_rs %>% 
  glmmulti("group_fu", 1)

# Denominator model
fit_dem = colon_rs %>% 
  glmmulti("group_fu", explanatory)

# Check denominator model
fit_dem %>% 
  fit2df()
```


```{r}
colon_rs = colon_rs %>% 
  mutate(
    p = predict(fit_num, type = "response"),
    pi = predict(fit_dem, type = "response"),
    weights = case_when(
      group_fu == 1 ~ p / pi,
      group_fu == 0 ~ (1-p) / (1-pi) 
    )
  )

mean(colon_rs$weights)
# Whoop whoop if ~1.0
```

A weighted proportion can then be generated, which gives a very similar result to the multiple imputation. 

```{r message=FALSE, warning=FALSE}
colon_rs %>%  
  mutate(
    mort_5yr_wgt = (as.numeric(mort_5yr) - 1) * weights
  ) %>% 
  summary_factorlist(explanatory = c(" mort_5yr", " mort_5yr_wgt"), digits = c(4, 4, 3, 1, 1) )

```