LRT for analysing GP

We assume the situation where one treatment factor with I levels exists and N_i repetitions (e.g. Petri dishes, trays, plates, etc.) are done for each treatment level i. We use i for denoting the treatment level and j for the replication. Let n_ij and g_ij be the number of tested seeds and the number of germinated seeds for the jth replication of the ith treatment level, respectively. Finally, p_i is the GP of the ith treatment level. In this situation, we are interested in testing whether there is a significant difference in GP between the treatment levels (H ₀: p ₁ = ⋅⋅⋅ = p_I vs. H ₁: not H ₀) and which levels differ significantly from one another in case there is a difference. Our GP data are (g_ij, n_ij), i = 1, … ,I, j = 1, … ,N_i.

The LRT is a statistical test used to compare the fit of two nested statistical models. A nested model can be obtained by imposing restrictions on the parameters of a more general model. In our case, the model specified under H ₁ is g_ij ~ binomial(n_ij,p_i) for i = 1, … ,I, while the model under H ₀ is g_ij ~ binomial(n_ij,p_i) for i = 1, … ,I, subject to the constraint p ₁ = ⋅⋅⋅ = p_I. Consequently, under H ₀, there is only one parameter. The likelihood function measures the probability of the observed data, given the model and its parameters. Therefore, the likelihood ratio compares the support provided by the data for one model relative to the other.

In our case of analysing GP, the estimated parameters under H ₀ and H ₁ through the maximum likelihood method are

(1)$$\matrix{ {{\hat{\,p}}_{i, H_0} = \displaystyle{{\mathop \sum \nolimits_{i = 1}^I \mathop \sum \nolimits_{\,j = 1}^{N_i} g_{ij}} \over {\mathop \sum \nolimits_{i = 1}^I \mathop \sum \nolimits_{\,j = 1}^{N_i} n_{ij}}}\;\;\;{\rm and}\;\;\;{\hat{\,p}}_{i, H_1} = \displaystyle{{\mathop \sum \nolimits_{\,j = 1}^{N_i} g_{ij}} \over {\mathop \sum \nolimits_{\,j = 1}^{N_i} n_{ij}}}, \;} \cr } $$

respectively. Note that $\hat{p}_{i, H_0}$ does not depend on the treatment level i, but $\hat{p}_{i, H_1}$ does depend. $\hat{p}_{i, H_0}$ is the pooled estimate across all treatment groups. Furthermore, the log-likelihoods of the two models are given as follows:

(2)$$\matrix{ {l( {H_0} ) = l( {{\hat{\,p}}_{i, H_0}; \;{\rm data}} ) = \log \left({\mathop \prod \limits_{i = 1}^I \mathop \prod \limits_{\,j = 1}^{N_i} {\hat{\,p}}_{i, H_0}^{g_{ij}} {( {1-{\hat{\,p}}_{i, H_0}} ) }^{n_{ij}-g_{ij}}} \right)} \cr } $$

(3)$$\matrix{ {l( {H_1} ) = l( {{\hat{\,p}}_{i, H_1}; \;{\rm data}} ) = \log \left({\mathop \prod \limits_{i = 1}^I \mathop \prod \limits_{\,j = 1}^{N_i} {\hat{\,p}}_{i, H_1}^{g_{ij}} {( {1-{\hat{\,p}}_{i, H_1}} ) }^{n_{ij}-g_{ij}}} \right)} \cr } $$

The LRT typically calculates two times the difference of the two log-likelihoods Λ = 2 (l(H ₁) − l(H ₀))and decides which model to support depending on the size of Λ. When Λ is large, we support the model under H ₁ and when Λ is close to zero, we support the model H ₀. Note that Λ ≥ 0 since the model under H ₁ nests the model under H ₀. Thanks to Wilks’ theorem (Wilks, Reference Wilks1938), we know that under H ₀ the test statistic Λ is asymptotically chi-squared distributed with degrees of freedom I − 1, which is the difference in the number of parameters of the two models. Simply speaking, when Λ is greater than the (1 − α)-quantile of the chi-square distribution with degrees of freedom I − 1, then we conclude that there is a significant difference in germination rates at significance level α.

One of the significant benefits of the LRT framework we are exploring is that it permits testing even when the GPs are all zeros or ones at particular treatment levels. Moreover, unlike ANOVA, it considers the difference in confidence in the germination rate that arises from variations in the number of seeds used in each repetition when testing the difference in GP. Figure 1 shows that the LRT test does consider the number of tested seeds utilized to compare GP across levels. Once it is confirmed that the GP is not equal as a whole, a post-hoc test can be conducted to pinpoint the specific treatment levels contributing to the GP difference. This can be achieved by conducting our LRT test pairwisely and adjusting the P-values for multiple comparisons. Alternatively, other tests designed for comparing two proportions can be utilized for pairwise comparisons. Based on our simulation results (not presented in the paper), it appears that the pairwise comparison with the LRT yields slightly better outcomes compared to the two-sample proportion test. Our LRT method, which we have implemented in our developed R package called gerstat, can be executed by utilizing the function ger_test. For more information about our R package, please refer to the Appendix.

