# resin.R library(car) library(lsmeans) library(MASS) library(nlme) resdata <- read.table(file="resin.dat",header=T) # read the data in head(resdata) resdata\$Temp <- factor(resdata\$temp) #side-by-side boxplots. plot(failtime~Temp,data=resdata) title(main="Boxplots of fail time for each temp") m1 <- aov(failtime~Temp,resdata) summary(m1) # get ls means lsmeans(m1, specs = ~ Temp) # Conduct modified Levene's test: leveneTest(m1) # agrees with SAS # get resids vs fitteds plot and normal Q-Q plot plot(m1, which=c(1,2)) # Determine a power transformation of y that will be more # normal and homoscedastic # First use regression approach: # compute log(mean) and log(SD) of failtime within each Temp sumStats <- matrix(0,nrow=5,ncol=2) j <- 0 for(i in unique(resdata\$temp)){ j <- j+1 sumStats[j,] <- log( c(mean(resdata\$failtime[resdata\$temp==i]), sd(resdata\$failtime[resdata\$temp==i])) ) } sumData <- data.frame(logmn = sumStats[,1], logsd = sumStats[,2]) head(sumData) #regress log(SD) on log(mean) to estimate power transformation lm(logsd ~ logmn,data=sumData) # suggests taking y^lambda where lambda=(1-0.863) which is close to 0 # so use log transformation instead # Alternatively, can use the Box-Cox approach, which is implemented in # the boxcox() function in the MASS package boxcox(failtime~Temp, data=resdata) # also suggests a log transformation resdata\$logfail <- log(resdata\$failtime) m2 <- aov(logfail~Temp,resdata) summary(m2) # get ls means lsmeans(m2, specs = ~ Temp) # Conduct modified Levene's test: leveneTest(m2) # agrees with SAS # get resids vs fitteds plot and normal Q-Q plot plot(m2, which=c(1,2)) # Alternatively, could analyze failtime (no transformation) allowing for # heterogeneous variances. This can be done in the gls() function in # the nlme package m3 <- gls(failtime~Temp,data=resdata,weights=varIdent(form= ~ 1 | Temp)) summary(m3) # note that the estimated variances in each treatment are as follows: # Temp=175: (resid se)^2= 13.029^2 = 169.75 # Temp=194: [(resid se)*(0.7340)]^2 = 91.449 # Temp=213: [(resid se)*(0.4919)]^2 = 41.074 # Temp=231: [(resid se)*(0.1330)]^2 = 3.001 # Temp=250: [(resid se)*(0.2840)]^2 = 13.69 # These agree with SAS plot(m3) # gives standardized resids vs fitteds plot qqnorm(m3) # gives normal Q-Q plot of standardized resids