In certain situations, such as when a test chemical is anticipated to have low toxicity or is poorly soluble, a limit test may be conducted. This test aims to show that the No Observed Effect Concentration (NOEC) is greater than or equal to the tested limit concentration, and that the LD50 or EC50 is also greater than the tested limit dose, provided no effects are observed during the study.
Below is an example for Limit Test analysis.
data("test_cases_data")
testdata <- test_cases_data%>% filter(Design =="Limit")
metainfo <- testdata[1,]
design <- ifelse(metainfo$Design=="Limit","limit test", "full dose response study")
nconc <- length(unique(testdata$Dose))The test is a limit test for test organism Myriophyllum. There are 2 test concentrations. The interested endpoint is Growth Rate for Total shoot length at 14 d.
Preliminary Assessment
ctr <- testdata %>% filter(Dose == "0")
ctr0 <- mean(ctr$Response)
sres <- testdata %>% group_by(Dose)%>% dplyr::summarize(Mean=mean(Response),SD=sd(Response)) %>% mutate(`% Inhibition` = - ((Mean-ctr0)/ctr0)*100)
sres %>% knitr::kable(.,digits = 3)| Dose | Mean | SD | % Inhibition |
|---|---|---|---|
| 0 | 0.136 | 0.010 | 0.000 |
| 100 | 0.126 | 0.005 | 7.361 |
Outlier Check
dixon.test <- outliers::dixon.test
DixonQ[DixonQ$n==6,"Q_critical"]
#> [1] 0.56
outRes <- testdata %>% group_by(Dose) %>% nest() %>% mutate(normtest=map(data,~dixon.test(.x$Response)),
tidies= map(normtest,broom::tidy))
outRes <- outRes %>% dplyr::select(-c(data,normtest)) %>% unnest(c(tidies))
outRes %>% knitr::kable(.,digits = 3)| Dose | statistic | p.value | method | alternative |
|---|---|---|---|---|
| 0 | 0.377 | 0.423 | Dixon test for outliers | lowest value 0.120915680098082 is an outlier |
| 100 | 0.342 | 0.524 | Dixon test for outliers | lowest value 0.117455425985384 is an outlier |
## by specifying the opposite: a logical indicating whether you want to check not the value with largest difference from the mean, but opposite (lowest, if most suspicious is highest etc.)
outRes <- testdata %>% group_by(Dose) %>% nest() %>% mutate(normtest=map(data,~dixon.test(.x$Response,opposite = TRUE)),
tidies= map(normtest,broom::tidy))
outRes <- outRes %>% dplyr::select(-c(data,normtest)) %>% unnest(c(tidies))
outRes %>% knitr::kable(.,digits = 3)| Dose | statistic | p.value | method | alternative |
|---|---|---|---|---|
| 0 | 0.146 | 0.716 | Dixon test for outliers | highest value 0.149892124046466 is an outlier |
| 100 | 0.048 | 0.252 | Dixon test for outliers | highest value 0.131517109035102 is an outlier |
Normality Check
Welch’s t-test only requires normality for each group.
Normality for residuals
mod <- lm(Response~factor(Dose),data=testdata)
normalres <- shapiro.test(residuals(mod))Normality for each group
normRes <- testdata %>% group_by(Dose) %>% nest() %>% mutate(normtest=map(data,~shapiro.test(.x$Response)),
tidies= map(normtest,broom::tidy))
normRes <- normRes %>% dplyr::select(-c(data,normtest)) %>% unnest(c(tidies))
normRes %>% knitr::kable(.,digits = 3)| Dose | statistic | p.value | method |
|---|---|---|---|
| 0 | 0.968 | 0.880 | Shapiro-Wilk normality test |
| 100 | 0.892 | 0.326 | Shapiro-Wilk normality test |
Variance Homogeneity
(varTest <- car::leveneTest(Response~factor(Dose),data=testdata,center=mean))
#> Levene's Test for Homogeneity of Variance (center = mean)
#> Df F value Pr(>F)
#> group 1 1.4879 0.2505
#> 10
## lawstat::levene.test(testdata$Response,testdata$Dose,location="mean")
car::leveneTest(Response~factor(Dose),data=testdata,center=median)
#> Levene's Test for Homogeneity of Variance (center = median)
#> Df F value Pr(>F)
#> group 1 1.5126 0.2469
#> 10Student’s t-test
two-sided
(res <- t.test(testdata$Response~testdata$Dose,var.equal = TRUE))
#>
#> Two Sample t-test
#>
#> data: testdata$Response by testdata$Dose
#> t = 2.0993, df = 10, p-value = 0.06215
#> alternative hypothesis: true difference in means between group 0 and group 100 is not equal to 0
#> 95 percent confidence interval:
#> -0.0006161316 0.0207027221
#> sample estimates:
#> mean in group 0 mean in group 100
#> 0.1364410 0.1263977
tres <- broom::tidy(res)
tres1 <- c(tres$estimate1,tres$estimate2, tres$parameter,sres$`% Inhibition`[2], -tres$statistic,tres$p.value)
names(tres1) <- c("Mean Control", "Mean Test Item", "df", "%Inhibition", "T-statistic", "p-value")
tres1
#> Mean Control Mean Test Item df %Inhibition T-statistic
#> 0.13644102 0.12639773 10.00000000 7.36090584 -2.09934893
#> p-value
#> 0.06214521one-sided greater
Note that if using formula, the greater means test group is greater than control
(res <- t.test(testdata$Response~testdata$Dose,var.equal = TRUE,alternative="greater"))
#>
#> Two Sample t-test
#>
#> data: testdata$Response by testdata$Dose
#> t = 2.0993, df = 10, p-value = 0.03107
#> alternative hypothesis: true difference in means between group 0 and group 100 is greater than 0
#> 95 percent confidence interval:
#> 0.001372473 Inf
#> sample estimates:
#> mean in group 0 mean in group 100
#> 0.1364410 0.1263977
tres <- broom::tidy(res)
tres1 <- c(tres$estimate1,tres$estimate2, tres$parameter,sres$`% Inhibition`[2], -tres$statistic,tres$p.value)
names(tres1) <- c("Mean Control", "Mean Test Item", "df", "%Inhibition", "T-statistic", "p-value")
tres1
#> Mean Control Mean Test Item df %Inhibition T-statistic
#> 0.13644102 0.12639773 10.00000000 7.36090584 -2.09934893
#> p-value
#> 0.03107261one-sided smaller
Note that if using formula, the greater means test group is greater than control
(res <- t.test(testdata$Response~testdata$Dose,var.equal = TRUE,alternative="less"))
#>
#> Two Sample t-test
#>
#> data: testdata$Response by testdata$Dose
#> t = 2.0993, df = 10, p-value = 0.9689
#> alternative hypothesis: true difference in means between group 0 and group 100 is less than 0
#> 95 percent confidence interval:
#> -Inf 0.01871412
#> sample estimates:
#> mean in group 0 mean in group 100
#> 0.1364410 0.1263977
tres <- broom::tidy(res)
tres1 <- c(tres$estimate1,tres$estimate2, tres$parameter,sres$`% Inhibition`[2], -tres$statistic,tres$p.value)
names(tres1) <- c("Mean Control", "Mean Test Item", "df", "%Inhibition", "T-statistic", "p-value")
tres1
#> Mean Control Mean Test Item df %Inhibition T-statistic
#> 0.1364410 0.1263977 10.0000000 7.3609058 -2.0993489
#> p-value
#> 0.9689274Welch’s t-test
Welch Two Sample t-test
(res <- t.test(testdata$Response~testdata$Dose))
#>
#> Welch Two Sample t-test
#>
#> data: testdata$Response by testdata$Dose
#> t = 2.0993, df = 7.5933, p-value = 0.07085
#> alternative hypothesis: true difference in means between group 0 and group 100 is not equal to 0
#> 95 percent confidence interval:
#> -0.001092339 0.021178930
#> sample estimates:
#> mean in group 0 mean in group 100
#> 0.1364410 0.1263977
tres <- broom::tidy(res)
tres1 <- c(tres$estimate1,tres$estimate2, tres$parameter,sres$`% Inhibition`[2], -tres$statistic,tres$p.value)
names(tres1) <- c("Mean Control", "Mean Test Item", "df", "%Inhibition", "T-statistic", "p-value")
tres1
#> Mean Control Mean Test Item df %Inhibition T-statistic
#> 0.13644102 0.12639773 7.59330452 7.36090584 -2.09934893
#> p-value
#> 0.07085487one-sided greagter
(res <- t.test(testdata$Response~testdata$Dose,alternative ="less"))
#>
#> Welch Two Sample t-test
#>
#> data: testdata$Response by testdata$Dose
#> t = 2.0993, df = 7.5933, p-value = 0.9646
#> alternative hypothesis: true difference in means between group 0 and group 100 is less than 0
#> 95 percent confidence interval:
#> -Inf 0.01900155
#> sample estimates:
#> mean in group 0 mean in group 100
#> 0.1364410 0.1263977
tres <- broom::tidy(res)
tres1 <- c(tres$estimate1,tres$estimate2, tres$parameter,sres$`% Inhibition`[2], -tres$statistic,tres$p.value)
names(tres1) <- c("Mean Control", "Mean Test Item", "df", "%Inhibition", "T-statistic", "p-value")
tres1
#> Mean Control Mean Test Item df %Inhibition T-statistic
#> 0.1364410 0.1263977 7.5933045 7.3609058 -2.0993489
#> p-value
#> 0.9645726one-sided smaller
(res <- t.test(testdata$Response~testdata$Dose,alternative ="greater"))
#>
#> Welch Two Sample t-test
#>
#> data: testdata$Response by testdata$Dose
#> t = 2.0993, df = 7.5933, p-value = 0.03543
#> alternative hypothesis: true difference in means between group 0 and group 100 is greater than 0
#> 95 percent confidence interval:
#> 0.001085044 Inf
#> sample estimates:
#> mean in group 0 mean in group 100
#> 0.1364410 0.1263977
tres <- broom::tidy(res)
tres1 <- c(tres$estimate1,tres$estimate2, tres$parameter,sres$`% Inhibition`[2], -tres$statistic,tres$p.value)
names(tres1) <- c("Mean Control", "Mean Test Item", "df", "%Inhibition", "T-statistic", "p-value")
tres1
#> Mean Control Mean Test Item df %Inhibition T-statistic
#> 0.13644102 0.12639773 7.59330452 7.36090584 -2.09934893
#> p-value
#> 0.03542743Nonparametric test
Kruskal–Wallis test
(kwres <- kruskal.test(Response~Dose,data=testdata))
#>
#> Kruskal-Wallis rank sum test
#>
#> data: Response by Dose
#> Kruskal-Wallis chi-squared = 4.3333, df = 1, p-value = 0.03737Dunn’s test (many-to-one)
In essence, Dunn’s test can be viewed as a series of Wilcoxon rank sum tests applied to pairs of groups, with adjustments made for multiple comparisons.
The link between Dunn’s test and the Wilcoxon rank sum test lies in their non-parametric nature and rank-based approach. While the Wilcoxon rank sum test is used for comparing two groups, Dunn’s test extends the concept of rank comparisons to multiple groups following the Kruskal-Wallis test.
## (dunnres <- PMCMRplus::kwManyOneDunnTest(Response~factor(Dose),data=testdata))Wilcoxon test
Wilcoxon rank-sum test (also called Mann–Whitney U-test) when the data deviates significantly from normal.
(Ures <- wilcox.test(Response~factor(Dose),data=testdata))
#>
#> Wilcoxon rank sum exact test
#>
#> data: Response by factor(Dose)
#> W = 31, p-value = 0.04113
#> alternative hypothesis: true location shift is not equal to 0
## same results would be obtained by specifying x and y instead of a formula
## wilcox.test(testdata$Response[testdata$Dose=="0"],testdata$Response[testdata$Dose=="100"],data=testdata)
(Ures <- wilcox.test(Response~factor(Dose),data=testdata,alternative="greater"))
#>
#> Wilcoxon rank sum exact test
#>
#> data: Response by factor(Dose)
#> W = 31, p-value = 0.02056
#> alternative hypothesis: true location shift is greater than 0
(Ures <- wilcox.test(Response~factor(Dose),data=testdata,alternative="less"))
#>
#> Wilcoxon rank sum exact test
#>
#> data: Response by factor(Dose)
#> W = 31, p-value = 0.987
#> alternative hypothesis: true location shift is less than 0