STAT 417 - Nathaniel Cinnamon Alex Hellkamp Colin Barrett

---
title: "League Data"
author: "Alex Hellkamp, Nathaniel Cinnamon, Colin Barrett"
date: "3/19/2019"
output: html_document
---

```{r}
library(tidyverse)
library(survival)
library(rebus)
library(survminer)
```

```{r}
leagueRaw <- read.csv("D:/Desktop/League_Data_2018.csv", header = TRUE)

str(leagueRaw)
head(leagueRaw, 40)
```

```{r}
#Manipulating data set into new data frame leagueWork

leagueWork <- leagueRaw %>%

  #Inlcude only winning teams
  filter(result == 1) %>%

  #Select only wanted variables
  select(gameid, position, side, herald, heraldtime, goldat15, csat15, xpat10) %>%

  #Change herald varibale from 3 values (team got it, enemy got it, no one got it) to 2 values (team got it,   team didn't get it)
  #Also replaced blank herald times with the maximum observation time (20 min)
  mutate(herald = case_when(herald == 1 ~ 1, TRUE ~ 0), heraldtime = ifelse(is.na(heraldtime), 20,  heraldtime)) %>%

  #Select only the team observations, getting rid of player observations
  filter(position == "Team") %>%

  #Remove unneccesary variables
  select(-position)

#Check to make sure data table looks fine
str(leagueWork)
head(leagueWork, 50)
```

```{r}
#Finding Descriptive Statistics of Explanatory Variables

  #Max, Min Median, and IQR for Gold, XP, and CS

  #Gold
  print(max(leagueWork$goldat15))
  print(min(leagueWork$goldat15))
  print(median_iqr(leagueWork$goldat15))

  #XP
  print(max(leagueWork$xpat10))
  print(min(leagueWork$xpat10))
  print(median_iqr(leagueWork$xpat10))

  #CS
  print(max(leagueWork$csat15))
  print(min(leagueWork$csat15))
  print(median_iqr(leagueWork$csat15))

```

```{r}
#Making Kaplan-Meier Curves for Explanatory Variables

#Kapla-Meier Curves for Map Side
km.obj.mapside <- survfit(Surv(heraldtime, herald)~side, data = leagueWork)
plot(km.obj.mapside,lty=1:2,xlab="Time Until Rift Herald Acquired (min)", ylab = "Survival Probability", main = "Estimated Survival Probabilities of Blue Side vs Red Side")
legend(1, .8, c("Blue","Red"),lty=1:2)

#Kaplan-Meier Curves for Gold at 15 min
km.obj.goldat15 <- survfit(Surv(heraldtime, herald)~1, data = leagueWork)
cr.object.goldat15 <- coxph(Surv(heraldtime, herald)~goldat15, data = leagueWork)

pred.max.goldat15 <- data.frame(goldat15=max(leagueWork$goldat15), data = leagueWork)
adj.surv.max <- survfit(cr.object.goldat15, newdata = pred.max.goldat15)

print(pred.max.goldat15)

pred.min.goldat15 <- data.frame(goldat15=min(leagueWork$goldat15), data = leagueWork)
adj.surv.min <- survfit(cr.object.goldat15, newdata = pred.min.goldat15)

print(pred.min.goldat15)

pred.median.goldat15 <- data.frame(goldat15=median(leagueWork$goldat15), data = leagueWork)
adj.surv.median <- survfit(cr.object.goldat15, newdata = pred.median.goldat15)

print(pred.median.goldat15)

plot(adj.surv.max, main = "Survival Curves for Observed Max, Min, and Median Gold", xlab = "Time until acquiring Rift Herald (min)", ylab = "Survival Probability", lty=1, ylim=c(0,1))
lines(adj.surv.min, lty = 2)
lines(adj.surv.median, lty = 3)
legend(1,.8,c("Max (28728G)","Min (22114G)", "Median (24521G)"),lty = 1:3)


#Kaplan-Meier Curves for Experience at 10 min

km.obj.xpat10 <- survfit(Surv(heraldtime, herald)~1, data = leagueWork)
cr.object.xpat10 <- coxph(Surv(heraldtime, herald)~xpat10, data = leagueWork)

pred.max.xpat10 <- data.frame(xpat10=max(leagueWork$xpat10), data = leagueWork)
adj.surv.max <- survfit(cr.object.xpat10, newdata = pred.max.xpat10)

print(pred.max.xpat10)

pred.min.xpat10 <- data.frame(xpat10=min(leagueWork$xpat10), data = leagueWork)
adj.surv.min <- survfit(cr.object.xpat10, newdata = pred.min.xpat10)

print(pred.min.xpat10)

pred.median.xpat10 <- data.frame(xpat10=median(leagueWork$xpat10), data = leagueWork)
adj.surv.median <- survfit(cr.object.xpat10, newdata = pred.median.xpat10)

print(pred.median.xpat10)

plot(adj.surv.max, main = "Survival Curves for Observed Max, Min, and Median Experience", xlab = "Time until acquiring Rift Herald (min)", ylab = "Survival Probability", lty=1, ylim=c(0,1))
lines(adj.surv.min, lty = 2)
lines(adj.surv.median, lty = 3)
legend(1,.8,c("Max (20306xp)","Min (16780xp)", "Median (18924xp)"),lty = 1:3)


#Kaplan-Meier Curves for CS at 10 min

km.obj.csat15 <- survfit(Surv(heraldtime, herald)~1, data = leagueWork)
cr.object.csat15 <- coxph(Surv(heraldtime, herald)~csat15, data = leagueWork)

pred.max.csat15 <- data.frame(csat15=max(leagueWork$csat15), data = leagueWork)
adj.surv.max <- survfit(cr.object.csat15, newdata = pred.max.csat15)

print(pred.max.csat15)

pred.min.csat15 <- data.frame(csat15=min(leagueWork$csat15), data = leagueWork)
adj.surv.min <- survfit(cr.object.csat15, newdata = pred.min.csat15)

print(pred.min.csat15)

pred.median.csat15 <- data.frame(csat15=median(leagueWork$csat15), data = leagueWork)
adj.surv.median <- survfit(cr.object.csat15, newdata = pred.median.csat15)

print(pred.median.csat15)

plot(adj.surv.max, main = "Survival Curves for Observed Max, Min, and Median CS", xlab = "Time until acquiring Rift Herald (min)", ylab = "Survival Probability", lty=1, ylim=c(0,1))
lines(adj.surv.min, lty = 2)
lines(adj.surv.median, lty = 3)
legend(1,.8,c("Max (580 CS)","Min (446 CS)", "Median (521 CS)"),lty = 1:3)
```

