Appendix 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 *t _{0 }*and

*C*are taken to be some initial time and concentration respectively.

_{0}For fixed *θ, k, t _{0}*, and

*C*

_{0},

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, t _{0}*, and

*C*as

_{0},*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(*t _{c}*) = 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 *t _{c}*. Constraining the model so that it is real valued
for all time greater than

*t*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,

_{c}*t*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

_{c}*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 *t _{P}*, we re-parameterized the model to
include

*C*as a parameter rather than

_{P}*C*. It is possible to show that

_{0}

as long as *t _{P}* is in the interval where the model is real
valued. Although we chose

*t*to be 1975, we found that in practice, we had better convergence when we used the second form of the equation with

_{0}*t*taken to be 1990. We speculate that this could be because 1975 was relatively close to

_{P}*t*and the unbounded behavior of the model as it approaches

_{c}*t*could effect the convergence of the optimization routines.

_{c}