kyjov<- read.delim("c:\\Users\\lubop\\Documents\\kurs\\kyjov.txt", header=T)
xk = kyjov[which(kyjov[,3]=="1034"),]
x = xk[,6]

library(mclust)
xs<- sort(x)
dmc <- Mclust(xs,G=NULL)
dmcBIC <- mclustBIC(x); dmcBIC
dmcSummary <- summary(dmcBIC, data = xs); dmcSummary

nn = 3
dmc <- Mclust(xs,G=nn)
dmcz <-round(dmc$z,1); dmcz
dmcc <-dmc$classification; dmcc
dmcp <-dmc$parameters; dmcp
dmcpp <-dmc$parameters$pro; dmcpp
dmcpm <-dmc$parameters$mean; dmcpm
dmcps <-sqrt(as.vector(dmcp$variance$sigmasq));dmcps
xss = seq(min(x), max(x), 0.01)
PDXX = dmcp$pro[1]*dnorm(xss,dmcpm[1],dmcps[1])+dmcp$pro[2]*dnorm(xss,dmcpm[2],dmcps[2])+dmcp$pro[3]*dnorm(xss,dmcpm[3],dmcps[3])
plot(xss,PDXX,type="l",col=1, xlab="koncentrace P2O5 (ug/g)",ylab="Density")
rug(xs)
PDX1 = dmcp$pro[1]* dnorm(xss,dmcpm[1], dmcps[1]); points(xss,PDX1,type="l",col=3)
PDX2 = dmcp$pro[2]* dnorm(xss,dmcpm[2], dmcps[2]); points(xss,PDX2,type="l",col=2)
PDX3 = dmcp$pro[3]* dnorm(xss,dmcpm[3],dmcps[3]); points(xss,PDX3,type="l",col=4)

hw<- bw.SJ(xs, method = "ste"); hw
# method = c("ste", "dpi")
dx<-density(xs, bw=hw, adjust = 1,kernel = "gaussian", weights = NULL)  
# kernel = c("gaussian", "epanechnikov", "rectangular","triangular", "biweight","cosine", "optcosine")
xx = dx$x
yy = dx$y
plot(dx$x,dx$y,type="l",xlab="koncentrace P2O5 (ug/g)",ylab="Density")
rug(xs)

a1 = dmcp$pro[1];a1
a2 = dmcpm[1];a2   
a3 = dmcps[1];a3
a4 = dmcp$pro[2];a4
a5 = dmcpm[2];a5   
a6 = dmcps[2];a6
a7 = dmcp$pro[3];a7
a8 = dmcpm[3];a8   
a9 = dmcps[3];a9

library(nls2)
fitG1 <- nls2(yy~(C1*exp(-(xx-mean1)**2/(2*sigma1**2))+ C2*exp(-(xx-mean2)**2/(2*sigma2**2))+ C3*exp(-(xx-mean3)**2/(2*sigma3**2))), 
              start=list(C1=a1,mean1=a2,sigma1=a3,C2=a4,mean2=a5,sigma2=a6,C3=a7,mean3=a8,sigma3=a9),data=as.data.frame(cbind(xx,yy)))
#,lower=c(0,-Inf,0,0,-Inf,0,0,-Inf,0),upper=c(Inf,0,Inf,Inf,0,Inf,Inf,0,Inf))
summary(fitG1)
coef(fitG1)
a1=coef(fitG1)[1]; a2=coef(fitG1)[2];a3=coef(fitG1)[3];
a4=coef(fitG1)[4]; a5=coef(fitG1)[5]; a6=coef(fitG1)[6];
a7=coef(fitG1)[7]; a8=coef(fitG1)[8]; a9=coef(fitG1)[9]

gauss <- a1*exp(-(xx-a2)**2/(2*a3**2))+a4*exp(-(xx-a5)**2/(2*a6**2)+a7*exp(-(xx-a8)**2/(2*a9**2)))
gauss1 <- a1*exp(-(xx-a2)**2/(2*a3**2))
gauss2 <- a4*exp(-(xx-a5)**2/(2*a6**2))
gauss3 <- a7*exp(-(xx-a8)**2/(2*a9**2))

plot(xx,yy,type="l",xlab="koncentrace P2O5 (ug/g)",ylab="Density")
rug(xs)
points(xx,gauss,type="l",col=6,lwd=2)
points(xx,gauss1,type="l",col=3,lwd=2)
points(xx,gauss2,type="l",col=2,lwd=2)
points(xx,gauss3,type="l",col=4,lwd=2)


### Parove testy ##################################################################

# Weight of the mice before treatment
before <-c(200.1, 190.9, 192.7, 213, 241.4, 196.9, 172.2, 185.5, 205.2, 193.7)
summary(before)
# Weight of the mice after treatment
after <-c(392.9, 393.2, 345.1, 393, 434, 427.9, 422, 383.9, 392.3, 352.2)
summary(after)
my_data <- data.frame(weight=c(before,after),group=rep(c("before","after"),each=10))


my_data$group <- factor(my_data$group, levels=c("before","after"))
boxplot(my_data$weight~my_data$group)

library(PairedData)
pd <- paired(before, after)
plot(pd, type = "profile")

d <- before-after
shapiro.test(d) # Shapiro-Wilk normality test for the differences

t.test(before, after, paired = TRUE)
t.test(weight ~ group, data = my_data, paired = TRUE)
# average weight before treatment is less than the average weight after treatment
t.test(weight ~ group, data = my_data, paired = TRUE, alternative = "less")
# average weight before treatment is greater than the average weight after treatment
t.test(weight ~ group, data = my_data, paired = TRUE, alternative = "greater")


library(fmsb)
library(parviol)


### Analyza rozptylu #####################################################################################

# Byl sledován vliv druhu konzervacního cinidla A1 az A4 (faktor A) na trvanlivost potravinárského výrobku. 
# Vysetrete, zda na hladine významnosti 0.05 závisí trvanlivost potravin na druhu konzervacního cinidla. 
# Který druh konzervacního cinidla dosahuje silne odlisných výsledku vuci ostatním?

E503 <- read.table("C:/Users/lubop/Documents/kurs/Data1/ascii05/E/E503.txt")[-1,]
E503
data = stack(E503)
data
data = as.data.frame(data); colnames(data)=c("elem","lom")
elem = data[,1]
lom = data[,2]

### Distribucni grafy

boxplot(E503,ylab="conc",xlab="lab",las=1,cex=1,outline=TRUE) 
stripchart(E503, pch=1, col=1,vert=TRUE,cex=1,add=TRUE)

boxplot(elem~lom,ylab="conc",xlab="lab",las=1,cex=1,outline=TRUE) 
stripchart(elem~as.factor(lom), pch=1, col=1,vert=TRUE,cex=1,add=TRUE)


library(beeswarm)
boxplot(E503,ylab="",xlab="",las=1,cex=1,outline=TRUE,cex.lab=1) 
beeswarm(E503, col = 1, pch = 1, add = TRUE,horizontal=FALSE,pwcol = NULL,cex=1)

boxplot(elem~lom,ylab="",xlab="",las=1,cex=1,outline=TRUE,cex.lab=1) 
beeswarm(elem~as.factor(lom), col = 1, pch = 1, add = TRUE,horizontal=FALSE,pwcol = NULL,cex=1)


library(beanplot)
beanplot(E503, col = "bisque", ylim = "")
beanplot(elem ~ lom, data = data, col = "bisque", ylim = "")



#### Global test of model assumptions
library(gvlma)
gvmodel <- gvlma(lm(elem~lom)) 
summary(gvmodel)
gvmodel <- gvlma(lm(elem~lom-1)) 
summary(gvmodel)


