## pdf("E:\\blue_drive\\sta4211\\hops.pdf") hops <- read.table("http://www.stat.ufl.edu/~winner/data/hop_alphaacid.dat",header=F, col.names=c("variety","repnum","a_acid")) attach(hops) (ybar <- mean(a_acid)) (ybar.trt <- as.vector(tapply(a_acid,variety,mean))) (var.trt <- as.vector(tapply(a_acid,variety,var))) (rep.trt <- as.vector(tapply(a_acid,variety,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))) #### 95% 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(.975,df),df*s2.mu/qchisq(.025,df))) stripchart(a_acid ~ variety, xlab="alpha_acid", ylab="Variety") #### linear mixed effects model library(lmerTest) variety <- factor(variety) hops.mod1 <- lmer(a_acid ~ 1 + (1|variety)) summary(hops.mod1) rand(hops.mod1) confint(hops.mod1) ## dev.off()