# R code for chapter 63 Factorial ANOVA
# Book details:
# Health Science Statistics using R and R Commander
# by Robin Beaumont
# Paperback - ISBN 9781907904318
# http://www.amazon.co.uk/Health-Science-Statistics-using-Commander/dp/190790431X
# or direct from publisher at: 
# http://www.scionpublishing.com/
# book website: http://www.robin-beaumont.co.uk/rbook
###########################################################
#
headache <- 
  read.table("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/headache_stevens.txt",
   header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
### for me:
headache <- 
  read.table("D:/web_sites_mine/HIcourseweb new/book2data/headache_stevens.txt",
   header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
###########
   
headache$factor_headache_type <- as.factor(headache$headache_type)
headache$factor_treat_grp <- as.factor(headache$treat_grp)
headache
names(headache)

boxplot(headache$uncomf ~ headache$factor_treat_grp + headache$factor_headache_type, data=headache)  
 # or
 boxplot(uncomf ~ factor_headache_type *factor_treat_grp, data=headache)
### grouped boxplot
boxplot(uncomf ~interaction(factor_headache_type, factor_treat_grp), data=headache, col=c(2,15))
### but not gap between groups this can be done
boxplot(uncomf ~ headache$factor_treat_grp + headache$factor_headache_type, 
data=headache, at =c(1,2,3,4, 9,10,11,12))  

# interaction plots do not show individual points
interaction.plot(headache$factor_treat_grp, headache$factor_headache_type,
 headache$uncomf)
################### to develop
# also can't use jitter=0.2 with them
# for individual points use either boxplot or stripchart or boxplots (if you have sensible replication in each cell)
# below needs revising 
# stripchart(headache$uncomf ~ headache$factor_treat_grp + headache$factor_headache_type, vertical=TRUE, method="jitter", jitter=0.2)
# stripchart(headache$uncomf ~ headache$factor_treat_grp + headache$factor_headache_type, vertical=TRUE, method="jitter", pch="")
# stripchart(headache$uncomf[headache$factor_headache_type=="1"] ~ A[headache$factor_headache_type=="1"], method="jitter", pch=1, add=TRUE)
# stripchart(headache$uncomf[headache$factor_headache_type=="2"] ~ A[headache$factor_headache_type=="2"], method="jitter", pch=16, add=TRUE)
##################
#R commander generated interaction plot with ci's note the capital M in plotMeans
plotMeans(headache$uncomf, headache$factor_treat_grp, 
+   headache$factor_headache_type, error.bars="conf.int", level=0.95)

Alternative in the gplots package
# Plot Means with Error Bars
library(gplots)

plotmeans(headache$uncomf~ interaction(headache$factor_treat_grp, headache$factor_headache_type),
xlab="group", ylab="sensitivity score", main="Mean Plot\nwith 95% CI")
 
######################
######## model - R Commander generated
headache1 <- (lm(uncomf ~ factor_headache_type*factor_treat_grp, data=headache))
Anova(headache1)  #Anova Table (Type II tests)
tapply(headache$uncomf, 
   list(factor_headache_type=headache$factor_headache_type, 
   factor_treat_grp=headache$factor_treat_grp), mean, na.rm=TRUE) 
   # means
tapply(headache$uncomf, 
   list(factor_headache_type=headache$factor_headache_type, 
   factor_treat_grp=headache$factor_treat_grp), sd, na.rm=TRUE) 
   # std. deviations
tapply(headache$uncomf, 
   list(factor_headache_type=headache$factor_headache_type, 
   factor_treat_grp=headache$factor_treat_grp), function(x) sum(!is.na(x))) 
   # counts
#############
############################## spare development stuff
install.packages("phia")
library("phia")
interactionMeans(headache1)

###################### post hoc

factor_anova<- aov(uncomf ~ factor_headache_type*factor_treat_grp, data=headache)
TukeyHSD(factor_anova)
plot(TukeyHSD(factor_anova))  
# you can also specify a specific factor
TukeyHSD(factor_anova, "factor_headache_type")

 ######## if you want to use the previous lm model:
 library(multcomp)
summary(glht(headache1, linfct=mcp(factor_headache_type="Tukey", factor_treat_grp="Tukey")))
# or for all
summary(glht(headache1))
# to see raw code type stats::TukeyHSD.aov()
######################### pairwise comparisons
pairwise.t.test(uncomf, factor_treat_grp, p.adj="none", data="headache")
pairwise.t.test(uncomf, factor_treat_grp, p.adj="bonf", data="headache")
pairwise.t.test(uncomf, factor_treat_grp, p.adj="holm", data="headache")
#  to use individal sd's
pairwise.t.test(uncomf, factor_treat_grp, p.adj="none", data="headache" , pool.sd = FALSE)
pairwise.t.test(uncomf, factor_treat_grp, p.adj="bonf", data="headache" , pool.sd = FALSE)
pairwise.t.test(uncomf, factor_treat_grp, p.adj="holm", data="headache" , pool.sd = FALSE)

###########
factor_anovab<- aov(uncomf ~ factor_headache_type + factor_treat_grp, data=headache)
#### Dunnett's test
library(multcomp)           
test_dunnett=glht(factor_anovab,linfct=mcp(factor_treat_grp="Dunnett"))
confint(test_dunnett)
plot(test_dunnett)
###
# to change base level = 4:
factorlevels <- table(headache$factor_treat_grp) 
library(multcomp) 
contrasts <- contrMat(factorlevels, type = "Dunnett", base=4) 
dunnett_base4a <- glht(factor_anovab, linfct=mcp(factor_treat_grp=contrasts)) 
dunnett_base4a
confint(dunnett_base4a)
summary(dunnett_base4a)
plot(dunnett_base4a)
##  ignore below just providing example of a one sided ci
# dunnett_base4b <- glht(factor_anova, linfct=mcp(factor_treat_grp=contrasts),alternative='less') 
# dunnett_base4b 
# confint(dunnett_base4b)
# plot(dunnett_base4b)
##  
###### considering model without interaction
# to change base level = 4:
factorlevels <- table(headache$factor_treat_grp) 
contrasts <- contrMat(factorlevels, base=4) 
dunnett_base4_noint <- glht(factor_anovab, linfct=mcp(factor_treat_grp=contrasts)) 
dunnett_base4_noint 
confint(dunnett_base4_noint)
plot(dunnett_base4_noint)
summary(dunnett_base4_noint)
############## one sided
dunnett_base4a_noint <- glht(factor_anovab, linfct=mcp(factor_treat_grp=contrasts), alternative='less') 
dunnett_base4a_noint
confint(dunnett_base4a_noint)
plot(dunnett_base4a_noint)

##################### effect size measures
myanova <- anova.lm(factor_anovab) # calculate type I ANOVA table
myanova_ss <- myanova$"Sum Sq" # sums of squares
myanova_pes <- myanova_ss/(myanova_ss + myanova_ss[length(myanova_ss)]) # partial eta squared
myanova_pes[length(myanova_pes)] <- "" # clear effect size for residual
myanova$"Partial Eta Sq" <- myanova_pes # add effect size data to ANOVA table
myanova # display Type I ANOVA table with effect size column
r_stats <- summary.lm(factor_anovab) # summary of linear model
r_stats$"r.squared" # display R squared
r_stats$"adj.r.squared" # display Adjusted R squared
##############
##### original code from website
# datafilename="path/to/file.csv" # path to CSV data file
# mydata=read.csv(datafilename,header=TRUE,colClasses=c("numeric","factor","factor","factor")) # read CSV data file, specify data type
# mydata.aov = aov(DV~IV1*IV2*IV3,data=mydata)  # do the analysis of variance
# mydata.anova <- anova.lm(mydata.aov) # calculate type I ANOVA table
# mydata.ss <- mydata.anova$"Sum Sq" # sums of squares
# mydata.pes <- mydata.ss/(mydata.ss+mydata.ss[length(mydata.ss)]) # partial eta squared
# mydata.pes[length(mydata.pes)] <- "" # clear effect size for residual
# mydata.anova$"Partial Eta Sq" <- mydata.pes # add effect size data to ANOVA table
# mydata.anova # display Type I ANOVA table with effect size column
# r <- summary.lm(mydata.aov) # summary of linear model
# r$"r.squared" # display R squared
# r$"adj.r.squared" # display Adjusted R squared
##############################
### compare with the eta value produced by the etaSquared() function in the lsr  package
install.packages("lsr")
library("lsr") 
etaSquared(factor_anovab)  
# you can also specify what type of sums of squares you want it to use 
# type=1 
# type=2 # default
# type=3
# anova = FALSE default if set to true adds sums of squares etc
etaSquared(factor_anovab, type=1, anova = TRUE) 
etaSquared(factor_anovab, type=2, anova = TRUE) 
etaSquared(factor_anovab, type=3, anova = TRUE)
# you can see the code that forms the etaSquared function
# by typing 
# lsr::etaSquared 
###########################################

##  Omega squared taken from
# http://www.estudiosfonicos.cchs.csic.es/metodolo/1/
# http://www.estudiosfonicos.cchs.csic.es/metodolo/1/.Rprofile
##### OMEGA-SQUARED: ANOVA (a corrected effect-size index)
omega.squared <- function (model) {
	if (class(model)[1]=="aov") {
		df.factors=0
		for (i in 2:dim(model$model)[2]) {
			df.factors=df.factors+length(levels(model$model[,i]))-1
		}
	}
	if (class(model)[1]=="lm") {
		df.factors=summary(model)$fstatistic[2]
	}
	
	nvariables = dim(model$model)[2]
	nobservations = dim(model$model)[1]
	ncontrasts = 0
	for (i in 2:nvariables) {
		ncontrasts = ncontrasts + length(as.factor(levels(as.factor(model$model[,i]))))-1
	}
	endeffects = ncontrasts+1
	startresiduals = ncontrasts+2
	sseffects = sum(model$effects[2:endeffects]^2)
	ssresiduals = sum(model$effects[startresiduals:nobservations]^2)
	
	omega.sq = (sseffects-(df.factors*(ssresiduals/model$df.residual)))/(sseffects+ssresiduals+(ssresiduals/model$df.residual))
	attr(omega.sq,"names") <- NULL
	return(omega.sq)
}
#################
factor_anova<- aov(uncomf ~ factor_headache_type*factor_treat_grp, data=headache)
factor_anovab<- aov(uncomf ~ factor_headache_type + factor_treat_grp, data=headache)
omega.squared(factor_anova) # with interaction
omega.squared(factor_anovab) # additive terms
#  above only produces one value?
######  andy field omega squared for two way anova with interaction
omega_factorial<-function(n,a,b, SSa, SSb, SSab, SSr)
{
	MSa<-SSa/(a-1)
	MSb<-SSb/(b-1)
	MSab<-SSab/((a-1)*(b-1))
	MSr<-SSr/(a*b*(n-1))
	varA<-((a-1)*(MSa-MSr))/(n*a*b)
	varB<-((b-1)*(MSb-MSr))/(n*a*b)
	varAB<-((a-1)*(b-1)*(MSab-MSr))/(n*a*b)
	varTotal<-varA + varB + varAB + MSr
	
	print(paste("Omega-Squared A: ", varA/varTotal))
	print(paste("Omega-Squared B: ", varB/varTotal))
	print(paste("Omega-Squared AB: ", varAB/varTotal))
	}	
################ parameter details
# n    = number in cell; assumms balanced design
# a    = number of levels in factor a
# b    = number of levels in factor b
# SSa  = ssq for factor a
# SSb  = ssq for factor b
# SSab = ssq for interaction 
# SSr  = ssq for residual
# example call#	
#	omega_factorial(8,2,3,169,3332,1978,3488)
#######################  source option
# this onlt shows a fell lines:
# source("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/r_code/omega_squared_andy_field.r", echo=TRUE)
# to see all the data
thecode <- source("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/r_code/omega_squared_andy_field.r")

omega_factorial(11, 2, 4, 14.52, 73.38, 0 ,879.69)	
##################
# Olejnik & Algina, 2003, suggest alternatives
# Generalised Eta value aov.car() function in the afex package. 
install.packages("afex", repos="http://R-Forge.R-project.org")
library("afex")
# need to add a subject id to dataframe
headache$rowid <- rownum(headache)
# ####
aov.car(uncomf ~ factor_headache_type + factor_treat_grp + Error(rowid), return="Anova", data=headache)
aov.car(uncomf ~ factor_headache_type + factor_treat_grp + Error(rowid), return="nice", data=headache)
# next option produces ln model with all output!!
aov.car(uncomf ~ factor_headache_type + factor_treat_grp + Error(rowid), return="lm", data=headache)	
### this approach also can divid the factors up into observed and manipulated factor/variables
# > names(headache)
#[1] "headache_type"        "treat_grp"           
#  "pre_uncomf"           "pre_fun"              "uncomf"               "dfunpl"               
# "factor_headache_type" "factor_treat_grp"  [9] "rowid" 
can use pre_uncomf as a covariate
aov.car(uncomf ~ factor_headache_type + factor_treat_grp +pre_uncomf 
 +  Error(rowid), observed = c("pre_uncomf"), return="Anova", data=headache)
aov.car(uncomf ~ factor_headache_type + factor_treat_grp +pre_uncomf 
 +  Error(rowid), observed = c("pre_uncomf"), return="nice", data=headache)
 # or possibly
 aov.car(uncomf ~ factor_headache_type + factor_treat_grp  +  
 Error(rowid), observed = c("factor_headache_type"),  return="Anova",   data=headache)
 aov.car(uncomf ~ factor_headache_type + factor_treat_grp  +  
 Error(rowid), observed = c("factor_headache_type"),  return="nice",   data=headache)
 ## above provides identical result to treating factor_headache-type as a 
 ## manipulated factor.

##################################################
install.packages("compute.es")
library("compute.es")
# the function takes the following parameters 
# mes(mean1, mean2, sd1, sd2, n1, n2)
#######cell details
# Uncomf =outcome variable, , Treat_grp factor (4 levels), marginal , Level:, 1, 2, 3, 4, 
# Treat_grp level 1 ,2, 3, 4; headache_type factor level 1  = migraine (=1)  
# x¯11=2.22, sd=1.148,n=11
# x¯12=3.1,  sd=3.18, n=11
# x¯13=4.78, sd=3.85, n=12  # treatment group
# x¯14=2.19, sd=3.17, n=11  # control group
# Treat_grp level 1 ,2, 3, 4; headache_type factor level 2 =tension
# x¯21=3.36, sd=3.50, n=14
# x¯22=4.53, sd=2.91, n=11
# x¯23=4.69, sd=3.81, n=16
# x¯24=2.87, sd=1.69, n=12
############
# compare control group to treatment group for migraine suffers (=1):
mes(4.78, 2.19, 3.85, 3.17, 12, 11)

# compare control group to treatment group for tension headache suffers (=2):
mes(4.69, 2.87, 3.81, 1.69, 16, 12)

############################# adding observations to the boxplots
 # basic method
stripchart(uncomf ~ factor_headache_type + factor_treat_grp, pch =16,cex = 1, vertical=TRUE, data=headache)
 
boxplot(uncomf ~ headache_type+ treat_grp, data = headache, 
outline = FALSE,     ## avoid double-plotting outliers, if any
main = 'boxplot + stripchart', add = TRUE)
 
# Boxplots with jitter method 1 using beeswarm package
install.packages('beeswarm')
library(beeswarm)

boxplot(uncomf ~ headache_type+ treat_grp, data = headache, 
outline = FALSE,     ## avoid double-plotting outliers, if any
main = 'boxplot + beeswarm')
beeswarm(uncomf ~ headache_type+ treat_grp, method = "square", data = headache, 
col = 4, pch = 16, add = TRUE)
###########################	
##### using GGplot2 functionsfor points + boxplots 
##### ignore neew work
# First creat temorary dataframe  with just grouping var and outcome
#### get variable names
names(headache)
# [1] "headache_type"        "treat_grp"           
# "pre_uncomf"           "pre_fun"             
# "uncomf"               "dfunpl"               
# "factor_headache_type" "factor_treat_grp"[9] "rowid" 
 
headache_reduced1<- headache[,c("headache_type", "treat_grp", "uncomf")]
headache_reduced1
boxplot(uncomf ~ headache_type + treat_grp, data=headache_reduced1)
library(reshape2) 
longdf <- melt(mydataframe)
library(ggplot2)
ggplot(data=longdf) + geom_boxplot(aes(x=variable,y=value))
or
qplot(variable, value,  geom="boxplot", data=longdf)
Showing individual points as well.
qplot(variable, value, geom=c("boxplot", "jitter"), data=longdf)
########################

##############  ignore below
data(Headache)
str(Headache)
# basic MLM, specifying between-S effects
# headache.mod <- lm(cbind(u1, du1, u2, du2) ~ type * treatment, data=Headache)
##############################
## between-S tests
##############################
Anova(headache.mod, test="Roy")
# test each contrast separately

print(linearHypothesis(headache.mod, hypothesis="treatment1", test="Roy"), SSP=FALSE)
print(linearHypothesis(headache.mod, hypothesis="treatment2", test="Roy"), SSP=FALSE)
print(linearHypothesis(headache.mod, hypothesis="treatment3", test="Roy"), SSP=FALSE)
heplot(headache.mod, variables=c(1,3),
hypotheses=paste("treatment", 1:3, sep=""),
hyp.labels=c("extra.exp", "no.inst", "explicit.inst"),
xlab="u1: Initial sensitivity", ylab="u2: Final sensitivity",
main="Headache data: Unpleasant levels")
abline(0, 1, lty=5, col="gray")

heplot(headache.mod, variables=c(2,4),
hypotheses=paste("treatment", 1:3, sep=""),
xlab="du1: Initial sensitivity", ylab="du2: Final sensitivity",
main="Headache data: Definitely Unpleasant levels")
abline(0, 1, lty=5, col="gray")
pairs(headache.mod)
##############################
# between-S and within-S tests
##############################
idata = expand.grid(level=factor(c("U", "DU")), phase=factor(1:2))
Anova(headache.mod, idata=idata, idesign=~level*phase)


############## testing structure of lm compared to aov()
# factor_anovab$"Sum Sq" # does not work with a aov() object
lmfactor_anovab <- anova.lm(factor_anovab)
lmfactor_anovab$"Sum Sq"
length(as.factor(levels(as.factor(factor_anova$model[,1]))))
length(as.factor(levels(as.factor(factor_anova$model[,2]))))
length(as.factor(levels(as.factor(factor_anovab$model[,1]))))  
length(as.factor(levels(as.factor(factor_anovab$model[,2])))) 

################################################
read.table