### outliers

# Bonferroni Outlier Test
library(car)
ou = outlierTest(lm(elem~lom), cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))

# Lund Outlier Test
# modified https://stackoverflow.com/questions/1444306/how-to-use-outlier-tests-in-r-code/1444548#1444548
lundout<-function(x1,x2,a) {
  testoutlm<-lm(x1~x2)
  standardresid<-rstandard(testoutlm)
  nl<-length(x1)
  q<-length(testoutlm$coefficients)
  F<-qf(c(1-(a/nl)),df1=1,df2=nl-q-1,lower.tail=TRUE)
  crit<-((nl-q)*F/(nl-q-1+F))^0.5
  oul = which(abs(standardresid)>crit)
  return(oul)
}
lundout(elem,lom,0.05)

library(referenceIntervals)
vanderLoo.outliers(stresid)
cook.outliers(stresid)
dixon.outliers(stresid)
horn.outliers(stresid)



#### shoda rozptylu 

## Bartlett
bartlett.test(elem~as.factor(lom))
library(RVAideMemoire)
pairwise.var.test(elem,as.factor(lom))

library(onewaytests)
homog.test(elem~lom, data = data, method = "Bartlett", alpha = 0.05, na.rm = TRUE, verbose = TRUE)

library(RVAideMemoire)
perm.bartlett.test(elem~lom, nperm = 999, progress = TRUE)
pairwise.perm.var.test(elem,as.factor(lom),nperm=999)


## Fligner Kileen
fligner.test(elem~lom)

library(onewaytests)
homog.test(elem~lom, data = data, method = "Fligner", alpha = 0.05, na.rm = TRUE, verbose = TRUE)

library(coin)
fligner_test(elem~lom, data=data, subset = NULL, weights = NULL,  conf.int = TRUE, conf.level = 0.95, ties.method = "average-scores")
fligner_test(elem~lom, data=data, subset = NULL, weights = NULL,  conf.int = TRUE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000))
# ties.method = c("mid-ranks", "average-scores"),


## Conover, Johnson, Johnson
library(coin)
conover_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores")
conover_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores", distribution = approximate(nresample = 1000))


## Levene
library(lawstat)
lawstat::levene.test(elem, group = as.factor(lom), location="mean", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
# location=c("median", "mean", "trim.mean")
# correction.method=c("none","correction.factor","zero.removal","zero.correction")

library(car)
car::leveneTest(elem~as.factor(lom),center=mean)
car::leveneTest(elem~as.factor(lom),center=median)
# center = c("median", "mean")

library(DescTools)
DescTools::LeveneTest(elem~as.factor(lom),center=median)
DescTools::LeveneTest(elem~as.factor(lom),center=mean)

library(onewaytests)
homog.test(elem ~ lom, data = data, method = "Levene", alpha = 0.05, na.rm = TRUE, verbose = TRUE)

library(s20x)
s20x::levene.test(elem ~ lom, data=data, digit = 5, show.table = TRUE)

library(asbio)
modlevene.test(elem,lom)


## Brown Forsythe
library(HH)
HH::hov(elem~lom,method = "bf")

library(onewaytests)
bf.test(elem ~ lom, data = data, alpha = 0.05, na.rm = TRUE, verbose = TRUE)


# Hartley test
library(PMCMRplus) 
hartleyTest(elem~lom)


# Cochran C test
library(outliers)
cochran.test(elem~lom)

library(GAD)
C.test(lm(elem ~ lom))

library(PMCMRplus)
cochranTest(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "greater")
# alternative = c("greater", "less")

library(outliers) # Cochran
cochran.test(elem~lom, inlying = FALSE)
cochran.test(elem~lom, inlying = TRUE)


library(PMCMRplus) 
GSTTest(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided")
# dist = c("Chisquare", "KruskalWallis")

ans <- aov(elem ~ as.factor(lom)) # spravnejsi 
GSTTest(residuals(ans) ~ as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided")


### ANOM for variances
library(ANOM)
# Absolute deviations from the median (robust Levene residuals)
mmedian <- tapply(elem, lom, median)[lom]
absdev <- abs(elem - mmedian)
dataf = data.frame(cbind(absdev,lom))
dataf[,1] = as.numeric(absdev)
dataf[,2] = as.factor(lom)
library(multcomp)
mod <- lm(absdev~lom, dataf)
test <- glht(mod, mcp(lom="GrandMean"))
ANOM(test,xlab="")

# variance
mmean <- tapply(elem, lom, mean)[lom]
sqdev <- (elem - mmean)^2
dataf = data.frame(cbind(sqdev,lom))
dataf[,1] = as.numeric(sqdev)
dataf[,2] = as.factor(lom)
library(multcomp)
mod <- lm(sqdev~lom, dataf)
test2 <- glht(mod, mcp(lom="GrandMean"))
ANOM(test2,xlab="")


#### Trend v rozptylech

# Mudholkar-McDermott-Aumont Test
library(lawstat)
mma.test(elem, lom, tail = "both")
# tail = c("right", "left", "both")

# robustni Mudholkar-McDermott-Aumont Test
library(lawstat)
robust.mmm.test(elem, lom, tail = "both")
# tail = c("right", "left", "both")


# Ltrend
library(lawstat)
ltrend.test(elem,lom,score=NULL,location="mean", tail="both", trim.alpha=0.25,bootstrap=FALSE,num.bootstrap=1000,
            correction.method = "none", correlation.method = "pearson")
# location = c("median", "mean", "trim.mean")
# tail = c("right", "left", "both")
# correction.method = c("none", "correction.factor", "zero.removal", "zero.correction")
# correlation.method = c("pearson", "kendall", "spearman")

# Lntrend
library(lawstat)
lnested.test(elem,lom,location="mean", tail="both",trim.alpha=0.25, bootstrap=FALSE, num.bootstrap=1000, correction.method="none", correlation.method="pearson")
# location = c("median", "mean", "trim.mean")
# tail = c("right", "left", "both")
# correction.method = c("none", "correction.factor", "zero.removal", "zero.correction")
# correlation.method = c("pearson", "kendall", "spearman")


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

#### ANOVA homoskedastic

fit <- aov(elem~lom)
summary(fit) # display Type I ANOVA table
summary.lm(fit)
anova(fit)
model.tables(fit,se=T)
model.tables(fit,"means",se=T)
library(car)
Anova(fit, type="II")
### If you use type="III", you need the following line before the analysis
### options(contrasts = c("contr.sum", "contr.poly"))

library(FSA)   
Summarize(elem ~ lom, data=data, digits=3)

library(multcompView)
multcompBoxplot(elem ~ lom, data = data, horizontal = FALSE, compFn = "TukeyHSD",sortFn = "mean", decreasing = TRUE, plotList = list(boxplot = list(fig = c(0, 0.75, 0, 1)),                                                                                                                                multcompTs = list(fig = c(0.7, 0.85, 0, 1)),multcompLetters = list(fig = c(0.87, 0.97, 0.03, 0.98), fontsize = 20,fontface = "bold")))


# Student-Newman-Keuls (SNK) test
library(PMCMRplus)  
out = snkTest(elem,lom, subset=NULL, na.action=NULL,alternative = "two.sided")
summary(out)


# Tukey Honestly Significant Differences test, Tukey - Kramer
library(agricolae)
out <- HSD.test(aov(elem~lom),"lom",alpha = 0.95, group=TRUE,console=TRUE,main="")
plot(out,las=2)

library(PMCMRplus)
out = tukeyTest(elem,lom, subset=NULL, na.action=NULL,alternative = "two.sided")
summary(out)
summaryGroup(out)


## Scheffe test
library(agricolae)
out <- scheffe.test(aov(elem~lom),"lom",alpha = 0.95, group=TRUE,console=TRUE,main="")
plot(out,las=2)

library(PMCMRplus)
out = scheffeTest(elem~lom, subset=NULL, na.action=NULL)
summary(out)
summaryGroup(out)


## Waller test
library(agricolae)
out <- waller.test(aov(elem~lom),"lom", K = 100, group=TRUE, main = NULL,console=TRUE)
plot(out,las=2)


# least significant difference all-pairs comparisons test for normally distributed data with equal group variances.
library(PMCMRplus)
out = lsdTest(elem,lom, subset=NULL, na.action=NULL); out
summary(out)
summaryGroup(out)

library(agricolae)
out <- LSD.test(aov(elem~lom),"lom", alpha=0.95, group=TRUE,main=NULL,console=TRUE,p.adj="bonferroni")
# p.adj=c("none","holm","hommel","hochberg", "bonferroni", "BH", "BY", "fdr")
plot(out,las=2)

# Ryan, Einot and Gabriel and Welsch multiple range test
library(agricolae)
out<- REGW.test(aov(elem~lom),"lom",alpha = 0.95, group=TRUE, main = NULL,console=TRUE)
plot(out,las=2)
summary(out)


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

#### ANOVA heteroskedastic

# trimmed means
library(WRS2)
t1way(elem~lom, data, tr = 0.2, alpha = 0.05, nboot = 500)
vx=lincon(elem~lom, data, tr = 0.2, alpha = 0.05)
cbind(vx$comp[,c(1,2)],round((vx$comp[,6]*100),0))


## Alexander-Govern test.
library(onewaytests)
out = ag.test(elem ~ lom, data = data, alpha = 0.05, na.rm = TRUE, verbose = TRUE)
out
out =paircomp(out, adjust.method = "BH")
# p.adjust.methods: "holm", "hochberg","hommel","bonferroni","BH","BY","fdr","none","tukey","games-howell"


## Welch test

library(car)
mod = aov(elem~lom, data=data, na.action=na.omit)
# mod = lm(elem~lom, data=as.data.frame(cbind(as.numeric(elem),lom)), na.action=na.omit)
Anova(mod, white.adjust=TRUE)
# white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4")

oneway.test(elem~as.factor(lom), subset=NULL, na.action=NULL, var.equal = FALSE)

library(welchADF)
summary(welchADF.test(elem~lom))

library(asbio)
pairw.oneway(elem,as.factor(lom), conf = 0.95, digits = 5, method = "fdr")

library(userfriendlyscience)
oneway(elem,lom, fullDescribe = TRUE, posthoc = 'fdr',corrections=TRUE,plot=TRUE)
# p.adjust.methods: "holm", "hochberg","hommel","bonferroni","BH","BY","fdr","none","tukey","games-howell"

library(coin)
oneway_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores")
oneway_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000))


## Games-Howell
library(userfriendlyscience)
posthocTGH(elem,lom, method= "games-howell", conf.level = 0.95, digits=2, p.adjust="holm",formatPvalue = TRUE)
# method=c("games-howell", "tukey")
# p.adjust= c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")

library(PMCMRplus)
out = gamesHowellTest(elem~lom,alternative = "two.sided",subset=NULL, p.adjust="BH")
summary(out)
summaryGroup(out)


## Tamhane T2
library(PMCMRplus)
out = tamhaneT2Test(elem,lom, subset=NULL, na.action=NULL,welch = TRUE,alternative = "two.sided")
summary(out)
summaryGroup(out)


## Ury Wiggins Hochberg
library(PMCMRplus)
out = uryWigginsHochbergTest(elem, lom, subset=NULL, na.action=NULL,alternative = "two.sided", p.adjust.method ="holm")
# p.adjust.method = p.adjust.methods, ...)
summary(out)
summaryGroup(out)


library(robustbase)
r.mod = lmrob(elem~lom,data=data, na.action=na.omit,weights=NULL)
summary(r.mod)


#### shoda distribuce #########################################################################################

library(PMCMRplus)
adKSampleTest(elem~lom, subset=NULL, na.action=NULL)
out = adAllPairsTest(elem~lom, subset=NULL, na.action=NULL, p.adjust.method = "BH"); out
# p.adjust.methods = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")
summary(out)
summaryGroup(out)
barPlot(out, alpha = 0.05)
plot(out, alpha = 0.05)

library(PMCMRplus) # Murakami k-sample Baumgartner Weiss Schindler Test of Equal Distributions
bwsKSampleTest(elem,lom, subset=NULL, na.action,nperm = 1000)
out <- bwsAllPairsTest(elem, lom, p.adjust="fdr"); out
summary(out)
summaryGroup(out)


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

#### Non-parametric

## Kruskal-Wallis test  # homoscedastic

kruskal.test(elem~lom)

library(agricolae)
kruskal(elem,lom, alpha = 0.05, p.adj="holm",group=TRUE, main = NULL,console=TRUE)
# p.adj=c("none","holm","hommel", "hochberg", "bonferroni", "BH", "BY", "fdr")

library(asbio)
pairw.kw(elem, lom, conf=0.05)
plot(pairw.kw(elem, lom, conf=0.05),las=2)


library(NSM3)
pKW(elem,lom,method="Monte Carlo",n.mc=500) # Kruskal-Wallis
# method=c("Exact", "Monte Carlo", or "Asymptotic")

library(coin)
kruskal_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores") 
kruskal_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000)) 

