### Practicals 10
library(ggplot2)
library(readxl)

## Task 1
hair<-read_excel("10_data.xlsx", sheet="1")
summary(hair)

ggplot(hair, aes(faculty, hair.length))+
  geom_boxplot()+theme_classic() # It is evident from the plot that 1. variance is 
# strongly heterogeneous and 2. residual distribution in strongly skewed


ggplot(hair, aes(faculty, hair.length))+
  geom_boxplot()+theme_classic()+
  scale_y_log10()# This looks much better. Thus, we 
# will use the log-transformation in analysis

t.test(log(hair.length)~faculty, data=hair)
# 	Welch Two Sample t-test
# 
# data:  log(hair.length) by faculty
# t = -5.2944, df = 21.075, p-value = 2.967e-05
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#  -1.8933169 -0.8255749
# sample estimates:
# mean in group Law mean in group Sci 
#          1.623870          2.983316 

library(ggplot2)
library(Hmisc)

ggplot(hair, aes(faculty, hair.length))+
  stat_summary(fun.data=mean_cl_normal, geom="pointrange")+
  scale_y_log10()+theme_classic()


## Task 2
## This task covers the topic of species-area relationship (SAR) - one of the principal ecological laws.
# See the theory here: https://en.wikipedia.org/wiki/Species%E2%80%93area_relationship

# The relationship is described by the non-linear equation S=c * A^z
# Our task is to find the c and z parameters by fitting a statistical model
# this can be done by fitting a linear regression on log-transformed data because, 
# after log-transformation, we change the SAR to a linear function:
# S=c * A^z -> log(S) = log(c) + z * log(A), with log(c) being the intercept and z the slope


sar<-read_excel("10_data.xlsx", sheet="2")
summary(sar)


sar.1<-lm(log(sp)~log(area), data=sar)
summary(sar.1)
# lm(formula = log(sp) ~ log(area), data = sar)
# 
# Residuals:
#       Min        1Q    Median        3Q       Max 
# -0.271520 -0.033880  0.005651  0.072867  0.250488 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  5.35276    0.09983   53.62 1.17e-14 ***
# log(area)    0.23335    0.01076   21.68 2.24e-10 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.1425 on 11 degrees of freedom
# Multiple R-squared:  0.9771,	Adjusted R-squared:  0.9751 
# F-statistic: 470.1 on 1 and 11 DF,  p-value: 2.244e-10

# The value of z is directly displayed as the coefficient estimate of log(area)
# The value of c must be computed by back-transforming the intercept like this:
exp(5.35276)
# [1] 211.1904

ggplot(sar, aes(area, sp))+geom_point()+theme_classic()#Plotting the SAR data
# on the linear scale (no easy way how to display the fitted curve in ggplot)

ggplot(sar, aes(area, sp))+geom_point()+theme_classic()+
  geom_smooth(formula=y~x, method="lm")+
  scale_x_log10()+scale_y_log10()#Plotting the SAR on the log scale 

### Task 3
football<-read_excel("10_data.xlsx", sheet="3")
# In this case we should perform a correlation analysis but the data are clearly not
# suitable for Pearson correlation because like many semi-quantitative measures, the variables cannot have normal distribution.
# Therefore, we can use the Spearman correlation here

cor.test(football$league, football$aggressiveness, method="spearman")
# 	Spearman's rank correlation rho
# 
# data:  football$league and football$aggressiveness
# S = 1789.1, p-value = 0.136
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#        rho 
# -0.3452061 
# 
# Warning message:
# In cor.test.default(football$league, football$aggressiveness, method = "spearman") :
#   Cannot compute exact p-value with ties

# Note the warning message here - it says that there are ties (same values repeated 
# several times, which prevents perfect ordering of the values); as a result, the
# p-values are approximated.

# Conclusion: there is no significant association between league level and fan 
# aggressiveness. 


# Task 4
pizza<-read_excel("10_data.xlsx", sheet="4")
summary(pizza)
## In this case, the arrangement of the data is paired (by the evaluator identity), 
# which leads us to a paired-test. The data are again semiquantitative and thus the 
# parametric test is not suitable. We should use Wilcoxon paired test test instead:

