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


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


library(EBImage) # Bioc
Eg1 <- readImage("c:\\Users\\lubop\\Documents\\image\\Egypt1.jpg")
display(Eg1)
print(Eg1)
img <- Eg1

# Brightness
img1 = img+0.4
img1 = img-0.4
display(img1)

# Contrast
img2 = 2*img
img2 = .5*img
display(img2)

# Gamma correction
img3 = img^3
img3 = img^.7
img3 = (0.2+img)^3
display(img3)

# Treshold
img4 = img>0.5
display(img4)

imgc=img
colorMode(imgc) = Grayscale
display(imgc)

colorMode(imgc) = Color
display(imgc)

imgb = rgbImage(red=img, green=img*0, blue=img*0)
display(imgb)
imgb = rgbImage(red=img*0, green=img, blue=img*0)
display(imgb)
imgb = rgbImage(red=img*0, green=img*0, blue=img)
display(imgb)
imgb = rgbImage(red=img>0.6, green=img*0.5, blue=img*1.5)
display(imgb)

hist(img)


## Low-pass disc-shaped filter
flo = makeBrush(11,shape='disc', step=FALSE)^2
flo = flo/sum(flo)
imgflo = filter2(imgc, flo)
display(imgflo)

## High-pass Laplacian filter
fhi = matrix(1, nc=3, nr=3)
fhi[2,2] = -8
imgfhi = filter2(imgc, fhi)
display(imgfhi)
display(1-imgfhi)
display(.8-imgfhi)

## Low-pass Gaussian filter
sigma=6
imgflg = gblur(img, sigma, radius = 2 * ceiling(3 * sigma) + 1)
display(imgflg)


## Median filter
imgm = medianFilter(img, 3)
display(imgm)


red=as.vector(Eg1[,,1])
green=as.vector(Eg1[,,2])
blue=as.vector(Eg1[,,3])
plot(density(red),col=2,ylim=c(0,6.5),lwd=2,main="",sub="",cex.axis=1.5,cex.lab=1.5)
points(density(green),col=3,lty=1,type="l",lwd=2)
points(density(blue),col=4,lty=1,type="l",lwd=2)


library(squash)
Mtr1 <- xyzmat2xyz(img[,,1])
Mtr2 <- xyzmat2xyz(img[,,2])
Mtr3 <- xyzmat2xyz(img[,,3])
Mtr <- cbind(Mtr1$z,Mtr2$z,Mtr3$z)


library(e1071)
N <- 4
fc<-cmeans(Mtr,centers=N,iter.max=100,verbose=FALSE,dist="euclidean",method="cmeans",m=1.1,rate.par=NULL,weights=1)
fcl<-fc$cluster
xyz <- cbind(Mtr1$x,Mtr1$y[length(Mtr1$y):1],fcl)


library(raster)
Eg1tr <- rasterFromXYZ(xyz)

image(Eg1tr,col=c(1:N),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")

image(Eg1tr,col=c(1,0,0,0,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,2,0,0,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,4,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,3,0,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,0,5,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,0,0,6,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,0,0,0,7),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")

fclm = fcl*0
fclm[which(fcl==1)]= "black"
fclm[which(fcl==2)]= "red"
fclm[which(fcl==3)]= "green3"
fclm[which(fcl==4)]= "blue"
fclm[which(fcl==5)]= "cyan"
fclm[which(fcl==6)]= "magenta"
fclm[which(fcl==7)]= "yellow"
qwfc = Image(fclm,dim=c(dim(img)[1],dim(img)[2]))
display(qwfc)           

fclm = rep("white",length(fcl))
fclm[which(fcl==1)]= "black"
qwfc1 = Image(fclm,dim=c(dim(img)[1],dim(img)[2]))
display(qwfc1)

CL = 1
fcmem<-fc$membership
fcm<-fcmem[,CL]
xyz <- cbind(Mtr1$x,Mtr1$y[length(Mtr1$y):1],fcm)
library(raster)
Eg1tM <- rasterFromXYZ(xyz)
image(Eg1tM,col=gray.colors(12, start = 0.1, end = 0.9, gamma = 2.2, alpha = NULL),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
qwfcm = Image(fcm,dim=c(dim(img)[1],dim(img)[2]))
display(qwfcm)



#####################################################################################################################################]

rgb(0, 1, 0)
rgb(0.3, 0.7, 0.9)
col2rgb(paste0("gold", 1:4))
col2rgb("peachpuff")
col2rgb(c(blu = "royalblue", reddish = "tomato")) 


library(plotrix)
color.id("#cc00cc")
sapply(rainbow(4), color.id)


library(ggtern)
rgb2hex(crgb[,1])
rgb2hex(255,0,0)

library(gplots)
col2hex(c("red","yellow","lightgrey"))


library(jjb)
rgb_to_hex(255,255,255) # Hexadecimal with pound sign
rgb_to_hex(255,255,255,FALSE) # Heaxadecimal without pound sign
shade(c(22, 150, 230), shade_factor = 0.5)
tint(rgb_value, tint_factor = 0.2)


library(aqp)
# Napthol Red PR112 (8R-5-14)
munsell2rgb(the_hue="8R", the_value=5, the_chroma=14, return_triplets=TRUE)
munsell2rgb(the_hue=c('10YR', '2.5YR'), the_value=c(3, 5), the_chroma=c(5, 6), return_triplets=TRUE)


library(munsellinterpol)
# Napthol Red PR112 (8R-5-14)
# Ultramarine Blue Reddish PB29 (7PB-2-12)
hns = HueNumberFromString( '8R' )
(mrgb = MunsellToRGB( c(hns,2,12), space='sRGB', maxSignal=255, adapt='Bradford'))
rgbc = as.numeric(unlist(mrgb)[5:7])
library(jjb)
col= rgb_to_hex(rgbc[1],rgbc[2],rgbc[3])
barplot(5,col=col, axes = FALSE)


# Transmitance to absorbance
library(photobiology)
T2A(0.99, action = NULL, byref = FALSE, clean = TRUE) # transmittance into absorbance

(crgb = col2rgb(w_length2rgb(c(400, 500, 600, 700), color.name=c("a","b","c","d"))))

col2rgb(w_length_range2rgb(500:600))
col2rgb(w_length_range2rgb(c(500,600)))