library(onewaytests)
out = kw.test(elem ~ lom, data = data, alpha = 0.05, na.rm = TRUE, verbose = TRUE)
# paircomp(out, adjust.method = "bonferroni")
# adjust.method = c("bonferroni", "holm", "hochberg", "hommel", "BH","BY", "fdr", "none")

library(PMCMRplus)
kruskalTest(elem, lom, subset=NULL, na.action=NULL, dist = "KruskalWallis")
# dist = c("Chisquare", "KruskalWallis", "FDist")


## Nemenyi
library(PMCMR)
posthoc.kruskal.nemenyi.test(elem,lom, method="Tukey")
# method="Tukey", "Chisq"

library(PMCMRplus)
out = kwAllPairsNemenyiTest(elem, lom, subset=NULL, na.action=NULL,dist = "Tukey"); out
# dist = c("Tukey", "Chisquare")
summary(out)
summaryGroup(out)

library(PMCMRplus)
out=kwAllPairsDunnTest(elem, lom, subset=NULL, na.action=NULL, p.adjust.method ="fdr"); out
summary(out)
summaryGroup(out)

library(PMCMR)
out = posthoc.kruskal.dunn.test(elem,lom,p.adjust.method ="holm"); out
summary(out)


## Conover's non-parametric all-pairs comparison test for Kruskal-type ranked data.
library(conover.test)
conover.test(elem,lom, method="holm", kw=TRUE, label=TRUE,wrap=FALSE, table=TRUE, list=FALSE, rmc=FALSE, alpha=0.95)
# method=c("none", "bonferroni", "sidak", "holm", "hs", "hochberg", "bh", "by")

library(PMCMRplus)
out = kwAllPairsConoverTest(elem, lom, subset=NULL, na.action=NULL,p.adjust.method = "holm"); out
# p.adjust.method = c("single-step", p.adjust.methods); single-step = Tukey's p-adjustment
summary(out)
summaryGroup(out)


library(PMCMR)
posthoc.kruskal.conover.test(elem, lom,p.adjust.method ="holm")


## Dwass, Steel, Critchlow and Fligner
library(PMCMRplus)
out = dscfAllPairsTest(elem, lom, subset=NULL, na.action=NULL, p.adjust.method = "holm", alternative = "two.sided")
summary(out)
summaryGroup(out)


