---
title: "Enrolled/NotEnrolled"
author: "Nathan Leon Pace, MD, MStat"
fontsize: 12pt
output:
  pdf_document: default
  word_document: default
  latex_engine: xelatex
geometry: margin = 0.5in
classoption:
- landscape
- legalpaper
---

This analysis was performed on `r date()`.  

This analysis used the R platform (`r citation()`).  

```{r setup, include = F}

library(knitr)
use.knitr = T
knitr::opts_chunk$set(echo = T)
opts_knit$set(root.dir = '..')
opts_chunk$set(cache=T)
options(max.print = 99999)
options(dplyr.print_min = Inf)
options(tibble.width = Inf)
options(width = 150)
options(tinytex.verbose = T)

```

```{r snapshot, include = F}

library(readit)
library(tinytex)
library(magrittr)
library(tidyverse)
library(MASS)
library(skimr)
library(ggplot2)
library(ggmosaic)
library(plotly)
library(jtools)

```

# Set skimmer arguments
```{r, include = F}

skim_with(numeric = list(hist = NULL))

skim_with(integer = list(hist = NULL))

skim_with(factor = list(ordered = NULL))

```

# Enrolled data
```{r, include = F}

RawEnrolled.df <- readit('Data/ForNathanFig1enrolled_20190418.xlsx', sheet = 1)

names(RawEnrolled.df) <- c('InstitutionName', 'THA', 'TKA', 'Spine', 'Chest',
                        'Breast', 'Abdomen')

Enrolled.df <- gather(RawEnrolled.df, 'THA', 'TKA', 'Spine', 'Chest',
                        'Breast', 'Abdomen', key = 'SurgeryType',
       value = 'cases')

Enrolled.df <- data.frame(Enrolled.df, Status = 'Enrolled')

Enrolled.df <- Enrolled.df %>%
  filter(InstitutionName != 'Utrecht' &
           InstitutionName != 'Vermont' &
           InstitutionName != 'Virginia' &
           InstitutionName != 'WashingtonStLouis') %>%
  mutate(InstitutionName = fct_drop(InstitutionName))

Enrolled.df <- arrange(Enrolled.df, InstitutionName, SurgeryType)

```

# Eligible data
```{r, include = F}

RawEligible.df <- readit('Data/ForNathanFig1eligible_20190418.xlsx', sheet = 1)

names(RawEligible.df) <- c('InstitutionName', 'THA', 'TKA', 'Spine', 'Chest',
                        'Breast', 'Abdomen')

Eligible.df <- gather(RawEligible.df, 'THA', 'TKA', 'Spine', 'Chest',
                        'Breast', 'Abdomen', key = 'SurgeryType',
       value = 'cases')

Eligible.df <- data.frame(Eligible.df, Status = 'Eligible')

Eligible.df <- Eligible.df %>%
  filter(InstitutionName != 'Utrecht' &
           InstitutionName != 'Vermont' &
           InstitutionName != 'Virginia' &
           InstitutionName != 'WashingtonStLouis') %>%
  mutate(InstitutionName = fct_drop(InstitutionName))

Eligible.df <- arrange(Eligible.df, InstitutionName, SurgeryType)

```

# NotEnrolled
```{r, include = F}

NotEnrolled.df <-
  left_join(Eligible.df,
            Enrolled.df,
            by = c('InstitutionName',
                   'SurgeryType')) %>%
  mutate(cases = cases.x - cases.y)

NotEnrolled.df <- data.frame(NotEnrolled.df$InstitutionName,
                             NotEnrolled.df$SurgeryType,
                      NotEnrolled.df$cases) %>%
  mutate(Status = 'NotEnrolled')

names(NotEnrolled.df) <- c('InstitutionName', 'SurgeryType', 'cases', 'Status')

```

# Combined data
```{r, echo = F}

Total.df <- bind_rows(NotEnrolled.df, Enrolled.df) %>%
  mutate(InstitutionName = as_factor(InstitutionName),
         SurgeryType = as_factor(SurgeryType),
         Status = as_factor(Status)) 

TallTotal.df <- uncount(Total.df, cases)

skim(Total.df)

```

\newpage
# Data Counts
##  Mosaic plot
##  The area of each box is proportional to the cell frequency.  
```{r}

TallTotal.gg <- ggplot(data = TallTotal.df) +
  geom_mosaic(aes(x = product(InstitutionName, Status), fill = SurgeryType), conds = product(SurgeryType), na.rm = T) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
  xlab('SurgeryType:Status') + ylab('InstitutionName')

TallTotal.gg

ggsave('Plots/TallTotal.pdf', width = 11, height = 8.5)

```

\newpage
# Loglinear model
```{r, include = F}

Total.loglm <- loglm(cases ~ Status + InstitutionName + SurgeryType, Total.df)

```

# Loglinear coefficients
```{r}

Total.loglm.coef <- coef(Total.loglm)

exp(Total.loglm.coef$Status)

exp(Total.loglm.coef$InstitutionName)

exp(Total.loglm.coef$SurgeryType)

rm(Total.loglm.coef)

```

# Poisson regression
```{r, echo = F, eval = T}

Total.glm <- glm(cases ~  Status + InstitutionName + SurgeryType, Total.df,
                 family = poisson(link = 'log'),
                 contrasts = list(Status = contr.sum,
                                  InstitutionName = contr.sum,
                                  SurgeryType = contr.sum))

summ(Total.glm)

```

\newpage
# Poisson regression
## Effect plot
##  Status
```{r, eval = T}

effect_plot(Total.glm, pred = Status, interval = T, plot.points = T, jitter = 0.2)

```

\newpage
# Poisson regression
## Effect plot
##  SurgeryType
```{r, eval = T}

effect_plot(Total.glm, pred = SurgeryType, interval = T, plot.points = T, jitter = 0.2)

```

\newpage
# Poisson regression  
## Effect plot  
##  InstitutionName  
```{r, eval = T}

effect_plot(Total.glm, pred = InstitutionName, interval = T, plot.points = T, jitter = 0.2)

```


```{r, include = F, eval = F}

save.image('Rdata/Enrolled.Rdata')

```