pdf("E:\\sta4211\\KNNL_TA25_01.pdf") apex <- read.table("http://www.stat.ufl.edu/~winner/sta4211/datasets/CH25TA01.txt",header=F, col.names=c("rating","officer","candidate")) attach(apex) (ybar <- mean(rating)) (ybar.trt <- as.vector(tapply(rating,officer,mean))) (var.trt <- as.vector(tapply(rating,officer,var))) (rep.trt <- as.vector(tapply(rating,officer,length))) (num.trt <- length(ybar.trt)) sse <- 0; sstr <- 0 for (i1 in 1:num.trt) { sse <- sse + (rep.trt[i1]-1)*var.trt[i1] sstr <- sstr + rep.trt[i1]*((ybar.trt[i1]-ybar)^2) } df.trt <- num.trt-1 df.err <- sum(rep.trt)-num.trt (mse <- sse/df.err) (mstr <- sstr/df.trt) (f_star <- mstr/mse) (p_f_star <- 1-pf(f_star,df.trt,df.err)) ### 95% CI for mu_dot (ybar + c(qt(.025,df.trt)*sqrt(mstr/sum(rep.trt)),qt(.975,df.trt)*sqrt(mstr/sum(rep.trt)))) #### 95% CI for rho_I L <- (1/rep.trt[1])*((mstr/mse)*(1/qf(.975,df.trt,df.err))-1) U <- (1/rep.trt[1])*((mstr/mse)*(1/qf(.025,df.trt,df.err))-1) (c(L/(1+L),U/(1+U))) #### 95% CI for sigma^2 (c(sse/qchisq(.975,df.err),sse/qchisq(.025,df.err))) #### 90% CI for sigma_mu^2 s2.mu <- (mstr-mse)/rep.trt[1] df.num <- (rep.trt[1]*s2.mu)^2 df.den <- (mstr^2/df.trt) + (mse^2/df.err) (df <- round(df.num/df.den)) (c(df*s2.mu/qchisq(.95,df),df*s2.mu/qchisq(.05,df))) stripchart(rating ~ officer, xlab="Rating", ylab="Officer") dev.off()