E518 <- read.table("C:/Users/lubop/Documents/kurs/Data1/ascii05/E/E518.txt")
data = stack(E518[-1,])
data = as.data.frame(data); colnames(data)=c("elem","lom")
elem = data[,1]
lom = data[,2]


# Jonckheere-Terpstra (testuje trend v poradi souboru)
library(DescTools)
lomc =  ordered(sapply(lom,function(x) sub("V", "", x)))   
JonckheereTerpstraTest(elem, lomc, alternative = "increasing",nperm = 500)
# alternative = c("two.sided", "increasing", "decreasing")


#### normal scores

## Waerden test
library(agricolae)
out<-waerden.test(elem,lom,alpha=0.95,group=TRUE,console=TRUE,main="")
plot(out,las=2)


# van-der-Waerden normal scores test
library(PMCMRplus)
vanWaerdenTest(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided")
out = vanWaerdenAllPairsTest(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided",p.adjust.method = "fdr")
#  p.adjust.method = c("single-step", p.adjust.methods)
summary(out)
summaryGroup(out)

library(PMCMR)
vanWaerden.test(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided")
out = posthoc.vanWaerden.test(elem,as.factor(lom),p.adjust.method ="holm")
summary(out)

library(coin)
normal_test(as.numeric(elem)~as.factor(lom), ties.method = "mid-ranks", conf.int = FALSE, conf.level = 0.95)
normal_test(as.numeric(elem)~as.factor(lom), ties.method = "mid-ranks", conf.int = FALSE, conf.level = 0.95, distribution = approximate(B = 1000))
# method = c("mid-ranks", "average-scores")


# Lu-Smith normal scores test, transformace dat na normalitu
library(PMCMRplus)
normalScoresTest(elem,lom, subset=NULL, na.action=NULL,alternative = "two.sided")

out = normalScoresAllPairsTest(elem,lom, subset=NULL, na.action=NULL,alternative = "two.sided", p.adjust.method = "holm"); out
# p.adjust.method = c("single-step", p.adjust.methods)
summary(out)
summaryGroup(out)


library(ARTool)
oart = art(elem ~ lom, data = data)
summary(oart)
anova(oart, response = "art",type = "I", factor.contrasts = "contr.sum", test = "F", all.rows = FALSE)
# response = c("art", "aligned")
# type = c("III","II", "I", 3, 2, 1)
# factor.contrasts = "contr.sum"
# test = c("F", "Chisq")
library(ART)
aligned.rank.transform(elem ~ lom, data = data, perform.aov = TRUE, SS.type = "I")
# SS.type = c("III","II", "1I")


## ordinal logistic regression
library(rms)
olr.mod = orm(elem~lom)
olr.mod



## Median test
library(agricolae)
out<-Median.test(elem,lom,correct=TRUE,simulate.p.value = TRUE,alpha=0.95,group=TRUE,console=TRUE,main="")
plot(out,las=2)

library(coin)
median_test(elem~lom, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95)
median_test(elem~lom, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95, distribution = approximate(nresample = 1000))

library(RVAideMemoire)
mood.medtest(elem, lom, exact = TRUE)
pairwise.mood.medtest(elem, lom, exact = TRUE, p.method = "fdr")


### Standard ANOM (Gaussian, homoscedastic)
library(multcomp)
model <- lm(elem ~ lom, data)
hom <- glht(model, mcp(lom="GrandMean"), alternative="two.side")
ANOM(hom,ylab="",xlab="")


### Heteroscedastic ANOM
# With sandwich covariance matrix estimate (Herberich et al. 2010)
library(ANOM)
library(multcomp)
library(sandwich)
model <- lm(elem ~ lom, data)
het1 <- glht(model, mcp(lom="GrandMean"), alternative="two.side", vcov=vcovHC)
ANOM(het1,ylab="",xlab="")


### Nonparametric ANOM
library(ANOM)
ss <- tapply(elem, lom, length)
library(multcomp)
Mat <- contrMat(ss, "GrandMean")
## Using a multivariate t approximation
library(nparcomp)
mult <- mctp(elem~lom, data=data, type="UserDefined", contrast.matrix=Mat, alternative="two.side", info=FALSE, correlation=TRUE, asy.method="mult.t")
mult
ANOM(mult)



##### 2-way ANOVA ####################################################################################################################################

library(openxlsx)
shn = getSheetNames("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx")
shn
sal <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "2wANOVA", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
data = as.data.frame(sal)
nitr = data[,1]
subst = data[,2]
hnoj = data[,3]

### Global test of model assumptions
library(gvlma)
gvmodel <- with(data,gvlma(lm(nitr~subst*hnoj)))
summary(gvmodel)

fit = lm(nitr~subst*hnoj,data=data)
fit = lm(nitr~subst*hnoj-1,data=data)
summary(fit)
anova(fit)

fit2 = aov(nitr~subst*hnoj,data=data)
anova(fit2)
summary(fit2)


library(Hmisc)
with(data,summary(nitr~subst*hnoj))
with(data,summary(nitr~subst*hnoj-1))

library(car)
(my_anova <- aov(nitr~subst*hnoj,data=data))
anova(my_anova)
Anova(my_anova, type = "II") # preferable
# Anova(my_anova, type = "III")

with(data,interaction.plot(subst,hnoj,nitr,ylab="",type ="l"))

summary(aov(nitr~subst*hnoj,data=data))
with(data,interaction.plot(subst,hnoj,nitr,ylab=""))
anova(aov(nitr~subst*hnoj,data=data))

library(ExpDes)
with(data, fat2.crd(subst, hnoj, nitr, quali = c(TRUE, TRUE), mcomp = "tukey",
                    fac.names = c("substrat", "hnojeni"), sigT = 0.05, sigF = 0.05))

library(welchADF)
summary(welchADF.test(lm(nitr~subst*hnoj,data=data)))

library(ARTool)
oart = art(nitr~as.factor(subst)*as.factor(hnoj),data=data)
summary(oart)
anova(oart, response = "art",type = "I", factor.contrasts = "contr.sum", test = "F", all.rows = FALSE)
# response = c("art", "aligned")
# type = c("III","II", "I", 3, 2, 1)
# factor.contrasts = "contr.sum"
# test = c("F", "Chisq")


library(ART)
aligned.rank.transform(nitr~subst*hnoj,data=data,perform.aov = TRUE, SS.type = "I")
# SS.type = c("III","II", "1I")
aligned.rank.transform(nitr~subst*hnoj,data=data,perform.aov = TRUE, SS.type = "II")


library(rankFD)
model <- rankFD(nitr~subst*hnoj,data=data, CI.method = "Normal", effect = "unweighted", hypothesis = "H0F")
model


## ordinal logistic regression
library(rms)
orm(nitr~subst*hnoj,data=data)


library(lmPerm)
# perm "Exact", "Prob", "SPR" will produce permutation probabilities. Anything else, such as "", will produce F-test probabilities.
out = aovp(nitr~subst*hnoj,data=data,perm="Prob",nCycle = 5000); out
summary(out)
anova(lmp(nitr~subst*hnoj,data=data))

library(permuco)
out = aovperm(nitr~subst*hnoj,data=data, np = 5000, method = "freedman_lane"); out
summary(out)

library(robustbase)
model = lmrob(nitr~subst*hnoj-1,data=data) 
summary(model)
plot(fitted(model), residuals(model))
abline(h=0, col=2, lty=2)


library(robustbase)
model = lmrob(nitr~subst*hnoj,data=data) 
summary(model)
plot(fitted(model), residuals(model))
abline(h=0, col=2, lty=2)


##### KORELACE ################################################################################


library(smovie)
hscale = 2
correlation(n = 30, rho = 0, panel_plot = TRUE, hscale = 2, vscale = hscale, delta_n = 1, delta_rho = 0.1)


## Anscombovo kvarteto
library(asbio)
data(anscombe)
cor(anscombe[,"x1"],anscombe[,"y1"])
cor(anscombe[,"x2"],anscombe[,"y2"])
cor(anscombe[,"x3"],anscombe[,"y3"])
cor(anscombe[,"x4"],anscombe[,"y4"])
par(mfrow=c(2,2))
plot(anscombe[,"x1"],anscombe[,"y1"])
plot(anscombe[,"x2"],anscombe[,"y2"])
plot(anscombe[,"x3"],anscombe[,"y3"])
plot(anscombe[,"x4"],anscombe[,"y4"])
par(mfrow=c(1,1))


library(datasauRus)
# Anscombe
par(mfrow = c(5, 3), mar=c(1,3,3,1))
nms <- names(datasaurus_dozen_wide)
for (i in seq(1, 25, by = 2)){
  nm <- substr(nms[i], 1, nchar(nms[i]) - 2)
  plot(datasaurus_dozen_wide[[nms[i]]],
       datasaurus_dozen_wide[[nms[i+1]]],
       xlab = "", ylab = "", main = nm)
}
# Simpson
par(mfrow = c(1, 2))
plot(simpsons_paradox_wide[["simpson_1_x"]],
     simpsons_paradox_wide[["simpson_1_y"]],
     xlab = "x", ylab = "y", main = "Simpson's Paradox 1")
plot(simpsons_paradox_wide[["simpson_2_x"]],
     simpsons_paradox_wide[["simpson_2_y"]],
     xlab = "x", ylab = "y", main = "Simpson's Paradox 2")
par(mfrow = c(1, 1))



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

library(ACSWR)

data(viscos)
data=viscos
x = data[,1]
y = data[,2]

plot(x,y)

# kovariance
cov(x,y)

# korelace
cov(x,y)/(sd(x)*sd(y))
cor(x,y)


library(openxlsx)
pisky <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "Pisky", startRow = 1, colNames = TRUE, rowNames = TRUE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
summary(pisky)
pisky = pisky[,-1]
colnames(pisky)
x =  pisky[,"Sr"]
y =  pisky[,"Ba"]


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


cor(x,y, method="pearson")
cor(x,y, method="spearman")
cor(x,y, method="kendall")


library(StatDA)
cor.sign(cbind(x,y), method="pearson")
cor.sign(cbind(x,y), method="kendall")
cor.sign(cbind(x,y), method="spearman")

library(Hmisc)
rs<-rcorr(x, y, type="spearman"); rs
# type=c("pearson","spearman")


library(StatDA)
RobCor.plot(x, y, quan=1/2, alpha=0.025, colC=3, colR=4)

library(robustbase)
mr<-covMcd(cbind(x,y), alpha=0.975, cor=TRUE)$cor; mr

library(fmsb)
spearman.ci.sas(x, y, adj.bias=TRUE, conf.level=0.95)


plot(x,y,log="")

library(aplpack)
bagplot(x,y,show.baghull=TRUE,show.loophull=TRUE,precision=1,dkmethod=2)


### test vyznamnosti
cor.test(x, y, alternative = "two.sided", method = "pearson", exact = NULL, conf.level = 0.95, continuity = FALSE)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "kendall", "spearman")



library(UsingR)
scatter.with.hist(x, y, hist.col = gray(0.95),trend.line = "lm")
simple.scatterplot(x, y)



# rank scatter plot
xr <-rank(x)
yr <-rank(y)
plot(xr,yr, col=1, xlab=" x", ylab="y",cex=1.2)
# text(xr,yr, labels=rownames(mb))
abline(0,1,lty=2)


library(asbio)
chi.plot(x,y)


library(lcopula)
K.plot(cbind(x,y), add = F)

set.seed(123)
data=cbind(x.y)
library(boot)
b3 <- boot(data, statistic = function(data, i) {
         cor(data[i, "x"], data[i, "y"], method='pearson')},  R = 1000)
b3

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

library(openair)
corPlot(pisky, pollutants = NULL, type = "default",cluster = TRUE, cols = "hue", r.thresh = 0.8,text.col = c("black", "black"), auto.text = FALSE)
cols =c("default","increment","heat","spectral","hue", "brewer1","greyscale") 


m = cor(pisky, method="pearson")

library(ellipse)
plotcorr(m)


library(ellipse)
colors <- c("#A50F15","#DE2D26","#FB6A4A","#FCAE91","#FEE5D9","white",
            "#EFF3FF","#BDD7E7","#6BAED6","#3182BD","#08519C")
plotcorr(m, col=colors[5*m + 6])


library(corrgram)
corrgram(m, order=TRUE)


library(StatDA)
mp<-cor(pisky, method="pearson")
ms<-cor(pisky, method="spearman")
CorCompare(mp, ms, labels1="pearson", labels2="spearman", ndigits=3, cex.cor=0.8, cex.label=0.8)



# TESTY KORELACE

library(fBasics)
correlationTest(x, y, method = "pearson", title = NULL, description = NULL)
# method = c("pearson", "kendall", "spearman")
pearsonTest(x, y, title = NULL, description = NULL)
kendallTest(x, y, title = NULL, description = NULL)
spearmanTest(x, y, title = NULL, description = NULL)


library(psych)
corr.test(x, y, use = "pairwise",method="pearson",adjust="holm")
# use= c("pairwise", "complete")
# method= c("pearson", "spearman", "kendall")
# adjustment for multiple tests: adjust = c("holm", "hochberg", "hommel","bonferroni", "BH", "BY", "fdr", "none")


library(psych)
# correlation matrix is an identity matrix
cortest.bartlett(m, n = NULL)
# Steiger test
cortest.normal(m, R2 = NULL, n1 = NULL, n2 = NULL, fisher = TRUE)
cortest(pisky,R2=NULL,n1=NULL,n2 = NULL, fisher = TRUE,cor=TRUE) 
# Jennrich test
cortest.jennrich(m,R2,n1=NULL, n2=NULL)
# alternative test
cortest.mat(m,R2=NULL,n1=NULL,n2 = NULL)


rb = corCi(pisky, keys = NULL, n.iter = 500,  p = 0.05,overlap = FALSE,poly = FALSE, method = "pearson", plot=F,minlength=5,n=NULL)
rb$rho
rb$ci
corCi(pisky, keys = NULL, n.iter = 500,  p = 0.05,overlap = FALSE,poly = FALSE, method = "pearson", plot=T,minlength=5,n=NULL)



### Dunnettuv test ######################################################################


data <- data.frame(value = c(76, 77, 77, 81, 82, 82, 83, 84, 85, 89,
                             81, 82, 83, 83, 83, 84, 87, 90, 92, 93,
                             77, 78, 79, 88, 89, 90, 91, 95, 95, 98),
                   Group = rep(c("control", "V1", "V2"), each = 10))
data$Group<-as.factor(data$Group)
head(data)

boxplot(value~Group,data = data,main = "",xlab = "Groups",ylab = "Value",col = "gray",border = "black")

model <- aov(value ~ Group, data = data)
summary(model)

library(DescTools)
DunnettTest(x=data$value, g=data$Group)

library(PMCMRplus)
dunnettTest(x=data$value, g=data$Group, alternative = "two.sided")
# alternative = c("two.sided", "greater", "less")


boxplot(Ozone ~ Month, data = airquality)
DunnettTest(Ozone ~ Month, data = airquality)
DunnettTest(Ozone ~ Month, data = airquality, control="8", conf.level=0.9)


## Mucociliary efficiency from the rate of removal of dust in normal
##  subjects, subjects with obstructive airway disease, and subjects
##  with asbestosis.
reference <- c(2.9, 3.0, 2.5, 2.6, 3.2) # normal subjects
OAD <- c(3.8, 2.7, 4.0, 2.4)      # with obstructive airway disease
asbestosis <- c(2.8, 3.4, 3.7, 2.2, 2.0) # with asbestosis
data2 = list(reference, OAD, asbestosis)
names(data2) = c("normal", "OAD", "asbestosis")
boxplot(data2)
DunnettTest(data2)


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


### Analýza cenzorovaných dat 

library(NADA)
v = c('<1', 1, '<1', 1, 2)
splitQual(v)

obs = c(0.5, 0.5, 1.0, 1.5, 5.0, 10, 100)
censored = c(TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE)

pctCen(obs, censored)

# Cadmium concentrations in fish for two regions of the Rocky Mountains.
data(Cadmium)
head(Cadmium)
obs = Cadmium$Cd
censored = Cadmium$CdCen
groups = Cadmium$Region

# Cd differences between regions?
cendiff(obs, censored, groups)
cenfit(Cen(obs, censored)~groups)
cenboxplot(obs, censored, groups, log=FALSE, range=0)
cenboxplot(obs, censored, groups, log=TRUE, range=0)
library(NADA2)
cboxplot(obs, censored, groups)

#SRKYMT
censtats(obs[1:9],censored[1:9])

#COLOPLT
censtats(obs[10:19],censored[10:19])


# Hodnoty 17-alfa-hydroxyprogesteronu v krvi novorozenců
library(openxlsx)
ohpx <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\OHP.xlsx",
                  sheet = "OHP2", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
ohp2 <- rep(0, length(ohpx[,1]))
ohp2[which(ohpx[,1]=="<DL")]<- 1
ohp1 <- ohpx[,1]
ohp1[which(ohpx[,1]=="<DL")]<- 30
ohp1  <- as.numeric(ohp1)
ohp  <- cbind(ohp1,ohp2); ohp
colnames(ohp) = c("OHP","CENZ")
head(ohp)

as.logical(ohp[,"CENZ"])
pctCen(obs=ohp[,"OHP"], censored=as.logical(ohp[,"CENZ"]))
cenfit(Cen(obs=ohp[,"OHP"], censored=as.logical(ohp[,"CENZ"]))~1)


# podil hodnot pod DL (DL = 30)
t<-table(ohp[,"CENZ"])
names(t)<-c("necenzorováno","cenzorováno")
100*t/sum(t) # procenta cenzorovaných

# prumer vypocteny pouze z necenzorovanych hodnot
mean(ohp[ohp[,"CENZ"]==0,"OHP"]) 
median(ohp[ohp[,"CENZ"]==0,"OHP"]) 

# prumer vypocteny ze vsech hodnot pokud misto cenzorovani nechame detekcni limit
mean(ohp[,"OHP"])
median(ohp[,"OHP"]) 

# prumer vypocteny ze vsech hodnot pokud cenzorovane hodnoty nahradime nulou
xp<-ohp
xp[xp[,"CENZ"]==1,"OHP"]<-0
mean(xp[,"OHP"])
median(xp[,"OHP"]) 

# prumer vypocteny ze vsech hodnot pomoci survival analysis
library(survival)
# za predpokladu normalniho rozlozeni.
m<-survreg(formula=Surv(ohp[,"OHP"], ohp[,"CENZ"], type="left")~1, dist="gaussian")   
summary(m)
m[[1]]
# za predpokladu logaritmicko normalniho rozdeleni
m<- survreg(Surv(ohp[,"OHP"],ohp[,"CENZ"],type="left")~1,dist='loggaussian')
summary(m)
exp(m[[1]])

boxplot(ohp[,"OHP"])
hist(ohp[,"OHP"],freq = F)

cenboxplot(ohp[,"OHP"], as.logical(ohp[,"CENZ"]), groups=NULL, log=FALSE, range=0)

censtats(ohp[,"OHP"], as.logical(ohp[,"CENZ"]))


library(NADA2)
cboxplot(ohp[,"OHP"], as.logical(ohp[,"CENZ"]))

# Dissolved iron concentrations over several years in the Brazos River, Texas. Summer concentrationswere used.
# Objective is to determine if there is a trend over time.  Iron has two detection limits, at 3 and 10ug/L. 
data(DFe)
DFe
cenxyplot(x=DFe$Year,xcen=DFe$YearCen,y=DFe$Summer,ycen=DFe$SummerCen, log="", lty="dashed")
with(DFe, cenxyplot(Year, YearCen, Summer, SummerCen))

with(DFe, censummary(Summer, SummerCen))


data(CuZn)
with(CuZn, cenxyplot(Cu, CuCen, Zn, ZnCen))


# Lead concentrations in the blood and several organs of herons in Virginia.
data(PbHeron)
ppw.test(PbHeron$Liver,PbHeron$LiverCen,PbHeron$Bone,PbHeron$BoneCen)

# with groups
with(PbHeron, cen_ecdf(Liver, LiverCen, DosageGroup))
# all data
with(PbHeron, cen_ecdf(Liver, LiverCen))

# test of differences in means between two groups
cen2means(PbHeron$Liver,PbHeron$LiverCen,PbHeron$DosageGroup, LOG = TRUE, printstat = TRUE)
cenperm2(PbHeron$Liver,PbHeron$LiverCen,PbHeron$DosageGroup,alternative="t",R = 9999, alternative = "two.sided", printstat = TRUE)
cen1way(PbHeron$Liver,PbHeron$LiverCen,PbHeron$DosageGroup)

# test of differences in means between groups
cen1way(PbHeron$Liver,PbHeron$LiverCen,PbHeron$Group, mcomp.method = "BH", printstat = TRUE)
cenanova(PbHeron$Liver,PbHeron$LiverCen,PbHeron$Group, LOG = FALSE, printstat = TRUE)
cenpermanova(PbHeron$Liver,PbHeron$LiverCen,PbHeron$Group, R = 9999, printstat = TRUE)


# Atrazine concentrations in a series of Nebraska wells before (June) and after (September) the grow-ing season.
data(atrazine)
# Comparing standard/guideline value
cen_paired(atrazine$June, atrazine$JuneCen, 0.01, alternative = "greater")

cen_paired(atrazine$June,atrazine$JuneCen,atrazine$Sept,atrazine$SeptCen, alternative = "two.sided", printstat = TRUE)

cen_signtest(atrazine$June,atrazine$JuneCen,atrazine$Sept,atrazine$SeptCen, alternative = "two.sided", printstat = TRUE)
cen_signedranktest(atrazine$June,atrazine$JuneCen,atrazine$Sept,atrazine$SeptCen, alternative = "two.sided", printstat = TRUE)


# Mercury concentrations in fish across the United States.
data(Brumbaugh)
cenCompareQQ(Brumbaugh$Hg,Brumbaugh$HgCen)
cenCompareCdfs(Brumbaugh$Hg,Brumbaugh$HgCen)


library(EnvStats)
boxcoxCensored(x=ohp[,1], censored=ohp[,2], censoring.side = "left",lambda = seq(-2, 2, by = 0.5), optimize = FALSE,
               objective.name = "PPCC", eps = .Machine$double.eps,include.x.and.censored = TRUE, prob.method = "michael-schucany", plot.pos.con = 0.375)

qqPlotCensored(x=ohp[,1], censored=ohp[,2])

distChooseCensored(x=ohp[,1], censored=ohp[,2], censoring.side = "left", alpha = 0.05,
                   method = "sf", choices = c("norm", "gamma", "lnorm"),
                   est.arg.list = NULL, prob.method = "hirsch-stedinger",
                   plot.pos.con = 0.375, warn = TRUE, keep.data = TRUE,
                   data.name = NULL, censoring.name = NULL)

gofTestCensored(x=ohp[,1], censored=ohp[,2], test = "sf")


###### REGRESE ##################################################################################################

library(animation)
least.squares(x, y, n = 15, ani.type = c("slope", "intercept"), a, b, a.range, b.range,
              ab.col = c("gray", "black"), est.pch = 19, v.col = "red", v.lty = 2, rss.pch = 19,
              rss.type = "o", mfrow = c(1, 2))

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

library(investr)
data(arsenic)

library(ACSWR)
data(viscos)
data=viscos
x = data[,1]
y = data[,2]

library(openxlsx)
sal <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "regrese2", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
data = as.data.frame(sal)
x = data[,1]
y = data[,2]
z = data[,3]

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

plot(x,y)

lmMod <- lm(y~x, data=data, weights=NULL, na.action = na.pass)
summary(lmMod)

plot(x,y)
abline(lm(y~x, data=data), col=2, lty=2)
abline(lm(y~x-1, data=data), col=4, lty=2)

cp = predict(lmMod, newdata = 0, interval = 'confidence')
lines(x,cp[,2],col=6,lty=3)
lines(x,cp[,3],col=6,lty=3)

pp = predict(lmMod, newdata = 0, interval = 'prediction')
lines(x,pp[,2],col=3,lty=3)
lines(x,pp[,3],col=3,lty=3)

confint(lmMod,level = 0.95)
confint(lmMod,level = 0.99)

lm.out <- lm(y ~ x)
newx = seq(min(x),max(x),by = 0.05)

conf_interval <- predict(lm.out, newdata=data.frame(x=newx), interval="confidence", level = 0.95)
plot(x, y, xlab="x", ylab="y", main=" ")
abline(lm.out, col="lightblue")
lines(newx, conf_interval[,2], col="blue", lty=2)
lines(newx, conf_interval[,3], col="blue", lty=2)

pred_interval <- predict(lm.out, newdata=data.frame(x=newx), interval="prediction", level = 0.95)
lines(newx, pred_interval[,2], col="green", lty=2)
lines(newx, pred_interval[,3], col="green", lty=2)


# test koeficientu
library(lmtest)
coeftest(lmMod, vcov. = NULL, df = NULL)
coefci(lmMod, parm = NULL, level = 0.95, vcov. = NULL, df = NULL)

library(car)
confidenceEllipse(lmMod)


# fitted
fitted.values(lmMod)

library(visreg)
visreg(lmMod, "x")


# residuals
resid(lmMod) # klasicka
plot(fitted.values(lmMod), resid(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), resid(lmMod)), col = 4)

rstandard(lmMod) # standardizovana
plot(fitted.values(lmMod), rstandard(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), rstandard(lmMod)), col = 4)

