nba2016 <- read.csv("http://www.stat.ufl.edu/~winner/data/nba_player_201617a.csv") attach(nba2016); names(nba2016) pts36 <- 36*PTS/Minutes # pts36 <- log(36*PTS/Minutes+1) (n.play <- max(Player_ID1)) (n.play.game <- length(PTS)) game.play <- as.vector(tapply(PTS,Player_ID1,length)) mean.play <- as.vector(tapply(pts36,Player_ID1,mean)) (mu <- mean(mean.play)) alpha.play <- mean.play-mu hist(alpha.play,breaks=30) var.e <- rep(0,n.play) eps.play.game <- rep(0, n.play.game) count <- 0 for (i1 in 1:n.play) { count1 <- count + 1 count2 <- count + game.play[i1] for(i2 in 1:game.play[i1]) { eps.play.game[count1:count2] <- pts36[Player_ID1==i1] - mean.play[i1] } count <- count2 var.e[i1] <- var(pts36[Player_ID1==i1] - mean.play[i1]) } sum(eps.play.game[1:game.play[1]]) hist(eps.play.game,breaks=100) plot(var.e ~ alpha.play) abline(lm(var.e ~ alpha.play)) summary(lm(var.e ~ alpha.play)) (sigma2a <- var(alpha.play)) (sigma2e <- mean(var.e)) (rho_I <- sigma2a/(sigma2a+sigma2e)) n.sample <- 10000 n.play.s <- 10 n.game.play.s <- 8 set.seed(7531) mean.out <- matrix(rep(0,n.sample*n.play.s),ncol=n.play.s) var.out <- matrix(rep(0,n.sample*n.play.s),ncol=n.play.s) MS_TRT.out <- rep(0,n.sample) for (i in 1:n.sample) { play.s <- sample(1:n.play,n.play.s) for (i1 in 1:n.play.s) { play.s1 <- pts36[Player_ID1 == play.s[i1]] play.g.s <- sample(1:game.play[play.s[i1]],n.game.play.s) mean.out[i,i1] <- mean(play.s1[play.g.s]) var.out[i,i1] <- var(play.s1[play.g.s]) } MS_TRT.out[i] <- n.game.play.s* sum((mean.out[i,] - mean(mean.out[i,]))^2)/(n.play.s-1) } (sigma2_muhat <- sigma2a/n.play.s + sigma2e/(n.play.s*n.game.play.s)) muhat <- rowMeans(mean.out) summary(muhat) var(muhat) hist(muhat,breaks=50) s2e <- rowSums(var.out)/n.play.s summary(s2e) hist(s2e,breaks=50) s2a <- (MS_TRT.out-s2e)/n.game.play.s summary(s2a) hist(s2a,breaks=50) rho_I_hat <- s2a / (s2a + s2e) summary(rho_I_hat) hist(rho_I_hat, breaks=50) F_L <- qf(.025,n.play.s-1,n.play.s*(n.game.play.s-1)) F_U <- qf(.975,n.play.s-1,n.play.s*(n.game.play.s-1)) F_obs <- MS_TRT.out / s2e rho_L <- 1/(n.game.play.s*(1/(F_obs/F_U-1))+1) rho_U <- 1/(n.game.play.s*(1/(F_obs/F_L-1))+1) sum((rho_L <= rho_I) & (rho_U >= rho_I)) / n.sample