isaactpetersen

lme(aav ~ 1 + numTeams * I(numTeams^2) * auctionValue * positionRank * position, random = ~ 1 + numTeams|player, method="REML", data=aavData)

isaactpetersen

#Calculate SD of fantasy points per week
sdWeeklyPts <- matrix(nrow=NROW(fantasyPts[[1]]), ncol=simulations)

for(i in 1:simulations){  
  sdWeeklyPts[,i] <- apply(fantasyPts[[i]], 1, function(x) sd(x, na.rm=TRUE))
}

#Calculate robust average of weekly SD
projections$weeklySD <- apply(sdWeeklyPts, 1, function(x) tryCatch(wilcox.test(x, conf.int=TRUE, na.action="na.exclude")$estimate, error=function(e) median(x, na.rm=TRUE)))

isaactpetersen

#Calculate fantasy points per week
passYdsPts <- list()
passTdsPts <- list()
passIntPts <- list()
rushYdsPts <- list()
rushTdsPts <- list()
recPts <- list()
recYdsPts <- list()
recTdsPts <- list()
twoPtsPts <- list()

isaactpetersen

#Simulation
simulations <- 100
games <- 16

passYds <- list()
passTds <- list()
passInt <- list()
rushYds <- list()
rushTds <- list()
rec <- list()

isaactpetersen

simulateIntegers <- function(n, sum, sd, pos.only = TRUE){
  if(sum == 0 & pos.only == TRUE){
    vec <- rep(0, n)
  } else{
    vec <- rnorm(n, sum/n, sd)
    if (abs(sum(vec)) < 0.01) vec <- vec + 1
    vec <- round(vec / sum(vec) * sum)
    deviation <- sum - sum(vec)
    for (. in seq_len(abs(deviation))){
      vec[i] <- vec[i <- sample(n, 1)] + sign(deviation)

isaactpetersen

#Calculate week-to-week SD for each statistical category across same player/year combinations
sdPassYds <- apply(weeklyDataWide[,grep("passYds", names(weeklyDataWide))], 1, sd)
sdPassTds <- apply(weeklyDataWide[,grep("passTds", names(weeklyDataWide))], 1, sd)
sdPassInt <- apply(weeklyDataWide[,grep("passInt", names(weeklyDataWide))], 1, sd)
sdRushYds <- apply(weeklyDataWide[,grep("rushYds", names(weeklyDataWide))], 1, sd)
sdRushTds <- apply(weeklyDataWide[,grep("rushTds", names(weeklyDataWide))], 1, sd)
sdRec <- apply(weeklyDataWide[,grep("rec", names(weeklyDataWide))], 1, sd)
sdRecYds <- apply(weeklyDataWide[,grep("recYds", names(weeklyDataWide))], 1, sd)
sdRecTds <- apply(weeklyDataWide[,grep("recTds", names(weeklyDataWide))], 1, sd)

isaactpetersen

#Long to Wide
weeklyDataWide <- reshape(weeklyDataLong,
                          timevar = c("week"),
                          idvar = c("name","year"),
                          direction = "wide",
                          sep="")

isaactpetersen

#Libraries
library("XML")

#Specify info to scrape
years <- 2011:2013
weeks <- 17

#Scrape data
qb <- list()
rb1 <- list()

isaactpetersen

#Number of players at each position drafted in Top 100 (adjust for your league)
qbReplacements <- 17
rbReplacements <- 35
wrReplacements <- 35
teReplacements <- 13

qbValueOfReplacement <- mean(c(qb$projections[qb$positionRank==qbReplacements],qb$projections[qb$positionRank==(qbReplacements-1)],qb$projections[qb$positionRank==(qbReplacements+1)]))
rbValueOfReplacement <- mean(c(rb$projections[rb$positionRank==rbReplacements],rb$projections[rb$positionRank==(rbReplacements-1)],rb$projections[rb$positionRank==(rbReplacements+1)]))
wrValueOfReplacement <- mean(c(wr$projections[wr$positionRank==wrReplacements],wr$projections[wr$positionRank==(wrReplacements-1)],wr$projections[wr$positionRank==(wrReplacements+1)]))
teValueOfReplacement <- mean(c(te$projections[te$positionRank==teReplacements],te$projections[te$positionRank==(teReplacements-1)],te$projections[te$positionRank==(teReplacements+1)]))

isaactpetersen

#Calculate shape and rate parameters of average projections for Weibull distribution
weibullShape <- fitdistr(projections$projectedPts, 'weibull')$estimate[[1]]
weibullScale <- fitdistr(projections$projectedPts, 'weibull')$estimate[[2]]

projectedPtsLatentWeibull <- qweibull(pnorm(projectedPtsLatent), shape=weibullShape, scale=weibullScale)

#Recale distribution to have same range as average projections
rescaleRange function(variable, minOutput, maxOutput){
  minObserved min(variable)
  maxObserved max(variable)