# Mann whitney U example
# mydataframe <- read.delim("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/mwu_example1.dat", header=TRUE)

# local:
mydataframe <- read.delim("D:\\web_sites_mine\\HIcourseweb new\\book2data\\mwu_example1.dat", header=TRUE)
names(mydataframe)
#######
# get a overall histogram of the result variable
hist(mydataframe$result, col="lightgray")
############
#now create 2 new columns to be able to produce separate histograms
gpconsults <- mydataframe[mydataframe$grouping == 1,"result"]
gpconsults
### same thing for the nurse results
nurseconsults <- mydataframe[mydataframe$grouping ==2,"result"]
nurseconsults
hist(gpconsults, col="darkgray")
hist(nurseconsults, col="red")
############## alternative way uses the density plot
plot(density(nurseconsults))
lines(density(gpconsults), lty=2)
#####################
# boxplot
boxplot(result ~grouping, data = mydataframe)
#  now do the Mannwhitney U test
wilcox.test(result ~ grouping, alternative="two.sided", data=mydataframe, conf.int=TRUE)
### effect size
effectsize <- median(nurseconsults) - median(gpconsults)
effectsize

  ##################### 
  
# alternative method:

plot(density(mydataframe$age[mydataframe$gender == "Female"]), ylim = c(0, 0.07), 
main = "densityplots for males/females[dotted] for age", xlab= "age (years)" )
# You now need to overlay the female density line.
lines(density(mydataframe$age[mydataframe$gender == "Male"]), lty = 2)

########################

attach(mydataframe)
length(result)
# check number in each group
table(grouping)
# obtain ranks
ranked_result <- rank(result)
# obtain frequency table of ranks
table(ranked_result)
# obtain histogram of ranks
hist(ranked_result, breaks=length(ranked_result))

#####################
# mean ranks for each group
# copy the data to a newdataframe
newdataframe <- mydataframe
# add the ranked score column to the new dataframe
newdataframe["ranking"] <- rank(mydataframe$result)
# attach it
attach(newdataframe)
# print it out to check
newdataframe
# now get the means for each column divided by group
lapply((split(newdataframe$ranking, newdataframe$grouping)), mean)

########################
# exact permutation test
# the exactRankTests package is no 
# longer being developed
# when you use it it works but recommends using
# coin
############
install.packages("exactRankTests", dependencies=TRUE)
library(exactRankTests)
wilcox.exact(result ~ grouping, alternative="two.sided", data=mydataframe,conf.int=TRUE) 

################# using coin instead
install.packages("coin", dependencies=TRUE)
library(coin)

wilcox_test(result ~ grouping, alternative="two.sided",distribution = "exact", conf.int = TRUE, 
 data=mydataframe)
######################### small sample example 
################  this is shown at start of chapter
control <- c(2,3,1,4,3)
treatment <- c(5,6,4,7,6)
jogging_reading <- data.frame(control, treatment)
## boxplot = shown at start of chapter
boxplot(control, treatment, data=jogging_reading, 
names = c("control", "treatment"),
main="High ratings: jogging v. reading",
id.method="y", col = "gray90", medlwd=6)
# medlwd option increases thickness of medean line
# gray90 = 90% grey; you can use gray or grey 

### stripchart
stripchart(list(jogging_reading$control,jogging_reading$treatment), 
vertical = TRUE,
group.names= c("control","treatment"), 
main = "stripchart of jogging versus reading results")

wilcox.test(jogging_reading$control, jogging_reading$treatment, conf.int=TRUE)



newdataframe <- mydataframe
attach(newdataframe)
newdataframe["ranking"] <- rank(mydataframe$result)
lapply((split(newdataframe$ranking, newdataframe$grouping)), mean)

# drawing and calculating various values
n1<-32  # group 1 size
n2<-33   # group 2 size
U_result <- 506.506
# using U1=n1n2 - U2
U_comp <- n1*n2 - U_result
U_comp
#expected value is (n1 * n2)/2
expect_u <- (n1*n2)/2
expect_u
# need lines each side of the mean value
# now for the plot
efsize<-(n1*n2) + 1
efsize
x <- seq(0,efsize)
plot(x, dwilcox(x,n1,n2), type="h", xlab="U1 value", ylab="p", col="pink", 
main=paste("U  pdf (n =",n1, ",", n2 ,"), ","expected U =",expect_u ))
abline(v=U_result)
abline(v=U_comp)

############### man whitney U pdf
#  n=4, 4
n1<-4  # group 1 size
n2<-4   # group 2 size
U_result <- 13
U_comp <- n1*n2 - U_result
U_comp
expect_u <- (n1*n2)/2
expect_u
efsize<-(n1*n2) + 1
efsize
x <- seq(0,efsize)
plot(x, dwilcox(x,n1,n2), type="h", xlab="U1 value", ylab="p", col="black", 
main=paste("U  distribution (n =",n1, ",", n2 ,"), ","expected U =",expect_u ))
abline(v=U_result) # one U value
abline(v=U_comp)    # the other U value, they are symmetrical

#  n=5, 5
n1<-5  # group 1 size
n2<-5   # group 2 size
U_result <- 506.506
U_comp <- n1*n2 - U_result
U_comp
expect_u <- (n1*n2)/2
expect_u
efsize<-(n1*n2) + 1
efsize
x <- seq(0,efsize)
plot(x, dwilcox(x,n1,n2), type="h", xlab="U1 value", ylab="p", col="black", 
main=paste("U  pdf (n =",n1, ",", n2 ,"), ","expected U =",expect_u ))
abline(v=U_result) # one U value
abline(v=U_comp)    # the other U value, they are symmetrical



#  n=30, 30
n1<-30  # group 1 size
n2<-30   # group 2 size
#this optionally allow you to show where an observed U value lies
U_result <- 506.506
U_comp <- n1*n2 - U_result
U_comp
expect_u <- (n1*n2)/2
expect_u
efsize<-(n1*n2) + 1
efsize
x <- seq(0,efsize)
plot(x, dwilcox(x,n1,n2), type="h", xlab="U1 value", ylab="p", col="black", 
main=paste("U distribution (n =",n1, ",", n2 ,"), ","expected U =",expect_u ))
abline(v=U_result) # one U value
abline(v=U_comp)    # the other U value, they are symmetrical