```{r}

#determining statistical significance in survival curves with log-rank tests

survdiff(Surv(heraldtime,herald)~side, data = leagueWork)


```

```{r}
#Fitting full and reduced models
leagueCoxReduced <- coxph(Surv(heraldtime, herald) ~ side + goldat15 + csat15 + xpat10, data = leagueWork)

leagueCoxFull <- coxph(Surv(heraldtime, herald) ~ side + goldat15 + csat15 + xpat10 + csat15*xpat10 + csat15*goldat15, data = leagueWork)

#Model outputs
summary(leagueCoxReduced)
summary(leagueCoxFull)

#Doesn't violate PH assumptions
test.ph <- cox.zph(leagueCoxReduced)
test.ph
ggcoxzph(test.ph)

#No obvious outliers in deviance
ggcoxdiagnostics(leagueCoxReduced, type = "deviance", linear.predictions = FALSE, ggtheme = theme_bw())

#Martingale residuals might be showing a problem with csat15
mart <- residuals(leagueCoxReduced, type = "martingale")

par(mfrow = c(1,3))
plot(leagueWork$goldat15, mart, xlab = "Gold at 15 min", ylab = "Martingale Residuals", main = "Predictor: Gold at 15 min")
smooth.sres <- lowess(leagueWork$goldat15, mart)
lines(smooth.sres$x, smooth.sres$y, lty = 1)

plot(leagueWork$csat15, mart, xlab = "CS at 15 min", ylab = "Martingale Residuals", main = "Predictor: CS at 15 min")
smooth.sres <- lowess(leagueWork$csat15, mart)
lines(smooth.sres$x, smooth.sres$y, lty = 1)

plot(leagueWork$xpat10, mart, xlab = "XP at 10 min", ylab = "Martingale Residuals", main = "Predictor: XP at 10 min")
smooth.sres <- lowess(leagueWork$xpat10, mart)
lines(smooth.sres$x, smooth.sres$y, lty = 1)
```


Inclusion of interaction terms do not improve model based on PLRT and AIC comparison
```{r}
leagueCoxFull$loglik
leagueCoxReduced$loglik

partiallikelihoodTestStat <- (2 * (-287.5465 - -287.7593))
print(partiallikelihoodTestStat)
1-pchisq(partiallikelihoodTestStat, 2)

AIC(leagueCoxFull, leagueCoxReduced)
```