Response feature analyses

One possible line of attack is to take a summary measure for each dish, e.g. the total number of germinated seeds. Taking a single value for each dish brings us back to more common methods of data analysis: for example, we can fit some sort of GLM to test the significance of effects (storage, temperature and their interaction), within a fully factorial design.

Although the above method is not wrong, undoubtedly, it may be sub-optimal. Indeed, dishes may contain the same total number of germinated seeds, but, nonetheless, they may differ for some other germination traits, such as velocity or uniformity. Indeed, we do not want to express a judgment about one specific characteristic of the seed lot, we would like to express a judgment about the whole seed lot. In other words, we are not specifically asking: “do the seed lots differ for their germination capability?”. We are, more generally, asking “are the seed lots different?”.

In order to get a general assessment, a different method of analysis should be sought, which considers the entire time series (from 1 to 15 days) and not only one single summary measure. This method exists and it is available within the time-to-event platform, which has shown very useful and appropriate for seed germination studies (Onofri et al., 2011; Ritz et al., 2013; Onofri et al., 2019).

The germination time-course

It is necessary to re-organise the dataset in a more useful way. A good format can be obtained by using the ‘makeDrm()’ function in the ‘drcSeedGerm’ package, which can be installed from gitHub (see the code at: this link). The function needs to receive a dataframe storing the counts (dataset[,4:18]), a dataframe storing the factor variables (dataset[,2:3]), a vector with the number of seeds in each Petri dish (rep(25, 16)) and a vector of monitoring times.

library(drcSeedGerm)
library(dplyr)
datasetR <- makeDrm(dataset[,4:18], dataset[,2:3], rep(25, 16), 1:15)
head(datasetR, 16)
##      Storage Temp Dish timeBef timeAf count nCum propCum
## 1        Low  Fix    1       0      1     0    0    0.00
## 1.1      Low  Fix    1       1      2     0    0    0.00
## 1.2      Low  Fix    1       2      3     0    0    0.00
## 1.3      Low  Fix    1       3      4     0    0    0.00
## 1.4      Low  Fix    1       4      5     0    0    0.00
## 1.5      Low  Fix    1       5      6     0    0    0.00
## 1.6      Low  Fix    1       6      7     0    0    0.00
## 1.7      Low  Fix    1       7      8     0    0    0.00
## 1.8      Low  Fix    1       8      9     3    3    0.12
## 1.9      Low  Fix    1       9     10     4    7    0.28
## 1.10     Low  Fix    1      10     11     6   13    0.52
## 1.11     Low  Fix    1      11     12     0   13    0.52
## 1.12     Low  Fix    1      12     13     1   14    0.56
## 1.13     Low  Fix    1      13     14     0   14    0.56
## 1.14     Low  Fix    1      14     15     3   17    0.68
## 1.15     Low  Fix    1      15    Inf     8   NA      NA
datasetR <- datasetR %>% 
  mutate(across(c(Storage, Temp, Dish), .fns = factor))

The snippet above shows the first dish. Roughly speaking, we have gone from a WIDE format to a LONG format. The column ‘timeAf’ contains the time when the inspection was made and the column ‘count’ contains the number of germinated seeds (e.g. 9 seeds were counted at day 9). These seeds did not germinate exactly at day 9; they germinated within the interval between two inspections, that is between day 8 and day 9. The beginning of the interval is given in the variable ‘timeBef’. Apart from these columns, we have additional columns, which we are not going to use for our analyses. The cumulative counts of germinated seeds are in the column ‘nCum’; these cumulative counts have been converted into cumulative proportions by dividing by 25 (i.e., the total number of seeds in a dish; see the column ‘propCum’).

We can use a time-to-event model to parameterise the germination time-course for this dish. This is easily done by using the ‘drm()’ function in the ‘drc’ package (Ritz et al., 2013):

modPre <- drm(count ~ timeBef + timeAf, fct = LL.3(), 
              data = datasetR, 
              type = "event", subset = c(Dish == 1))
plot(modPre, log = "", xlab = "Time", 
     ylab = "Proportion of germinated seeds",
     xlim = c(0, 15))

Please, note the following:

1. we are using the counts (‘count’) as the dependent variable
2. as the independent variable: we are using the extremes of the inspection interval, within which germinations were observed (count ~ timeBef + time Af)
3. we have assumed a log-logistic distribution of germination times (fct = LL.3()). A three parameter model is necessary, because there is a final fraction of ungerminated seeds (truncated distribution)
4. we have set the argument ‘type = “event”’. Indeed, we are fitting a time-to-event model, not a nonlinear regression model, which would be incorrect, in this setting (see this link here ).

As we have determined the germination time-course for dish 1, we can do the same for all dishes. However, we have to instruct ‘drm()’ to define a different curve for each combination of storage and temperature. It is necessary to make an appropriate use of the ‘curveid’ argument. Please, see below.

