API tools faq
paste
Guest User

STAT 417 - Nathaniel Cinnamon Alex Hellkamp Colin Barrett

a guest
Mar 18th, 2019
130
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 7.11 KB | None | 0 0
raw download clone embed print report
  1. ---
  2. title: "League Data"
  3. author: "Alex Hellkamp, Nathaniel Cinnamon, Colin Barrett"
  4. date: "3/19/2019"
  5. output: html_document
  6. ---
  7.  
  8. ```{r}
  9. library(tidyverse)
  10. library(survival)
  11. library(rebus)
  12. library(survminer)
  13. ```
  14.  
  15. ```{r}
  16. leagueRaw <- read.csv("D:/Desktop/League_Data_2018.csv", header = TRUE)
  17.  
  18. str(leagueRaw)
  19. head(leagueRaw, 40)
  20. ```
  21.  
  22. ```{r}
  23. #Manipulating data set into new data frame leagueWork
  24.  
  25. leagueWork <- leagueRaw %>%
  26.  
  27.   #Inlcude only winning teams
  28.   filter(result == 1) %>%
  29.  
  30.   #Select only wanted variables
  31.   select(gameid, position, side, herald, heraldtime, goldat15, csat15, xpat10) %>%
  32.  
  33.   #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)
  34.   #Also replaced blank herald times with the maximum observation time (20 min)
  35.   mutate(herald = case_when(herald == 1 ~ 1, TRUE ~ 0), heraldtime = ifelse(is.na(heraldtime), 20,  heraldtime)) %>%
  36.  
  37.   #Select only the team observations, getting rid of player observations
  38.   filter(position == "Team") %>%
  39.  
  40.   #Remove unneccesary variables
  41.   select(-position)
  42.  
  43. #Check to make sure data table looks fine
  44. str(leagueWork)
  45. head(leagueWork, 50)
  46. ```
  47.  
  48. ```{r}
  49. #Finding Descriptive Statistics of Explanatory Variables
  50.  
  51.   #Max, Min Median, and IQR for Gold, XP, and CS
  52.  
  53.   #Gold
  54.   print(max(leagueWork$goldat15))
  55.   print(min(leagueWork$goldat15))
  56.   print(median_iqr(leagueWork$goldat15))
  57.  
  58.   #XP
  59.   print(max(leagueWork$xpat10))
  60.   print(min(leagueWork$xpat10))
  61.   print(median_iqr(leagueWork$xpat10))
  62.  
  63.   #CS
  64.   print(max(leagueWork$csat15))
  65.   print(min(leagueWork$csat15))
  66.   print(median_iqr(leagueWork$csat15))
  67.  
  68. ```      
  69.  
  70. ```{r}
  71. #Making Kaplan-Meier Curves for Explanatory Variables
  72.  
  73. #Kapla-Meier Curves for Map Side
  74. km.obj.mapside <- survfit(Surv(heraldtime, herald)~side, data = leagueWork)
  75. 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")
  76. legend(1, .8, c("Blue","Red"),lty=1:2)
  77.  
  78. #Kaplan-Meier Curves for Gold at 15 min
  79. km.obj.goldat15 <- survfit(Surv(heraldtime, herald)~1, data = leagueWork)
  80. cr.object.goldat15 <- coxph(Surv(heraldtime, herald)~goldat15, data = leagueWork)
  81.  
  82. pred.max.goldat15 <- data.frame(goldat15=max(leagueWork$goldat15), data = leagueWork)
  83. adj.surv.max <- survfit(cr.object.goldat15, newdata = pred.max.goldat15)
  84.  
  85. print(pred.max.goldat15)
  86.  
  87. pred.min.goldat15 <- data.frame(goldat15=min(leagueWork$goldat15), data = leagueWork)
  88. adj.surv.min <- survfit(cr.object.goldat15, newdata = pred.min.goldat15)
  89.  
  90. print(pred.min.goldat15)
  91.  
  92. pred.median.goldat15 <- data.frame(goldat15=median(leagueWork$goldat15), data = leagueWork)
  93. adj.surv.median <- survfit(cr.object.goldat15, newdata = pred.median.goldat15)
  94.  
  95. print(pred.median.goldat15)
  96.  
  97. 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))
  98. lines(adj.surv.min, lty = 2)
  99. lines(adj.surv.median, lty = 3)
  100. legend(1,.8,c("Max (28728G)","Min (22114G)", "Median (24521G)"),lty = 1:3)
  101.  
  102.  
  103. #Kaplan-Meier Curves for Experience at 10 min
  104.  
  105. km.obj.xpat10 <- survfit(Surv(heraldtime, herald)~1, data = leagueWork)
  106. cr.object.xpat10 <- coxph(Surv(heraldtime, herald)~xpat10, data = leagueWork)
  107.  
  108. pred.max.xpat10 <- data.frame(xpat10=max(leagueWork$xpat10), data = leagueWork)
  109. adj.surv.max <- survfit(cr.object.xpat10, newdata = pred.max.xpat10)
  110.  
  111. print(pred.max.xpat10)
  112.  
  113. pred.min.xpat10 <- data.frame(xpat10=min(leagueWork$xpat10), data = leagueWork)
  114. adj.surv.min <- survfit(cr.object.xpat10, newdata = pred.min.xpat10)
  115.  
  116. print(pred.min.xpat10)
  117.  
  118. pred.median.xpat10 <- data.frame(xpat10=median(leagueWork$xpat10), data = leagueWork)
  119. adj.surv.median <- survfit(cr.object.xpat10, newdata = pred.median.xpat10)
  120.  
  121. print(pred.median.xpat10)
  122.  
  123. 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))
  124. lines(adj.surv.min, lty = 2)
  125. lines(adj.surv.median, lty = 3)
  126. legend(1,.8,c("Max (20306xp)","Min (16780xp)", "Median (18924xp)"),lty = 1:3)
  127.  
  128.  
  129. #Kaplan-Meier Curves for CS at 10 min
  130.  
  131. km.obj.csat15 <- survfit(Surv(heraldtime, herald)~1, data = leagueWork)
  132. cr.object.csat15 <- coxph(Surv(heraldtime, herald)~csat15, data = leagueWork)
  133.  
  134. pred.max.csat15 <- data.frame(csat15=max(leagueWork$csat15), data = leagueWork)
  135. adj.surv.max <- survfit(cr.object.csat15, newdata = pred.max.csat15)
  136.  
  137. print(pred.max.csat15)
  138.  
  139. pred.min.csat15 <- data.frame(csat15=min(leagueWork$csat15), data = leagueWork)
  140. adj.surv.min <- survfit(cr.object.csat15, newdata = pred.min.csat15)
  141.  
  142. print(pred.min.csat15)
  143.  
  144. pred.median.csat15 <- data.frame(csat15=median(leagueWork$csat15), data = leagueWork)
  145. adj.surv.median <- survfit(cr.object.csat15, newdata = pred.median.csat15)
  146.  
  147. print(pred.median.csat15)
  148.  
  149. 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))
  150. lines(adj.surv.min, lty = 2)
  151. lines(adj.surv.median, lty = 3)
  152. legend(1,.8,c("Max (580 CS)","Min (446 CS)", "Median (521 CS)"),lty = 1:3)
  153. ```
  154.  
  155. ```{r}
  156.  
  157. #determining statistical significance in survival curves with log-rank tests
  158.  
  159. survdiff(Surv(heraldtime,herald)~side, data = leagueWork)
  160.  
  161.  
  162. ```
  163.  
  164. ```{r}
  165. #Fitting full and reduced models
  166. leagueCoxReduced <- coxph(Surv(heraldtime, herald) ~ side + goldat15 + csat15 + xpat10, data = leagueWork)
  167.  
  168. leagueCoxFull <- coxph(Surv(heraldtime, herald) ~ side + goldat15 + csat15 + xpat10 + csat15*xpat10 + csat15*goldat15, data = leagueWork)
  169.  
  170. #Model outputs
  171. summary(leagueCoxReduced)
  172. summary(leagueCoxFull)
  173.  
  174. #Doesn't violate PH assumptions
  175. test.ph <- cox.zph(leagueCoxReduced)
  176. test.ph
  177. ggcoxzph(test.ph)
  178.  
  179. #No obvious outliers in deviance
  180. ggcoxdiagnostics(leagueCoxReduced, type = "deviance", linear.predictions = FALSE, ggtheme = theme_bw())
  181.  
  182. #Martingale residuals might be showing a problem with csat15
  183. mart <- residuals(leagueCoxReduced, type = "martingale")
  184.  
  185. par(mfrow = c(1,3))
  186. plot(leagueWork$goldat15, mart, xlab = "Gold at 15 min", ylab = "Martingale Residuals", main = "Predictor: Gold at 15 min")
  187. smooth.sres <- lowess(leagueWork$goldat15, mart)
  188. lines(smooth.sres$x, smooth.sres$y, lty = 1)
  189.  
  190. plot(leagueWork$csat15, mart, xlab = "CS at 15 min", ylab = "Martingale Residuals", main = "Predictor: CS at 15 min")
  191. smooth.sres <- lowess(leagueWork$csat15, mart)
  192. lines(smooth.sres$x, smooth.sres$y, lty = 1)
  193.  
  194. plot(leagueWork$xpat10, mart, xlab = "XP at 10 min", ylab = "Martingale Residuals", main = "Predictor: XP at 10 min")
  195. smooth.sres <- lowess(leagueWork$xpat10, mart)
  196. lines(smooth.sres$x, smooth.sres$y, lty = 1)
  197. ```
  198.  
  199.  
  200. Inclusion of interaction terms do not improve model based on PLRT and AIC comparison
  201. ```{r}
  202. leagueCoxFull$loglik
  203. leagueCoxReduced$loglik
  204.  
  205. partiallikelihoodTestStat <- (2 * (-287.5465 - -287.7593))
  206. print(partiallikelihoodTestStat)
  207. 1-pchisq(partiallikelihoodTestStat, 2)
  208.  
  209. AIC(leagueCoxFull, leagueCoxReduced)
  210. ```
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!