# During the research of the content of various elements in soils,
# the concentration of arsenic in soils from four localities
# (A, B, C and D) was examined. 
# As the arsenic content was assumed to be very closely
# related to the aluminum content, the aluminum concentration
# in the samples was also determined. 
# The task of the analysis is to determine whether the content
# of arsenic in soils is different in individual localities
# and whether it is also related to the concentration of aluminum.
# A logarithmic transformation was performed in the data
# because of linearity.


#data input
soil=read.table("http://user.mendelu.cz/drapela/Forest_Biometry/Data/arsenic.txt",header=TRUE)
soil
soil=data.frame(soil)
attach(soil)

# additional libraries
library(ggpubr)
library(gplots)
library(HH)
library(rstatix)
library(emmeans)


# 1-F ANOVA 
aov=aov(LN_AS~Locality)
summary(aov)
plotmeans(LN_AS~Locality)
as.al=lm(LN_AS~LN_AL)
summary(as.al)
# dependance of AS on AL for all data
plot(LN_AS~LN_AL)
abline(as.al)
as.al1=data.frame(x=LN_AL,y=LN_AS)
lm1=lm(y ~ x, data=as.al1)
ci.plot(lm1)


# categorised scatter plot and regression models
ggscatter(
  soil, x = "LN_AL", y = "LN_AS",
  color = "Locality", add = "reg.line"
)+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = Locality)
  )

#test of parallelism of slopes
paralel=aov(LN_AS~LN_AL*Locality)
summary(paralel)


#ANCOVA
ancova_soil=soil %>%anova_test(LN_AS~LN_AL+Locality)
ancova_soil


# test od significance of adjusted means (estimated marginal means)
library(emmeans)
pwc <- soil %>% 
  emmeans_test(
    LN_AS ~ Locality, covariate = LN_AL,
    p.adjust.method = "bonferroni"
  )
pwc


# graphical comparison of adjusted means
pwc <- pwc %>% add_xy_position(x = "Locality", fun = "mean_se")
ggline(get_emmeans(pwc), x = "Locality", y = "emmean") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) + 
  stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) +
  labs(
    subtitle = get_test_label(ancova_soil, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )




detach(soil)