###################################################################
#           Statistics and Explaratoroty Data Analysis            #
#                        Laboratory 9                             #
#                     ANOVA & related                             #
###################################################################

if(!require(Rcpp)){install.packages("Rcpp")}
library(Rcpp) # rcompanion - R / C++ interface
#?Rcpp

if(!require(rcompanion)){install.packages("rcompanion")}
library(rcompanion) # plotNormalHistogram(x), groupWiseMean(x)
# ?`rcompanion-package`

if(!require(psych)){install.packages("psych")}
library(psych) #describe()
# ?psych


if(!require(exactRankTests)){install.packages("exactRankTests")}
library(exactRankTests) #Ansari-Bradley and Wilcoxon exact test with ties

if(!require(DescTools)){install.packages("DescTools")}
library(DescTools) #sign test, gtest


if(!require(normtest)){install.packages("normtest")}
library(normtest) #jb.norm.test & others

if(!require(lattice)){install.packages("lattice")}
library(lattice) #drawings tunning

if(!require(RVAideMemoire)){install.packages("RVAideMemoire")}
library(RVAideMemoire) #mood.medtest

if(!require(FSA)){install.packages("FSA")}
library(FSA) #Summarize

if(!require(ggplot2)){install.packages("ggplot2")}
library(ggplot2) #drawings tuning


#ANOVA

(!require(car)){install.packages("car")}
library(car) #Anova

if(!require(lsmeans)){install.packages("lsmeans")}
library(lsmeans) #lsmeans, cld

#KRUSKAL - WALLIS

if(!require(dunn.test)){install.packages("dunn.test")}
library(dunn.test) #dunn.test

if(!require(conover.test)){install.packages("conover.test")}
library(conover.test) #conover.test


### Mixed models

if(!require(lme4)){install.packages("lme4")}
library(lme4)

if(!require(lmerTest)){install.packages("lmerTest")}
library(lmerTest) #anova after lmer

if(!require(multcompView)){install.packages("multcompView")}
library(multcompView) #multcompLetters

if(!require(phia)){install.packages("phia")}
library(phia) #plots of interactions
?phia


if(!require(nlme)){install.packages("nlme")}
library(nlme)

if(!require(AICcmodavg)){install.packages("AICcmodavg")}
library(AICcmodavg)


#Set English in R
Sys.setenv(LANG = "en")

options(scipen=999) #avoiding e10 notation

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


# ANOVA with blocks


