R program for one-way ANOVA with the tensile strength data

# There are two easy ways to run these commands in R.  One is to merely 'copy and paste' 
# these commands into R's 'Console window'.  They will be executed automatically, and 
# you can then go back and see which commands resulted in what output.

# A better way, in general, is to (within R) click on 'File -> New script', and then 
# copy and paste the commands into the script file that has opened up.  With that 
# window highlighted, click on 'Edit -> Run all'.  
# You can save this script file, and run it again later, perhaps after making changes to it.

# In doing assigned questions it is a good idea to put your commands into a script file, and
# just keep running/changing/running the script file, until you get things right.  That way, 
# once you get some of the commands right, you won't have to keep re-entering them.



A <- c(7, 7, 15, 11, 9)
B <- c(12, 17, 12, 18, 18)
C <- c(14, 18, 18, 19, 19)
D <- c(19, 25, 22, 19, 23)
E <- c(7, 10, 11, 15, 11)

data <- rbind(A, B, C, D, E)

data #This time data is a matrix, rather than a frame, 
       # with the treatments as rows

# We suppose that the measurements were made in the following, randomly chosen, order:
times <- matrix(sample(1:25), nrow=5)
times

#The sample mean for each treatment
apply(data, 1, mean)

#The sample mean for each replicate
apply(data, 2, mean)

# The overall mean
mean(data)

#Use 'sum' rather than 'mean' for totals


# Draw a scatterplot of 'data' as a data frame with treatments as columns
stripchart(data.frame(t(data)), vertical = TRUE)

#Arrange responses by column
strength <- c(data)
strength

# Total SS:
(25-1)*var(strength)

#Set the factor at 5 levels

d <- rep(c("A", "B", "C", "D", "E"), times=5)
content <- as.factor(d)
content

#Perform one-way ANOVA (lm stands for 'linear model').
g <- lm(strength~content)
anova(g)

names(g)  # Here are the components of g; 

windows()
par(mfrow=c(2,2))  # This arranges the plots in a 2 by 2 matrix, 
				#row by row.
qqnorm(g$resid)
qqline(g$resid)
plot(g$fitted,g$resid)
abline(h=0)  #Add the horizontal axis
time.order <- c(times)
plot(time.order, g$resid)
abline(h=0) 

# Bartlett's test for equality of variances between treatment groups
bartlett.test(g$resid, content)

# Carry out Levene's test:  First get the absolute values of the 
# deviations from the medians:

m <- apply(data, 1, median) # apply "median" to each row of data
diff <- abs(strength-rep(m, times=5))

m
diff

g.d <- lm(diff~content)
anova(g.d)       

# Do a nonparametric test.
# First just do ANOVA on the ranks:
R <- rank(strength)
strength
R
g.rank <- lm(R~content)
anova(g.rank)

# Now do Kruskal-Wallis:
kruskal.test(strength,content)