Section 7 Exploring the results of a time-to-event fit: model parameters
Once we have fit a time-to-event model, we are usually interested in exploring the results, to get all possible information from the fitted model. If we have fitted a parametric model, the value of the estimated parameters is usually of extreme interest, as it gives information on the main traits of germination/emergence (e.g., capability, speed and uniformity).
For example, let’s consider the hydro-time model we have fitted in our previous post at this link:
P(t,Ψ)=Φ{Ψ−(θH/t)−ΨbσΨb}
where P is the cumulative proportion of germinated seeds at time t and water potential Ψ, Φ is a gaussian cumulative distribution function for base water potential, Ψb is the median base water potential in the seed lot (in MPa), θH is the hydro-time constant (in MPa day or MPa hour) and σΨb represents the variability of Ψb within the population. Clearly, these parameters have a clear biological meaning and getting to know about their value represents the reason why we have fitted such a model. The box below shows the code we used in our previous post:
library(drcSeedGerm)
library(drcte)
data(rape)
<- drmte(nSeeds ~ timeBef + timeAf + Psi,
modHTE data = rape, fct = HTnorm())
coef(modHTE)
## ThetaH:(Intercept) Psib50:(Intercept) sigmaPsib:(Intercept)
## 0.7510691 -0.9069810 0.2369954
Indeed, we got parameter estimates, but we are not happy with this. We also need standard errors, to present along with estimates in our papers. The easiest way to obtain parameters and their standard errors altogether is to use the summary()
method for ‘drcte’ objects:
summary(modHTE)
##
## Model fitted: Hydrotime model with normal distribution of Psib (Bradford et al., 2002)
##
## Robust estimation: no
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## ThetaH:(Intercept) 0.7510691 0.0394604 19.034 < 2.2e-16 ***
## Psib50:(Intercept) -0.9069810 0.0118081 -76.810 < 2.2e-16 ***
## sigmaPsib:(Intercept) 0.2369954 0.0071406 33.190 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Unfortunately, the above standard errors are not correct: indeed, they are obtained assuming that the observational units (i.e., the seeds) are independent, while they are clustered within randomisation units (Petri dishes, in this case). Consequently, seeds in the same Petri dish are more similar than seeds in different Petri dishes (there is intra-class correlation, we say). How can we obtain standard errors that account for such lack of independence?
If we look at the literature about survival analysis (that is where we borrowed our methods from), we can see that cluster robust sandwich estimators of standard errors have proven useful and reliable (Yu and Peng, 2008). Therefore, we have implemented them in ‘drcte’; the ‘units’ argument in the summary()
method can be used to provide a variable for the Petri dishes and calculate cluster-robust SEs, by way of the facilities provided in the ‘sandwich’ package (Zeileis et al. 2020).
summary(modHTE, robust = T, units = Dish)
##
## Model fitted: Hydrotime model with normal distribution of Psib (Bradford et al., 2002)
##
## Robust estimation: Cluster robust sandwich SEs
##
## Parameter estimates:
##
## Estimate Std. Error t value Pr(>|t|)
## ThetaH:(Intercept) 0.751069 0.131968 5.6913 3.075e-08 ***
## Psib50:(Intercept) -0.906981 0.039530 -22.9444 < 2.2e-16 ***
## sigmaPsib:(Intercept) 0.236995 0.031309 7.5696 4.974e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We see that the difference between ‘naive’ and cluster robust SEs is remarkable.
There might be other methods to obtain cluster robust standard errors (e.g., jackknife and bootstrap methods) and we are looking for ways to implement them in a reliable way. So far, we recommend that you make sure that your standard errors for model parameters do not neglect the clustering of seeds within randomisation units (petri dishes, pots, boxes or plots).