###########################################################
Input <- ("
Year       Student  Wrkday
          '1st year'  a        1200
          '1st year'  b        1400
          '1st year'  c        1350
          '1st year'  d         950
          '1st year'  e        1400
          '1st year'  f        1150
          '1st year'  g        1300
          '1st year'  h        1325
          '1st year'  i        1425
          '1st year'  j        1500
          '1st year'  k        1250
          '1st year'  l        1150
          '1st year'  m         950
          '1st year'  n        1150
          '1st year'  o        1600
          '1st year'  p        1300
          '1st year'  q        1050
          '1st year'  r        1300
          '1st year'  s        1700
          '1st year'  t        1300
          '2nd year'  u        1100
          '2nd year'  v        1200
          '2nd year'  w        1250
          '2nd year'  x        1050
          '2nd year'  y        1200
          '2nd year'  z        1250
          '2nd year'  aa       1350
          '2nd year'  ab       1350
          '2nd year'  ac       1325
          '2nd year'  ad       1525
          '2nd year'  ae       1225
          '2nd year'  af       1125
          '2nd year'  ag       1000
          '2nd year'  ah       1125
          '2nd year'  ai       1400
          '2nd year'  aj       1200
          '2nd year'  ak       1150
          '2nd year'  al       1400
          '2nd year'  am       1500
          '2nd year'  an       1200
          '3rd year'  u        1600
          '3rd year'  v        1700
          '3rd year'  w        1450
          '3rd year'  x        1650
          '3rd year'  y        1800
          '3rd year'  z        1550
          '3rd year'  aa       1950
          '3rd year'  ab       1750
          '3rd year'  ac       1925
          '3rd year'  ad       1825
          '3rd year'  ae       1625
          '3rd year'  af       1525
          '3rd year'  ag       1800
          '3rd year'  ah       1725
          '3rd year'  ai       1200
          '3rd year'  aj       1600
          '3rd year'  ak       1950
          '3rd year'  al       1100
          '3rd year'  am       1400
          '3rd year'  an       1600
          ")


Data <- read.table(textConnection(Input),header=TRUE)
rm(Input)

#Create blocking variable
los<-rbinom(nrow(Data), size = 1, prob=0.7)
Data$Avg<-ifelse(los==1,"Above 4.0","Below 4.0")

#descriptive statistics
Summarize(Wrkday ~ Year + Avg, data=Data, digits=3)

### ANOVA with blocks  - we put blocks directly into model to separate its influence

#linear model  
model = lm(Wrkday ~ Year + Avg, data = Data)
summary(model)   ### Will show overall p-value and r-squared

#anova analysis 
Anova(model, type = "II")   
### Type II sum of squares (A|B) & (B|A)

#diagnostics
#normality
res<- residuals(model)
plotNormalHistogram(res)
shapiro.test(res)

# equality of variance tests for more then 1 IV
Data$YearAvg <- with(Data, interaction(Year,  Avg))
bartlett.test(Wrkday ~ YearAvg, data = Data) #assuming normality
# bartlett.test(Wrkday ~ interaction(Year,Avg), data=Data)
leveneTest(Wrkday ~ YearAvg, data = Data)
fligner.test(Wrkday ~ YearAvg, data = Data)

#post-hoc analysis
detach("package:lmerTest", unload=TRUE) #we need to unload lmerTest package as it's changing lsmeans function for mixed models
library(lsmeans)   
# We are not interesed in blocks so we only compare means for Year
leastsquare = lsmeans(model, pairwise ~ Year, adjust = "tukey")

leastsquare$contrasts

# Results are averaged over the levels of: Avg 
#plot  
CLD = cld(leastsquare, alpha   = 0.05, Letters = letters, adjust  = "tukey")
CLD

### Order the levels for printing
CLD$Year = factor(CLD$Year, levels=c("1st year", "2nd year", "3rd year"))

###  Remove spaces in .group  
CLD$.group=gsub(" ", "", CLD$.group)


### Plot

ggplot(CLD,
       aes(x     = Year, y     = lsmean, label = .group)) +
  geom_point(shape  = 15, size   = 4) +
  geom_errorbar(aes(ymin  =  lower.CL, ymax  =  upper.CL), width =  0.2, size  =  0.7) +
  theme_bw() +
  theme(axis.title   = element_text(face = "bold"), axis.text    = element_text(face = "bold"),  plot.caption = element_text(hjust = 0)) +
  ylab("Least square mean\n Wrkday cost") +
  geom_text(nudge_x = c(0, 0, 0),
            nudge_y = c(150, 150, 150),
            color   = "black")


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


# ANOVA with random blocks   


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

#Mixed-effect model

model = lmer(Wrkday ~ Year + (1|Avg), data=Data, REML=TRUE)
# (1|Avg) - random intercept for each Avg
# REML - restricted maximum likelihood estimator

anova(model)

#Diagnostics
#are residuals normal?

res<- residuals(model)
plotNormalHistogram(res)
shapiro.test(res)

#are variances equal to each other?
# equality of variance tests for more then 2 samples
# H0:all samples have the same variance
# H1:variance is not the same for all samples
bartlett.test(Wrkday ~ Year, data = Data) #assuming normality
# bartlett.test(len ~ interaction(supp,dose), data=ToothGrowth)
leveneTest(Wrkday ~ Year, data = Data)
fligner.test(Wrkday ~ Year, data = Data)

# comparison of a model with null modell
# null model including random effect
model.null = lmer(Wrkday ~ 1 + (1|Avg), data=Data, REML=TRUE)

#2 models comparisons
#H0: General model is statistitcally indifferent to specific model
#H1: Genereal model is better than specific model
anova(model,model.null)

#R2
nagelkerke(fit  = model, null = model.null)

# null model not including random effects
model.null.2 = lm(Wrkday ~ Year, data = Data)
anova(model, model.null.2)
nagelkerke(fit  = model, null = model.null.2)

detach("package:lmerTest", unload=TRUE) #we need to unload lmerTest package as it's changing lsmeans function for mixed models
library(lsmeans)

#POST_HOC
leastsquare = lsmeans(model, pairwise ~ Year,   adjust="tukey")         ###  Tukey-adjusted comparisons
leastsquare$contrasts


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


# 2-way ANOVA     


##################################################################
#Specification as A+B+A*B

Input = ("
         Year        Avg  Wrkday
         '1st year'   A           1200
         '1st year'   A           1400
         '1st year'   A           1350
         '1st year'   A            950
         '1st year'   A           1400
         '1st year'   B           1150
         '1st year'   B           1300
         '1st year'   B           1325
         '1st year'   B           1425
         '1st year'   B           1500
         '1st year'   C           1250
         '1st year'   C           1150
         '1st year'   C            950
         '1st year'   C           1150
         '1st year'   C           1600
         '1st year'   D           1300
         '1st year'   D           1050
         '1st year'   D           1300
         '1st year'   D           1700
         '1st year'   D           1300
         '2nd year'   A           1100
         '2nd year'   A           1200
         '2nd year'   A           1250
         '2nd year'   A           1050
         '2nd year'   A           1200
         '2nd year'   B           1250
         '2nd year'   B           1350
         '2nd year'   B           1350
         '2nd year'   B           1325
         '2nd year'   B           1525
         '2nd year'   C           1225
         '2nd year'   C           1125
         '2nd year'   C           1000
         '2nd year'   C           1125
         '2nd year'   C           1400
         '2nd year'   D           1200
         '2nd year'   D           1150
         '2nd year'   D           1400
         '2nd year'   D           1500
         '2nd year'   D           1200
         '3rd year'  A            900
         '3rd year'  A           1100
         '3rd year'  A           1150
         '3rd year'  A            950
         '3rd year'  A           1100
         '3rd year'  B           1150
         '3rd year'  B           1250
         '3rd year'  B           1250
         '3rd year'  B           1225
         '3rd year'  B           1325
         '3rd year'  C           1125
         '3rd year'  C           1025
         '3rd year'  C            950
         '3rd year'  C            925
         '3rd year'  C           1200
         '3rd year'  D           1100
         '3rd year'  D            950
         '3rd year'  D           1300
         '3rd year'  D           1400
         '3rd year'  D           1100
         ")


Data <- read.table(textConnection(Input),header=TRUE)


#Preanalysis
#descriptive stats
Summarize(Wrkday ~ Year+Avg, data=Data, digits=3)

#2-way ANOVA
#define a model Year:Avg is interaction
model <- lm(Wrkday ~ Year + Avg + Year:Avg, data = Data)
summary(model)
anova(model)
#anova analysis type 2 - why?
Anova(model,type = "II")
Anova(model,type = "III")
#interactions not significant

model <- lm(Wrkday ~ Year + Avg, data = Data)

#Diagnostics
res<- residuals(model)
plotNormalHistogram(res)
shapiro.test(res)

# equality of variance tests for more than 2 samples
# H0:all samples have the same variance
# H1:variance is not the same for all samples

bartlett.test(Wrkday ~ interaction(Year,Avg), data=Data)
leveneTest(Wrkday ~ interaction(Year,Avg), data = Data)
fligner.test(Wrkday ~ interaction(Year,Avg), data = Data)



#### POST-HOC LS MEANS
detach("package:lmerTest", unload=TRUE) #we need to unload lmerTest package as it's changing lsmeans function for mixed models
library(lsmeans)

#means for years
lsYear <- lsmeans(model, pairwise ~ Year, adjust = "tukey")
lsYear$contrasts

CLDYear = cld(lsYear, alpha   = 0.05,   Letters = letters,   adjust  = "tukey")         
CLDYear
CLDYear$.group=gsub(" ", "", CLDYear$.group)

ggplot(CLDYear,
       aes(x     = Year, y     = lsmean, label = .group)) +
  geom_point(shape  = 15,   size   = 4) +
  geom_errorbar(aes(ymin  =  lower.CL,     ymax  =  upper.CL),  width =  0.2,  size  =  0.7) +
  theme_bw() +
  theme(axis.title   = element_text(face = "bold"),  axis.text    = element_text(face = "bold"),  plot.caption = element_text(hjust = 0)) +
  ylab("Least square means Years") +
  geom_text(nudge_x = c(0, 0, 0), nudge_y = c(120, 120, 120),   color   = "black")


#means for AVG  
lsAvg <- lsmeans(model, pairwise ~ Avg, adjust = "tukey")
lsAvg$contrasts

CLDAvg = cld(lsAvg,  alpha   = 0.05,    Letters = letters,  adjust  = "tukey")   
CLDAvg
CLDAvg$.group=gsub(" ", "", CLDAvg$.group)  

ggplot(CLDAvg,
       aes(x     = Avg, y     = lsmean, label = .group)) +
  geom_point(shape  = 15,   size   = 4) +
  geom_errorbar(aes(ymin  =  lower.CL,     ymax  =  upper.CL),  width =  0.2,  size  =  0.7) +
  theme_bw() +
  theme(axis.title   = element_text(face = "bold"),  axis.text    = element_text(face = "bold"),  plot.caption = element_text(hjust = 0)) +
  ylab("Least square means AVG") +
  geom_text(nudge_x = c(0, 0, 0), nudge_y = c(120, 120, 120),   color   = "black")


### What to do if interactions are significant?
# present an interaction plot of the means for both main effects together on one plot. 

leastsquare = lsmeans(model, pairwise ~ Year + Avg, adjust = "tukey")
CLD = cld(leastsquare,  alpha   = 0.05,  Letters = letters,  adjust  = "tukey")   
###  Remove spaces in .group  
CLD$.group=gsub(" ", "", CLD$.group)
CLD

### Plot

pd = position_dodge(0.6)    ### How much to jitter the points on the plot
ggplot(CLD,
       aes(x     = Year, y     = lsmean, color = Avg,  label = .group)) +
  geom_point(shape  = 15, size   = 4,  position = pd) +
  geom_errorbar(aes(ymin  =  lower.CL,  ymax  =  upper.CL),  width =  0.2,  size  =  0.7,  position = pd) +
  theme_bw() +
  theme(axis.title   = element_text(face = "bold"),  axis.text    = element_text(face = "bold"),  plot.caption = element_text(hjust = 0)) +
  ylab("Least square means Year*Avg") +
  geom_text(nudge_x = c(-0.22, 0.08, -0.22, 0.22, 0.08, 0.08, -0.08, -0.22, 0.22, 0.22, -0.08, -0.08),nudge_y = rep(270, 12),
            color   = "black") +
  scale_color_manual(values = c("steelblue", "steelblue2","springgreen4", "springgreen3"))

### Produce interaction plot
IM = interactionMeans(model)
IM

### Produce interaction plot
plot(IM)

### Return the graphics device to its default 1-plot-per-window state
par(mfrow=c(1,1))




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



# Reapeated measures ANOVA      



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


Input = ("
         Wrk_type        Student  Month   Salary
         'Journalist'     a        1       2000
         'Journalist'     a        2       1978
         'Journalist'     a        3       1962
         'Journalist'     a        4       1873
         'Journalist'     a        5       1782
         'Journalist'     a        6       1737
         'Journalist'     b        1       1900
         'Journalist'     b        2       1826
         'Journalist'     b        3       1782
         'Journalist'     b        4       1718
         'Journalist'     b        5       1639
         'Journalist'     b        6       1644
         'Journalist'     c        1       2100
         'Journalist'     c        2       2067
         'Journalist'     c        3       2065
         'Journalist'     c        4       2015
         'Journalist'     c        5       1994
         'Journalist'     c        6       1919
         'Journalist'     d        1       2000
         'Journalist'     d        2       1981
         'Journalist'     d        3       1987
         'Journalist'     d        4       2016
         'Journalist'     d        5       2010
         'Journalist'     d        6       1946
         'Data Scientist'     e        1       2100
         'Data Scientist'     e        2       2004
         'Data Scientist'     e        3       2027
         'Data Scientist'     e        4       2109
         'Data Scientist'     e        5       2197
         'Data Scientist'     e        6       2294
         'Data Scientist'     f        1       2000
         'Data Scientist'     f        2       2011
         'Data Scientist'     f        3       2089
         'Data Scientist'     f        4       2124
         'Data Scientist'     f        5       2199
         'Data Scientist'     f        6       2234
         'Data Scientist'     g        1       2000
         'Data Scientist'     g        2       2074
         'Data Scientist'     g        3       2141
         'Data Scientist'     g        4       2199
         'Data Scientist'     g        5       2265
         'Data Scientist'     g        6       2254
         'Data Scientist'     h        1       2000
         'Data Scientist'     h        2       1970
         'Data Scientist'     h        3       1951
         'Data Scientist'     h        4       1981
         'Data Scientist'     h        5       1987
         'Data Scientist'     h        6       1969
         'Other'     i        1       1950
         'Other'     i        2       2007
         'Other'     i        3       1978
         'Other'     i        4       1965
         'Other'     i        5       1984
         'Other'     i        6       2020
         'Other'     j        1       2000
         'Other'     j        2       2029
         'Other'     j        3       2033
         'Other'     j        4       2050
         'Other'     j        5       2001
         'Other'     j        6       1988
         'Other'     k        1       2000
         'Other'     k        2       1976
         'Other'     k        3       2025
         'Other'     k        4       2047
         'Other'     k        5       2033
         'Other'     k        6       1984
         'Other'     l        1       2000
         'Other'     l        2       2020
         'Other'     l        3       2009
         'Other'     l        4       2017
         'Other'     l        5       1989
         'Other'     l        6       2020
         ")      


Data = read.table(textConnection(Input),header=TRUE)

#test whether there are significant differences between wrk_type/months/their interaction
#Note: We could run the model using gls as:
#Generalized Least Squares estimation 
# without random effects

# model = gls(Salary ~ Wrk_type + Month + Wrk_type*Month,
#             correlation = corAR1(form = ~ Month | Student,
#                                  value = 0.8990),
#             data=Data,
#             method="REML")
# 
# library(car)
# 
# Anova(model)

#Instead we may also use mixed effects model
#As opposed to gls in the mixed model below
#Students are treated as random variables - this takes out any differences between students

#Mixed effects model:
model = lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~1|Student, correlation = corAR1(form = ~ Month | Student), data=Data, method="REML")
summary(model)

# correlation = Structure(form  = ~ time | subjvar)
#   Structure -  autocorrelation structure.
#     corAR1	- autoregressive process of order 1.
#     corARMA	- autoregressive moving average process, with arbitrary orders for the autoregressive and moving average components.
#     corCompSymm - compound symmetry structure corresponding to a constant correlation.
#     corSymm	-  general correlation matrix, with no additional structure.
# time - time variable 
# subjvar - unit id Autocorrelation is modeled within levels of the subjvar, and not between them.
# 

#to test the significance of the random effects we would compare the models with and without the random effects:
model.fixed = gls(Salary ~ Wrk_type + Month + Wrk_type*Month, data=Data, method="REML")
anova(model, model.fixed)


# AIC selection of correlation structure
# Model selection for fixed effects is only appropriate with method=ML
# REML (default) should only be used to select random effects for a constant set of fixed effects

Cand.models <- list( )

Cand.models[[1]] <- lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~1|Student, data=Data, method="ML")
Cand.models[[2]] <- lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~Month|Student, data=Data, method="ML")
Cand.models[[3]] <- lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~1|Student, correlation = corCompSymm(form = ~ Month | Student), data=Data, method="ML")
Cand.models[[4]] <- lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~Month|Student, correlation = corCompSymm(form = ~ Month | Student), data=Data, method="ML")
Cand.models[[5]] <- lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~1|Student, correlation = corAR1(form = ~ Month | Student), data=Data, method="ML")
Cand.models[[6]] <- lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~1|Student, correlation = corSymm(form = ~ Month | Student), data=Data, method="ML")