rstudent(lmMod) # studentizovana
plot(fitted.values(lmMod), rstudent(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), rstudent(lmMod)), col = 4)

library(s20x)
ciReg(lmMod, conf.level = 0.95, print.out = TRUE)

cooks20x(lmMod, main = "Cook's Distance plot", xlab = "observation number",ylab = "Cook's distance", line = c(0.5, 1.2, 2), cex.labels = 1,axisOpts = list(xAxis = TRUE, yAxisTight = FALSE))

eovcheck(lmMod, xlab = "Fitted values",ylab = "Residuals", col = NULL, smoother = FALSE, twosd = FALSE,levene = FALSE)

modcheck(lmMod, plotOrder = 1:4, args = list(eovcheck = list(smoother = FALSE, twosd = FALSE, levene = FALSE), normcheck = list(xlab = c("Theoretical Quantiles", ""), ylab = c("Sample Quantiles",""), main = c("", ""), col = "light blue", bootstrap = FALSE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = FALSE, whichPlot = 1:2, usePar =   TRUE), cooks20x = list(main = "Cook's Distance plot", xlab = "observation number", ylab = "Cook's distance", line = c(0.5, 0.1, 2), cex.labels = 1, axisOpts = list(xAxis = TRUE))), parVals = list(mfrow = c(2, 2), xaxs = "r", yaxs = "r", pty = "s", mai= c(0.2, 0.2, 0.05, 0.05)))

