`nlstimedist`

PackageThis vignette presents the `nlstimedist`

package, a method to fit a new distribution model to the time distribution of a biological phenomenon (Franco 2016). The model differentiates between three essential aspects of a time distribution: the rate at which the process is expected to occur (parameter \(r\)), the rate of change of \(r\) with time, which is reflected in the time concentration of the distribution (parameter \(c\)), and a measure of the overall distribution time lag (parameter \(t\)). The fitting method incorporates the `minpack.lm::nlsLM`

function (Elzhov et al. 2016) to estimate these three parameters and to plot the estimated time distribution. The `nlstimedist`

package, however, also estimates the standard distribution moments. The method is being proposed to analyse the time distribution of biological events such as germination, phenology, invasion, conclusion of a race, etc. Because parameter values have clear, unique effects on three different aspects of the distribution’s shape (and are correlated but not identical to specific moments), they have clear biological interpretation. This allows the user to further investigate the effect that biological (e.g., species, gender, health, etc.) and environmental factors (e.g., temperature) have on a biological time course. For example, are differences between the sexes in the completion of a marathon race reflected in a particular parameter? If so, what do these differences mean in terms of their size, musculature, aerobic capacity, etc.? If the parameters have a biological interpretation, how are they affected by ambient temperature, hydration, sugar levels, etc.?

In the model, time is represented by variable \(x\) and the biological phenomenon is represented by variable \(y\). The values in each \(y\) column should be proportions and should be calculated from the cumulative number of events. This must be completed for each column in a dataset. If data have been set up in this manner, skip ahead to the modelling section. If the data have not been set up in this format and it is in a raw format of counts vs. time, they must first be cleaned using the `tdData`

function.

`tdData`

The `nlstimedist`

package comes with several example datasets, one being the `lobelia`

dataset.

```
head(lobelia)
# A tibble: 6 x 3
Day Temperature Germination
<int> <dbl> <int>
1 0 9.8 0
2 1 9.8 0
3 2 9.8 0
4 3 9.8 0
5 4 9.8 0
6 5 9.8 0
```

We can clean and prepare the data for modelling using the `tdData`

function.

```
tdLobelia <- tdData(lobelia, x = "Day", y = "Germination",
group = "Temperature")
tdLobelia
# A tibble: 82 x 5
Day Temperature Germination cumN propMax
<int> <dbl> <dbl> <dbl> <dbl>
1 18 9.8 3 3 0.188
2 19 9.8 3 6 0.375
3 20 9.8 1 7 0.438
4 21 9.8 1 8 0.5
5 22 9.8 1 9 0.562
6 23 9.8 2 11 0.688
7 25 9.8 1 12 0.75
8 28 9.8 2 14 0.875
9 29 9.8 2 16 1
10 9 12.5 1 1 0.0435
# … with 72 more rows
```

The model is fitted by nonlinear regression employing the Levenberg-Marquardt algorithm. This requires three starting values for \(r\), \(c\) and \(t\), respectively.

Suggestions for appropriate starting values for each parameter are as follows:

Parameter | Recommendation |
---|---|

\(r\) | \(\frac{1}{\text{the period of the time course}}\), e.g., if completion of the process (all individual events) occurred in 25 days, an appropriate starting value for \(r\) would be around \(\frac{1}{25}\) = 0.04. |

\(c\) | This requires some trial and error with your particular dataset. We suggest you start with 0.5 and increase (or decrease) it along a logarithmic scale to get a feel of how it is changing. Increasing values of \(c\) reduce the spread of the distribution: \(c\) is a measure of concentration of the distribution. |

\(t\) | This tends to be close to the mid-point of the monitoring period, but it varies with the skew produced by the combination of parameter values. Nonetheless, as a rule of thumb choose a number near the middle of your time range – if completion of a process (e.g., a marathon race) was closed after 10 hours, choose \(t = 5\). |

`timedist`

FunctionThe model is fitted to the data using the `timedist`

function.

```
# Fitting the model to data already in the format x = time and y = proportion
# of cumulative number of events.
lobelia12_5 <- tdLobelia$clean[tdLobelia$clean$Temperature == 12.5, ]
model12_5 <- timedist(lobelia12_5, x = "Day", y = "propMax",
r = 0.03, c = 0.5, t = 14.5)
model12_5
Nonlinear regression model
model: propMax ~ 1 - (1 - (r/(1 + exp(-c * (Day - t)))))^Day
data: data
r c t
0.08339 0.62678 12.09364
residual sum-of-squares: 0.03901
Number of iterations to convergence: 16
Achieved convergence tolerance: 1.49e-08
```

On rare occasions the model may fail to converge within 50 iterations. This may occur if a very small dataset is used. It is possible to overcome this issue by fixing or setting upper and lower bounds for one of the starting values. The parameter \(r\) is the most appropriate parameter to do this with. It is suggested that you calculate the starting value for \(r\) as in the starting values section and set the upper and lower bounds around this figure (see below).

```
modelFix <- timedist(data = lobelia12_5, x = "Day", y = "propMax",
r = 0.03, c = 0.5, t = 14.5,
upper = c(0.1, Inf, Inf), lower = c(0.01, -Inf, -Inf))
modelFix
Nonlinear regression model
model: propMax ~ 1 - (1 - (r/(1 + exp(-c * (Day - t)))))^Day
data: data
r c t
0.08339 0.62678 12.09364
residual sum-of-squares: 0.03901
Number of iterations to convergence: 16
Achieved convergence tolerance: 1.49e-08
```

To assess how well the model has fit the data, and the reliability of parameter estimates, it is suggested that the standard errors, correlations of the estimates, and confidence intervals are obtained. In each example we have used the model `model12_5`

.

