# R Examples of Analysis of Random Effects Model. # We analyze the Apex Enterprises data from Chapter 25. # Save the data file into a directory and # use the full path name: apex.data <- read.table(file = "z:/My Documents/teaching/stat_705/apexdata.txt", header=FALSE, col.names = c('rating', 'officer', 'candidate')) # attaching the data frame: attach(apex.data) # Making "officer" a factor: officer <- factor(officer) # Fitting the no-interaction ANOVA model and producing the ANOVA table: apex.fit <- lm(rating ~ officer); anova(apex.fit) # There is no equivalent of the RANDOM statement in the lm() function, # but keep in mind that "officer" has random levels. # The F* test statistic is calculated in the same way as in the fixed effects case. # We see there is significant variation in rating among the population # of officers. (F* = 5.39, P-value = 0.0068.) ###############################################################################* # 90% CI for mu-dot and 90% CI for sigma^2 # The 90% CI for mu-dot: # Run this function in R: CI.for.mu.dot <- function(response, conf.level, r, MSTR, n){ alpha <- 1 - conf.level LL <- mean(response) - qt(1 - alpha/2, df=r-1)*sqrt(MSTR/(r*n)) UL <- mean(response) + qt(1 - alpha/2, df=r-1)*sqrt(MSTR/(r*n)) CI <- list(paste(conf.level*100, "percent CI for mu.dot:", round(LL,4), round(UL,4))) return(CI) } # And then call the function with the correct values of # response, conf.level, r, MSTR, n: # From the anova() output, MSTR = 394.93. # Also: conf.level = 0.90, r = 5, n = 4. CI.for.mu.dot(response=rating, conf.level = 0.90, r = 5, n = 4, MSTR=394.93) # The 90% CI for mu-dot is (61.98, 80.92). (Find this on the R output.) # The 90% CI for sigma^2: # Run this function in R: CI.for.sigma.sq <- function(conf.level, r, MSE, n){ alpha <- 1 - conf.level LL <- r*(n-1)*MSE / qchisq(1 - alpha/2, df=r*(n-1)) UL <- r*(n-1)*MSE / qchisq(alpha/2, df=r*(n-1)) CI <- list(paste(conf.level*100, "percent CI for sigma^2:", round(LL,4), round(UL,4))) return(CI) } # And then call the function with the correct values of # conf.level, r, MSE, n: # From the anova() output, MSE = 73.28. # Also: conf.level = 0.90, r = 5, n = 4. CI.for.sigma.sq(conf.level = 0.90, r = 5, n = 4, MSE = 73.28) # The 90% CI for sigma^2 is (43.98, 151.39). (Find this on the R output.) ###############################################################################* # CI for intraclass correlation coefficient and CI for sigma_mu^2 ## For these data: # n = 4; # input number of observations per cell # r = 5; # input number of levels of the random factor # MSTR = 394.93; # input MSTR from output above # MSE = 73.28; # input MSE from output above # Copy this function into R: otherCIs <- function(conf.level, MSTR, MSE, r, n){ # lots of calculations alpha <- 1-conf.level fl<-qf(1-alpha/2,r-1,r*(n-1)); fu<-qf(alpha/2,r-1,r*(n-1)); l<-1/n*((MSTR/MSE)*(1/fl)-1); u<-1/n*((MSTR/MSE)*(1/fu)-1); LCL_ICC<-l/(1+l); UCL_ICC<-u/(1+u); c1<-1/n; c2<- -1/n; lhat<-(MSTR-MSE)/n; f1<-qf(1-alpha/2,r-1,9999); f2<-qf(1-alpha/2,r*(n-1),9999); f3<-qf(1-alpha/2,9999,r-1); f4<-qf(1-alpha/2,9999,r*(n-1)); f5<-qf(1-alpha/2,r-1,r*(n-1)); f6<-qf(1-alpha/2,r*(n-1),r-1); g1<-1-(1/f1); g2<- 1-(1/f2); g3<-((f5-1)^2-(g1*f5)^2-(f4-1)^2)/f5; g4<-f6*(((f6-1)/f6)^2-(((f3-1)/f6)^2)-g2^2); hl<-(((g1*c1*MSTR)^2)+(((f4-1)*c2*MSE)^2)-g3*c1*c2*MSTR*MSE)^.5; hu<-(((f3-1)*c1*MSTR)^2+((g2*c2*MSE)^2)-(g4*c1*c2*MSTR*MSE))^.5; LCL_sigma_mu_sq<-lhat-hl; UCL_sigma_mu_sq<-lhat+hu; CIs <- list(paste(100*conf.level, "% CI for ICC:", round(LCL_ICC,3), round(UCL_ICC,3), "and", 100*conf.level, "% CI for sigma_mu^2:", round(LCL_sigma_mu_sq,3), round(UCL_sigma_mu_sq,3)) ) return(CIs) } # And then call the function with the correct values of # conf.level, MSTR, MSE, r, n: otherCIs(conf.level=0.90, MSTR=394.93, MSE=73.28, r=5, n=4) # A 90% CI for the intraclass correlation coefficient is (0.160, 0.884). # Compare to book's result, at top of page 1041. # A 90% CI for sigma_mu^2 is (20.22, 536.32). # Compare to book's result, at bottom of page 1046. ###############################################################################*