HomeAppendix 1: Analytical and Numerical Issues in Fitting the Mixed-Order Model
In Stow et.al.(1999), the mixed-order model was given as a generalization of first-order exponential decay (growth). Our initial assumption was that the mixed-order model had qualitative behavior similar to that of exponential decay. It seemed reasonable to think that for different parameter regimes the behavior of the model would be such that it was decreasing (increasing) for all time. Instead, we found that the model has a variety of behaviors that can give nonsensical predictions.
Looking at the model given by Stow et al., we have:
where t0 and C0 are taken to be some initial time and concentration respectively.
For fixed θ, k, t0, and C0,
must be non-negative for a given t in order for a modeled value at that time to be real (except for special cases of θ). If θ≠1 and k≠0, this cannot be true for all time. To see this, note that in order for it to be true, it must be the case that
(1)
With fixed θ, k, t0, and C0,
as t→∞, the right hand side of (1) is unbounded in either
the positive or negative direction depending on the sign of
kּ(1-θ). As t→-∞, the right hand
side of (1) is unbounded in the opposite direction, so it follows that (1)
holds on only one of the two intervals
and
, with
.
C(tc) = 0 if θ<1 and is unbounded if θ>1.
In order to fit the model and use it for prediction, it was necessary to constrain the parameters in such a way that predictions have real values. The simplest constraint is to make kּ(1-θ)≤0. This constraint guarantees the model is real valued for all time greater than tc. Constraining the model so that it is real valued for all time greater than tc is more than necessary. We generally found a significantly better model (as measured by the maximum likelihood) when we only constrained the model to be real valued on a specified finite interval (that is, tc was forced to be outside the finite interval of interest, but the type of interval was not forced as with the previous constraint). A consequence of this relaxation is that it is not generally possible to make predictions for all time. However, given the ad hoc nature of this model and the time scale on which we would expect it to be accurate, we think this limitation is acceptable. The key to generating these constraints is the fact that the right hand side of (1) is monotone with respect to t. This monotonicity implies that constraining the endpoints to satisfy (1) will guarantee that the entire interval satisfies (1).
There were numerical issues that complicated the process of fitting the model to the data. A Newton-Raphson based approach to maximize the likelihood function was not reliable because the second derivative matrix was ill-conditioned with our data. Instead, we used a much slower derivative free algorithm called the downhill simplex method developed by Nelder and Mead (1965). Even this algorithm failed to converge to the optimal solution when maximizing simultaneously over all the parameters (though it was usually not far off). This fact was discovered when the profiling approach we used in prediction converged to better fitting models.
In order to find a profile-likelihood based confidence interval for the mean in a given year tP, we re-parameterized the model to include CP as a parameter rather than C0. It is possible to show that
as long as tP is in the interval where the model is real valued. Although we chose t0 to be 1975, we found that in practice, we had better convergence when we used the second form of the equation with tP taken to be 1990. We speculate that this could be because 1975 was relatively close to tc and the unbounded behavior of the model as it approaches tc could effect the convergence of the optimization routines.