wilcox.test(pizza$score[pizza$cook=="Fra"], 
            pizza$score[pizza$cook=="Gia"],paired=T)
# 	Wilcoxon signed rank test with continuity correction
# 
# data:  score by cook
# V = 9.5, p-value = 0.4821
# alternative hypothesis: true location shift is not equal to 0
# 
# Warning messages:
# 1: In wilcox.test.default(x = c(1L, 2L, 1L, 1L, 3L, 1L, 1L, 2L, 1L,  :
#   cannot compute exact p-value with ties
# 2: In wilcox.test.default(x = c(1L, 2L, 1L, 1L, 3L, 1L, 1L, 2L, 1L,  :
#   cannot compute exact p-value with zeroes

# There is no significant difference between the cook. 
# But who wins the competition?
tapply(pizza$score, pizza$cook, mean)
# Fra Gia 
# 1.5 1.8 

# Francesco is the winner

# But actually he is not significantly better than Giacomo. For you as a customer, 
# there is no difference in visiting the restaurants of these two guys - in both
# you should get the same (good) quality. In case of Giacomo, it may actually be cheaper,
# because Francesco may charge a premium for the winner of the pizza competition :-).


## Task 5 Lettuce

lettuce<-read_excel("10_data.xlsx", sheet="5")
summary(lettuce)

wilcox.test(lettuce$taste~lettuce$leaf.col, paired=F)
# 	Wilcoxon rank sum test with continuity correction
# 
# data:  lettuce$taste by lettuce$leaf.col
# W = 2378, p-value = 0.001308
# alternative hypothesis: true location shift is not equal to 0

#Despite called Wilcoxon text, it is actually the Mann-Whitney (unpaired) test. 

# Alternatively, a permutation test may be performed
install.packages("coin")
library(coin)

oneway_test(lettuce$taste~as.factor(lettuce$leaf.col))
# 	Asymptotic Two-Sample Fisher-Pitman Permutation Test
# 
# data:  lettuce$taste by lettuce$leaf.col (green, red)
# Z = 3.2398, p-value = 0.001196
# alternative hypothesis: true mu is not equal to 0

# As you can the results of both tests are very similar and indicate that there is 
# a significant association between leaf color and taste.

tapply(lettuce$taste, lettuce$leaf.col, mean)
#    green      red 
# 2.300000 1.783333 

# Red lettuce is better.

## Task 6
authors<-read_excel("10_data.xlsx", sheet="6")
summary(authors)

ggplot(authors,aes(Author, grade))+geom_boxplot()+theme_classic()

kt.1<-kruskal.test(grade~Author, data=authors)
kt.1
# 	Kruskal-Wallis rank sum test
# 
# data:  grade by Author
# Kruskal-Wallis chi-squared = 9.1957, df = 3, p-value = 0.0268

# There is a significant difference among the authors.

# Pair-wise post-hoc comparison Dunn test
install.packages('FSA')
library(FSA)
dunnTest(authors$grade, authors$Author)
# Dunn (1964) Kruskal-Wallis multiple comparison
#   p-values adjusted with the Holm method.

#            Comparison           Z     P.unadj      P.adj
# 1 Dickens - Hemingway -0.50658151 0.612448485 1.00000000
# 2   Dickens - Pushkin  0.08443025 0.932714356 0.93271436
# 3 Hemingway - Pushkin  0.59101176 0.554512541 1.00000000
# 4   Dickens - Tolstoy  2.27961679 0.022630426 0.11315213
# 5 Hemingway - Tolstoy  2.78619830 0.005333024 0.03199814
# 6   Pushkin - Tolstoy  2.19518654 0.028150219 0.11260088

# The only significant difference in pairwise comparison is Hemingway vs. Tolstoy.

# Alternatively, the main test may be computed using a permutation test
library(coin)
oneway_test(grade~as.factor(Author), data=authors, distribution=approximate(nresample=9999))
# 	Approximative K-Sample Fisher-Pitman Permutation Test
# 
# data:  grade by
# 	 Author (Dickens, Hemingway, Pushkin, Tolstoy)
# chi-squared = 9.4135, p-value = 0.02002