##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")

##generate AICc table
aictab<-aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)
#which model is the best?
aictab

#the best is model 5
model = lme(Salary ~ Wrk_type + Month + Wrk_type:Month, random = ~1|Student, correlation = corAR1(form = ~ Month | Student), data=Data, method="REML")
summary(model)

# diagnostics

#model significance check
model.null.2 = gls(Salary ~ 1, data = Data)
nagelkerke(model, model.null.2)

#residuals normality check
x = residuals(model)
plotNormalHistogram(x)
shapiro.test(x)


#anova analysis

#if interactions are significant
#2 is disputable - because we omit interactions which are important
Anova(model, type="II") # (A|B),(B|A), (AB|A,B)

#3 type is appropriate
Anova(model, type="III") # (A|B,AB),(B|A,AB), (AB|A,B)

#Interpretation
#Interaction Wrk_type:Month - difference in month trends in salaries for different Wrk_type
#Linear trend in Salaries for base Wrk_type (Data Scientist)
#Different Salaries for Wrk_type at Month 1

#post-hoc analysis
detach("package:lmerTest", unload=TRUE) #we need to unload lmerTest package as it's changing lsmeans function for mixed models
library(lsmeans)

lsmeans(model,  pairwise ~ Wrk_type:Month,  adjust="holm")

