###################################################################################################################
######### Testy shody ###############################################################################################
##########################################################################################################

b314<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B314.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b314)
x = unlist(b314)


#### t-test ############################################################################################
# one sample t-test: t.test(x,mu=3, alternative = "two.sided") 

t.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.one <- pwr.t.test(d=(abs(mean(x)-2.4)), power=0.9, type="one.sample", alternative="two.sided")
plot(p.t.one)


#### medianovy test ############################################################################################

library(BSDA)
SIGN.test(x, y = NULL, md = 2.4, alternative = "two.sided", conf.level = 0.95)


#### Wilcoxonuv (Mann-Whitney) test ############################################################################################

wilcox.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, exact = NULL, correct = TRUE, conf.int = TRUE, conf.level = 0.95) 

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


##### Dvouvyberove testy  ###########################################################################################################################################


###### Scale #######################################################################################

### F test
var.test(x, y, ratio = 1, alternative = "two.sided", conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

### Ansari Bradley two sample test
ansari.test(x, y, alternative = "two.sided", exact = FALSE, conf.int = FALSE, conf.level = 0.95)


###### Location ###################################################################################

b325ab<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B325ab.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b325ab)
X = as.data.frame(b325ab)
# x = unlist(b325ab)
x= unlist(X[,1])
y= unlist(X[,2])


#### t-test ############################################################################################

# t.test(x,mu=3, alternative = "two.sided") 

sd(x)
sd(y)
# srovnani pomoci F-testu nebo jineho testu shody rozptylu

boxplot(x,y)

t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95)
t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) # Welchuv test
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.two <- pwr.t.test(d=(abs(mean(x)-mean(y))), power=0.9, type="two.sample", alternative="two.sided")
plot(p.t.two)
plot(p.t.two, xlab="sample size per group")

p.t.two2 <- pwr.t2n.test(n1 = length(x), n2= length(y), d = (abs(mean(x)-mean(y))), sig.level = 0.05, power = NULL,
                         alternative ="two.sided")
# alternative = c("two.sided","less","greater")
plot(p.t.two2)
plot(p.t.two2, xlab="sample size per group")


#### 2 sample medianovy test ############################################################################################

# one sample sign test 
library(BSDA)
SIGN.test(x, y, md = 0, alternative = "two.sided", conf.level = 0.95)


#### 2 sample Wilcoxonuv (Mann-Whitney) test ############################################################################################


wilcox.test(x, y,alternative = "two.sided", mu = 0, paired = FALSE, exact = FALSE, correct = TRUE,conf.int = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(exactRankTests)
wilcox.exact(x, y, mu = 0,alternative = "two.sided",paired = FALSE, exact = FALSE,conf.int = TRUE, conf.level = 0.95)

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

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

### Srovnani 2 distribuci ##########################################################################################


### Kolmogorov Smirnov test
ks.test(x,y, alternative = "two.sided", exact = FALSE)
# alternative = c("two.sided", "less", "greater")

library(KScorrect)
ks_test_stat(x, y)

library(dgof)
ks.test(x, y,alternative = "two.sided", exact = NULL, tol=1e-8, simulate.p.value=TRUE, B=2000)
# alternative = c("two.sided", "less", "greater")

library(twosamples)
ks_test(x, y, nboots = 2000, p = 0.025)


### Cramer von Mises
library(CDFt)
CramerVonMisesTwoSamples(x, y)

library(cramer)
cramer.test(x,y,conf.level=0.95,replicates=1000,sim="ordinary",just.statistic=FALSE,kernel="phiCramer", maxM=2^14, K=160)



lxy = list(x,y); names(lxy)= c("x","y")
data = stack(lxy)
data = as.data.frame(data); colnames(data)=c("xy","z")
xy = data[,1]
z = data[,2]
stripchart(xy~z, method="jitter",pch=1,col=1,cex=1.5) 
boxplot(xy~z, pch=1,col=5,cex=1.5) 


### Anderson Darling
library(kSamples)
ad.test(xy~z, data = data, method = "simulated", dist = FALSE, Nsim = 10000)
# method = c("asymptotic", "simulated", "exact")


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


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


### Distribucni grafy

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(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(elem ~ lom, data = data, col = "bisque", ylim = "")



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



# 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)


### outliers
stresid<-rstandard(lm(elem~lom))
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(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(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")


#### 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")))


## 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)

# 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)


#### ANOVA heteroskedastic


## 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 = "fdr")
# 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)


#### 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="holm"); out
summary(out)
summaryGroup(out)


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

#### Non-parametic

## 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

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

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)



## 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")



##### 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(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)


############################################################################################################
###### Regrese #####################################################################
############################################################################################################

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)


### regrese s LOD
library(lodr)
x2 = 
  fitl <- lod_lm(data=lod_data, frmla=y~x, lod=0.1,var_LOD="x2", nSamples=100, boots=500)
summary(fitl)
coef(fit)
residuals(fitl)


library(truncreg)
truncreg(y~x, data=data, point = 2, direction = "left",model = TRUE, y = FALSE, x = FALSE, scaled = FALSE)


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

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(bfsl)
library(MethComp)



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)


# 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)

data(ferritin)
temp <- ferritin[ferritin$period <4,]
plot(temp$old.lot, temp$new.lot, type='n', log='xy', xlab="Old lot", ylab="New Lot")
text(temp$old.lot, temp$new.lot, temp$period,col=temp$period)


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.PEFR  # data
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(MethComp)
library(blandr)


# Linuv koeficient konkordance
library(DescTools)
CCC(arsenate$aes, arsenate$aas, ci = "z-transform", conf.level = 0.95, na.rm = FALSE)
library(epiR)
epi.ccc(arsenate$aes, arsenate$aas, ci = "z-transform", conf.level = 0.95, rep.measure = FALSE, subjectid=NULL)
library(agRee)
ratings = cbind(arsenate$aes, arsenate$aas)
agree.ccc(ratings, conf.level=0.95,method="bootstrap",nboot=999, nmcmc=10000,mvt.para=list(prior=list(lower.v=4, upper.v=25,Mu0=rep(0, ncol(ratings)),Sigma0=diag(10000, ncol(ratings)),p=ncol(ratings),V=diag(1, ncol(ratings))),initial=list(v=NULL, Sigma=NULL)),NAaction="omit")
# method=c("jackknifeZ", "jackknife","bootstrap","bootstrapBC","mvn.jeffreys", "mvn.conjugate","mvt", "lognormalNormal", "mvsn", "mvst")
# NAaction=c("fail", "omit")


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