normcheck(lmMod, xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = FALSE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = FALSE, whichPlot = 1:2, usePar = TRUE)

residPlot(lmMod, f = 1)


library(car)
residualPlots(lmMod, terms = ~., layout = NULL, main = "", fitted = TRUE, AsIs=TRUE, plot = TRUE,tests = TRUE)
influencePlot(lmMod)
infIndexPlot(lmMod)


library(UsingR)
simple.lm(x, y, show.residuals=TRUE, show.ci=TRUE, conf.level=0.95,pred=NULL) 

par(mfrow=c(2,2))
plot(lmMod)
par(mfrow=c(1,1))

library(lindia)
gg_boxcox(lmMod, showlambda = TRUE, lambdaSF = 3, scale.factor = 1)
gg_cooksd(lmMod, label = TRUE, show.threshold = TRUE,threshold = "convention", scale.factor = 1)
gg_diagnose(lmMod, theme = NULL, ncol = NA, plot.all = TRUE,scale.factor = 1, boxcox = FALSE, max.per.page = NA)
gg_qqplot(lmMod, scale.factor = 1)
gg_resfitted(lmMod, scale.factor = 1)
gg_reshist(lmMod, bins = 20)
gg_resleverage(lmMod, method = "loess", se = FALSE, scale.factor = 1)
gg_resX(lmMod, plot.all = TRUE, scale.factor = 1, max.per.page = NA, ncol = NA)
gg_scalelocation(lmMod, method = "loess", scale.factor = 1, se = FALSE)