#Because Month is an integer variable, not a factor variable, it is listed in the lsmeans cld table as its average only.

#plot 
Sum = groupwiseMean(Salary ~ Wrk_type + Month, data   = Data,  conf   = 0.95, digits = 3,   traditional = FALSE,    percentile  = TRUE)
Sum

pd = position_dodge(.2)

ggplot(Sum, aes(x =    Month, y =Mean, color = Wrk_type)) +
  geom_errorbar(aes(ymin=Percentile.lower,  ymax=Percentile.upper), width=.2, size=0.7, position=pd) +
  geom_point(shape=15, size=4, position=pd) +
  theme_bw() +
  theme(axis.title = element_text(face = "bold")) +
  ylab("Mean salaries per work type")


#EXTRA
# # What to do if we want to check wheter time has an influence on each wrk_type
# # and when the significant difference between wrk_type starts?
# 
# #ordered factor
# Data$Mnt<-factor(Data$Month,ordered=T)
# #model with ordered factor (note: we do no assume linear influence month)
# model = lme(Salary ~ Wrk_type + Mnt + Wrk_type:Mnt, random = ~1|Student, correlation = corAR1(form = ~ Month | Student), data=Data, method="REML")
# summary(model)
# 
# 
# #adjustment for multiple comparisons
# lsmeans(model,  pairwise ~ Wrk_type|Mnt, adjust="tukey")
# lsmeans(model,  pairwise ~ Mnt|Wrk_type,  adjust="tukey")

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


