Fitting and selecting models
Base functions
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:
- Linear regression models that are fitted to a specific transformation of the disease data according to each of the models.
- Non-linear regression models fitted to original data
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.
Basic usage
Dataset
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.01000000
## 2 1 10 0.02672668 0.01969015
## 3 1 20 0.06946279 0.06893039
## 4 1 30 0.16868385 0.17026089
## 5 1 40 0.35549349 0.36379455
## 6 1 50 0.59989486 0.57301510The 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()
Linear regression using 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.9980 0.9961 0.1997
## Gompertz 0.9765 0.9540 0.5006
## Monomolecular 0.9340 0.8762 0.6633
## Exponential 0.9144 0.8423 0.6171
##
## Infection rate:
## Estimate Std.error Lower Upper
## Logistic 0.09957284 0.000719583 0.09813936 0.09813936
## Gompertz 0.07115465 0.001803955 0.06756098 0.06756098
## Monomolecular 0.05506327 0.002390469 0.05030121 0.05030121
## Exponential 0.04450957 0.002223781 0.04007957 0.04007957
##
## Initial inoculum:
## Estimate Linearized lin.SE Lower Upper
## Logistic 0.010178364 -4.577260 0.04257111 9.358512e-03 0.011069237
## Gompertz 0.000031155 -2.339547 0.10672345 2.666663e-06 0.000227313
## Monomolecular -1.795273306 -1.027930 0.14142203 -2.704903e+00 -1.108976501
## Exponential 0.028743875 -3.549331 0.13156068 2.211691e-02 0.037356504Model fit stats
The Stats element of the list shows how each of the four
models predicted the observations based on three measures:
- Lin’s concordance correlation coefficient
CCC(Lin 2000), a measure of agreement that takes both bias and precision into account - Coefficient of determination
r_squared(R2), a measure of precision - Residual standard deviation
RSEfor 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.
Model coefficients
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.
Global stats
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.20e-4 0.0981 0.101 -4.58 0.0426 0.996 0.200
## 2 2 Gomp… 0.0712 1.80e-3 0.0676 0.0747 -2.34 0.107 0.954 0.501
## 3 3 Mono… 0.0551 2.39e-3 0.0503 0.0598 -1.03 0.141 0.876 0.663
## 4 4 Expo… 0.0445 2.22e-3 0.0401 0.0489 -3.55 0.132 0.842 0.617
## # ℹ 4 more variables: CCC <dbl>, y0 <dbl>, y0_ci_lwr <dbl>, y0_ci_upr <dbl>Model predictions
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.01000000 Exponential -4.60517019 0.028743875 -0.0187438753
## 2 0 0.01000000 Monomolecular 0.01005034 -1.795273306 1.8052733065
## 3 0 0.01000000 Logistic -4.59511985 0.010178364 -0.0001783642
## 4 0 0.01000000 Gompertz -1.52717963 0.000031155 0.0099688450
## 5 10 0.01969015 Exponential -3.92763666 0.044858828 -0.0251686757
## 6 10 0.01969015 Monomolecular 0.01988659 -0.611712361 0.6314025133Plot of predictions
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.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the epifitter package.
## Please report the issue at <https://github.com/AlvesKS/epifitter/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Publication-ready plots
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 "
)
Non-linear regression
Two-parameters
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.9977 0.9955 0.0272
## Gompertz 0.9954 0.9921 0.0391
## Monomolecular 0.9116 0.8621 0.1568
## Exponential 0.8884 0.8266 0.1743
##
## Infection rate:
## Estimate Std.error Lower Upper
## Logistic 0.10077756 0.002250322 0.09629469 0.09629469
## Gompertz 0.07072098 0.002194214 0.06634988 0.06634988
## Monomolecular 0.02252155 0.001507002 0.01951945 0.01951945
## Exponential 0.01916342 0.001396564 0.01638132 0.01638132
##
## Initial inoculum:
## Estimate Std.error Lower Upper
## Logistic 9.207821e-03 9.807453e-04 7.254076e-03 1.116157e-02
## Gompertz 4.906495e-08 8.184169e-08 -1.139721e-07 2.121020e-07
## Monomolecular -1.912683e-01 4.681954e-02 -2.845376e-01 -9.799898e-02
## Exponential 1.784771e-01 2.074455e-02 1.371518e-01 2.198023e-01We can check the results using plot_fit.
Estimating 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.9977 0.9955 0.0219
## Gompertz 0.9716 0.9709 0.0748
## Monomolecular 0.9443 0.9016 0.1036
##
## Infection rate:
## Estimate Std.error Lower Upper
## Logistic 0.10078900 0.002739828 0.095329778 0.095329778
## Gompertz 0.10330774 0.012970636 0.077463189 0.077463189
## Monomolecular 0.01610849 0.003158304 0.009815437 0.009815437
##
## Initial inoculum:
## Estimate Std.error Lower Upper
## Logistic 7.363085e-03 8.961341e-04 5.577499e-03 9.148671e-03
## Gompertz 1.012327e-21 2.436108e-20 -4.752818e-20 4.955283e-20
## Monomolecular -1.453687e-01 3.297468e-02 -2.110721e-01 -7.966517e-02
##
## Maximum disease intensity:
## Estimate Std.error Lower Upper
## Logistic 0.7999605 0.005440932 0.7891192 0.8108017
## Gompertz 0.6794420 0.016001832 0.6475576 0.7113263
## Monomolecular 1.0000000 0.111379506 0.7780715 1.2219285
NOTE: The exponential model is not included because it doesn’t have a maximum asymptote. The estimated value of
Kis the expected 0.8.
Fit models to multiple DPCs
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.
Data
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.00100000
## 2 1 1 5 0.01515505 0.01291719
## 3 1 1 10 0.07878459 0.07994754
## 4 1 1 15 0.21411521 0.18039002
## 5 1 1 20 0.39266393 0.16276550
## 6 1 1 25 0.56723412 0.58871013We 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.
Using 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.10595945 0.003026184 0.09988118 0.11203772
## 2 1 2 Logistic 0.16534902 0.006962635 0.15136416 0.17933388
## 3 1 3 Monomolecular 0.07731844 0.003820213 0.06964532 0.08499156
## 4 1 4 Exponential 0.08803058 0.008677199 0.07060191 0.10545925
## 5 2 1 Gompertz 0.12149314 0.002107657 0.11725979 0.12572650
## 6 2 2 Monomolecular 0.09507345 0.003357728 0.08832926 0.10181765
## v0 v0_se r_squared RSE CCC y0
## 1 -2.0562726 0.10699177 0.9608149 0.4082546 0.9800159 0.0004029172
## 2 -4.7390650 0.24616632 0.9185628 0.9393108 0.9575530 0.0086709772
## 3 -0.7049836 0.13506492 0.8912167 0.5153749 0.9424797 -1.0238135591
## 4 -4.0340814 0.30678531 0.6730358 1.1706181 0.8045683 0.0177019344
## 5 -2.0058465 0.07451693 0.9851755 0.2843385 0.9925324 0.0005917757
## 6 -0.7917888 0.11871360 0.9412959 0.4529822 0.9697603 -1.2073414433
## y0_ci_lwr y0_ci_upr
## 1 6.185293e-05 0.001826796
## 2 5.306487e-03 0.014138346
## 3 -1.654535e+00 -0.542952451
## 4 9.558946e-03 0.032781697
## 5 1.782431e-04 0.001662867
## 6 -1.801716e+00 -0.739061415Similarly, all data can be returned.
## DPC time y model linearized predicted residual
## 1 1 0 0.00100000 Exponential -6.90775528 0.0177019344 -0.0167019344
## 2 1 0 0.00100000 Monomolecular 0.00100050 -1.0238135591 1.0248135591
## 3 1 0 0.00100000 Logistic -6.90675478 0.0086709772 -0.0076709772
## 4 1 0 0.00100000 Gompertz -1.93264473 0.0004029172 0.0005970828
## 5 1 5 0.01291719 Exponential -4.34919633 0.0274901243 -0.0145729347
## 6 1 5 0.01291719 Monomolecular 0.01300134 -0.3749143380 0.3878315275If 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## DPC model y0 y0_se r r_se df
## 1 1 Gompertz 0.000155693 0.0001560439 0.10641384 0.004736988 50
## 2 1 Logistic 0.021808953 0.0041414919 0.15388274 0.007524517 50
## 3 1 Monomolecular -0.191074957 0.0497424790 0.04192926 0.002890979 50
## 4 1 Exponential 0.218008103 0.0291040197 0.02847442 0.002750494 50
## 5 2 Gompertz 0.001316706 0.0009040714 0.11240869 0.005046403 50
## 6 2 Logistic 0.031505686 0.0054861648 0.16310258 0.008121478 50
## CCC r_squared RSE y0_ci_lwr y0_ci_upr r_ci_lwr
## 1 0.9922281 0.9845813 0.04879904 -0.0001577303 0.0004691164 0.09689931
## 2 0.9913334 0.9836448 0.05139205 0.0134905213 0.0301273837 0.13876930
## 3 0.9379506 0.8994216 0.12961549 -0.2909856665 -0.0911642477 0.03612256
## 4 0.8644663 0.7906910 0.18418430 0.1595509592 0.2764652471 0.02294989
## 5 0.9922107 0.9847190 0.04801676 -0.0004991751 0.0031325867 0.10227269
## 6 0.9912929 0.9835343 0.05065014 0.0204864001 0.0425249727 0.14679011
## r_ci_upr best_model
## 1 0.11592836 1
## 2 0.16899618 2
## 3 0.04773596 3
## 4 0.03399895 4
## 5 0.12254469 1
## 6 0.17941505 2Want to estimate K?
If you want to estimate K, set nlin = TRUE
and estimate_K = TRUE.
NOTE: If you do not set both arguments
TRUE,Kwill not be estimated, becausenlindefaut 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
## 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.0001556867 0.0001966899 0.10641400 0.007042914 1
## 2 1 Logistic 0.0218089645 0.0046998139 0.15388272 0.009243869 1
## 3 1 Monomolecular -0.1910765702 0.0551678550 0.04192940 0.007694815 1
## 4 2 Gompertz 0.0013167114 0.0010751024 0.11240865 0.006978295 1
## 5 2 Logistic 0.0315056770 0.0060336522 0.16310259 0.009573314 1
## 6 2 Monomolecular -0.1721258929 0.0494344750 0.04946617 0.007073582 1
## K_se df CCC r_squared RSE y0_ci_lwr y0_ci_upr
## 1 0.01754255 49 0.9922281 0.9845813 0.04929448 -0.0002395765 0.0005509499
## 2 0.01515799 49 0.9913334 0.9836448 0.05191381 0.0123643348 0.0312535942
## 3 0.08040264 49 0.9379507 0.8994212 0.13093141 -0.3019405255 -0.0802126150
## 4 0.01462132 49 0.9922107 0.9847190 0.04850425 -0.0008437879 0.0034772106
## 5 0.01303618 49 0.9912929 0.9835343 0.05116437 0.0193805989 0.0436307551
## 6 0.05461727 49 0.9519795 0.9220498 0.11489666 -0.2714681897 -0.0727835962
## r_ci_lwr r_ci_upr K_ci_lwr K_ci_upr best_model
## 1 0.09226074 0.12056727 0.9647469 1.035253 1
## 2 0.13530647 0.17245897 0.9695389 1.030461 2
## 3 0.02646609 0.05739271 0.8384249 1.161575 3
## 4 0.09838525 0.12643206 0.9706174 1.029383 1
## 5 0.14386429 0.18234088 0.9738028 1.026197 2
## 6 0.03525127 0.06368106 0.8902425 1.109758 3Graphical outputs
Use ggplot2
to produce elegant data visualizations of models curves and the
estimated parameters.
DPCs and fitted curves
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))
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"
)
Apparent infection rate
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()
References
Lin L (2000). A note on the concordance correlation coefficient. Biometrics 56: 324 - 325.