Numerical simulation

This section assesses the efficacy of our proposed LRT method in analysing GP data relative to conventional methods under various scenarios. The methods we compare include ANOVA, ANOVA with arcsine transformed response, logistic regression, and Kruskal–Wallis test. Furthermore, we consider the test suggested by Sill (Reference Sill2007) and Schaarschmidt et al. (Reference Schaarschmidt, Sill and Hothorn2008) and will henceforth refer to this method as the MCPAN method. For our simulations, we employ a fully randomized design with three replications at each level and four treatment levels.

Three primary simulation scenarios were taken into consideration: (A) a situation where specific levels exhibit seed dormancy, (B) a scenario where there is no seed dormancy, but a difference in GP exists among the levels and (C) a situation where there is no seed dormancy, and there is no difference in GP. For Scenario A, the first two treatment levels exhibit seed dormancy (p ₁ = p ₂ = 0), but the third and fourth treatment levels have GP of p ₃ = 0.15 and p ₄ = 0.5, respectively. For Scenario B, we have (p ₁,p ₂,p ₃,p ₄) = (0.1,0.2,0.3,0.4). Finally, we have the same GP (p ₁ = p ₂ = p ₃ = p ₄ = 0.3) all across the levels. We generate all the data from the model g_ij ~ binomial(n_ij,p_i), i = 1, … ,4, with n_ij being 15 or 30 (the number of seeds per replication). The number of simulations is 1000 and the significance level α is 0.05.

Once it is established that the GP values are not uniformly equal, a post-hoc test is performed to identify the specific treatment levels contributing to the GP difference. This post-hoc test entails conducting multiple pairwise comparison tests. To control the overall rate of false positives, multiple testing methods are implemented during the post-hoc test to control either the false positive rate (FPR) or the family-wise error rate (FWER). As multi-step methods are more efficient in controlling the FPR or the FWER than single-step methods, we utilized the Holm procedure, a multi-step method, for the entire simulation. To conduct pairwise comparison tests in both ANOVA and logistic regression, we utilized the glht function from the R package multcomp. This package streamlines pairwise comparisons across diverse models by furnishing suitable inference procedures for testing general linear hypothesis tests. In the context of Kruskal–Wallis, we conducted pairwise comparisons similar to Dunn's test (Dunn, Reference Dunn1964), a frequently employed post-hoc test in nonparametric scenarios. As for the LRT method, we executed the pairwise comparison as outlined in the section ‘LRT for analysing GP’. To apply the MCPAN method, we conducted pairwise comparisons and adjusted P-values following the guidelines provided by Sill (Reference Sill2007). The MCPAN method was implemented using the function binomRDtest in the R package MCPAN.

In our comparison of the methodologies under consideration, we evaluate their performance on two key dimensions. Firstly, we calculate the proportion of overall tests that have detected differences between groups, if any. Secondly, we assess the accuracy of the follow-up tests in detecting treatment pairs with distinct GPs. Although various GP vectors could be examined for a given scenario, such simulations (which are not included in the paper) did not yield significantly different conclusions compared to those presented in the paper.

Scenario A ((p ₁,p ₂,p ₃,p ₄) = (0,0,0.15,0.5))

In this subsection, we present the results obtained when certain treatment levels exhibit seed dormancy. As shown in Table 1, the performance differences between the methodologies are insignificant when it comes to testing for differences in GP. However, as highlighted in Table 2, the methods exhibit noticeable performance differences in identifying the specific treatment levels that contribute to the difference in GP. Logistic regression struggles with detecting GP differences due to the issue of complete separation (Mansournia et al., Reference Mansournia, Geroldinger, Greenland and Heinze2017). On the other hand, the Kruskal–Wallis test has a lower power, which is a characteristic of nonparametric methods, and can only detect pairs with relatively large differences in germination rates (p ₁ = 0 vs. p ₄ = 0.5 or p ₂ = 0 vs. p ₄ = 0.5). In most cases, the proposed LRT demonstrates superior performance compared to ANOVA, arcsine-transformed response-based ANOVA and the MCPAN method in detecting treatment pairs with distinct germination rates.

Table 1. Detection rate of whether there is a difference in germination rates when (p ₁,p ₂,p ₃,p ₄) = (0,0,0.15,0.5) (Scenario A)

Table 2 Percentage that identified a difference in germination rate for a given treatment pair when (p ₁,p ₂,p ₃,p ₄) = (0,0,0.15,0.5) (Scenario A)

Note.The treatment pair (p₁ vs. p₂) has no difference in germination rate.

Scenario B ((p ₁,p ₂,p ₃,p ₄) = (0.1,0.2,0.3,0.4))

