---
title: "R Notebook"
output:
html_document:
df_print: paged
---
V podniku byly instalovány 4 čističe oleje. Z každého byly odebráno 6 vzorků a stanoveny nečistoty, které filtr nezachytil. Rozhodněte, zda čističe propouštějí stejné množství nečistot (na hladině významnosti 5 %).
```{r}
xc1 = c(8, 12, 7, 2, 19, 7)
xc2 = c(18, 11, 17, 16, 10, 15)
xc3 = c(14, 6, 5, 12, 4, 13)
xc4 = c(15, 6, 10, 9, 14, 5)
data = as.data.frame(cbind(xc1,xc2, xc3,xc4))
data2 = stack(data)
boxplot(values~as.factor(ind), data = data2,ylab="",names=colnames(data),las=2,cex=1) # ,ylim=c(0,30)
stripchart(values~as.factor(ind), data = data2, pch=20, col=1,vert=TRUE,cex=1.5,add=TRUE)
library(beeswarm)
boxplot(values~as.factor(ind), data = data2, ylab="",xlab="",las=1,cex=1,outline=TRUE,cex.lab=1)
beeswarm(values~as.factor(ind), data = data2, col = 1, pch = 1, add = TRUE,horizontal=FALSE,pwcol = NULL,cex=1)
library(beanplot)
beanplot(values~as.factor(ind), data = data2, col = "bisque", ylim = "")
library(PASWR2)
oneway.plots(data2$values,as.factor(data2$ind))
modelr <- lm(values~ind, data = data2)
summary(modelr)
modela <- aov(values~as.factor(ind), data = data2)
summary(modela)
library(bayesanova)
assumption.check(x1=xc1,x2=xc2,x3=xc3,x4=xc4,x5=NULL,x6=NULL,conf.level=0.95)
#### shoda rozptylu
## Bartlett
library(onewaytests)
homog.test(values~as.factor(ind), data = data2, method = "Bartlett", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(RVAideMemoire)
perm.bartlett.test(values~as.factor(ind), data = data2, nperm = 999, progress = TRUE)
## Fligner Kileen
library(coin)
fligner_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores")
fligner_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(B = 1000))
# ties.method = c("mid-ranks", "average-scores"),
## Conover, Johnson, Johnson
library(coin)
conover_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores")
conover_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores", distribution = approximate(B = 1000))
## Levene
library(onewaytests)
homog.test(values~as.factor(ind), data = data2, method = "Levene", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(car)
car::leveneTest(values~as.factor(ind), data = data2, center=mean)
car::leveneTest(values~as.factor(ind), data = data2,center=median)
# center = c("median", "mean")
library(lawstat)
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="mean", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="median", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="trim.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")
## Brown Forsythe
library(onewaytests)
out = bf.test(values~as.factor(ind), data = data2, alpha = 0.05, na.rm = TRUE, verbose = TRUE)
out
# Hartley test
library(PMCMRplus)
hartleyTest(values~as.factor(ind), data = data2)
### outliers
# Bonferroni Outlier Test
library(car)
ou = outlierTest(values~as.factor(ind), data = data2, cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))
# Lund Outlier Test
library(rapportools)
rp.outlier(rstandard(values~as.factor(ind), data = data2))
library(referenceIntervals)
stresid<-rstandard(lm(values~as.factor(ind), data = data2))
vanderLoo.outliers(stresid)
cook.outliers(stresid)
### normalita
library(onewaytests)
nor.test(values~as.factor(ind), data = data2, method = "SW", alpha = 0.05, plot = "qqplot", mfrow = NULL, na.rm = TRUE, verbose = TRUE)
# method = c("SW", "SF", "LT", "AD", "CVM", "PT"),
# plot = c("qqplot-histogram", "qqplot", "histogram")
library(s20x)
normcheck(lm(values~as.factor(ind), data = data2), xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = TRUE, whichPlot = 1:2, usePar = TRUE)
fit <- aov(values~as.factor(ind), data = data2)
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")
library(FSA)
Summarize(values~as.factor(ind), data = data2, digits=3)
# Permutacni ANOVA
library(permuco)
modp <- aovperm(values~as.factor(ind), data = data2, np = 10000)
summary(modp)
# Bootstrap ANOVA
library(lmboot)
modb <- ANOVA.boot(values~as.factor(ind), data = data2, B = 1000, type = "residual", wild.dist = "normal", seed = NULL,keep.boot.resp = FALSE)
modb$"p-values"
# Box-Cox ANOVA
library(AID)
modbc <- boxcoxfr(data2$values, data2$ind, option = "both", lambda = seq(-3, 3, 0.01), lambda2 = NULL, tau = 0.05, alpha = 0.05, verbose = TRUE)
model_bc <- lm(modbc$tf.data ~ data2$ind)
confInt(modbc)
# Bayes ANOVA
result=bayes.anova(n=1000,first = xc1, second=xc2, third=xc3,fourth=xc4,ci=0.95)
anovaplot(result)
anovaplot(result, type="effect")
post.pred.check(result, ngroups = 4, out = c(xc1,xc2,xc3,xc4), reps = 25, eta = c(1/4,1/4,1/4,1/4))
# ANOVA zalozena na kvantilech
library(WRS2)
Qanova(values~as.factor(ind), data = data2, q = 0.5, nboot = 500)
library(coin)
median_test(values~as.factor(ind), data = data2, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95)
median_test(values~as.factor(ind), data = data2, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95, distribution = approximate(nresample = 1000))
```
Byl sledován vliv kombinace minerálních hnojiv (draselná, dusíkatá a fosfátová) na produkci rostlinné hmoty (v metrických centech na jednotku plochy). Každá varianta (varianty kombinace hnojiv a kontrolní vzorek bez hnojení) byla testována čtyřikrát.Ovlivňují hnojiva přírůstek produkce rostlinné hmoty? Zhodnoťte vliv jednotlivých kombinací hnojiv.
```{r}
bez = c(67, 67, 55, 42)
KN = c(98, 96, 91, 66)
KP = c(60, 69, 50, 35)
NP = c(79, 64, 81, 70)
KPN = c(90, 70, 79, 88)
data = as.data.frame(cbind(bez, KN, KP, NP, KPN))
data2 = stack(data)
boxplot(values~as.factor(ind), data = data2,ylab="",names=colnames(data),las=2,cex=1)
stripchart(values~as.factor(ind), data = data2, pch=20, col=1,vert=TRUE,cex=1.5,add=TRUE)
library(beeswarm)
boxplot(values~as.factor(ind), data = data2, ylab="",xlab="",las=1,cex=1,outline=TRUE,cex.lab=1)
beeswarm(values~as.factor(ind), data = data2, col = 1, pch = 1, add = TRUE,horizontal=FALSE,pwcol = NULL,cex=1)
library(beanplot)
beanplot(values~as.factor(ind), data = data2, col = "bisque", ylim = "")
library(ggplot2)
ggplot(data2, aes(x=ind, y=values)) +
geom_violin(trim=FALSE)+
geom_boxplot(width=0.1)
library(PASWR2)
oneway.plots(data2$values,as.factor(data2$ind))
library(bayesanova)
assumption.check(x1=bez,x2=KN,x3=KP,x4=NP,x5=KPN,x6=NULL,conf.level=0.95)
#### shoda rozptylu
## Bartlett
library(onewaytests)
homog.test(values~as.factor(ind), data = data2, method = "Bartlett", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(RVAideMemoire)
perm.bartlett.test(values~as.factor(ind), data = data2, nperm = 999, progress = TRUE)
## Fligner Kileen
library(coin)
fligner_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores")
fligner_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(B = 1000))
# ties.method = c("mid-ranks", "average-scores"),
## Conover, Johnson, Johnson
library(coin)
conover_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores")
conover_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores", distribution = approximate(B = 1000))
## Levene
library(onewaytests)
homog.test(values~as.factor(ind), data = data2, method = "Levene", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(car)
car::leveneTest(values~as.factor(ind), data = data2, center=mean)
car::leveneTest(values~as.factor(ind), data = data2,center=median)
# center = c("median", "mean")
library(lawstat)
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="mean", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="median", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="trim.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")
## Brown Forsythe
library(onewaytests)
out = bf.test(values~as.factor(ind), data = data2, alpha = 0.05, na.rm = TRUE, verbose = TRUE)
out
# Hartley test
library(PMCMRplus)
hartleyTest(values~as.factor(ind), data = data2)
### outliers
# Bonferroni Outlier Test
library(car)
ou = outlierTest(values~ind, data = data2, cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))
# Lund Outlier Test
library(rapportools)
rp.outlier(rstandard(values~as.factor(ind), data = data2))
library(referenceIntervals)
stresid<-rstandard(lm(values~as.factor(ind), data = data2))
vanderLoo.outliers(stresid)
cook.outliers(stresid)
### normalita
library(onewaytests)
nor.test(values~as.factor(ind), data = data2, method = "SW", alpha = 0.05, plot = "qqplot", mfrow = NULL, na.rm = TRUE, verbose = TRUE)
# method = c("SW", "SF", "LT", "AD", "CVM", "PT"),
# plot = c("qqplot-histogram", "qqplot", "histogram")
library(s20x)
normcheck(lm(values~as.factor(ind), data = data2), xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = TRUE, whichPlot = 1:2, usePar = TRUE)
fit <- aov(values~as.factor(ind), data = data2)
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")
library(FSA)
Summarize(values~as.factor(ind), data = data2, digits=3)
# Permutacni ANOVA
library(permuco)
modp <- aovperm(values~as.factor(ind), data = data2, np = 10000)
summary(modp)
# Bootstrap ANOVA
library(lmboot)
modb <- ANOVA.boot(values~as.factor(ind), data = data2, B = 1000, type = "residual", wild.dist = "normal", seed = NULL,keep.boot.resp = FALSE)
modb$"p-values"
# Box-Cox ANOVA
library(AID)
modbc <- boxcoxfr(data2$values, data2$ind, option = "both", lambda = seq(-3, 3, 0.01), lambda2 = NULL, tau = 0.05, alpha = 0.05, verbose = TRUE)
model_bc <- lm(modbc$tf.data ~ data2$ind)
summary(model_bc)
confInt(modbc)
# Bayes ANOVA
result=bayes.anova(n=1000,first = xc1, second=xc2, third=xc3,fourth=xc4,ci=0.95)
anovaplot(result)
anovaplot(result, type="effect")
post.pred.check(result, ngroups = 4, out = c(xc1,xc2,xc3,xc4), reps = 25, eta = c(1/4,1/4,1/4,1/4))
# ANOVA zalozena na kvantilech
library(coin)
median_test(values~as.factor(ind), data = data2, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95)
median_test(values~as.factor(ind), data = data2, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95, distribution = approximate(nresample = 1000))
# Post-hoc testy
# Student-Newman-Keuls (SNK) test
library(PMCMRplus)
out = snkTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,alternative = "two.sided")
summary(out)
summaryGroup(out)
# Tukey Honestly Significant Differences test, Tukey - Kramer
library(PMCMRplus)
out = tukeyTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,alternative = "two.sided")
summary(out)
summaryGroup(out)
## Scheffe test
library(PMCMRplus)
out = scheffeTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL)
summary(out)
summaryGroup(out)
# LSD test
library(PMCMRplus)
out = lsdTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL); out
summary(out)
summaryGroup(out)
library(rcompanion)
PT = pairwisePermutationTest(x=data2$values,g = as.factor(data2$ind),method = "fdr"); PT
cldList(p.adjust ~ Comparison,data = PT,threshold = 0.05)
# Dunnett test
library(PMCMRplus)
out = dunnettT3Test(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL); out
summary(out)
summaryGroup(out)
library(DescTools)
DunnettTest(x=data2$values, g=data2$ind)
## Kruskal-Wallis test # homoscedastic
library(NSM3)
pKW(data2$values,as.factor(data2$ind), method="Monte Carlo",n.mc=500) # Kruskal-Wallis
# method=c("Exact", "Monte Carlo", or "Asymptotic")
library(coin)
kruskal_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores")
kruskal_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000))
# Jonckheere Terpstra test
library(DescTools)
lomc = ordered(sapply(data2$ind,function(x) sub("V", "", x)))
JonckheereTerpstraTest(data2$values,lomc, alternative = "increasing", nperm = 999, exact = NULL)
# alternative = c("two.sided", "increasing", "decreasing")
# Post-hoc testy ke Kruskall Wallisove testu
# Nemenyi test
library(PMCMRplus)
out = kwAllPairsNemenyiTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,dist = "Tukey"); out
# dist = c("Tukey", "Chisquare")
summary(out)
summaryGroup(out)
# Conover's non-parametric all-pairs comparison test for Kruskal-type ranked data.
library(conover.test)
conover.test(data2$values,as.factor(data2$ind), 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(values~as.factor(ind), data = data2, 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)
# Dwass, Steel, Critchlow and Fligner
library(PMCMRplus)
out = dscfAllPairsTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL, p.adjust.method = "holm", alternative = "two.sided")
summary(out)
summaryGroup(out)
# Transformace na normalitu
# van-der-Waerden normal scores test
library(PMCMRplus)
vanWaerdenTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,alternative = "two.sided")
out = vanWaerdenAllPairsTest(values~as.factor(ind), data = data2, 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(coin)
normal_test(values~as.factor(ind), data = data2, ties.method = "mid-ranks", conf.int = FALSE, conf.level = 0.95)
normal_test(values~as.factor(ind), data = data2, 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(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,alternative = "two.sided")
out = normalScoresAllPairsTest(values~as.factor(ind), data = data2, 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)
```
ANOVA ze sumarnich dat
```{r}
Mean <- c(90,85,92,100,102,106)
SD <- c(9.035613,11.479667,9.760268,7.662572,9.830258,9.111457)
SampleSize <- c(9,9,9,9,9,9)
library(rpsychi)
fm1 <- ind.oneway.second(Mean, SD, SampleSize)
fm1
TukeyHSD(fm1$anova.table)
plot(TukeyHSD(fm1$anova.table, conf.level=.95), las = 2)
```
Porovnejte obsah mědi (v %) při třech různých technologiích lití bronzu (hladina významnosti 0.05).
```{r}
xt1 = c(80.1, 78.7, 79.1, 76.1, 78.1)
xt2 = c(82.8, 80.5, 83.6)
xt3 = c(79.8, 81.2, 82.7, 80.4)
data = as.data.frame(cbind(xt1,xt2,xt3))
data2 = stack(data)
boxplot(values~as.factor(ind), data = data2,ylab="",names=colnames(data),las=2,cex=1)
stripchart(values~as.factor(ind), data = data2, pch=20, col=1,vert=TRUE,cex=1.5,add=TRUE)
library(beeswarm)
boxplot(values~as.factor(ind), data = data2, ylab="",xlab="",las=1,cex=1,outline=TRUE,cex.lab=1)
beeswarm(values~as.factor(ind), data = data2, col = 1, pch = 1, add = TRUE,horizontal=FALSE,pwcol = NULL,cex=1)
library(beanplot)
beanplot(values~as.factor(ind), data = data2, col = "bisque", ylim = "")
library(ggplot2)
ggplot(data2, aes(x=ind, y=values)) +
geom_violin(trim=FALSE)+
geom_boxplot(width=0.1)
library(PASWR2)
oneway.plots(data2$values,as.factor(data2$ind))
library(bayesanova)
assumption.check(x1=xt1,x2=xt2,x3=xt3,x4=NULL,x5=NULL,x6=NULL,conf.level=0.95)
#### shoda rozptylu
## Bartlett
library(onewaytests)
homog.test(values~as.factor(ind), data = data2, method = "Bartlett", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(RVAideMemoire)
perm.bartlett.test(values~as.factor(ind), data = data2, nperm = 999, progress = TRUE)
## Fligner Kileen
library(coin)
fligner_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores")
fligner_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(B = 1000))
# ties.method = c("mid-ranks", "average-scores"),
## Conover, Johnson, Johnson
library(coin)
conover_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores")
conover_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores", distribution = approximate(B = 1000))
## Levene
library(onewaytests)
homog.test(values~as.factor(ind), data = data2, method = "Levene", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(car)
car::leveneTest(values~as.factor(ind), data = data2, center=mean)
car::leveneTest(values~as.factor(ind), data = data2,center=median)
# center = c("median", "mean")
library(lawstat)
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="mean", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="median", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(data2$values, group = as.factor(data2$ind), location="trim.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")
## Brown Forsythe
library(onewaytests)
out = bf.test(values~as.factor(ind), data = data2, alpha = 0.05, na.rm = TRUE, verbose = TRUE)
out
# Hartley test
library(PMCMRplus)
hartleyTest(values~as.factor(ind), data = data2)
### outliers
# Bonferroni Outlier Test
library(car)
ou = outlierTest(values~ind, data = data2, cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))
# Lund Outlier Test
library(rapportools)
rp.outlier(rstandard(values~as.factor(ind), data = data2))
library(referenceIntervals)
stresid<-rstandard(lm(values~as.factor(ind), data = data2))
vanderLoo.outliers(stresid)
cook.outliers(stresid)
### normalita
library(onewaytests)
nor.test(values~as.factor(ind), data = data2, method = "SW", alpha = 0.05, plot = "qqplot", mfrow = NULL, na.rm = TRUE, verbose = TRUE)
# method = c("SW", "SF", "LT", "AD", "CVM", "PT"),
# plot = c("qqplot-histogram", "qqplot", "histogram")
library(s20x)
normcheck(lm(values~as.factor(ind), data = data2), xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = TRUE, whichPlot = 1:2, usePar = TRUE)
fit <- aov(values~as.factor(ind), data = data2)
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")
library(FSA)
Summarize(values~as.factor(ind), data = data2, digits=3)
# Permutacni ANOVA
library(permuco)
modp <- aovperm(values~as.factor(ind), data = data2, np = 10000)
summary(modp)
# Bootstrap ANOVA
library(lmboot)
modb <- ANOVA.boot(values~as.factor(ind), data = data2, B = 1000, type = "residual", wild.dist = "normal", seed = NULL,keep.boot.resp = FALSE)
modb$"p-values"
# Box-Cox ANOVA
library(AID)
modbc <- boxcoxfr(data2$values, data2$ind, option = "both", lambda = seq(-3, 3, 0.01), lambda2 = NULL, tau = 0.05, alpha = 0.05, verbose = TRUE)
model_bc <- lm(modbc$tf.data ~ data2$ind)
summary(model_bc)
confInt(modbc)
# Bayes ANOVA
result=bayes.anova(n=1000,first = xc1, second=xc2, third=xc3,fourth=xc4,ci=0.95)
anovaplot(result)
anovaplot(result, type="effect")
post.pred.check(result, ngroups = 4, out = c(xc1,xc2,xc3,xc4), reps = 25, eta = c(1/4,1/4,1/4,1/4))
# ANOVA zalozena na kvantilech
library(coin)
median_test(values~as.factor(ind), data = data2, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95)
median_test(values~as.factor(ind), data = data2, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95, distribution = approximate(nresample = 1000))
# Post-hoc testy
# Student-Newman-Keuls (SNK) test
library(PMCMRplus)
out = snkTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,alternative = "two.sided")
summary(out)
summaryGroup(out)
# Tukey Honestly Significant Differences test, Tukey - Kramer
library(PMCMRplus)
out = tukeyTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,alternative = "two.sided")
summary(out)
summaryGroup(out)
## Scheffe test
library(PMCMRplus)
out = scheffeTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL)
summary(out)
summaryGroup(out)
# LSD test
library(PMCMRplus)
out = lsdTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL); out
summary(out)
summaryGroup(out)
# Parovy permutacni test
library(rcompanion)
PT = pairwisePermutationTest(x=data2$values,g = as.factor(data2$ind),method = "fdr"); PT
cldList(p.adjust ~ Comparison,data = PT,threshold = 0.05)
# trimmed means
library(WRS2)
t1way(values~as.factor(ind), data = data2, tr = 0.2, alpha = 0.05, nboot = 500)
vx=lincon(values~as.factor(ind), data = data2, tr = 0.2, alpha = 0.05)
## ANOVA heteroskedastic
## Alexander-Govern test.
library(onewaytests)
out = ag.test(values~as.factor(ind), data = data2, 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(welchADF)
summary(welchADF.test(values~as.factor(ind), data = data2))
library(coin)
oneway_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores")
oneway_test(values~as.factor(ind), data = data2, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000))
# Games-Howell
library(PMCMRplus)
out = gamesHowellTest(values~as.factor(ind), data = data2,alternative = "two.sided",subset=NULL, p.adjust="BH")
summary(out)
summaryGroup(out)
# Tamhane T2
library(PMCMRplus)
out = tamhaneT2Test(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,welch = TRUE,alternative = "two.sided")
summary(out)
summaryGroup(out)
# Ury, Wiggins a Hochberg
library(PMCMRplus)
out = uryWigginsHochbergTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL,alternative = "two.sided", p.adjust.method ="holm")
# p.adjust.method = p.adjust.methods, ...)
summary(out)
summaryGroup(out)
#### shoda distribuce
# Anderson Darling
library(PMCMRplus)
adKSampleTest(values~as.factor(ind), data = data2, subset=NULL, na.action=NULL)
out = adAllPairsTest(values~as.factor(ind), data = data2, 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)
# Murakami k-sample Baumgartner Weiss Schindler Test of Equal Distributions
library(PMCMRplus)
bwsKSampleTest(values~as.factor(ind), data = data2, subset=NULL, na.action,nperm = 1000)
out <- bwsAllPairsTest(values~as.factor(ind), data = data2, p.adjust="holm"); out
summary(out)
summaryGroup(out)
```
Bylo provedeno 12 měření průměru strojní součástky: aritmetický průměr byl 14.306 mm a směrodatná odchylka 0.572 mm. Odhadněte střední hodnotu základního souboru pomocí 95% intervalu spolehlivosti. Určete přibližný rozsah výběru, který umožní odhadnout skutečnou hodnotu průměru s přeností 0.2 mm.
```{r}
m = 14.306
s = 0.572
n = 12
m + c(-s, s)*qt(0.975,df=n-1)/sqrt(n)
library(LearningStats)
Mean.CI(m, sigma = NULL, sc = NULL, s = s, n = n, conf.level= 0.95)$CI
n1 = 15
qt(0.975,(n1-1))*s/sqrt(n1)
n2 = 25
qt(0.975,(n2-1))*s/sqrt(n2)
n3 = 35
qt(0.975,(n3-1))*s/sqrt(n3)
n4 = 34
qt(0.975,(n4-1))*s/sqrt(n4)
sigma = 0.572
delta = 0.2
qt = 3
nt = ((sigma*qt)/delta)^2; nt
# iterace
nt = nt-1
nt = ((sigma*qt(0.975,(nt-1)))/delta)^2; nt
```
25 ml 0.1 M NaOH bylo titrováno 0.1 M HCl (v ml), a to jednak bez zahřátí a jednak po zahřátí NaOH. Celkem bylo provedeno 29 dvojic stanovení.
```{r}
xpr = c(26.55, 26.60, 26.80, 26.45, 26.40,
26.70, 26.80, 26.75, 26.70, 26.50,
26.40, 26.40, 26.30, 26.40, 26.60,
26.45, 26.70, 26.65, 26.55, 26.50,
26.60, 26.70, 26.75, 26.60, 26.60,
26.75, 26.70, 26.75, 26.6)
xpo = c(26.80, 26.80, 26.95, 26.70, 26.75,
26.90, 26.95, 26.90, 26.85, 26.80,
26.55, 26.60, 26.60, 26.60, 26.75,
26.65, 26.90, 26.90, 26.85, 26.65,
26.75, 26.85, 26.85, 26.80, 26.85,
26.90, 26.90,26.95, 26.80)
XP = as.data.frame(cbind(xpr,xpo))
xp = stack(XP)
colnames(xp) = c("Spotreba HCl","Zahrati")
## Deming regrese
library(deming)
demr = deming(xpo~xpr,conf=.95,jackknife=TRUE,model=TRUE)
demr
(cor(xpr,xpo))^2 # R^2
ggplot(XP, aes(x=xpr, y=xpo))+
geom_point(shape= 16, size = 2)+
scale_x_continuous(limits = c(26.2, 27))+scale_y_continuous(limits = c(26.2, 27))+
labs(title="Deming regression",x="Spotreba pred zahratim (ml)",y="Spotreba po zahrati (ml)",size=5)+
geom_abline(mapping = NULL,data = NULL,slope= demr$coefficients[2], intercept= demr$coefficients[1], na.rm = FALSE, show.legend = NA,col = "#C42126",size = 0.5)+
geom_abline(mapping = NULL,data = NULL,slope= 1, intercept= 0, na.rm = FALSE, show.legend = NA, size = 0.5,linetype = 2,col="black")+
theme(plot.title = element_text(hjust = 0.5))
# Passing Bablok
library(deming)
pbr <- pbreg(xpo~xpr, conf=.95,nboot = 1000, method=1, x = FALSE, y = FALSE, model = TRUE)
pbr
ggplot(XP, aes(x=xpr, y=xpo))+
geom_point(shape= 16, size = 2)+
scale_x_continuous(limits = c(26.2, 27))+scale_y_continuous(limits = c(26.2, 27))+
labs(title="Passing-Bablok regression",x="Spotreba pred zahratim (ml)",y="Spotreba po zahrati (ml)",size=5)+
geom_abline(mapping = NULL,data = NULL,slope= pbr$coefficients[2], intercept= pbr$coefficients[1], na.rm = FALSE, show.legend = NA,col = "#C42126",size = 0.5)+
geom_abline(mapping = NULL,data = NULL,slope= 1, intercept= 0, na.rm = FALSE, show.legend = NA, size = 0.5,linetype = 2,col="black")+
theme(plot.title = element_text(hjust = 0.5))
library(PairedData)
summary(paired(xpr, xpo),tr=0.1)
effect.size(paired(xpr, xpo),tr=0.1)
slidingchart(paired(xpr, xpo))
library(blandr)
bas2 = blandr.statistics (xpr,xpo, sig.level=0.95)
bas2$bias
bas2$biasUpperCI
bas2$biasLowerCI
bas2$regression.equation
cor(bas2$means,bas2$differences)^2 # R^2
blandr.plot.qq(bas2)
blandr.display.and.draw(method1=xpr, method2=xpo, plotter = "ggplot",method1name = "Spotreba pred zahratim (ml)", method2name = "Spotreba po zahrati (ml)",plotTitle = "Bland-Altmanuv graf",sig.level = 0.95, annotate = TRUE, ciDisplay = TRUE,ciShading = TRUE, normalLow = FALSE, normalHigh = FALSE,lowest_y_axis = FALSE, highest_y_axis = FALSE, point_size = 2)
library(s20x)
normcheck((xpo-xpr))
library(PairedData)
pd <- paired(xpr, xpo)
plot(pd, type = "profile") + theme_bw()
library(parviol)
parviol(XP, violinplot=TRUE, main=NULL, sub=NULL)
parviol(XP, violinplot=FALSE, main=NULL, sub=NULL)
# t-test
t.test(xpr, xpo, paired = TRUE, alternative = "two.sided")
# t.test(elem ~ lom, data = data, paired = TRUE, alternative = "two.sided")
# Bootstrap t-Test
library(MKinfer)
boot.t.test(xpr, xpo, alternative = "two.sided",mu = 0, paired = TRUE, var.equal = FALSE, conf.level = 0.95, R = 9999, symmetric = FALSE)
# alternative = c("two.sided", "less", "greater")
# Yuen's trimmed mean test
library(PairedData)
yuen.t.test(xpr, xpo, tr = 0.2, alternative = "two.sided", mu = 0, paired = TRUE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")
# Bayesovsky t-test
library(Bolstad)
bayes.t.test(xpr, xpo,alternative = "two.sided", mu = 0, paired = TRUE,var.equal = TRUE, conf.level = 0.95,prior = "jeffreys", m = NULL,n0 = NULL,sig.med = NULL, kappa = 1,sigmaPrior = "chisq",nIter = 10000, nBurn = 1000)
# alternative = c("two.sided", "less", "greater")
# prior = c("jeffreys", "joint.conj"),
## Wilcoxonuv parovy test
wilcox.test(xpr, xpo, paired = TRUE, alternative = "two.sided", conf.int = TRUE)
library(PairedData)
wilcox.test(paired(xpr, xpo))
## medianovy (znamenkovy) parovy test
library(EnvStats)
signTest(xpr, xpo, alternative = "two.sided", mu = 0, paired = TRUE, conf.level = 0.95)
# Permutation test for paired data
library(CarletonStats)
permTestPaired(xpr, xpo, B = 9999,alternative = "two.sided", plot.hist = TRUE,legend.loc = "topright", plot.qq = TRUE,x.name = deparse(substitute(x)), y.name = deparse(substitute(y)))
```
Vědci chtějí zjistit, zda čtyři různé léky vedou k rozdílným reakčním časům. Aby to ověřili, změřili reakční dobu pěti pacientů užívajících čtyři různé léky. Vzhledem k tomu, že u každého pacienta je měřena reakční doba na každém ze čtyř léků, použijeme metodu ANOVA s opakovanými měřeními, abychom zjistili, zda se průměrná reakční doba u jednotlivých léků liší.
```{r}
dmed <- data.frame(patient=rep(1:5, each=4),
drug=rep(1:4, times=5),
response=c(30, 28, 16, 34,
14, 18, 10, 22,
24, 20, 18, 30,
38, 34, 20, 44,
26, 28, 14, 30))
library(onewaytests)
homog.test(response~as.factor(drug), data = dmed, method = "Bartlett", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(RVAideMemoire)
perm.bartlett.test(response~as.factor(drug), data = dmed, nperm = 999, progress = TRUE)
## Levene
library(onewaytests)
homog.test(response~as.factor(drug), data = dmed, method = "Levene", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(car)
car::leveneTest(concentration~as.factor(method), data = cyc, center=mean)
car::leveneTest(response~as.factor(drug), data = dmed,center=median)
# center = c("median", "mean")
library(lawstat)
lawstat::levene.test(dmed$response, group=dmed$drug, location="mean", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(dmed$response, group=dmed$drug, location="median", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(dmed$response, group=dmed$drug, location="trim.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")
### normalita
library(onewaytests)
nor.test(response~as.factor(drug), data = dmed, method = "SW", alpha = 0.05, plot = "qqplot", mfrow = NULL, na.rm = TRUE, verbose = TRUE)
# method = c("SW", "SF", "LT", "AD", "CVM", "PT"),
# plot = c("qqplot-histogram", "qqplot", "histogram")
library(s20x)
normcheck(lm(response~as.factor(drug), data = dmed), xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = TRUE, whichPlot = 1:2, usePar = TRUE)
model <- aov(response~factor(drug)+Error(factor(patient)), data = dmed)
summary(model)
library(MKinfer)
rm.oneway.test(x=dmed$response, g=dmed$drug, id=dmed$patient, method = "aov")
# method = "aov","lme","friedman","quade"
rm.oneway.test(x=dmed$response, g=dmed$drug, id=dmed$patient, method = "lme")
# method = "aov","lme","friedman","quade"
rm.oneway.test(x=dmed$response, g=dmed$drug, id=dmed$patient, method = "friedman")
# method = "aov","lme","friedman","quade"
rm.oneway.test(x=dmed$response, g=dmed$drug, id=dmed$patient, method = "quade")
# method = "aov","lme","friedman","quade"
```
K porovnání tří laboratorních metod stanovení procentního podílu cyklamátu sodného bylo Všemi metodami analyzováno 12 vzorků komerčně vyráběného pomerančového nápoje. Pomocí Friedmanova testu zjistěte, zda mezi těmito třemi metodami existuje statisticky významný rozdíl.
```{r}
cycl <- read.delim("c:\\Users\\lubop\\Dropbox\\kurs\\cyclamate.txt", header = TRUE)
head(cycl)
cycl = as.data.frame(cycl)
rownames(cycl) = cycl[,1]
cycl1 = as.matrix(cycl[,-1])
library(reshape2)
cyc = setNames(melt(cycl1), c('sample', 'method', 'concentration'))
# cyc$sample = as.character(cyc$sample)
boxplot(cyc$concentration~cyc$method)
library(parviol)
parviol(cycl1, violinplot=TRUE, main=NULL, sub=NULL)
parviol(cycl1, violinplot=FALSE, main=NULL, sub=NULL)
library(onewaytests)
homog.test(concentration~as.factor(method), data = cyc, method = "Bartlett", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(RVAideMemoire)
perm.bartlett.test(concentration~as.factor(method), data = cyc, nperm = 999, progress = TRUE)
## Levene
library(onewaytests)
homog.test(concentration~as.factor(method), data = cyc, method = "Levene", alpha = 0.05, na.rm = TRUE, verbose = TRUE)
library(car)
car::leveneTest(concentration~as.factor(method), data = cyc, center=mean)
car::leveneTest(concentration~as.factor(method), data = cyc,center=median)
# center = c("median", "mean")
library(lawstat)
lawstat::levene.test(cyc$concentration, group = as.factor(cyc$method), location="mean", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(cyc$concentration, group = as.factor(cyc$method), location="median", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
lawstat::levene.test(cyc$concentration, group = as.factor(cyc$method), location="trim.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")
### normalita
library(onewaytests)
nor.test(concentration~as.factor(method), data = cyc, method = "SW", alpha = 0.05, plot = "qqplot", mfrow = NULL, na.rm = TRUE, verbose = TRUE)
# method = c("SW", "SF", "LT", "AD", "CVM", "PT"),
# plot = c("qqplot-histogram", "qqplot", "histogram")
library(s20x)
normcheck(lm(concentration~as.factor(method), data = cyc), xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = TRUE, whichPlot = 1:2, usePar = TRUE)
model <- aov(concentration~factor(method)+Error(factor(sample)), data = cyc)
summary(model)
library(MKinfer)
rm.oneway.test(x=cyc$concentration, g=cyc$method, id=dmed$sample, method = "aov")
# method = "aov","lme","friedman","quade"
rm.oneway.test(x=cyc$concentration, g=cyc$method, id=dmed$sample, method = "lme")
# method = "aov","lme","friedman","quade"
## Friedman test
friedman.test(y=cyc$concentration, groups=cyc$method, blocks=cyc$sample)
library(MKinfer)
rm.oneway.test(x=cyc$concentration, g=cyc$method, id=cyc$sample, method = "friedman")
# method = "aov","lme","friedman","quade"
# Quade test
rm.oneway.test(x=dmed$response, g=dmed$drug, id=dmed$patient, method = "quade")
# method = "aov","lme","friedman","quade"
library(PMCMR)
posthoc.friedman.nemenyi.test(y=cyc$concentration, groups=cyc$method, blocks=cyc$sample, p.adjust.method = "fdr")
library(PMCMRplus)
frdAllPairsExactTest(y=cyc$concentration, groups=cyc$method, blocks=cyc$sample, p.adjust.method = "fdr")
library(pgirmess)
friedmanmc(y=cyc$concentration, groups=cyc$method, blocks=cyc$sample,probs=0.05)
GT = xtabs(concentration ~ method + sample, data = cyc); GT
library(DescTools)
KendallW(GT, correct=TRUE, test=TRUE)
library(rstatix)
friedman_effsize(data = cyc,concentration ~ method | sample, ci = TRUE, conf.level = 0.95, ci.type = "perc", nboot = 500)
```
Pomoci korelace ověřte vztah mezi výškou rostliny a velikostí produkce u rajčat dané odrůdy.
```{r}
vys = c(21, 33, 18, 26, 28, 22, 27, 25, 26, 29, 31, 19, 20, 22, 26, 30)
prod = c(27, 42, 26, 35, 39, 31, 38, 22, 33, 36, 41, 23, 29, 25, 32, 37)
# scatter plot
plot(vys,prod)
library(UsingR)
scatter.with.hist(vys,prod, hist.col = gray(0.95),trend.line = "lm")
library(aplpack)
bagplot(vys,prod,show.baghull=TRUE,show.loophull=TRUE,precision=1,dkmethod=2)
# rank scatter plot
xr <-rank(vys)
yr <-rank(prod)
plot(xr,yr, col=1, xlab=" x", ylab="y",cex=1.2)
abline(0,1,lty=2)
# Chi-plot
library(asbio)
chi.plot(vys,prod)
# Kendalluv graf
library(lcopula)
K.plot(cbind(vys,prod), add = F)
# Pearsonuv korelacni koeficient
cp = cor(vys,prod, method="pearson"); cp
# bootstrap korelacni koeficient
library(CarletonStats)
bootCor(vys,prod, conf.level = 0.95, B = 10000,plot.hist = TRUE, hist.title = NULL, plot.qq = FALSE,legend.loc = "topright", x.name = deparse(substitute(x)),y.name = deparse(substitute(y)))
# permutacni korelacni koeficient
library(CarletonStats)
permTestCor(vys,prod, B = 999, alternative = "two.sided",plot.hist = TRUE, legend.loc = "topright", plot.qq = FALSE,x.name = deparse(substitute(x)), y.name = deparse(substitute(y)))
# robustni Pearsonuv korelacni koeficient
library(robustbase)
covMcd(cbind(vys,prod), alpha=0.95, cor=TRUE)$cor[2,1]
library(StatDA)
RobCor.plot(vys,prod, quan=1/2, alpha=0.025, colC=3, colR=4)
library(confintr)
ci_cor(vys,prod, probs = c(0.025, 0.975),method ="pearson", type = "normal")
# type = c("normal", "bootstrap")
ci_cor(vys,prod, probs = c(0.025, 0.975),method ="pearson", type = "bootstrap", boot_type = "norm",R = 9999,seed = NULL)
# boot_type = c("bca", "perc", "norm", "basic")
### test vyznamnosti
cor.test(vys,prod, alternative = "two.sided", method = "pearson", exact = FALSE, conf.level = 0.95, continuity = FALSE)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "kendall", "spearman")
library(jmuOutlier)
perm.cor.test(vys,prod, alternative = "two.sided", method = "pearson", num.sim = 20000)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "spearman")
## sila testu
library(WebPower)
wp.correlation(n = length(vys), r = cp, power = NULL, p = 0, rho0 = 0, alpha = 0.05, alternative = "two.sided")
# alternative = c("two.sided", "less", "greater")
## bayesovsky test
library(BFpack)
bpct = cor_test(as.data.frame(cbind(vys,prod)), formula = NULL, iter = 5000)
bpct$correstimates
library(BayesFactor)
samples = correlationBF(vys,prod,posterior = TRUE, iterations = 5000)
plot(samples[,"rho"])
## rs vs. n vs. confidence width
library(presize)
#n = length(sraz)
conf.width = cirs$interval[2]-cirs$interval[1]
conf.width = conf.width/2
pcsp = prec_cor(rs,n = NULL, conf.width = conf.width,conf.level = 0.95,method = "pearson")
pcsp$n
# pcsp$conf.width
# method = c("pearson", "kendall", "spearman")
# Spearmanuv koeficient
rs = cor(vys,prod, method="spearman"); rs
library(fmsb)
spearman.ci.sas(vys,prod, adj.bias=TRUE, conf.level=0.95)
library(confintr)
library(confintr)
ci_cor(vys,prod, probs = c(0.025, 0.975),method ="spearman", type = "normal")
# type = c("normal", "bootstrap")
ci_cor(vys,prod, probs = c(0.025, 0.975),method ="spearman", type = "bootstrap", boot_type = "norm",R = 9999,seed = NULL)
# boot_type = c("bca", "perc", "norm", "basic")
### test vyznamnosti
cor.test(vys,prod, alternative = "two.sided", method = "spearman", exact = FALSE, conf.level = 0.95, continuity = FALSE)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "kendall", "spearman")
library(pspearman)
spearman.test(vys,prod,alternative = "two.sided",approximation = "exact")
# alternative = c("two.sided", "less", "greater"),
# approximation = c("exact", "AS89", "t-distribution"))
library(jmuOutlier)
perm.cor.test(vys,prod, alternative = "two.sided", method = "spearman", num.sim = 20000)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "spearman")
## rs vs. n vs. confidence width
library(presize)
#n = length(sraz)
conf.width = cirs$interval[2]-cirs$interval[1]
conf.width = conf.width/2
pcsp = prec_cor(rs,n = NULL, conf.width = conf.width,conf.level = 0.95,method = "spearman")
pcsp$n
# pcsp$conf.width
# method = c("pearson", "kendall", "spearman")
```
Při výrobě plastu byl pozorován vliv aditiva (v g) na pevnost materiálu (kp/cm2). Vyjádřete tento vztah pomocí korelace.
```{r}
adit = c(1.50, 1.80, 1.22, 0.73, 1.55, 1.73, 2.30, 1.89, 1.20, 0.72, 1.95, 0.94, 2.16, 1.68, 2.02, 1.77, 0.65, 1.84, 2.40, 1.51, 1.20, 2.39, 2.12, 1.11, 1.60, 0.96, 2.23, 1.87, 1.00, 1.23)
pevn = c(20.7, 23.9, 22.7, 19.1, 24.8, 21.2, 25.6, 26.8, 20.2, 20.5, 23.3, 23.0, 26.6, 22.8, 24.2, 25.5, 18.4, 24.6, 26.6, 23.1, 24.5, 28.1, 25.3, 20.6, 22.0, 21.2, 23.9, 22.7, 18.9, 19.2)
# scatter plot
plot(adit,pevn)
library(UsingR)
scatter.with.hist(adit,pevn, hist.col = gray(0.95),trend.line = "lm")
library(aplpack)
bagplot(adit,pevn,show.baghull=TRUE,show.loophull=TRUE,precision=1,dkmethod=2)
# rank scatter plot
xr <-rank(adit)
yr <-rank(pevn)
plot(xr,yr, col=1, xlab=" x", ylab="y",cex=1.2)
abline(0,1,lty=2)
# Chi-plot
library(asbio)
chi.plot(adit,pevn)
# Kendalluv graf
library(lcopula)
K.plot(cbind(adit,pevn), add = F)
# korelacni koeficient
cp = cor(adit,pevn, method="pearson"); cp
# bootstrap korelacni koeficient
library(CarletonStats)
bootCor(adit,pevn, conf.level = 0.95, B = 10000,plot.hist = TRUE, hist.title = NULL, plot.qq = FALSE,legend.loc = "topright", x.name = deparse(substitute(x)),y.name = deparse(substitute(y)))
# permutacni korelacni koeficient
library(CarletonStats)
permTestCor(adit,pevn, B = 999, alternative = "two.sided",plot.hist = TRUE, legend.loc = "topright", plot.qq = FALSE,x.name = deparse(substitute(x)), y.name = deparse(substitute(y)))
# robustni korelacni koeficient
library(robustbase)
covMcd(cbind(adit,pevn), alpha=0.95, cor=TRUE)$cor[2,1]
library(StatDA)
RobCor.plot(adit,pevn, quan=1/2, alpha=0.025, colC=3, colR=4)
library(confintr)
ci_cor(adit,pevn, probs = c(0.025, 0.975),method ="pearson", type = "normal")
# type = c("normal", "bootstrap")
ci_cor(adit,pevn, probs = c(0.025, 0.975),method ="pearson", type = "bootstrap", boot_type = "norm",R = 9999,seed = NULL)
# boot_type = c("bca", "perc", "norm", "basic")
### test vyznamnosti
cor.test(adit,pevn, alternative = "two.sided", method = "pearson", exact = FALSE, conf.level = 0.95, continuity = FALSE)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "kendall", "spearman")
library(jmuOutlier)
perm.cor.test(adit,pevn, alternative = "two.sided", method = "pearson", num.sim = 20000)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "spearman")
## sila testu
library(WebPower)
wp.correlation(n = length(adit), r = cp, power = NULL, p = 0, rho0 = 0, alpha = 0.05, alternative = "two.sided")
# alternative = c("two.sided", "less", "greater")
## bayesovsky test
library(BFpack)
bpct = cor_test(as.data.frame(cbind(adit,pevn)), formula = NULL, iter = 5000)
bpct$correstimates
library(BayesFactor)
samples = correlationBF(adit,pevn,posterior = TRUE, iterations = 5000)
plot(samples[,"rho"])
## rs vs. n vs. confidence width
library(presize)
#n = length(sraz)
conf.width = cirs$interval[2]-cirs$interval[1]
conf.width = conf.width/2
pcsp = prec_cor(rs,n = NULL, conf.width = conf.width,conf.level = 0.95,method = "pearson")
pcsp$n
# pcsp$conf.width
# method = c("pearson", "kendall", "spearman")
#Spearmanuv koeficient
rs = cor(vys,prod, method="spearman"); rs
library(fmsb)
spearman.ci.sas(adit,pevn, adj.bias=TRUE, conf.level=0.95)
library(confintr)
library(confintr)
ci_cor(adit,pevn,probs = c(0.025, 0.975),method ="spearman", type = "normal")
# type = c("normal", "bootstrap")
ci_cor(adit,pevn, probs = c(0.025, 0.975),method ="spearman", type = "bootstrap", boot_type = "norm",R = 9999,seed = NULL)
# boot_type = c("bca", "perc", "norm", "basic")
### test vyznamnosti
cor.test(adit,pevn, alternative = "two.sided", method = "spearman", exact = FALSE, conf.level = 0.95, continuity = FALSE)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "kendall", "spearman")
library(pspearman)
spearman.test(adit,pevn,alternative = "two.sided",approximation = "exact")
# alternative = c("two.sided", "less", "greater"),
# approximation = c("exact", "AS89", "t-distribution"))
library(jmuOutlier)
perm.cor.test(adit,pevn, alternative = "two.sided", method = "spearman", num.sim = 20000)
# alternative = c("two.sided", "less", "greater")
# method = c("pearson", "spearman")
## rs vs. n vs. confidence width
library(presize)
#n = length(sraz)
conf.width = cirs$interval[2]-cirs$interval[1]
conf.width = conf.width/2
pcsp = prec_cor(rs,n = NULL, conf.width = conf.width,conf.level = 0.95,method = "spearman")
pcsp$n
# pcsp$conf.width
# method = c("pearson", "kendall", "spearman")
```
Byla měřena závislost součinitele tření na při toku kapalin potrubím na Reynoldsově kriteriu. Nalezněte rovnici vystihující závislost.
lam = a*Re^b
```{r}
xRe = c(1e4, 2e4, 5e4, 1e5, 1.5e5, 2e5)
xla = c(3.2e2, 2.6e2, 2.0e2, 1.9e2, 1.6e2, 1.4e2)
data = as.data.frame(cbind(xla,xRe))
# log(lam) = log(a) + b*log(Re)
plot(xRe,xla)
plot(xRe,xla,log="xy")
plot(log(xRe),log(xla))
model1 = lm(log(xla)~log(xRe))
summary(model1)
confint(model1)
koef = coef(model1)
# model1$coefficients
library(robustbase)
model11 = lmrob(log(xla)~log(xRe),data=data, na.action= na.omit,weights=NULL)
summary(model11)
confint(model11)
library(MASS)
model12 = lqs(log(xla)~log(xRe), data=data, method = "lts")
# method = c("lts", "lqs", "lms", "S", "model.frame")
model12
confint(model12)
model13 = rlm(log(xla)~log(xRe), data=data, weights=NULL, method = "M", wt.method = "inv.var")
# method = c("M", "MM", "model.frame"),
# wt.method = c("inv.var", "case")
model13
confint(model13)
# Nelinearni regrese
a = koef[1]
b = koef[2]
start <- list(a=coef(model1[1]), b=coef(model1)[2])
rnls <- nls(xla ~ a*xRe^b, data = data, start = start)
smmr = summary(rnls); smmr
confint(rnls)
library(investr)
plotFit(rnls, interval = "both", pch = 19, shade = TRUE, col.conf = "skyblue4", col.pred = "lightskyblue2")
library(ellipse)
plot(ellipse(rnls, which = c(1, 3), level = 0.95),type="l")
points(smmr$coefficients[1,1],smmr$coefficients[3,1])
```
Vyjádřete závislost hustoty tepelného toku (v W/m^2) na teplotním rozdílu mezi topnou parou a roztokem, vroucím při 117 °C v odpařováku s šikmými trubkami.
```{r}
xt = c(4.6, 5.6, 5.8, 7.1, 7.3, 8.6, 8.7, 10.0, 10.4, 10.5)
xq = c(28e2, 107e2, 58e2, 181e2, 170e2, 250e2, 305e2, 424e2, 411e2, 416e2)
data = as.data.frame(cbind(xt,xq))
plot(xt,xq)
model2 = lm(xt~poly(xq,2)) # jen pro testovani, vypocet pomoci orthogonalnich polynomu
summary(model2)
model2a = lm(xt~xq+I(xq^2)) # pro vypocet koeficientu
summary(model2a)
model1 = lm(xt~xq)
summary(model1)
AIC(model1,model2)
library(AICcmodavg)
model.set <- list(model1, model2)
model.names <- c("linearni", "kvadraticky")
aictab(model.set, modnames = model.names)
library(performance)
compare_performance(model1,model2, rank = TRUE)
library(investr)
plotFit(model1, data=data,interval = "both", pch = 19, shade = TRUE, col.conf = "skyblue4", col.pred = "lightskyblue2")
library(s20x)
ciReg(model1, conf.level = 0.95, print.out = TRUE)
library(car)
confidenceEllipse(model1)
library(car)
residualPlots(model1, terms = ~., layout = NULL, main = "", fitted = TRUE, AsIs=TRUE, plot = TRUE,tests = TRUE)
influencePlot(model1)
infIndexPlot(model1)
library(s20x)
normcheck(model1, 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)
library(s20x)
eovcheck(model1, xlab = "Fitted values",ylab = "Residuals", col = NULL, smoother = FALSE, twosd = FALSE,levene = FALSE)
library(s20x)
modcheck(model1, 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)))
xre = resid(model1)
xtp = predict(model1)
cbind(xV,xt,round(xtp,2))
# Breusch-Pagan test na heteroskedasticitu
library(lmtest)
lmtest::bptest(xq~xt, varformula = NULL, studentize = TRUE, data = data)
# Harrison-McCabe test for heteroskedasticity
library(lmtest)
lmtest::hmctest(xq~xt, 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(model1)
# Goldfeld-Quandt test against heteroskedasticity
library(lmtest)
lmtest::gqtest(xq~xt, 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(xq~xt, order.by = NULL, alternative = "two.sided", iterations = 15, exact = NULL, tol = 1e-10, data = data)
# alternative = c("greater", "two.sided", "less")
library(car)
# Durbin-Watson test
durbinWatsonTest(model1)
# Shapiro-Wilk test rezidui
shapiro.test(resid(model1))
shapiro.test(rstandard(model1))
shapiro.test(rstudent(model1))
# Bonferroni Outlier Test
library(car)
ou = outlierTest(model1, cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))
# Lund test
library(rapportools)
rp.outlier(rstandard(model1))
# Cook distance plot for outliers
cooksD <- cooks.distance(model2)
plot(cooksD, main = "Cooks Distance")
abline(h = 4/nrow(data), lty = 2, col = "steelblue")
library(s20x)
cooks20x(model2, 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))
abline(h = 4/nrow(data),col=2,lty=2)
library(fdm2id)
cookplot(model1)
# Rainbow test linearity
library(lmtest)
lmtest::raintest(xq~xt, fraction = 0.5, order.by = NULL, center = NULL, data=data)
# Harvey-Collier test for linearity
library(lmtest)
lmtest::harvtest(xq~xt, order.by = NULL, data = data)
# Robustni regrese
library(robustbase)
xvt2 = lmrob(xt~xV+I(xV^2),data=data, na.action=na.omit,weights=NULL)
summary(xvt2)
confint(xvt2)
library(MASS)
xvt3 = lqs(xt~xV+I(xV^2), data=data, method = "lts")
# method = c("lts", "lqs", "lms", "S", "model.frame")
xvt3
confint(xvt3)
xvt4 = rlm(xt~xV+I(xV^2), data=data, weights=NULL, method = "M", wt.method = "inv.var")
# method = c("M", "MM", "model.frame"),
# wt.method = c("inv.var", "case")
xvt4
confint(xvt4)
```
Při pokusné filtraci byla sledována závislost objemu filtrátu (v cm3) na čase (min). Určete konstanty (V0 a C) rovnice filtrace.
V^2 + V0*V - C*t = 0
```{r}
# V^2 + V0*V - C*t = 0
# a = 1/C, b = 2*V0/C
# V^2 + b*V - t = 0
xV = seq(20,200,20)
xt = c(0.53, 1.53, 2.63, 4.40, 6.03, 8.11, 10.33, 13.13, 16.03, 18.94)
data = as.data.frame(cbind(xV,xt))
plot(xV,xt)
xvt = lm(xt~xV+I(xV^2))
summary(xvt)
confint(xvt)
koef = coef(xvt)
C = 1/(koef[3]); C
V0 = koef[2]*C/2; V0
library(investr)
plotFit(xvt, data=data, interval="both",pch = 19, shade = TRUE, col.conf = "skyblue4", col.pred = "lightskyblue2")
library(car)
residualPlots(xvt, terms = ~., layout = NULL, main = "", fitted = TRUE, AsIs=TRUE, plot = TRUE,tests = TRUE)
influencePlot(xvt)
infIndexPlot(xvt)
library(s20x)
normcheck(xvt, xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = TRUE, whichPlot = 1:2, usePar = TRUE)
library(s20x)
eovcheck(xvt, xlab = "Fitted values",ylab = "Residuals", col = NULL, smoother = FALSE, twosd = FALSE,levene = TRUE)
library(s20x)
modcheck(xvt, plotOrder = 1:4, args = list(eovcheck = list(smoother = FALSE, twosd = FALSE, levene = TRUE), normcheck = list(xlab = c("Theoretical Quantiles", ""), ylab = c("Sample Quantiles",""), main = c("", ""), col = "light blue", bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = TRUE, 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)))
xre = resid(xvt)
xtp = predict(xvt)
cbind(xV,xt,round(xtp,2))
library(car)
# Durbin-Watson test
durbinWatsonTest(xvt)
# Score Test for Non-Constant Error Variance
library(car)
car::ncvTest(xvt)
# Shapiro-Wilk test rezidui
shapiro.test(resid(xvt))
shapiro.test(rstandard(xvt))
shapiro.test(rstudent(xvt))
# Bonferroni Outlier Test
library(car)
ou = outlierTest(xvt, cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))
# Lund test
library(rapportools)
rp.outlier(rstandard(xvt))
library(fdm2id)
cookplot(xvt)
# Robustni regrese
library(robustbase)
xvt2 = lmrob(xt~xV+I(xV^2),data=data, na.action=na.omit,weights=NULL)
summary(xvt2)
confint(xvt2)
koef2 = coef(xvt2)
C = 1/(koef2[3]); C
V0 = koef2[2]*C/2; V0
library(MASS)
xvt3 = lqs(xt~xV+I(xV^2), data=data, method = "lts")
# method = c("lts", "lqs", "lms", "S", "model.frame")
xvt3
confint(xvt3)
koef3 = coef(xvt3)
C = 1/(koef3[3]); C
V0 = koef3[2]*C/2; V0
xvt4 = rlm(xt~xV+I(xV^2), data=data, weights=NULL, method = "M", wt.method = "inv.var")
# method = c("M", "MM", "model.frame"),
# wt.method = c("inv.var", "case")
xvt4
confint(xvt4)
koef4 = xvt4$coefficients
C = 1/(koef4[3]); C
V0 = koef4[2]*C/2; V0
```
Vypočítejte molární absorpční koeficient pro stanovení Fe3+ 1,10-fenanthrolinem v kyanidovém prostředí při vlnové délce 598 nm, v kyvetě optické délky 1 cm.
A = e * b * c
```{r}
cFe = c(1e-5, 2e-5, 3e-5, 5e-5, 7e-5, 9e-5) # mol/l
Ab = c(0.102, 0.201, 0.302, 0.502, 0.703, 0.902)
data = as.data.frame(cbind(cFe, Ab))
model1 = lm(Ab~cFe)
summary(model1)
model2 = lm(Ab~cFe-1)
summary(model2)
AIC(model1,model2)
library(AICcmodavg)
model.set <- list(model1, model2)
model.names <- c("s usekem", "bez useku")
aictab(model.set, modnames = model.names)
library(performance)
compare_performance(model1,model2, rank = TRUE)
library(EnvStats)
cal = calibrate(Ab ~ cFe, data=data, 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)
cal
detectionLimitCalibrate(cal, coverage = 0.99, simultaneous = TRUE)
library(chemCal)
calplot(model2,xlim = c("auto", "auto"),ylim = c("auto", "auto"),xlab = "Concentration",ylab = "Response",legend_x = "auto",alpha = 0.05,varfunc = NULL)
library(investr)
plotFit(model2, data=data,interval = "both", pch = 19, shade = TRUE, col.conf = "skyblue4", col.pred = "lightskyblue2")
## Prediction of x with confidence interval
abn = 0.458
prediction <- inverse.predict(model2, abn, alpha = 0.05)
prediction
## Critical value (decision limit)
crit <- lod(model2, alpha = 0.05, beta = 0.05)
lod(model2)
lod(model2, alpha = 0.05, beta = 0.05)
loq(model2)
loq(model2, n = 3) # better by using replicate measurements
```
Linearni interpolace
```{r}
x1 <- c(0, 5)
y1 <- c(0, 10)
data_approx1 <- approx(x1, y1)
data_approx1
plot(data_approx1$x, data_approx1$y)
points(x1, y1, col = "black", pch = 16)
x2 <- c(0, 5, 10, 15)
y2 <- c(0, 10, 100, 1000)
data_approx2 <- approx(x2, y2)
data_approx2
plot(x2, y2,type="b")
points(data_approx2$x,data_approx2$y,pch = 19,col = "red")
my_approxfun <- approxfun(x2, y2)
plot(x2, y2)
curve(my_approxfun,add = TRUE)
df <- data.frame(x=c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20),
y=c(4, 7, 11, 16, 22, 29, 38, 49, 63, 80))
plot(df$x, df$y,type="b")
#interpolate y value based on x value of 13
y_new = approx(df$x, df$y, xout=13)
y_new
y_new$x
y_new$y
points(13, y_new$y, col='red', pch=19)
```
Interpolace splajnem
```{r}
x2 <- c(0, 5, 10, 15)
y2 <- c(0, 10, 100, 1000)
data_spl2 <- spline(x2, y2)
data_spl2
plot(data_spl2$x,data_spl2$y)
points(x2, y2, col = "black",pch = 16)
my_splfun <- splinefun(x2, y2)
plot(x2, y2)
curve(my_splfun,add = TRUE)
df <- data.frame(x=c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20),
y=c(4, 7, 11, 16, 22, 29, 38, 49, 63, 80))
plot(df$x, df$y,type="b")
#interpolate y value based on x value of 13
y_new = spline(df$x, df$y, xout=13)
y_new
y_new$x
y_new$y
points(13, y_new$y, col='red', pch=19)
library(pracma)
y_new = cubicspline(df$x, df$y, xi = 13, endp2nd = FALSE, der = c(0, 0))
y_new
```
Interpolace polynomem
```{r}
df <- data.frame(x=c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20),
y=c(4, 7, 11, 16, 22, 29, 38, 49, 63, 80))
plot(df$x, df$y,type="b")
library(pracma) # Neville
neville(df$x, df$y, xs = 13) # interpolovana hodnota
library(polynom) # Lagrange
fp = poly.calc(df$x, df$y)
fp # copy+paste
fpol <- function(x) {
return(44 - 57.38552*x + 32.35672*x^2 - 9.575034*x^3 + 1.701248*x^4 - 0.1892397*x^5 +
0.01327582*x^6 - 0.0005699198*x^7 + 1.366025e-05*x^8 - 1.399395e-07*x^9 )
}
fpol(13) # interpolovana hodnota
plot(df$x, df$y,type="p")
curve(fpol,add = TRUE,col="gray")
```
Interpolace splajnem
```{r}
df <- data.frame(x=c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20),
y=c(4, 7, 11, 16, 22, 29, 38, 49, 63, 80))
plot(df$x, df$y,type="b")
library(pracma)
interp1(df$x, df$y, xi = 13, method = "cubic")
# method = c("linear", "constant", "nearest", "spline","cubic")
# linear: linear interpolation
# constant: constant between points
# nearest: nearest neighbor interpolation
# spline: cubic spline interpolation
# cubic: cubic Hermite interpolation
```
Vyhlazování (smoothing)
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
lines(lowess_values, col = "gray")
lines(lowess(Minutes, Temperature, f = 0.1), col = "red")
library(ggplot2)
ggplot(beaver, aes(x=Minutes, y=Temperature)) +
geom_point() +
geom_smooth(se = FALSE, method = "loess", formula = y ~ x, color = "gray")
ggplot(beaver, aes(x=Minutes, y=Temperature)) +
geom_point() +
geom_smooth(se = FALSE, method = "loess", span = 0.1, formula = y ~ x, color = "red")
```
Loess regrese
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
loessMod <- loess(Temperature ~ Minutes, data=beaver)
loessMod10 <- loess(Temperature ~ Minutes, data=beaver, span=0.10) # 10% smoothing span
loessMod25 <- loess(Temperature ~ Minutes, data=beaver, span=0.25) # 25% smoothing span
smoothed <- predict(loessMod)
smoothed10 <- predict(loessMod10)
smoothed25 <- predict(loessMod25)
plot(beaver$Temperature, x=beaver$Minutes, type="p", main="Loess Smoothing and Prediction", xlab="Date", ylab="Unemployment (Median)")
lines(smoothed, x=beaver$Minutes, col="gray")
lines(smoothed10, x=beaver$Minutes, col="red")
lines(smoothed25, x=beaver$Minutes, col="green")
predict(loessMod10, newdata = 600, se = FALSE,na.action = na.pass)
```
Splajny
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
spl_mod <- smooth.spline(Minutes, Temperature, cv= TRUE)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
lines(spl_mod, col = "gray")
predict(spl_mod, x = 600, deriv = 0)
```
Friedman super smoother
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
supsm = supsmu(x=Minutes, y=Temperature, span = "cv", periodic = FALSE, bass = 0, trace = FALSE)
supsm10 = supsmu(x=Minutes, y=Temperature, span = 0.1, periodic = FALSE, bass = 0, trace = FALSE)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
lines(supsm, col = "gray")
lines(supsm10, col = "red")
```
Moving averages smoothing
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
k=3
smva = filter(Temperature, rep(1/k,k))
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
lines(x= Minutes, y=smva, col = "gray")
library(data.table)
frollmean(Temperature, 3)
library(pracma)
movavg(Temperature, n = 3, type="s")
# type=c("s", "t", "w", "m", "e", "r")
```
Tukey's (Running Median) Smoothing
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
smm = smooth(x=Temperature, kind = "3RS3R", twiceit = FALSE, endrule = "Tukey", do.ends = FALSE)
# kind = c("3RS3R", "3RSS", "3RSR", "3R", "3", "S")
# endrule = c("Tukey", "copy")
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
lines(x= Minutes, y=smm, col = "gray")
```
Savitzky - Golay
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
library(pracma)
ssg = savgol(Temperature, fl = 5, forder = 4, dorder = 0)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
lines(x= Minutes, y=ssg, col = "gray")
```
Hampel Filter (MAD outliers)
```{r}
data(beavers)
head(beaver1)
Minutes <- seq(10, nrow(beaver1) * 10, 10) # extract minutes
Temperature <- beaver1$temp
beaver = as.data.frame(cbind(Minutes,Temperature))
lowess_values <- lowess(Minutes, Temperature)
plot(Minutes, Temperature, type = "p", main = "Body Temperature of Beavers Over Time")
library(pracma)
hampel(Temperature, k=2, t0 = 3)
```
Derivace
```{r}
library(titrationCurves)
# titration curve for a monoprotic weak base analyte using amonoprotic strong acid as the titrant.
wb1 = wb_sa(pka = 8)
head(wb1)
derivative(wb1, plot = TRUE)
tk_splfun <- splinefun(wb1$volume, wb1$ph)
plot(wb1$volume, wb1$ph)
curve(tk_splfun,add = TRUE, col=2)
library(numDeriv) # pouze 1. derivace
diftk1 = grad(tk_splfun, wb1$volume)
plot(wb1$volume, diftk1, type = "l")
wb1$volume[which.min(diftk1)]
library(pracma)
diftk1 = fderiv(f=tk_splfun, x=wb1$volume, n = 1, h = 0,method = "central")
# method = c("central", "forward", "backward")
plot(wb1$volume, diftk1, type = "l")
wb1$volume[which.min(diftk1)]
```
Integrace
```{r}
library(pracma)
data(titanium)
head(titanium)
wt_splfun <- splinefun(titanium$x, titanium$y)
plot(titanium$x, titanium$y)
curve(wt_splfun,add = TRUE, col=2)
integrate(wt_splfun,700,1000)
library(pracma)
trapzfun(wt_splfun,700,1000, maxit = 25, tol = 1e-07)
simpadpt(wt_splfun,700,1000, tol = 1e-6)
romberg(wt_splfun,700,1000, maxit = 25, tol = 1e-12)
library(cmna)
simp(wt_splfun,700,1000, m = 100)
adaptint(wt_splfun,700,1000, n = 10, tol = 1e-06)
# trapezoid
library(pracma)
trapz(titanium$x, titanium$y)
# Simpson
library(Bolstad2)
sintegral(x=titanium$x, fx=titanium$y, n.pts = 256)
library(PROcess) # bioc
# součet součinů průměru sousedních hodnot y a rozdílu sousedních hodnot x
intg(titanium$y, titanium$x)
plot(titanium$x, titanium$y)
nbl = which(titanium$x<700 | titanium$x>1000)
lws <- loess(titanium$y[nbl]~titanium$x[nbl])
sps <- smooth.spline(titanium$x[nbl], titanium$y[nbl], cv= TRUE)
lines(lws$x,lws$y,col=4)
lines(sps$x,sps$y,col=6)
library(pracma)
trapz(titanium$x,(titanium$y-predict(lws,titanium$x)))
trapz(titanium$x,(titanium$y-predict(sps,titanium$x)$y))
library(OrgMassSpecR)
DrawChromatogram(titanium$x, titanium$y,range = list(start = 700, stop = 1000),xlab = "Temperature (°C)")
```
Piky
```{r}
library(VPdtw)
data(reference)
plot(reference,type="l")
chrom = as.data.frame(cbind(as.numeric(names(reference)),reference))
colnames(chrom)= c("time","intensity")
summary(chrom)
chrom = na.omit(chrom)
plot(chrom[,1],chrom[,2],type="l")
library(IDPmisc)
peaks(x=chrom[,1],y=chrom[,2], minPH=1e9)
library(pracma)
findpeaks(chrom[,2], nups = 5, ndowns = 5, zero = "0", peakpat = NULL, minpeakheight = -Inf, minpeakdistance = 1, threshold = 0, npeaks = 0, sortstr = FALSE)
library(quantmod)
fp = findPeaks(chrom[,2], thresh=1e8)
chrom[,2][fp]
chrom[,1][fp]
findValleys(chrom[,2], thresh=1e9)
library(splus2R)
fp = which(peaks(chrom[,2], span=25, strict=TRUE, endbehavior=0))
chrom[,2][fp]
chrom[,1][fp]
```