# Bonferroni Outlier Test
library(car)
ou = outlierTest(lm(y~x), cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))

# Lund Outlier Test
lundout<-function(x1,x2,a) {
  testoutlm<-lm(x1~x2)
  standardresid<-rstandard(testoutlm)
  nl<-length(x1)
  q<-length(testoutlm$coefficients)
  F<-qf(c(1-(a/nl)),df1=1,df2=nl-q-1,lower.tail=TRUE)
  crit<-((nl-q)*F/(nl-q-1+F))^0.5
  oul = which(abs(standardresid)>crit)
  return(oul)
}
lundout(y,x,0.05)

# Shapiro-Wilk test rezidui
shapiro.test(resid(lmMod))
shapiro.test(rstandard(lmMod))
shapiro.test(rstudent(lmMod))


# Breusch-Pagan test na heteroskedasticitu
library(lmtest)
lmtest::bptest(y~x, varformula = NULL, studentize = TRUE, data = data) 

# Harrison-McCabe test for heteroskedasticity
library(lmtest)
lmtest::hmctest(y~x, point = 0.5, order.by = NULL, simulate.p = TRUE, nsim = 1000, plot = TRUE, data = data)

# Score Test for Non-Constant Error Variance
library(car)
car::ncvTest(lmMod)

# Goldfeld-Quandt test against heteroskedasticity
library(lmtest)
lmtest::gqtest(y~x, point = 0.5, fraction = 0, alternative = "two.sided",order.by = NULL, data = data)
# alternative = c("greater", "two.sided", "less")