In this subsection, we investigate the scenario where there is no seed dormancy, but a difference in GP exists among the treatment levels. Table 3 shows that the proposed LRT, logistic regression and the MCPAN method outperform the other methods in detecting whether there are differences in germination rates. Concerning the identification of specific treatment levels contributing to the difference in GP (Table 4), two primary trends are evident. Firstly, the Kruskal–Wallis test continues to exhibit low detection power. Secondly, the proposed LRT, logistic regression and the MCPAN method consistently outperform ANOVA and acsine-transformed response-based ANOVA in the majority of cases. Notably, the proposed LRT demonstrates the best performance.

Table 3 Detection rate of whether there is a difference in germination rates when (p ₁,p ₂,p ₃,p ₄) = (0.1,0.2,0.3,0.4) (Scenario B)

Table 4 Percentage that identified a difference in germination rate for a given treatment pair when (p ₁,p ₂,p ₃,p ₄) = (0.1,0.2,0.3,0.4) (Scenario B)

Scenario C ((p ₁,p ₂,p ₃,p ₄) = (0.3,0.3,0.3,0.3))

In Scenario C, there are no differences in germination rates among the treatment levels. Consequently, it is crucial to control the power of the overall test close to the significance level. Similarly, when employing multiple testing methods, it is important to ensure that the percentage of detected differences in each treatment pair is controlled at a value near the significance level. Upon analysing the findings presented in Tables 5 and 6, we observe that all the methods exhibit these characteristics. However, it is notable that the Kruskal–Wallis test proves to be conservative, with a detection rate lower than the significance level.

Table 5 Detection rate of whether there is a difference in germination rates when (p ₁,p ₂,p ₃,p ₄) = (0.1,0.2,0.3,0.4) (Scenario C)

Table 6 Percentage that identified a difference in germination rate for a given treatment pair when (p ₁,p ₂,p ₃,p ₄) = (0.3,0.3,0.3,0.3) (Scenario C)

An important observation across all scenarios is that both logistic regression and the proposed LRT yield identical results concerning the test for overall differences. This alignment arises from the equality between the deviation in deviance calculated between the full model and the null model (the model with only the intercept) used by logistic regression and our proposed LRT statistic. Despite the potential instability in coefficient estimates when there are consistent zeros or ones across all replicates in logistic regression, it does not adversely affect the likelihood values of the estimated model. Consequently, the detection rate of overall differences in logistic regression remains robust in situations with perfect separation. However, concerning pairwise comparisons, reliance on coefficient estimates significantly hampers testing performance in situations of perfect separation.

Reanalysis of the data from Rhie et al. (Reference Rhie, Kim, Lee and Kim2016)

To demonstrate the advantages of our method, we conducted a reanalysis of the data from Rhie et al. (Reference Rhie, Kim, Lee and Kim2016)'s study that examines the effects of GA₃ on germination under low temperatures. Figure 2 presents the raw data and the corresponding testing results. In Rhie et al. (Reference Rhie, Kim, Lee and Kim2016)'s original analysis, ANOVA was employed to test the overall difference, followed by the utilization of Tukey's honestly significant difference (HSD) test as a post-hoc test with a specific significance level α = 0.05. In our reanalysis, we applied the LRT described in the section ‘Numerical simulation’ with the same significance level α = 0.05. It is important to note that both the Holm–Bonferroni correction, which we used for our LRT analysis, and Tukey's HSD method for ANOVA are designed to control the family-wise error rate. Upon comparing the two analyses, it was determined that both reached the same conclusion that there is a difference in GP. However, there were variations in the results for the treatment pair that exhibited a difference in GP. In Rhie et al. (Reference Rhie, Kim, Lee and Kim2016)'s analysis, it was observed that there was no significant difference in GP between the 100 and 10 mg L⁻¹ doses of GA₃. However, it is challenging to accept the conclusion that there was no significant difference in germination between the two treatments. This is due to the fact that no germination occurred at the 10 mg L⁻¹ dose of GA₃, while certain replications at the 100 mg L⁻¹ dose showed high germination rates of up to 43%. The substantial disparity in germination outcomes between the two treatments raises doubts about the conclusion of no significant difference. In contrast, the LRT-based analysis classifies the 0 and 10 mg L⁻¹ treatments, which displayed no germination, into the same germination rate group (a). Furthermore, it distinguishes the 100 and 1000 mg L⁻¹ treatments, which exhibited a significant difference in average germination, into distinct germination rate groups (b and c). Considering our reanalysis of Rhie et al. (Reference Rhie, Kim, Lee and Kim2016)'s data and the findings from our simulations, where the ANOVA method occasionally failed to detect differences between treatments with zero germination and those without, we can conclude that the proposed LRT method can provide more reliable conclusions than traditional methods in various GP data analysis situations.

Figure 2. Reanalysis of the data from Rhie et al. (Reference Rhie, Kim, Lee and Kim2016)'s study. (*** means that the germination rate difference is significant at P < 0.001. The use of different letters to mark the averages signifies that the corresponding treatments have been identified as having significantly different germination rates through multiple comparisons.)