aov_analysis.Rmd

---
title: "Diagnostic atlas of variation"
output:
  html_document: default
  html_notebook: default
  pdf_document: default
editor_options:
  chunk_output_type: inline
---

# Introduction

```{r}
knitr::opts_chunk$set(cache = TRUE, echo = FALSE, warning = FALSE, message = FALSE)

library(dplyr)
library(ggplot2)
library(govstyle)
if(!require("factoextra"))install.packages("factoextra")
library(factoextra)
if(!require("FactoMineR"))install.packages("FactoMineR")
library(FactoMineR)

```





```{r}

aov_analysis <- readr::read_csv("aov_ccg.csv")

aov_analysis %>%
  group_by(map_no, Period) %>%
  count() %>%
  tidyr::spread(map_no, n, fill = 0) %>%
  knitr::kable(caption = "Table1: Available time periods for each indicator; CCG counts")

### need to check map 24 and map 25


aov_analysis %>%
  filter(map_no == "Map 24", Period == "2015/16") %>%
  distinct()

  





```


## Comparison of endoscopic procedures

I am going to compare rates of utilisation of radiological, endoscopic and related procedures - these are included in maps 16, 17, 19, 20 and 22.



```{r}
 aov_endo_latest <- aov_analysis %>%
  filter(map_no %in% c("Map 1", "Map 2","Map 16", "Map 17", "Map 19", "Map 20", "Map 22" )) %>%
  group_by(map_no) %>%
  do(tail(., 1)) %>%
  mutate(index = stringr::str_c(map_no,"-", Period)) %>%
  select(index)

aov_endo <- aov_analysis %>%
  mutate(index = stringr::str_c(map_no,"-", Period)) 

aov_endo %>%
  filter( index %in% aov_endo_latest$index, !`CCG name` == "ENGLAND") %>%
  select(map_no,  `CCG name`, Rate) %>%
  spread(map_no, Rate) %>%
  select(-`CCG name`) %>%
  pairs(pch = 20, panel = panel.smooth, main= "Scatter plot of CCG level scanning and endoscopy rates")


```


### Plot trends



```{r fig.height=8, fig.width=6, fig.cap= "Trends in rates of endoscopy and imaging"}
aov_endo %>%
  filter(!stringr::str_detect(Period, "Q"))%>%
  ggplot(aes(Period, Rate)) +
  geom_boxplot(notch = TRUE, fill = "#2B8CC4") +
  facet_wrap(~map, scales = "free", labeller = label_wrap_gen(30)) +
  theme_gov() +
  theme(axis.text.x = element_text(angle= 45, hjust= 1, size = rel(.9)), 
        strip.text.x = element_text(size = 8))




```


```{r fig.height=6, fig.width=6, fig.cap= "Trends in rates of endoscopy and imaging - quarterly data"}

aov_endo %>%
  filter(stringr::str_detect(Period, "Q"))%>%
  ggplot(aes(Period, Rate)) +
  geom_boxplot(notch = TRUE, fill = "#2B8CC4") +
  facet_wrap(~map, scales = "free", labeller = label_wrap_gen(30)) +
  theme_gov() +
  theme(axis.text.x = element_text(angle= 45, hjust= 1, size = rel(.9)), 
        strip.text.x = element_text(size = 8))

```

## Correlations

```{r}
## latest data
aov_latest <- aov_endo %>%
  filter( index %in% aov_endo_latest$index, !`CCG name` == "ENGLAND") %>%
  select(map_no,  `CCG name`, Rate) %>%
  spread(map_no, Rate) %>%
  select(-`CCG name`)

rownames(aov_latest) <- aov_latest$`CCG name`

library(corrr)

correlate(aov_latest)%>% 
  rearrange(method = "HC", absolute = TRUE) %>% 
  fashion()

aov_scale_wide <- purrr::map(aov_latest, scale) %>% data.frame() %>% na.omit()

mean(is.na(aov_scale_wide))

set.seed(123)

fviz_nbclust(aov_scale_wide, kmeans, method = "silhouette" )

aov_k <- kmeans(aov_scale_wide, 6, nstart = 25, iter.max = 15)


fviz_cluster(aov_k, aov_latest, ellipse.type = "norm") 

aov_pca <- PCA(aov_latest)

fviz_contrib(aov_pca, choice = "var", axes = 1)

fviz_contrib(aov_pca, choice = "var", axes = 2)


```