# Durbin-Watson test for autocorrelation (zavislost rezidui na x)
library(lmtest)
lmtest::dwtest(y~x, order.by = NULL, alternative = "two.sided", iterations = 15, exact = NULL, tol = 1e-10, data = data)
# alternative = c("greater", "two.sided", "less")
library(car)
durbinWatsonTest(lmMod)

### Testy linearity
# Rainbow test linearity
library(lmtest)
lmtest::raintest(y~x, fraction = 0.5, order.by = NULL, center = NULL, data=data)

# Harvey-Collier test for linearity
library(lmtest)
lmtest::harvtest(y~x, order.by = NULL, data = data)



###### regrese pocatkem ############################

s.lm <- lm(y ~ x, data = data)
summary(s.lm)

s.lm2 <- lm(y ~ 0 + x, data = data) # Adding the 0 term tells the lm() to fit the line through the origin
summary(s.lm2)

s.lm3 <- lm(y ~ x-1, data = data) # Adding the 0 term tells the lm() to fit the line through the origin
summary(s.lm3)

plot(x,y) # regrese bez useku
abline(lm(y ~ x-1, data = data), col=2, lty=2)


### Box-Cox #############################################################
library(AID)
lmbc = boxcoxlm(as.matrix(x), y, method = "lse", lambda = seq(-3,3,0.01), lambda2 = NULL, plot = TRUE,alpha = 0.05, verbose = TRUE)
lmbc$lambda.hat
caret::BoxCoxTrans(y)
library(forecast)
ynew = BoxCox(y, lambda=2)
lmMod_bc <- lm(ynew ~ x, data=data)
bptest(lmMod_bc)
plot(lmMod_bc)


### regrese polynomem
s.lmq <- lm(y ~ x + poly(x,2), data = data)
summary(s.lmq)

plot(x,y) # regrese polynomem
lines(lm(y ~ poly(x^2), data = data), col=2, lty=2)
abline(lm(y ~ x + poly(x,2), data = data), col=4, lty=2)

lmMod2 = lm(z ~ x, data = data)

### Srovnani 2 modelu
anova(lmMod,lmMod2)

library(performance)
compare_performance(lmMod, lmMod2, rank = TRUE)

# Srovnani 2 primek
library(gap)
chow.test(y1=y,x1=x,y2=z,x2=x,x=NULL)


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

### Resistant Regression
library(MASS)
lmM1 = lqs(y ~ x, data,method = "lts")
lmM1
# method = c("lts", "lqs", "lms", "S", "model.frame")

lmM2 = rlm(y ~ x, data, weights=NULL, method = "M", wt.method = "inv.var")
lmM2 
# method = c("M", "MM", "model.frame"),
# wt.method = c("inv.var", "case")

library(LearnEDA) # Fits Tukey resistant line of form a + b (x - xC).
lmM3 = rline(x,y,iter=1)


### truncated/censored regression 

library(NADA)
# TCE concentrations (ug/L) in ground waters of Long Island, New York, along with several possible explanatory variables.
data(TCEReg)
# Using the formula interface
cre=with(TCEReg, cenreg(Cen(TCEConc, TCECen)~PopDensity,dist="lognormal"))
# "extreme", "logistic", "gaussian", "weibull", "exponential", "rayleigh", "loggaussian", "lognormal", "loglogistic", "t"
summary(cre)
coef=as.numeric(cre@survreg$coef)
plot(TCEReg$PopDensity,TCEReg$TCEConc)
abline(coef[1],coef[2])
plot(TCEReg$PopDensity,TCEReg$TCEConc,log="y")


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

library(chemCal)
data(din32645)
din32645
data(massart97ex1)
massart97ex1
data(massart97ex3)
massart97ex3
data(utstats14)


m <- lm(y ~ x, data = din32645)
calplot(m,xlim = c("auto", "auto"), ylim = c("auto", "auto"),xlab = "Concentration", ylab = "Response", alpha=0.05, varfunc = NULL)
plot(m, which=1)
plot(m, which=2)
plot(m, which=3)


## Prediction of x with confidence interval
prediction <- inverse.predict(m, 3500, alpha = 0.05)

## Critical value (decision limit)
crit <- lod(m, alpha = 0.05, beta = 0.05)

m <- lm(y ~ x, data = din32645)
lod(m)
lod(m, alpha = 0.05, beta = 0.05)

m <- lm(y ~ x, data = din32645)
loq(m)
loq(m, n = 3) # better by using replicate measurements



library(EnvStats)
calibrate(y ~ x, data=din32645, test.higher.orders = TRUE, max.order = 4, p.crit = 0.05, F.test = "partial", method = "qr", model = TRUE, x = FALSE, y = FALSE, contrasts = NULL, warn = TRUE)


### nonlinear regression ##################################################################################################

library(investr)
data(Puromycin, package = "datasets")
Puromycin2 <- Puromycin[Puromycin$state == "treated", ][, 1:2]


detach(openxlsx)
write.xlsx(Puromycin, "c:\\Users\\lubop\\Documents\\kurs\\DVpiskyXX.xlsx", col.names=TRUE, row.names=TRUE)

write.table(Puromycin, file = "c:\\Users\\lubop\\Documents\\kurs\\DVpiskyXX.txt", append = FALSE, quote = TRUE, sep = " ",
            eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)



Puro.nls <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin2,start = c(Vm = 200, K = 0.05))
plotFit(Puro.nls, interval = "both", pch = 19, shade = TRUE, col.conf = "skyblue4", col.pred = "lightskyblue2")


################################################################################################################################3

library(BlandAltmanLeh)
bland.altman.PEFR  # data

library(deming)
data(arsenate)
data(ferritin)

# Deming
fit <- deming(aes ~ aas, data=arsenate, xstd=se.aas, ystd=se.aes, conf=.95, jackknife=TRUE, dfbeta=FALSE,x=FALSE, y=FALSE, model=TRUE)
print(fit)
fit$coefficients
fit$coefficients[1]
fit$coefficients[2]
fit$residuals
mean(fit$residuals)

# York
library(bfsl)
fit1 <- bfsl(x=arsenate$aas, y = arsenate$aes, sd_x = arsenate$se.aas, sd_y = arsenate$se.aes, r = 0, control = bfsl_control())
fit1

# Passing Bablok
afit1 <- pbreg(aes ~ aas, data= arsenate,conf=.95,nboot = 0, method=1, eps=sqrt(.Machine$double.eps), x = FALSE, y = FALSE, model = TRUE)
print(afit1)
# Theil Sen
afit2 <- theilsen(aes ~ aas, data= arsenate,conf=.95, nboot = 0, symmetric=FALSE, eps=sqrt(.Machine$double.eps), x = FALSE, y = FALSE, model = TRUE)
print(afit2)

library(mcr)
data(creatinine)
fit.lr <- mcreg(as.matrix(creatinine), method.reg="LinReg", na.rm=TRUE)
fit.wlr <- mcreg(as.matrix(creatinine), method.reg="WLinReg", na.rm=TRUE)
compareFit( fit.lr, fit.wlr )


#### Bland-Altman plot

library(BlandAltmanLeh)
bland.altman.plot(group1=arsenate$aes, group2=arsenate$aas, two = 1.96, mode = 1, graph.sys = "base", conf.int = 0, silent = TRUE, sunflower = FALSE, geom_count = FALSE)
bland.altman.stats(group1=arsenate$aes, group2=arsenate$aas, two = 1.96, mode = 1, conf.int = 0.95)

library(blandr)
blandr.statistics (arsenate$aes, arsenate$aas, sig.level=0.95 )
blandr.draw(arsenate$aes, arsenate$aas)
blandr.draw(arsenate$aes, arsenate$aas,ciDisplay = FALSE, ciShading = FALSE)


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