###### Task 1
library(ggplot2)
library(readxl) # loading the readxl library for the read_excel function

termometers<-read_excel("08_data.xlsx", sheet="1")

summary(termometers)

#1a F-test of variances
var.test(temp~manufacturer, data=termometers)
# 	F test to compare two variances
# 
# data:  temp by manufacturer
# F = 6.0674, num df = 9, denom df = 9, p-value = 0.01302
# alternative hypothesis: true ratio of variances is not equal to 1
# 95 percent confidence interval:
#   1.50706 24.42738
# sample estimates:
# ratio of variances 
#           6.067416

tapply(termometers$temp, termometers$manufacturer, var)#Computing variances for 
# the groups defined by manufacturer

### Conclusion: Variances of indicated temperature significantly differ. Thermom-produced
# thermometers show higher variability, i.e. inaccuracy

#1b testing of mean difference
# We have two groups here so we may use a t-test or ANOVA.
# t-test (with assumption of equality of variance - this is unrealistic but important for
# comparison with ANOVA result).
t.test(temp~manufacturer, data=termometers, var.equal=T)
# 	Two Sample t-test
# 
# data:  temp by manufacturer
# t = 0.11962, df = 18, p-value = 0.9061
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#  -1.656362  1.856362
# sample estimates:
# mean in group celsimet   mean in group termom 
#                   18.0                   17.9 

summary(aov(temp~manufacturer, data=termometers)) 
#              Df Sum Sq Mean Sq F value Pr(>F)
# manufacturer  1   0.05   0.050   0.014  0.906
# Residuals    18  62.90   3.494 

# As you can see, the p-values are exactly the same, which is because the two approaches
# are fully equivalent here. They agree that there is no significant difference in 
# mean temperature.

# Overall conclusion: the two termometers do not show a significant bias when compared to each other but
# they Termom is significanlty more precise. 

## BTW - even Termom is still quite imprecise.

# Graphical illustration of these data may be done with a boxplot:
ggplot(termometers, aes(x = manufacturer, y = temp))+geom_boxplot(fill="grey")+
           theme_classic()+labs(x="Manufacturer", y = "Temperature (°C)")

# or maybe a dotchart with CIs
ggplot(termometers, aes(x = manufacturer, y = temp))+stat_summary(fun.data=mean_cl_normal, geom="pointrange")+
  theme_classic()+labs(x="Manufacturer", y = "Temperature (°C)")


####### Task 2
fert<-read_excel("08_data.xlsx", sheet="2")
summary(fert)

aov.f<-aov(Biomass~as.factor(Fertilizer), data=fert)

plot(aov.f) # Checking the ANOVA assumptions. Seem reasonably fulfilled here.

summary(aov.f)
#             Df Sum Sq Mean Sq F value  Pr(>F)    
# Fertilizer   4  26740    6685   45.96 3.6e-11 ***
# Residuals   25   3636     145                    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## Conclusion: There is a significant effect of fertilizer on biomass production. 

#Explained variability
26740/(26740+3636)
# [1] 0.8803002

# We may proceed to post-hoc comparisons to see which means are significantly different from which
# Post-hoc multiple comparisons by Tukey HSD test

TukeyHSD(aov.f)
#   Tukey multiple comparisons of means
#     95% family-wise confidence level
# 
# Fit: aov(formula = Biomass ~ Fertilizer, data = fert)
# 
# $Fertilizer
#                               diff        lwr        upr     p adj
# Mineral-Control           69.83333  49.383672  90.282994 0.0000000
# MinSlow-Control           86.00000  65.550339 106.449661 0.0000000
# NitrateAmm-Control        47.16667  26.717006  67.616328 0.0000040
# OrganicManure-Control     68.83333  48.383672  89.282994 0.0000000
# MinSlow-Mineral           16.16667  -4.282994  36.616328 0.1713159
# NitrateAmm-Mineral       -22.66667 -43.116328  -2.217006 0.0245596
# OrganicManure-Mineral     -1.00000 -21.449661  19.449661 0.9998971
# NitrateAmm-MinSlow       -38.83333 -59.282994 -18.383672 0.0000769
# OrganicManure-MinSlow    -17.16667 -37.616328   3.282994 0.1311420
# OrganicManure-NitrateAmm  21.66667   1.217006  42.116328 0.0340129

# Control: a
# Nitr. Amm: b
# OrganicManure, Mineral, Min. Slow: c


# Graphical illustrations
# One option is to use a boxplot
ggplot(fert, aes(x = Fertilizer, y = Biomass))+geom_boxplot(fill="grey")+
           theme_classic()+labs(x="Fertilizer type", y = "Biomass DW (g)")+
stat_summary(fun = "max", geom = "text", label= c("a", "c", "c", "b", "c"),
             hjust=-0.5)


# Another option is to use a dotchart 
# 
# We need top compute the upper CI first for placemen of the Tukey letters
CIs<-tapply(fert$Biomass, list(fert$Fertilizer), mean_cl_normal)
upperCI<-as.numeric(lapply(CIs, max))
#This is then used in annotate as the y-coordinate

ggplot(fert, aes(x = Fertilizer, y = Biomass))+
  stat_summary(fun.data=mean_cl_normal, geom="pointrange")+
           theme_classic()+labs(x="Fertilizer type", y = "Biomass DW (g)")+
annotate(x=1:5, geom="text", label= c("a", "c", "c", "b", "c"), y=upperCI, vjust=-0.5)

## Final conclusion: There is a significant effect of fertilizer on biomass production. 
# All fertilizers significantly increase biomass production compared to control and 
# mineral, mineral slow and organic manure are significantly better in this than nitrate ammonium.

# Task C
rats<-read_excel("08_data.xlsx", sheet="C")
summary(rats)

aov.r<-aov(time~nutrition, data=rats)

plot(aov.r)

summary(aov.r)
#           Df Sum Sq Mean Sq F value   Pr(>F)    
# nutrition    3   5607  1869.0   19.78 1.24e-05 ***
# Residuals   16   1512    94.5                     
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(aov.r)
#                diff         lwr       upr     p adj
# control-beer   16.6  -0.9900376 34.190038 0.0678509
# fat-beer       46.4  28.8099624 63.990038 0.0000064
# sugar-beer     25.6   8.0099624 43.190038 0.0036699
# fat-control    29.8  12.2099624 47.390038 0.0009242
# sugar-control   9.0  -8.5900376 26.590038 0.4805645
# sugar-fat     -20.8 -38.3900376 -3.209962 0.0178935

# beer: a
# control: ab
# sugar: b
# fat: c

# Graphical illustration by a dotchart
CIs<-tapply(rats$time, list(rats$nutrition), mean_cl_normal)
upperCI<-as.numeric(lapply(CIs, max))

ggplot(rats, aes(x=nutrition, y=time))+stat_summary(fun.data="mean_cl_normal",
                  geom="pointrange")+theme_classic()+
annotate(x=1:4, y=upperCI, geom="text", label=c("a", "ab", "c", "b"), vjust=-0.5)