```
summary(model12_5, correlation = TRUE, symbolic.cor = FALSE)
Formula: propMax ~ 1 - (1 - (r/(1 + exp(-c * (Day - t)))))^Day
Parameters:
Estimate Std. Error t value Pr(>|t|)
r 0.083393 0.006615 12.607 6.99e-08 ***
c 0.626782 0.156122 4.015 0.00203 **
t 12.093637 0.462246 26.163 2.95e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.05955 on 11 degrees of freedom
Correlation of Parameter Estimates:
r c
c -0.49
t 0.72 -0.52
Number of iterations to convergence: 16
Achieved convergence tolerance: 1.49e-08
```

If a higher level of precision is required, the correlation of parameter estimates can be obtained separately.

```
cpe <- vcov(model12_5)
cov2cor(cpe)
r c t
r 1.0000000 -0.4857904 0.7153296
c -0.4857904 1.0000000 -0.5247587
t 0.7153296 -0.5247587 1.0000000
```

To produce accurate confidence intervals for the parameters in a nonlinear regression model fit, we can use the `confint2`

function.

```
confint2(model12_5)
2.5 % 97.5 %
r 0.06883386 0.09795187
c 0.28316063 0.97040290
t 11.07624152 13.11103281
```

There is no direct R-squared for non-linear regression. However, an R-squared value calculated as \(1-\bigg(\frac{\text{Residual Sum of Squares}}{\text{Corrected Sum of Squares}}\bigg)\) defines a similar quantity for nonlinear regression, is able to describe the proportion of variance explained by the model, and provides a very good estimate of how well the model fits the data. We can extract this value from our model using the `tdRSS`

function.

```
tdRSS(model12_5)
[1] 0.9681957
```

The following statistical moments for the fitted distribution can be calculated: mean, variance, standard deviation, skew, kurtosis and entropy.

```
model12_5$m$getMoments()
mean variance sd skew kurtosis entropy
1 15.75401 83.02729 9.111931 2.897078 12.17524 4.491001
```

The percentiles of the distribution can also be calculated. This can be achieved for a single percentile or for a sequence of percentiles.

```
# Extracting a single percentile
tdPercentiles(model12_5, n = 0.01)
1%
5.913667
# Extracting a sequence of percentiles from 10% to 90% in steps of 10.
tdPercentiles(model12_5, n = seq(0.1, 0.9, 0.1))
10% 20% 30% 40% 50% 60% 70%
9.159305 10.382816 11.269057 12.073122 12.918504 13.952189 15.516796
80% 90%
18.776037 26.446720
```

The package has two built-in graphing functions for plotting the estimated distribution as both a probability density function and a cumulative distribution function.

The PDF is produced using the function `tdPdfPlot`

. This function takes one or more objects produced by the model, a scaling parameter `S`

and values for the x-axis `xVals`

(which includes a value for smoothing the curve), as arguments to produce the PDF plot.

`tdPdfPlot(model12_5, S = 1, xVals = seq(0, 30, 0.01))`

Multiple models can be plotted on the same graph by providing the function with multiple model objects.

```
# Extract the individual data
lobelia9_8 <- tdLobelia$clean[tdLobelia$clean$Temperature == 9.8, ]
lobelia16_7 <- tdLobelia$clean[tdLobelia$clean$Temperature == 16.7, ]
lobelia20_2 <- tdLobelia$clean[tdLobelia$clean$Temperature == 20.2, ]
lobelia24_3 <- tdLobelia$clean[tdLobelia$clean$Temperature == 24.3, ]
lobelia28_5 <- tdLobelia$clean[tdLobelia$clean$Temperature == 28.5, ]
lobelia32 <- tdLobelia$clean[tdLobelia$clean$Temperature == 32, ]
# Create the models
model9_8 <- timedist(lobelia9_8, x = "Day", y = "propMax",
r = 0.1, c = 0.5, t = 25)
model16_7 <- timedist(lobelia16_7, x = "Day", y = "propMax",
r = 0.1, c = 0.5, t = 10)
model20_2 <- timedist(lobelia20_2, x = "Day", y = "propMax",
r = 0.1, c = 0.5, t = 10)
model24_3 <- timedist(lobelia24_3, x = "Day", y = "propMax",
r = 0.1, c = 1, t = 5)
model28_5 <- timedist(lobelia28_5, x = "Day", y = "propMax",
r = 0.1, c = 1, t = 5)
model32 <- timedist(lobelia32, x = "Day", y = "propMax",
r = 0.1, c = 0.5, t = 5)
# Generate the plot
tdPdfPlot(model9_8, model12_5, model16_7, model20_2, model24_3, model28_5,
model32,
S = c(0.213, 0.307, 0.533, 0.707, 0.867, 0.907, 0.840),
xVals = seq(0, 30, 0.001))
```

The CDF is produced using the function `tdCdfPlot`

. This function takes one or more objects produced by the model, a scaling parameter `S`

and values for the x-axis `xVals`

(which includes a value for smoothing the curve), as arguments to produce the CDF plot.

`tdCdfPlot(model12_5, S = 1, xVals = seq(0, 30, 0.01))`

```
tdCdfPlot(model9_8, model12_5, model16_7, model20_2, model24_3, model28_5,
model32,
S = c(0.213, 0.307, 0.533, 0.707, 0.867, 0.907, 0.840),
xVals = seq(0, 30, 0.001))
```

Elzhov, Timur V., Katharine M. Mullen, Andrej-Nikolai Spiess, and Ben Bolker. 2016. *Minpack.lm: R Interface to the Levenberg-Marquardt Nonlinear Least-Squares Algorithm Found in Minpack, Plus Support for Bounds*. https://CRAN.R-project.org/package=minpack.lm.

Franco, Miguel. 2016. *The Time-Course of Biological Phenomena – Illustrated with the London Marathon*.