epifitter
provides functions for the fit of
(non-flexible) two-parameter population dynamics models to disease
progress curve (DPC) data: exponential, monomolecular, logistic and
Gompertz.
The goal of fitting these models to DPCs is to estimate numerically
meaningful parameters: y0
that represents primary inoculum
and r
, representing the apparent infection rate. Hence, the
importance of choosing the model that best describe the epidemics to
better understand its proterties and compare epidemics.
Two approaches can be used to obtain the parameters:
Both approaches are available in epifitter. The simplest way
to fit these models to single epidemic data is using the
fit_lin()
and fit_nlin()
functions. For the
latter, the alternative fit_nlin2()
allows estimating a
third parameter, the upper asymptote, when maximum disease intensity is
not close to 100%. The fit_multi()
function is the most
flexible and allows to fit all models to multiple epidemic datasets.
First, we need to load the packages we’ll need for this tutorial.
To use epifitter data at least two variables are needed, one representing the time of each assessment of disease intensity during the course of the epidemics and the other representing a disease intensity variable as proportion (e.g. incidence, severity, prevalence). For the case of designed experiments with replicates, a third variable is needed.
Let’s simulate a DPC dataset for one epidemics measured at replicated
plots. The simulated data resemble a polyciclic epidemics of sigmoid
shape. We can do that using the epifitter
sim_logistic()
function of epifitter (more about
?sim_logistic
here).
dpcL <- sim_logistic(
N = 100, # duration of the epidemics in days
y0 = 0.01, # disease intensity at time zero
dt = 10, # interval between assessments
r = 0.1, # apparent infection rate
alpha = 0.2, # level of noise
n = 7 # number of replicates
)
Let’s give a look at the simulated dataset.
## replicates time y random_y
## 1 1 0 0.01000000 0.01352602
## 2 1 10 0.02672668 0.03721111
## 3 1 20 0.06946279 0.07387658
## 4 1 30 0.16868385 0.16010530
## 5 1 40 0.35549349 0.35136344
## 6 1 50 0.59989486 0.68519594
The dpc_L
object generated using
sim_logistic()
is a dataframe with four columns. The
y
variable is a vector for the disease intensity as
proportion (0 < y < 1). To facilitate visualization, let’s make a
plot using the ggplot
function of the ggplot2
package.
ggplot(
dpcL,
aes(time, y,
group = replicates
)
) +
geom_point(aes(time, random_y), shape = 1) + # plot the replicate values
geom_point(color = "steelblue", size = 2) +
geom_line(color = "steelblue") +
labs(
title = "Simulated 'complete' epidemics of sigmoid shape",
subtitle = "Produced using sim_logistic()"
)+
theme_minimal_hgrid()
fit_lin()
The fit_lin()
requires at least the time
and y
arguments. In the example, we will call the
random_y
which represents the replicates. A quick way to
call these variables attached to the dataframe is shown below.
fit_lin()
outputs a list object which contains several
elements. Three elements of the list are shown by default: stats of
model fit, Infection rate and Initial Inoculum
## Results of fitting population models
##
## Stats:
## CCC r_squared RSE
## Logistic 0.9981 0.9962 0.1962
## Gompertz 0.9797 0.9602 0.4616
## Monomolecular 0.9395 0.8859 0.6276
## Exponential 0.9102 0.8353 0.6401
##
## Infection rate:
## Estimate Std.error Lower Upper
## Logistic 0.09955211 0.0007070471 0.09814360 0.09814360
## Gompertz 0.07076200 0.0016635421 0.06744805 0.06744805
## Monomolecular 0.05456648 0.0022616595 0.05006102 0.05006102
## Exponential 0.04498563 0.0023069068 0.04039004 0.04039004
##
## Initial inoculum:
## Estimate Linearized lin.SE Lower Upper
## Logistic 1.018511e-02 -4.576591 0.04182947 9.378438e-03 0.0110603992
## Gompertz 3.960448e-05 -2.316150 0.09841648 4.407977e-06 0.0002407096
## Monomolecular -1.718803e+00 -1.000192 0.13380158 -2.549257e+00 -1.0826588341
## Exponential 2.797626e-02 -3.576399 0.13647845 2.131641e-02 0.0367168333
The Stats
element of the list shows how each of the four
models predicted the observations based on three measures:
CCC
(Lin
2000), a measure of agreement that takes both bias and precision into
accountr_squared
(R2),
a measure of precisionRSE
for each model.The four models are sorted from the high to the low CCC
.
As expected because the sim_logistic
function was used to
create the synthetic epidemic data, the the Logistic model was
superior to the others.
The estimates, and respective standard error and upper and lower 95%
confidence interval, for the two coefficients of interest are shown in
the Infection rate
and Initial inoculum
elements. For the latter, both the back-transformed (estimate) and the
linearized estimate are shown.
The element f_lin$stats_all
provides a wide format
dataframe with all the stats for each model.
## # A tibble: 4 × 14
## best_model model r r_se r_ci_lwr r_ci_upr v0 v0_se r_squared RSE
## <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 Logi… 0.0996 7.07e-4 0.0981 0.101 -4.58 0.0418 0.996 0.196
## 2 2 Gomp… 0.0708 1.66e-3 0.0674 0.0741 -2.32 0.0984 0.960 0.462
## 3 3 Mono… 0.0546 2.26e-3 0.0501 0.0591 -1.00 0.134 0.886 0.628
## 4 4 Expo… 0.0450 2.31e-3 0.0404 0.0496 -3.58 0.136 0.835 0.640
## # ℹ 4 more variables: CCC <dbl>, y0 <dbl>, y0_ci_lwr <dbl>, y0_ci_upr <dbl>
The predicted values are stored as a dataframe in the
data
element called using the same $
operator
as above. Both the observed and (y
) and the
back-transformed predictions (predicted
) are shown for each
model. The linearized value and the residual are also shown.
## time y model linearized predicted residual
## 1 0 0.01352602 Exponential -4.30314004 2.797626e-02 -0.014450240
## 2 0 0.01352602 Monomolecular 0.01361833 -1.718803e+00 1.732329203
## 3 0 0.01352602 Logistic -4.28952171 1.018511e-02 0.003340907
## 4 0 0.01352602 Gompertz -1.45934500 3.960448e-05 0.013486416
## 5 10 0.03721111 Exponential -3.29114778 4.386921e-02 -0.006658092
## 6 10 0.03721111 Monomolecular 0.03792112 -5.754280e-01 0.612639151
The plot_fit()
produces, by default, a panel of plots
depicting the observed and predicted values by all fitted models. The
arguments pont_size
and line_size
that control
for the size of the dots for the observation and the size of the fitted
line, respectively.
The plots are ggplot2 objects which can be easily customized
by adding new layers that override plot paramaters as shown below. The
argument models
allows to select the models(s) to be shown
on the plot. The next plot was customized for the logistic model.
# Customized plots
plot_fit(f_lin,
point_size = 2,
line_size = 1,
models = "Logistic")+
theme_minimal_hgrid(font_size =18) +
scale_x_continuous(limits = c(0,100))+
scale_color_grey()+
theme(legend.position = "none")+
labs(
x = "Time",
y = "Proportion disease "
)
The fit_nlin()
function uses the Levenberg-Marquardt
algorithm for least-squares estimation of nonlinear parameters. In
addition to time and disease intensity, starting values for
y0
and r
should be given in the
starting_par
argument. The output format and interpretation
is analogous to the fit_lin()
.
NOTE: If you encounter error messages saying “matrix at initial parameter estimates”, try to modify the starting values for the parameters to solve the problem.
f_nlin <- fit_nlin(
time = dpcL$time,
y = dpcL$random_y,
starting_par = list(y0 = 0.01, r = 0.03)
)
f_nlin
## Results of fitting population models
##
## Stats:
## CCC r_squared RSE
## Logistic 0.9982 0.9963 0.0246
## Gompertz 0.9959 0.9931 0.0368
## Monomolecular 0.9123 0.8631 0.1568
## Exponential 0.8848 0.8214 0.1773
##
## Infection rate:
## Estimate Std.error Lower Upper
## Logistic 0.10259801 0.002083434 0.09844760 0.09844760
## Gompertz 0.07201356 0.002122129 0.06778606 0.06778606
## Monomolecular 0.02278152 0.001523597 0.01974636 0.01974636
## Exponential 0.01899407 0.001407175 0.01619083 0.01619083
##
## Initial inoculum:
## Estimate Std.error Lower Upper
## Logistic 8.843319e-03 8.639662e-04 7.122210e-03 1.056443e-02
## Gompertz 3.028919e-08 4.978786e-08 -6.889331e-08 1.294717e-07
## Monomolecular -1.924571e-01 4.692156e-02 -2.859297e-01 -9.898459e-02
## Exponential 1.817184e-01 2.124845e-02 1.393894e-01 2.240475e-01
We can check the results using plot_fit
.
K
(maximum disease)In many epidemics the last measure (final time) of a DPC does not
reach the maximum intensity and, for this reason, estimation of maximum
asymptote (carrying capacity K
) may be necessary. The
fin_lin2()
provides an estimation of K
in
addition to the estimates provided by fit_lin()
.
Before demonstrating the function, we can transform our simulated
data by creating another variable with y_random2
with
maximum about 0.8 (80%). Simplest way is to multiply the
y_random
by 0.8.
Then we run the fit_nlin2()
for the new dataset.
f_nlin2 <- fit_nlin2(
time = dpcL2$time,
y = dpcL2$random_y,
starting_par = list(y0 = 0.01, r = 0.2, K = 0.6)
)
f_nlin2
## Results of fitting population models
##
## Stats:
## CCC r_squared RSE
## Logistic 0.9982 0.9963 0.0198
## Gompertz 0.9714 0.9722 0.0750
## Monomolecular 0.9449 0.9024 0.1034
##
## Infection rate:
## Estimate Std.error Lower Upper
## Logistic 0.10290702 0.002528312 0.09786925 0.09786925
## Gompertz 0.10287913 0.012940724 0.07709417 0.07709417
## Monomolecular 0.01626713 0.003144204 0.01000217 0.01000217
##
## Initial inoculum:
## Estimate Std.error Lower Upper
## Logistic 6.993492e-03 7.801060e-04 5.439097e-03 8.547888e-03
## Gompertz 2.722258e-20 6.003036e-19 -1.168908e-18 1.223354e-18
## Monomolecular -1.459606e-01 3.295105e-02 -2.116170e-01 -8.030419e-02
##
## Maximum disease intensity:
## Estimate Std.error Lower Upper
## Logistic 0.7989443 0.004795787 0.7893885 0.8085001
## Gompertz 0.6781694 0.015959810 0.6463688 0.7099700
## Monomolecular 1.0000000 0.109221950 0.7823705 1.2176295
NOTE: The exponential model is not included because it doesn’t have a maximum asymptote. The estimated value of
K
is the expected 0.8.
Most commonly, there are more than one epidemics to analyse either
from observational or experimental studies. When the goal is to fit a
common model to all curves, the fit_multi()
function is in
hand. Each DPC needs an unique identified to further combined in a
single data frame.
Let’s use the sim_
family of functions to create three
epidemics and store the data in a single data.frame
. The
Gompertz model was used to simulate these data. Note that we allowed to
the y0
and r
parameter to differ the DPCs. We
should combine the three DPCs using the bind_rows()
function and name the identifier (.id
), automatically
created as a character vector, for each epidemics as ‘DPC’.
epi1 <- sim_gompertz(N = 60, y0 = 0.001, dt = 5, r = 0.1, alpha = 0.4, n = 4)
epi2 <- sim_gompertz(N = 60, y0 = 0.001, dt = 5, r = 0.12, alpha = 0.4, n = 4)
epi3 <- sim_gompertz(N = 60, y0 = 0.003, dt = 5, r = 0.14, alpha = 0.4, n = 4)
multi_epidemic <- bind_rows(epi1,
epi2,
epi3,
.id = "DPC"
)
head(multi_epidemic)
## DPC replicates time y random_y
## 1 1 1 0 0.00100000 0.00137268
## 2 1 1 5 0.01515505 0.01928029
## 3 1 1 10 0.07878459 0.08434323
## 4 1 1 15 0.21411521 0.17902470
## 5 1 1 20 0.39266393 0.41933570
## 6 1 1 25 0.56723412 0.41985063
We can visualize the three DPCs in a same plot
p_multi <- ggplot(multi_epidemic,
aes(time, y, shape = DPC, group = DPC))+
geom_point(size =2)+
geom_line()+
theme_minimal_grid(font_size =18) +
labs(
x = "Time",
y = "Proportion disease "
)
p_multi
Or use facet_wrap()
for ploting them separately.
fit_multi()
fit_multi()
requires at least four arguments: time,
disease intensity (as proportion), data and the curve identifier
(strata_cols
). The latter argument accepts one or more
strata include as c("strata1",strata2")
. In the example
below, the stratum name is DPC
, the name of the
variable.
By default, the linear regression is fitted to data but adding
another argument nlin = T
, the non linear regressions is
fitted instead.
multi_fit <- fit_multi(
time_col = "time",
intensity_col = "random_y",
data = multi_epidemic,
strata_cols = "DPC"
)
All parameters of the list can be returned using the $ operator as below.
## DPC best_model model r r_se r_ci_lwr r_ci_upr
## 1 1 1 Gompertz 0.10175467 0.002246154 0.09724314 0.10626620
## 2 1 2 Monomolecular 0.07344722 0.003036010 0.06734921 0.07954523
## 3 1 3 Logistic 0.16072637 0.007293917 0.14607611 0.17537664
## 4 1 4 Exponential 0.08727915 0.008995850 0.06921046 0.10534785
## 5 2 1 Gompertz 0.11640533 0.002046837 0.11229414 0.12051653
## 6 2 2 Monomolecular 0.09092611 0.002866436 0.08516871 0.09668352
## v0 v0_se r_squared RSE CCC y0
## 1 -1.9635096 0.07941355 0.9762159 0.3030228 0.9879648 0.0008053075
## 2 -0.6317502 0.10733918 0.9212913 0.4095802 0.9590334 -0.8808996811
## 3 -4.6234652 0.25787892 0.9066417 0.9840032 0.9510352 0.0097232437
## 4 -3.9917150 0.31805132 0.6530954 1.2136065 0.7901485 0.0184680151
## 5 -1.8291347 0.07236661 0.9847760 0.2761334 0.9923296 0.0019724179
## 6 -0.6656043 0.10134381 0.9526614 0.3867033 0.9757569 -0.9456660284
## y0_ci_lwr y0_ci_upr
## 1 0.0002349107 0.002302151
## 2 -1.3334487024 -0.516118013
## 3 0.0058153023 0.016214530
## 4 0.0097494945 0.034983104
## 5 0.0007444043 0.004580772
## 6 -1.3849052006 -0.587323594
Similarly, all data can be returned.
## DPC time y model linearized predicted residual
## 1 1 0 0.00137268 Exponential -6.590990488 0.0184680151 -0.0170953354
## 2 1 0 0.00137268 Monomolecular 0.001373623 -0.8808996811 0.8822723608
## 3 1 0 0.00137268 Logistic -6.589616865 0.0097232437 -0.0083505641
## 4 1 0 0.00137268 Gompertz -1.885703639 0.0008053075 0.0005673722
## 5 1 5 0.01928029 Exponential -3.948671980 0.0285722532 -0.0092919638
## 6 1 5 0.01928029 Monomolecular 0.019468578 -0.3027978168 0.3220781062
If nonlinear regression is preferred, the nlim
argument
should be set to TRUE
multi_fit2 <- fit_multi(
time_col = "time",
intensity_col = "random_y",
data = multi_epidemic,
strata_cols = "DPC",
nlin = TRUE)
## Warning in log(y0/1): NaNs produced
## Warning in log(y0/1): NaNs produced
## Warning in log(y0/1): NaNs produced
## Warning in log(y0/1): NaNs produced
## Warning in log(y0/1): NaNs produced
## Warning in log(y0/1): NaNs produced
## Warning in log(y0/1): NaNs produced
## Warning in log(y0/1): NaNs produced
## DPC model y0 y0_se r r_se df
## 1 1 Gompertz 0.0009172858 0.0005485373 0.10050834 0.003651136 50
## 2 1 Logistic 0.0307905943 0.0051255574 0.14408262 0.006811398 50
## 3 1 Monomolecular -0.1827722118 0.0443567461 0.04276458 0.002642699 50
## 4 1 Exponential 0.2280574876 0.0289500577 0.02760485 0.002630089 50
## 5 2 Gompertz 0.0024386351 0.0010399392 0.11491018 0.003677132 50
## 6 2 Logistic 0.0378761089 0.0053843460 0.16452413 0.007085072 50
## CCC r_squared RSE y0_ci_lwr y0_ci_upr r_ci_lwr
## 1 0.9942737 0.9886172 0.04093061 -0.0001844839 0.002019055 0.09317482
## 2 0.9910036 0.9838327 0.05111148 0.0204956092 0.041085579 0.13040153
## 3 0.9498043 0.9177118 0.11501624 -0.2718653583 -0.093679065 0.03745656
## 4 0.8658364 0.7920559 0.17958999 0.1699095855 0.286205390 0.02232216
## 5 0.9959988 0.9920613 0.03386931 0.0003498558 0.004527414 0.10752445
## 6 0.9933463 0.9881513 0.04348315 0.0270613318 0.048690886 0.15029334
## r_ci_upr best_model
## 1 0.10784186 1
## 2 0.15776372 2
## 3 0.04807260 3
## 4 0.03288754 4
## 5 0.12229592 1
## 6 0.17875491 2
If you want to estimate K
, set nlin = TRUE
and estimate_K = TRUE
.
NOTE: If you do not set both arguments
TRUE
,K
will not be estimated, becausenlin
defaut isFALSE
. Also remember that when estimating K, we don’t fit the Exponential model.
multi_fit_K <- fit_multi(
time_col = "time",
intensity_col = "random_y",
data = multi_epidemic,
strata_cols = "DPC",
nlin = T,
estimate_K = T
)
## Warning in log(y0/K): NaNs produced
## Warning in log(y0/K): NaNs produced
## Warning in log(y0/K): NaNs produced
## Warning in log(y0/K): NaNs produced
## Warning in log(y0/K): NaNs produced
## Warning in log(y0/K): NaNs produced
## Warning in log(y0/K): NaNs produced
## DPC model y0 y0_se r r_se K
## 1 1 Gompertz 0.0009172856 0.0006884191 0.10050834 0.005531408 1
## 2 1 Logistic 0.0307905704 0.0058290206 0.14408266 0.008502368 1
## 3 1 Monomolecular -0.1827771887 0.0490730111 0.04276504 0.006891756 1
## 4 2 Gompertz 0.0024386352 0.0012145563 0.11491018 0.004969495 1
## 5 2 Logistic 0.0378761713 0.0058725189 0.16452409 0.008276520 1
## 6 2 Monomolecular -0.1622789864 0.0439639627 0.05208062 0.006395657 1
## K_se df CCC r_squared RSE y0_ci_lwr y0_ci_upr
## 1 0.015016867 49 0.9942737 0.9886172 0.04134616 -4.661443e-04 0.002300716
## 2 0.015415324 49 0.9910036 0.9838327 0.05163039 1.907671e-02 0.042504426
## 3 0.069116640 49 0.9498045 0.9177105 0.11618395 -2.813931e-01 -0.084161281
## 4 0.009783848 49 0.9959988 0.9920613 0.03421317 -2.107024e-06 0.004879377
## 5 0.010794912 49 0.9933463 0.9881513 0.04392462 2.607490e-02 0.049677440
## 6 0.044773228 49 0.9616234 0.9371599 0.10149506 -2.506279e-01 -0.073930096
## r_ci_lwr r_ci_upr K_ci_lwr K_ci_upr best_model
## 1 0.08939256 0.11162412 0.9698225 1.030178 1
## 2 0.12699651 0.16116880 0.9690217 1.030978 2
## 3 0.02891554 0.05661455 0.8611049 1.138895 3
## 4 0.10492361 0.12489676 0.9803386 1.019661 1
## 5 0.14789180 0.18115638 0.9783068 1.021693 2
## 6 0.03922806 0.06493317 0.9100248 1.089975 3
Use ggplot2
to produce elegant data visualizations of models curves and the
estimated parameters.
The original data and the predicted values by each model are stored
in multi_fit$Data
. A nice plot can be produced as
follows:
multi_fit$Data %>%
ggplot(aes(time, predicted, color = DPC)) +
geom_point(aes(time, y), color = "gray") +
geom_line(size = 1) +
facet_grid(DPC ~ model, scales = "free_y") +
theme_minimal_hgrid()+
coord_cartesian(ylim = c(0, 1))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Using the dplyr function filter
only the model
of interest can be chosen for plotting.
multi_fit$Data %>%
filter(model == "Gompertz") %>%
ggplot(aes(time, predicted, color = DPC)) +
geom_point(aes(time, y),
color = "gray",
size = 2
) +
geom_line(size = 1.2) +
theme_minimal_hgrid() +
labs(
x = "Time",
y = "Disease Intensity"
)
The multi_fit$Parameters
element is where all stats and
parameters as stored. Let’s plot the estimates of the apparent infection
rate.
multi_fit$Parameters %>%
filter(model == "Gompertz") %>%
ggplot(aes(DPC, r)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = r_ci_lwr, ymax = r_ci_upr),
width = 0,
size = 1
) +
labs(
x = "Time",
y = "Apparent infection rate"
) +
theme_minimal_hgrid()
Lin L (2000). A note on the concordance correlation coefficient. Biometrics 56: 324 - 325.