# 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