mod1 <- drm(count ~ timeBef + timeAf, fct = LL.3(),
            data = datasetR, type = "event", 
            curveid = Temp:Storage)
plot(mod1, log = "", legendPos = c(2, 1))

It appears that there are visible differences between the curves (the legend considers the curves in alphabetical order, i.e. 1: Fix-Low, 2: Fix-High, 3: Alt-Low and 4: Alt-High). We can test that the curves are similar by coding a reduced model, where we have only one pooled curve for all treatment levels. It is enough to remove the ‘curveid’ argument.

modNull <- drm(count ~ timeBef + timeAf, fct = LL.3(),
               data = datasetR, 
               type = "event")
anova(mod1, modNull, test = "Chisq")
## 
## 1st model
##  fct:      LL.3()
##  pmodels: Temp:Storage (for all parameters)
## 2nd model
##  fct:      LL.3()
##  pmodels: 1 (for all parameters)
## ANOVA-like table
## 
##           ModelDf  Loglik Df LR value p value
## 1st model     244 -753.54                    
## 2nd model     253 -854.93  9   202.77       0

Now, we can compare the full model (four curves) with the reduced model (one common curve) by using a Likelihood Ratio Test, which is approximately distributed as a Chi-square. The test is highly significant. Of course, we can also test some other hypotheses. For example, we can code a model with different curves for storage times, assuming that the effect of temperature is irrelevant:

mod2 <- drm(count ~ timeBef + timeAf, fct = LL.3(), 
            data = datasetR, type = "event", 
            curveid = Storage)
anova(mod1, mod2, test = "Chisq")
## 
## 1st model
##  fct:      LL.3()
##  pmodels: Temp:Storage (for all parameters)
## 2nd model
##  fct:      LL.3()
##  pmodels: Storage (for all parameters)
## ANOVA-like table
## 
##           ModelDf  Loglik Df LR value p value
## 1st model     244 -753.54                    
## 2nd model     250 -797.26  6   87.436       0

We see that such an assumption (temperature effect is irrelevant) is not supported by the data: the temperature effect cannot be removed without causing a significant decrease in the likelihood of the model. Similarly, we can test the effect of storage:

mod3 <- drm(count ~ timeBef + timeAf, fct = LL.3(), 
            data = datasetR, type = "event", 
            curveid = Temp)
anova(mod1, mod3, test = "Chisq")
## 
## 1st model
##  fct:      LL.3()
##  pmodels: Temp:Storage (for all parameters)
## 2nd model
##  fct:      LL.3()
##  pmodels: Temp (for all parameters)
## ANOVA-like table
## 
##           ModelDf  Loglik Df LR value p value
## 1st model     244 -753.54                    
## 2nd model     250 -849.48  6   191.87       0

Again, we get significant results. So, we need to conclude that temperature and storage time caused a significant influence on the germination behavior for the species under study.

Before concluding, it is necessary to mention that, in general, the above LR tests should be taken with care: the results are only approximate and the observed data are not totally independent, as multiple observations are taken in each experimental unit (Petri dish). In order to restore independence, we would need to add to the model a random effect for the Petri dish, which is not an easy task in a time-to-event framework (Onofri et al., 2019). However, we got very low p-levels, which leave us rather confident about the significance of effects. It may be a good suggestion, in general, to avoid formal hypothesis testing and compare the models by using the Akaike Information Criterion (AIC: the lowest is the best), which confirms that the complete model with four curves is, indeed, the best one.

AIC(mod1, mod2, mod3, modNull)
##          df      AIC
## mod1    244 1995.088
## mod2    250 2094.524
## mod3    250 2198.961
## modNull 253 2215.862

For those who are familiar with linear model parameterisation, it is possible to reach even a higher degree of flexibility by using the ‘pmodels’ argument, within the ‘drm()’ function. However, this will require another post. Thanks for reading!

References

1. Catara, S., Cristaudo, A., Gualtieri, A., Galesi, R., Impelluso, C., Onofri, A., 2016. Threshold temperatures for seed germination in nine species of Verbascum (Scrophulariaceae). Seed Science Research 26, 30–46.
2. Onofri, A., Mesgaran, M.B., Tei, F., Cousens, R.D., 2011. The cure model: an improved way to describe seed germination? Weed Research 51, 516–524.
3. Onofri, A., Piepho, H.-P., Kozak, M., 2019. Analysing censored data in agricultural research: A review with examples and software tips. Annals of Applied Biology 174, 3–13. https://doi.org/10.1111/aab.12477
4. Ritz, C., Pipper, C.B., Streibig, J.C., 2013. Analysis of germination data from agricultural experiments. European Journal of Agronomy 45, 1–6.