# Friedman & Quade tests


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


Input <- ("
          Wrk_type	Month	Salary
          'Journalist'	    1	  2000
          'Journalist'	    2	  1963
          'Journalist'	    3	  1949
          'Journalist'		  4	  1905.5
          'Journalist'	    5	  1856.25
          'Journalist'		  6	  1811.5
          'Data Scientist'	1	  2025
          'Data Scientist'	2	  2014.75
          'Data Scientist'	3	  2052
          'Data Scientist'	4	  2103.25
          'Data Scientist'	5	  2162
          'Data Scientist'	6	  2187.75
          'Other'		        1	  1987.5
          'Other'		        2	  2008
          'Other'		        3	  2011.25
          'Other'		        4	  2019.75
          'Other'		        5	  2001.75
          'Other'		        6	  2003
          ")


Data <- read.table(textConnection(Input),header=TRUE)

# The test is based on the differences between paired data
# The function pairwiseDifferences will create a new data frame of differences for all pairs of differences.  
# Note that the data must be ordered by the blocking variable so that the first observation 
# for 1st year will be paired with the first observation for 2nd year, and so on.

Data.diff = pairwiseDifferences(Salary ~ Wrk_type, data= Data, factorize=TRUE, plotit=F)
# plotit	-   If TRUE, then produces bar plots of the differences.
# factorize	- If TRUE, then adds a column to the output data frame consisting of the differences as a factor variable. 

#are differences symmetric?
histogram(~ Difference | Comparison,data=Data.diff,  type = "count",layout=c(1,3))      #  columns and rows of individual plots

#descriptive statistics, are distirbutions similar in spread and shape?
Summarize(Salary ~ Wrk_type, data=Data, digits=3)
describeBy(Data$Salary, group=Data$Wrk_type)

leveneTest(Salary ~ Wrk_type, data = Data)
fligner.test(Salary ~ Wrk_type, data = Data)


#friedman test
#H0:Same medians/same distributions

friedman.test(Salary ~ Wrk_type | Month, data = Data)

#POST-HOC
PT = pairwiseSignTest(Salary ~ Wrk_type, data   = Data,  method = "fdr")
PT
### Compact letter display
cldList(p.adjust ~ Comparison,data = PT,threshold  = 0.05)


#Quade test
#H0:Same medians/same distributions

quade.test(Salary ~ Wrk_type | Month, data = Data)

#CLD

PT = pairwise.wilcox.test(Data$Salary, Data$Wrk_type,  p.adjust.method="fdr",  paired=TRUE)$p.value
PT

#Convert a lower triangle matrix to a full matrix
PT1 = fullPTable(PT)
PT1

#Letter summary of similarities and differences
multcompLetters(PT1,  compare="<", threshold=0.05, Letters=letters)

# X symmetric matrix with p-values
# comapre - function for comapre if "<" then threshold than different
# Letters=letters - Vector of distinct characters (or character strings) 
# used to connect levels that are not significantly different. 


