\section{Introduction}

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), first reported in Wuhan, China at the end of 2019, spread across China and the globe and was declared a pandemic on March 11, 2020. As of middle June 2020, it had caused more than 8 million infections. As reported by \href{https://www.who.int/docs/default-source/coronaviruse/situation-reports/20200708-covid-19-sitrep-170.pdf?sfvrsn=bca86036_2}{WHO (2020)},  in recent months, USA and Europe have been epicentres of COVID-19 and are experiencing the rapid increase of its outbreak size. Although the growth of daily confirmed infections in these regions is slowing down, there is still no sign that the epidemic has gained any control or “flattened”  in these regions; even worse, Brazil, India and African region (all affected countries) have become the new epicentres of COVID-19 with increasing number of confirmed infection every day \citep{lancet2020covid, pearson2020projected}.\\

Mathematical modelling, including statistical modelling, is an important tool to understand and predict the dynamics of new diseases. Since the first identification of COVID-19, different approaches haves been developed to simulate and characterize its dynamics and spread trend \citep{grasselli2020critical, yang2020modified, roosa2020real, petropoulos2020forecasting, hengjian2020nonlinear, perc2020forecasting, sajadi2020temperature,zhang2020predicting}. The classical one is the dynamic infectious disease modeling, with using deterministic ODE models or stochastic individual based models, and such approach allows to incorporate with underlying mechanisms of spread and various risk factors in the simulation of transmission \citep{peng2020epidemic, wynants2020prediction, kucharski2020early,fanelli2020analysis}. This approach is commonly used to identify the crucial transmission parameters and  assess the potential impact of public health interventions. Although prediction of transmission trend can be also achieved, setting up such models needs heavy information of local demography and praxiology that is difficult to obtain accurately. For the purpose of prediction, data-driven (or phenomenological) methods are preferred, i.e., machine learning and statistical modeling \citep{alimadadi2020artificial, ribeiro2020short, benvenuto2020application, ceylan2020estimation, ribeiro2020short,zheng2020predicting,kavadi2020partial}. 

Almost investigated modeling approaches have been developed to characterize the transmission and impact of COVID-19 in the context of a specific country or region, from which this information can be estimated, e.g., forecasting the confirmed cases and deaths in China \citep{gao2020forecasting, al2020optimization}. However, for many regions, such modelling studies are still unavailable now. Knowledge of the inflection (“flattening the curve”) and maximal outbreak size is crucial to reflect the evolving trend of an epidemic; in the case of COVID-19, this information remains unclear but does influence the dates of changes to policy restrictions and the recovery of the global economy. Furthermore, an accurate prediction of such information at the early stage is difficult due to the lack of detailed data on testing availability and reporting/infection processes as well as governmental restrictions. \\

To address the problems, non-linear mixed effect model is used to model the (transformed) daily reported number of cumulative confirmed cases. The data was grouped according to country or region. In the countries at the later stage of the epidemics, the reported number of cumulative confirmed cases all show a sigmoidal shape with respect to time, so we use the four parameters Logistic model (FPLM), a generalization of Logistic growth model, to model the growth patterns. By fitting to such non-linear mixed effect model we can predict the inflection point and final size of outbreak in modeled countries and regions. \\

\section{Data}

The data set was downloaded from European Centre for Disease Prevention and Control, which provides the geographic distribution of COVID-19 cases worldwide, e.g. daily incidence (newly confirmed cases),  cumulative number of confirmed cases, population of each country, etc., from  \href{https://opendata.ecdc.europa.eu/covid19/casedistribution/csv}{European Centre for Disease Prevention and Control}. Note that we select two groups of countries/regions, the first is for the countries in the late stage of outbreak, including Australia, China, France, Italy, Germany, Spain and UK; and countries/regions in the second group are still in the early stage of the outbreak, including USA, Africa (as a whole), Brazil, India and Russia.
The data set is about the officially confirmed and reported cases, which is inevitably inaccurate and under-reported due to the limited coverage of testing, especially in the early period of the outbreak. Also, this study is more interested in the future growth trend and final outbreak size. Hence, early observations were thrown away to enable the model to fit the early period more flexibly.

\section{Growth curves with random coefficients}

Let $y(t)$ be the cumulative numbers at time $t$ and $Y’(t)$ be the derivative function representing the growth rate. If $x_t$ represents the explanatory variables(such as temperature, or behavour changes due to government restrictions) believed to be related to the growth rate, we need to incorporate their effect in the growth model via a link function
\begin{equation}
y’(t) = f(y(t))g(x_ty;\theta) + \epsilon_t,
\end{equation}
where $f(y)$ specifies the growth rate as a function of its current size under the constant environmental condition while
$g$ models how the growth rate might change when the environmental conditions ($x_t$) changes. Here $\sigma(t)$ is the stochastic error of zero mean representing the environmental or measurement perturbation  possibly with heteroscedasticity. More details can be seen in \cite{wang1999estimating}.

The simple linear function of $f$ corresponds to the asymptotic regression model (also known as von Bertalanffy growth curve).
However, we are particularly interested in a sigmoidal curves and the inflection point is of great interest. The well known  curves of this type include logistic (autocatalytic),  Richards and Gompertz \citep{Nonlinear1989}.\\

In this study, we will apply the mixed-effects models assuming each country follows the same curve but different set of parameters. These parameters can be potentially modelled as functions of population size and other attributes. More details can be seen in \cite{pinheiro2006mixed}.
A model for nonlinear mixed-effects can be written as
\begin{equation} \label{eq:nlmem}
 g(y_{ij}) = f(A_i \boldsymbol{\beta} + B_i \boldsymbol{b}_i, t_{ij})+\epsilon_{ij}, b_i \sim \mathcal{N}(0,\Sigma), \epsilon_{ij} \sim \mathcal{N}(0,\sigma^2)
\end{equation}
for observation $j, j=1,\ldots,n_i$ in group $i,i=1,\ldots,M$. In model~\eqref{eq:nlmem}, $\boldsymbol{\beta}_i$ includes both fixed effect $\boldsymbol{\beta}$ and random effects $\boldsymbol{b}_i$. Specifically in our case, $y$ is the cumulative number of confirmed cases, $g$ is the transformation used, $i$ is the index of country/region, and $t$ is time index of the observation (day). Here $A_i$ and $B_i$ are design matrices for $i$th group to determine the fixed and random effects. The advantage of mixed-effect model is to produce more precise estimate of  $\boldsymbol{b}_i$ and $\boldsymbol{\beta}_i$ by borrowing strength/information from the rest of the sample from the population. See \cite{pinheiro2006mixed} for more detains about non-linear mixed effect model.

To model the confirmed cases over time we use a four Four Parameter Logistic Model (FPLM)
\begin{equation}\label{eq:fpl}
f(\beta_i,t_{ij}) = \phi_{i1} + \frac{\phi_{i2} - \phi_{i1}}{1+\exp\left(\frac{\phi_{i3} - t_{ij}}{\phi_{i4}} \right)},
\end{equation}
where $\beta_i=(\phi_{i1},\phi_{i2},\phi_{i3},\phi_{i4}) $, and the parameters are:
\begin{itemize}
\item $\phi_{i1}$, the minimum theoretical value of $y_{ij}$ as time $t_{ij}\rightarrow -\infty$;
\item $\phi_{i2}$, the maximum value  as $t_{ij} \rightarrow \infty$;
\item $\phi_{i3}$, the inflection point, and  $t_{ij}=\phi_{i3}$ the response $y_{ij}$ is midway between the $\phi_{i1}$ and  $\phi_{i2}$;
\item $\phi_{i4}$, is a scale parameter for time.
\end{itemize}

We are particularly interested in the following two parameters. $n_{max}$: the maximum number (asymptote), and $n_{Infl}$: the number of cases at the inflection point ($\phi_1+\phi_2)/2$.

We first tried this on the raw data (Model 1). We then used power and logarithmic transformation of cumulative number of confirmed cases as the fitting response as follows.

\begin{enumerate}
\item[a.] Power (square root) transformation (Model 2)
\begin{equation}\label{eq.fplpw}
[y_{ij}]^{1/2}  = f(\beta_i,t_{ij}) = \phi_{i1} + \frac{\phi_{i2} - \phi_{i1}}{1+\exp\left(\frac{\phi_{i3} - t_{ij}}{\phi_{i4}} \right)}
\end{equation}
with $n_{max} =\phi_2^2$. $n_{Infl} =\sqrt{n_0 n_{\max}}$ where $n_0=\phi_1^2$;
\item[b.] Logarithmic Transformation (Model 3)
\begin{equation}\label{eq.fpllog}
\log_{10}[y_{ij}] = f(\beta_i,t_{ij}) = \phi_{i1} + \frac{\phi_{i2} - \phi_{i1}}{1+\exp\left(\frac{\phi_{i3} - t_{ij}}{\phi_{i4}} \right)}.
\end{equation}
\end{enumerate}
Note that the parameters of FPLM model have different interpretations under transformations. The models are validated by comparing the first-order difference of modeled $y_{ij}$ to the corresponding reported daily new confirmed cases.

\section{Results}

Here we fit the cumulative cases after power transformation to the well-known growth curves to describe the epidemic trends in the countries and regions of interest under the mixed-effect model and using the four-parameter logistic model. The advantage of the mixed-effect model is to “borrow information” from the members with rich information \citep{pinheiro2006mixed}. The four-parameter logistic model has been proved to perform well in describing epidemic growths \citep{wu2020generalized, chen2020reconstructing}.\\

We included 12 countries and regions in our study. Nine countries (Australia, China, France, Germany, Italy, Russia, Spain, UK and USA) have almost experienced a full growth curve, and other regions (Brazil, India and the African region) at the early phase of the epidemic will have certain similarity with some of these nine countries in transmission and response strategies. This may result in similarity in the growth of confirmed cases. \\

The fitted results for all selected countries/regions under three different models are shown in Table~\ref{tab:results}. Note that  $n_0= 10^{\phi_1}$, $n_{max} =10^{\phi_2}$, and $n_{Infl}$ is the number of cases at the inflection point, which is $\sqrt{10^{\phi_1+\phi_2}}=\sqrt{n_0 n_{max}}$. Model 1, 2 and 3 correspond to no transformation, power transformation and log transformation of the cumulative number of confirmed cases in ~\eqref{eq:fpl},~\eqref{eq.fplpw} and~\eqref{eq.fpllog}, respectively. Additionally, Figure~\ref{fig:fitm1}, ~\ref{fig:fitm2} and ~\ref{fig:fitm3} show the fitted curves and data (circles  in the figures) used for fitting of each country or region. (Each color represents the growth curve of one country or region. The dots in the figure are the reported number of cumulative confirmed cases and the curves are the estimated growth curves. The bottom horizontal line segments are the estimate of maximal outbreak size, while top horizontal line segments are the upper boundary of 95\% confidence interval.) Moreover, The reported number of daily confirmed cases (daily incidence in Figure~\ref{fig:vld1} and \ref{fig:vld2}) is used to validate the fitting model.\\

Additionally, here it is statistically sensible to consider a power transformation, $y^\theta$. The well known Box-Cox transformation, $(y^\theta-1)/\theta$, would be a sensible choice (the limiting case of $\theta \rightarrow 0$ corresponds to the log-transformation. An 'optimal' or estimate of $\theta$ can also be obtained via likelihood or other robust statistical approaches. We have chosen three different $\theta$ values (1, 0.5 and 0) in this analysis. It is interesting to see how robust the results are in terms of ballpark estimates of the inflection point and maximum number of cases.\\

Our models have confirmed that the nine countries aforementioned had passed the inflection point (the squared marks in Figure ~\ref{fig:fitm1}). Specifically, China had the inflection point on February 9, 2020, with a final outbreak size of 87k \citep{gu2020inflection}. Australia, Italy, Spain, Germany and France have passed their inflection point in later March and early April; in these countries the current outbreak size almost reaches the estimates of maximal level (with invisible upper confidence boundaries in Figure ~\ref{fig:fitm1}).  Unfortunately, cases in other countries (i.e., USA, Russia, and UK) will continue to increase, possibly to 2,357k (up to 2,425k); 538k (up to 546k); and 309k (up to 313k), respectively. \\

In terms of other three regions (Brazil, India, and African region), they are still in the early outbreak phase (before inflection). Prediction of inflection point with current data for these regions is very difficult, based on the “shape” from other countries, the non-linear mixed effect model provides sensible predictions albeit large error intervals (in the left panel of Figure ~\ref{fig:fitm1}). African region is predicted to have the largest total outbreak size of 3.9 million cases (2.2 to 6 million), and the inflection will come around September 13, 2020. Note that African region here includes all the African countries affected by COVID-19. Brazil and India are predicted to have a similar final outbreak size of around 2.5 million cases (1.1 to 4.3 million), with the inflection points arriving June 23 and July 26, respectively.  The epidemic in Brazil has entered the rapid growth stage and is increasing quickly in the number of the cumulative confirmed cases. The epidemic in India has the similar situation with Brazil, but the growth is predicted to have one-month delay. 
  
\section{Conclusions}

In these developing areas, community transmission and spread of COVID-19 has been ongoing for a while, but large-scale testing is only available until recent weeks. The rapid increase in confirmed cases can be attributed to the increasing level of testing coverage. Currently African region has a relatively small number of confirmed infection ($<$150k), but it might continue to increase to the level estimated in this study if not controlled effectively \citep{pearson2020projected}.\\

Although USA passed the inflection point around April 18,2020, its growth rate (or the number of daily new confirmed infections) failed to achieve rapid decrease. Russia has a similar epidemic curve to USA. This indicates that there are no effective intervention strategies in these countries to further curb the ongoing transmission of COVID-19. Form the validation figures in supplementary material, our model tends to underpredict the daily new confirmed infections in USA and Russia during the later phase of outbreak after inflection points; thus, USA and Russia may have a larger outbreak size than our estimates. Likewise, if these new epicentres (i.e., Brazil, India and African region) could not conduct effective measures to mitigate the spread and transmission, there will be a large number of confirmed infections even after the inflection points. In summary, our model predicted the inflection point and maximal outbreak size in Brazil, India, and African region; these regions might take over USA in terms of outbreak size at the end. \\

This work only fits the growth curves to reported number of confirmed infections, and incorporating localized intervention policies and behavior parameters will improve the performance of fitting and prediction. Furthermore, in this work we test three different power transformation of cumulative number of cases, however, they are likely not the optimal choice for all countries and regions. Other power transformations are worthwhile to test and possibility of using different transformation to different countries/regions should be explored under the same framework of mixed-effect model.   
\\

\section*{CRediT authorship contribution statement}

{\bf Qibin Duan:} Conceptualization, Methodology, Formal analysis, Validation, Writing - original draft, Writing - review \& editing. {\bf Wu:} Methodology, Writing - review \& editing. {\bf Jinran Wu:} Writing - review \& editing. {\bf You-Gan Wang:} Conceptualization, Methodology, Writing - original draft, Writing - review \& editing.

\section*{Declaration of Competing Interest}

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

\section*{Acknowledgement}

This work was supported by supported by the Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS), under grant number CE140100049.





\bibliographystyle{apacite}

