Modeling Dinophysis in Western Andalucía using an autoregressive hidden Markov model

Dinophysis spp. can produce diarrhetic shellfish toxins (DST) including okadaic acid and dinophysistoxins, and some strains can also produce non-diarrheic pectenotoxins. Although DSTs are of human health concern and have motivated environmental monitoring programs in many locations, these monitoring programs often have temporal data gaps (e.g., days without measurements). This paper presents a model for the historical time-series, on a daily basis, of DST-producing toxigenic Dinophysis in 8 monitored locations in western Andalucía over 2015–2020, incorporating measurements of algae counts and DST levels. We fitted a bivariate hidden Markov Model (HMM) incorporating an autoregressive correlation among the observed DST measurements to account for environmental persistence of DST. We then reconstruct the maximum-likelihood profile of algae presence in the water column at daily intervals using the Viterbi algorithm. Using historical monitoring data from Andalucía, the model estimated that potentially toxigenic Dinophysis algae is present at greater than or equal to 250 cells/L between < 1% and >10% of the year depending on the site and year. The historical time-series reconstruction enabled by this method may facilitate future investigations into temporal dynamics of toxigenic Dinophysis blooms.


Introduction
Following the identification of Dinophysis fortii as the causative agent of shellfish poisoning outbreaks in 1976 and 1977 in northeastern Japan, there has been interest in understanding the dinoflagellates of the genus Dinophysis (Yasumoto et al. 1978(Yasumoto et al. , 1980FAO and WHO 2016). Dinophysis has been shown to be toxigenic in a variety of environments (Fux et al. 2011;Mafra et al. 2013;Gao et al. 2017), and a global distribution system for these toxins is generally recognized (FAO and WHO 2016) as they occur in a multitude of different habitats. To date, 10 species of Dinophysis have been shown to produce one or two major types of lipophilic toxins, okadaic acid (OA) and its dinophysistoxins derivatives, and pectenotoxins (Reguera et al. 2012(Reguera et al. , 2014. We collectively refer to the set OA, its derivatives, and pectenotoxins as diarrhetic shellfish toxins (DST). While all DST are known to be harmful, OA is particularly important to study due to the knowledge gap about its potential for understudied chronic disease associations (IOC et al. 2005) as well as its likely tumor promotion properties (Suganuma et al. 1988;Fujiki et al. 2018;Valdiglesias et al. 2013). The potential for low-dose chronic exposures of human populations to OA (along with pectenotoxins) is environmentally plausible because the toxins may persist in the water for extended periods of time (Blanco et al. 2018;Pizarro et al. 2009) and the toxin has been found in the absence of Dinophysis (Fernández et al. 2019).
In 1994 the autonomous government of Andalucía implemented a phytoplankton monitoring system (Fernández et al. 2019). Since the late 1990s, multiple species of Dinophysis (e.g., D. acuminata complex, D. caudata, D. acuta, and D. fortii) have been detected in most of the sampled areas (Fernández et al. 2019). Levels of toxins have been found in various species of edible shellfish (e.g. Callista chione and Venus verrucose), sometimes exceeding the legal limit of 160 μg/kg (e.g. Donax trunculus, Chamelea gallina, Mytilus galloprovincialis, and Cerastoderma edule) (Fernández et al. 2019). The concentration of DST in a shellfish at a given time is a function of shellfish-specific features such as the rate of uptake, biotransformation, and elimination of each particular toxin (Reguera et al. 2014), but also reflects the environmental dynamics of toxigenic algae. Previous efforts to model the dynamics of toxigenic Dinophysis have had limitations.
Artificial neural networks have been applied to the coast of Huelva in Andalucía (Velo-Suárez and Estrada 2007); however, that modeling effort used the last 5 weekly D. acuminata concentrations as the only input variables to predict the upcoming week's concentration. Kulawiak (2016) used GIS and advanced very-high-resolution radiometer data to detect algae blooms in the Baltic sea, but this did not leverage toxin measurements. To our knowledge, it has been an unaddressed modeling challenge to account for the fact that DST can persist in water for extended periods of time, while also giving interpretable model parameters. Hidden Markov Models (HMM) have proven to be an effective modeling tool (Zucchini and Macdonald 2009). In addition, they have specifically been used to study algal blooms in other contexts and provide a potential scaffold for an improved Dinophysis model. Rousseeuw et al. (2015) used a hybrid HMM to detect and understand the dynamics of phytoplankton blooms in France using data on nutrients and water characteristics, but lacking direct data on algae. In the freshwater harmful algal bloom modeling field, Jiang et al. (2016) employed a continuous HMM alongside principal component analysis of water quality parameters and nutrient to forecast microcystins. Kim et al. (2020) analyzed chlorophyll-a concentrations, a metric to understand eutrophication, in the Nakdong river of South Korea using a continuous HMM. All of these approaches model a multivariate outcome of water quality and chemical parameters using a HMM with an unknown number of states. Rousseeuw et al. (2015) and Jiang et al. (2016) reduce the dimension of the multivariate outcome by clustering and principal components, respectively. Kim et al. (2020) models the spatial distribution of chlorophyll-a conditional on the latent state. In our work, we use a 2-state HMM to directly infer whether there is algae in the water column over time. This is important because we want to reconstruct a daily assessment of this variable for health surveillance. Autoregressive HMMs on the other hand, were initially developed for speech recognition Rabiner 1986, 1985). They have since been applied to various issues in recent years (Urban et al. 2020;Bartolucci et al. 2014;Shannon and Byrne 2010), but have not been used to model algae.
This paper presents a first-order autoregressive HMM approach to modeling potentially toxigenic Dinophysis in western Andalucía with the purpose of reconstructing the maximum-likelihood profile of whether algae were absent or present in the water column above a threshold count (e.g., ≥ 500 cells / L) at daily intervals using incomplete time-series historical data on both DST levels and algal counts. We model the presence/absence of algae in the water column by using algae cell counts from water column samples and DST measurements (in μg of OA equivalents per kg) from shellfish gathered from the regional government's website. Then, using the estimated model parameters we can reconstruct indicators of algae in the water column at a daily interval, even when data is often missing. Since OA can remain in the water for extended periods of time, it is important that we allow for serial dependence in the model formulation. Specifically, we introduce autoregressive dependence in the observed DST measurements after accounting for the hidden Markov structure. The forward-backward algorithm needed to be adapted (Stanculescu et al. 2014) for computing E-step calculations in order to implement the EM algorithm for maximumlikelihood in this setting. We assume a first order autoregressive model as algae blooms erupts quickly and this assumption provides a useful framework over a longer time horizon to capture the quickly moving event. Section 2 presents an in-depth explanation of the model while Sect. 3 explains the estimation procedure, along with the adapted forward-backward algorithm that can account for both missing data and dependence on previously observed states. In Sect. 4 we talk about three simulations which compare estimation with different amounts of missing data. In the 5 th section we apply our model to the Andalucía data and discuss the estimated parameters and we consider our results in Sect. 6. below. Also let S = S 1 , …, S T , and similarly for X and Y. S, X and Y all have the same follow up with equally spaced daily observations. We assume that algae in the water column can be modelled by a Markov chain, where S t is either 0 or 1 depending on if algae is absent or present in the water column. We define the notation r S t to be the probability of starting in state S t at time t = 1 and P S t − 1 S t to be the probability of transitioning from state S t − 1 to state S t at time t. Specifically, S t ∼ MC p 01 , p 10 , r 1 where p 01 denotes the probability of initiating water column algae over a day, p 10 indicates the probability of ending the episode over a day, and r 1 denotes the probability of being in the algae state at time t = 1. For use in calculations, we also define p 00 as the probability of algae remaining out of the water column over a day, p 11 as the probability of algae continuing to remain in the water column over a day, and r 0 as the probability of not being in the algae state at time t = 0 To model the algae cell counts from the water sample, we chose a negative binomial as it can account for overdispersion in the count data. X t takes on the integer count of algal cells in a sample of water and is conditional on S t . Thus the possible values of X t are 0 ≤ X t ≤ ∞. Although X t is directly observed, it is also subject to measurement error. Due to the spatial heterogeneity of the algae and the water sampling technique used (Fernández et al. 2019), excess zeros are possible based on the specific latitude and longitude sampled. When X t = 0 we can't be sure if the sample missed the algae or if algae is truly absent from the water column, however when X t > 0 we know for sure that algae is present in the water column. Reversing this, when S t = 0 then X t = 0 but when S t = 1 then it is possible but not necessary that X t be greater than 0. We chose to model the absence and presence of algae (quantified as above or below the 50 cells/L detection limit), however other thresholds can easily be chosen. Appendix A discusses the implications when two additional thresholds are considered. When algae is present in the water column, we model X t with a negative binomial distribution with mean μ a and size k a where E X t | S t = 1 = μ a and V X t | S t = 1 = μ a + μ a 2 k a . This relationship can be concisely summed up as: We discretized the DST measurements in to four states as a large proportion of values are below the quantification limit and are non-normally distributed on all commonly considered transformed scales. Discrete measurements also help with the computational feasibility of the method. Therefore, let Y t * represents the continuous toxin measurements and let Y t represent the discretized measurements as follows: These specific cutoffs were chosen as the limit of detection is 40μg and fisheries must close at 160μg while 100μg lies halfway between the two other constraints. The DST states follow an ordinal logistic regression model with two regression parameters and three intercept parameters. The full model is, where 0 ≤ c ≤ 2, α c is the intercept parameter, β 1 is the regression parameter for the last toxin measurement, and β 2 is the regression parameter for the current Markov chain state.
The regression coefficients can be interpreted as: there is e β 1 times the odds of Y t = c + 1 compared to Y t = c with each increase by one in Y t − 1 and when S t increases from 0 to 1 there is e β 2 times odds of Y t = c + 1 compared to Y t = c. The dependence of Y t on Y t − 1 creates issues when Y t − 1 is missing, however these problems will be dealt with in Sect. 3.
The joint distribution of the latent water column algae state, the observed water sample algae count, and DST measurement can be calculated by multiplying the following three components: (1) the probability of being in a water column algae state, (2) conditional on the water column state, the probability of the algae cell count from the water sample, and (3) conditional on the water column state and the last DST state, the probability of the current DST state. The complete-data joint distribution can be written as: We assume that Y 0 = 0 because most values are of Y t are zero. Additionally as a sensitivity analysis we ran our analysis where Y 0 = 1, Y 0 = 2, and Y 0 = 3 and the results did not change. The estimation procedures are described in the next section. Parametric bootstrap standard errors were calculated by simulating data 500 times per site using the estimated parameters. Bootstrap standard errors were then calculated for each parameter by calculating the standard deviation of the 500 samples. Finally, we apply the Viterbi algorithm to reconstruct the highest likelihood hidden state path.

Estimation
We will now introduce the procedure for estimation assuming no missing data. Define the indicator variable Z such that This indicator function is critical for later computations and can be used with variables other than S t . We adopt the notation where S t refers to the random variable, while s t refers to a possible value of the random variable S t . We use this notation across all random variables.
This indicator function is equal to 1 when the random variable of choice is equal to a specific realization of the random variable. We can then rewrite the complete-data joint distribution as: As S is not observable, to maximize this likelihood directly we would have to iterate over every possible value, thus calculating Maximizing this directly becomes intractable as T increases and as T = 2177 days, it is not feasible for our application. By using the expectation-maximization (EM) algorithm we can maximize this likelihood in a timely manner by alternating between an expectation and a maximization step, converging at the estimated parameters. The expectation step calculates the following complete-data log likelihood: where the expectations can all be calculated using the forward-backward algorithm (Baum 1972).

Estimation with missing data
Often times some parts of the observable data are missing. As shown in Fig. 1 most observations are missing (other areas are shown in Fig. 2). This time frame was chosen as spring and summer are often when most algae blooms occur. In our example, when X t is missing we simply leave out the second line of the likelihood calculation, however when Y t is missing a more complicated method is required. As the current DST state depends on the last DST state, when Y t is missing we must account for it to calculate the probability of Y t + 1 . By conditioning over all possible DST states for Y t , we can calculate the probability of Y t + 1 . The complete data joint distribution, accounting for missing data, can be written as where the complete-data log likelihood is now

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript To account for the missing data and the dependency in the emissions distribution we use the adapted forward-backward algorithm from Stanculescu et al. (2014), however our application has a bivariate rather than univariate emissions distribution. For the maximization step, we maximize equation (9) given the E step calculations. The E step calculations are hard to calculate so we use the Forward-Backward algorithm described in the next section. We consider convergence to occur when the log likelihood increase between iterations is less than 0.01.

Forward-backward algorithm
To account for missing toxin values, we will keep track of every possible DST state value. Assume that Y t is missing. By calculating the probability of all possible DST states at time t, we can then calculate the probability of Y t + 1 . We redefine the indicator variable Z to account for the scenario of missing data. Let The forward quantity is: The indicator function, not present in the standard forward quantity, allows us to incorporate different possible values of Y t when Y t is missing. If Y t is observed the forward quantity is zero except when y t = ω. However, when Y t is missing, ω corresponds to a possible DST state value at time t. This quantity is calculated recursively by The backward quantity is defined as where, similarly to the forward quantity, if Y t is observed the backward quantity is zero except when y t = ω and if Y t is missing ω corresponds to a possible DST state value at time t. It is also calculated recursively:

Calculating expectations
The expectations from the complete-data log likelihood are calculated as follows:

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript In equation 15 and 16 the function g x t | μ a k a is the probability density function of negative binomial distribution with parameters μ a and k a , calculating the probability of x t .

Simulation
We examine the performance of our proposed method by analyzing three simulations with varying amounts of missing data. We simulated data sets with no missing data, one-third of the data missing, and 85% of the data missing. For each category, 500 data sets were generated with the same follow up length as the data (2177 days). The simulated data structure corresponds to the application presented in our application section. These three amounts of missing data were chosen as they account for a wide variety of scenarios while also testing this specific application, which has (depending on the site) at most 83% of the data missing. By varying the level of missingness, we can measure how well our method preforms at recovering the true parameters with different levels of information available. It should also be noted that this is especially important to test for the DST measurements. When there is no missing data the estimation for the DST is straightforward, however with missing data the adapted forward backward algorithm explained in the estimation section is needed. With more missing data there are longer times between observations, meaning more reliance on the proposed adaption to account for Y t when calculating the probability of Y t + 1 .
For the simulations that have no missing data, the estimates are extremely accurate with minimal standard errors. Table 1 contains the simulation estimates for the three different levels of missing data, and shows that our method is extremely accurate at recovering the true parameters regardless of missing data. Additional computation time is required when our method encounters missing data as all possible values of the last DST state are iterated over. Thus, the time needed for our method scales linearly with the amount of missing data.
Even with most of the data missing, our method accurately estimates the parameters. However, as more data is missing, the standard errors increase. While this increase is quite small when one-third of the data is missing, it is much larger when 85% of the data is missing. The standard errors for the 33% missing simulation roughly double when compared to the 0% missing simulation, however the standard errors increase by a factor ranging from 9 to 37 when comparing the 85% and the 0% missing simulations. This can most easily be seen in the third row of Table 1 for μ a . The estimate itself is accurate across the three levels, however the standard error when 85% of the data is missing is extremely large, even compared to the standard error when 33% of the data is missing. Despite this high variability, there is no relationship between the initial and estimated value in the simulations.

Dataset description
We illustrate our method on data gathered from the regional government of Andalucía's website (Zonas de producción. (n.d.) xxxx). The Andalucían government established a phytoplankton toxin monitoring program for shellfish in 1994 to help deal with the recurrent blooms of Dinophysis that are linked to DSP outbreaks (Bouza and Aboal 2008). We used data on toxin levels, measured in μg of OA equivalent/kg, sampled from the bivalve Donax trunculus in the time frame from January 2015 to December 2020. The follow-up length is 2,177 days. Toxin levels were calculated as specified by Yasumoto et al. (1984) and liquid chromatography-tandem mass spectrometry was used as the chemical analysis technique (Fernández et al.2019).Water column samples used to calculate algae cell counts were gathered using a 10-meter-long weighted plastic hose. 25mL water samples from sedimentation chambers were then used to extrapolate the number of cells/L (Velo-Suárez and Estrada 2007). Data from eight geographical sites (areas 1, 2, 3, 4, 5, 6, 7, and 8) were analyzed separately. Areas 7 and 8 were recorded as a single area until May 2018 and were then split into distinct areas; we analyzed each site in a separate model. Table 2 contains some summary statistics about the data.

Results
Our HMM has two states representing the presence or absence of potentially toxigenic Dinophysis algae in the water column at a concentration exceeding a threshold (e.g., ≥ 500 cells/L). Both the initiation (p 01 )and termination (p 10 ) probabilities are low as can be seen in Table 2, indicating a tendency for algae to stay in or remain out of the water column for a number of days (corresponding to the 1 or 0 state of the HMM). Despite the minor differences between each of the different sites, there is broad homogeneity among the sites with the initiation probabilities being slightly lower than 10% and most termination probabilities hovering just above 10%. Within each area, initiation probabilities were lower than their corresponding termination probabilities. Although the state path of the hidden Markov model is unobserved, it is an important metric to recover as it can be useful in determining long-term changes in algae and can have implications for the effects of climate change. We reconstructed the hidden state paths using the Viterbi algorithm, producing the path with the highest likelihood. Figure 3 shows a visualization of the Viterbi path for area 1 in 2016 (other areas are shown in Fig. 4). Using the Viterbi path we can then calculate different summary statistics for algae presence/absence across time. As shown in Fig. 5, distinct trends can be seen within each year and across years. For instance, for areas 1-7 the earlier and later years have a higher proportion of algae in the water column when compared to the middle years. The proportion of days with algae was estimated to be at least 54.52 % and 61.1% for 2015 and 2019, while in 2017 it was estimated to be at most 48.22%.
As noted previously, we modeled the algae from the water sampled with a negative binomial model when S t = 1, and assumed that there cannot be any algae in the water sample when S t = 0. When there is algae in the water sample we can assume that S t = 1 because otherwise it would not be possible for there to be algae present. On the contrary, we cannot draw any conclusions when there is not algae in the water sample. This is the case because the algae cell count from the water sample is serving as an observable representative of algae in the water column with measurement error. Because the algae is not distributed evenly by either latitude, longitude, or depth, the water sample may not accurately capture whether algae is in the water column. As noted in the methods section, for this application we consider algae to be present in the water sample when it can be detected (a threshold of 50 cells per liter). We examine the consequences of higher thresholds in appendix A. Areas 3 and 5 have a mean parameter (μ a ) around 230-240, while areas 1, 2, 4, and 6 have a higher mean parameter ranging from 265 to 290. Area 8 has a larger mean parameter of 320, and area 7 has a significantly larger mean around 440. Barring area 7, the size parameter (k a in Eq. 1) is above 1 indicating that the negative binomial model is essential to help account for over-dispersion. As area 7 has the largest mean parameter and smallest size parameter (see expression 1), this leads to larger, more variable predictions for area 7.
The continuous DST measurements were discretized to form four (0 to 3) DST states, which follow an ordinal logistic model. The continuous measurements are highly skewed as the limit of quantification is 40 μg of OA equ/kg. Binning the continuous measurements less than 40 μg of OA equ/kg reduced the number of distributional assumptions. Unlike the algae cell count from the water sample that only depended on the current Markov chain state, we assume that the DST states are dependent on the current Markov chain state as well as the last DST state. This additional dependency is necessary in the emission distribution as major components of DST have been found to be very stable in the water column after a Dinophysis bloom (Blanco et al. 2018;Pizarro et al. 2009). This dependency cannot be estimated using the standard forward-backward algorithm as the standard HMM assumes that the observed states are all conditionally independent given the current latent state. In our model, the current observed state is dependent on the previous observed state and the current latent state, violating this assumption. Instead, by using procedures developed for autoregressive HMMs (Stanculescu et al. 2014), we can incorporate this dependency into the estimation procedure. We believe that a first order autoregressive model is applicable as algae blooms are rapid events. By having a shorter time dependency we are better able to model these events.
Our ordinal logistic model has five parameters: β 1 is the effect of the DST state at time t − 1, β 2 is the effect of the current Markov state, and α 0 , α 1 , and α 2 are the intercept parameters.
The effect of the last DST state is additive in relation to the log odds of the probability of the current DST state such that the effect of the last DST state is β 1 when Y t − 1 = 1, 2 × β 1 when Y t − 1 = 2, and 3 × β 1 when Y t − 1 = 3. Table 3 contains all parameter estimates for the eight different areas along with their standard errors.
Despite the somewhat high standard errors for the regression coefficients, the probabilities themselves have a low standard error. The probabilities and standard errors for area 1, along with the other areas, are shown in Table 4. Interestingly, we can see the difference in predictive power between S t and Y t by looking this table. For each area, the difference between the left and right halves of the table is not nearly as drastic as the difference between the rows, indicating that the autoregressive effect on Y t is indispensable to this model.

Discussion
In this paper we have focused on the historical reconstruction of an incomplete time-series by developing a model that recreates the most likely pattern of Dinophysis spp. algae occurrence at each of the eight different sites on a daily timescale, using a HMM with extensions to account for challenges inherent to the data. DST measurements were highly autocorrelated, even after accounting for the hidden states of the HMM, violating one of the standard assumptions of HMMs ). However, by using an autoregressive HMM we are able to model this. The sampling frequency of the monitoring program resulted in large amounts of missing data (at most 83%). Furthermore, the distribution of DST was skewed and left-censored at the assay limits of quantification, 40μg of OA equivalents per kg (Zonas de producción.(n.d.) xxxx). We addressed these challenges with an advanced HMM that included a bivariate emissions distribution with a negative binomial distribution for algae counts and ordinal autoregressive model for the serial DST measurements. We showed with simulations that the approach accurately estimated the parameters even with extensive missing measurements.
This paper presents a modified forward-backward algorithm in an EM context from Stanculescu et al. (2014) with an additional observed variable applied to data from a phytoplankton toxin monitoring program in western Andalucía. This generalized form allows us to estimate a model with both dependence in the emissions distribution and missing data. In our application, DST states are dependent on both the current Markov state as well as the last DST state. The proposed method works by keeping track of all possible DST states (when the DST state is missing) in the forward-backward algorithm. We can then condition on and sum over the most recent DST state to calculate the probability of the current DST state. Although this does lead to additional computation complexity, the time needed scales linearly with the amount of missing data and is still feasible when nearly all of the data is missing.
We applied this method to 2,177 days of algae water samples and DST data from eight geographical sites with dates ranging from January 2015 to December 2020. Despite the long stretch of time covered in the study, most days had no recorded data. Although the data available varied by site, it ranged between 377 (17%) and 524 (24%) days with recorded measurements out of the total 2,177 days. Although HMMs have not been applied to this problem in this area before, our application of this method shows that HMMs are capable of modeling complex processes that don't necessarily conform to the standard assumptions in the presence of large amount of missing data. Running our method on western Andalucía phytoplankton monitoring data we see that accounting for the last DST state requires the additional complexity of an autoregressive HMM.
One of the major advantages from our method is that we are able to reconstruct paths of the latent variable on a daily interval using historical time-series data in the presence of intermittent measurements and measurement error. Rather than forecast the future, our method focuses on predicting whether algae were absent or present in the water column for every day in our data set. This is useful when trying to identify long term algae trends for the different areas across time as well as for health surveillance. The Viterbi algorithm is ideal for computing estimates of the entire sequence of latent states. These sequences can later be used in downstream analyses that examine the relationship between toxicity and diseased risk. By aggregating these sequences, termed Viterbi paths, we are able to identify long term trends across years.
The proposed hidden Markov model makes a number of parametric assumptions including that the unobserved states follow a first-order Markov model and that the observed DST data follow a first-order autoregressive process after accounting for the HMM structure. We believe that these are reasonable assumptions since DST remains in the water for extended periods of time and algae blooms are rapid events. Latent state estimation should not be sensitive to small departures from these underlying assumptions. Therefore, we presume that the AR(1)-HMM framework adequately describes the biological process.
In the future, our method can be applied to other types of monitoring program data as well. Because most monitoring program data contains missing values, accounting for the temporal autocorrelation that is often present is not straightforward. Our method can adequately handle both complications simultaneously while also creating historical timeseries reconstructions. Using our method, we are also able to relate two separate processes together while we impute the maximum-likelihood profile of the variable of interest.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.
Matthew O. Gribble is an Associate Professor of Epidemiology at the University of Alabama at Birmingham. His research program focuses on substantive issues in Oceans and Human Health and drinking water epidemiology, and methodological challenges thereto pertaining.

Appendix A: Varying the algae threshold
This appendix examines the ramifications of adjusting the algae detection threshold. Noted previously in prior sections, when any amount of algae is present in the water sample we consider algae present in the water column. Although this metric is useful for specific scenarios, such as quantifying chronic exposure, it cannot single out larger events like harmful algae blooms. Because algae is present year-round, when a lower threshold is chosen larger algae events are mixed in with smaller algae events. As this interpretation may not be sufficient for a study of harmful algae blooms, we consider three different thresholds. By raising the threshold for what we consider algae presence to be, we can study larger algae events beyond presence. The two additional thresholds examined are at 250 and 500 cells/L. 500 cells/L was established by the Andalucía HAB monitoring program as a critical threshold, and 250 cells/L was chosen as a halfway point.
As can be seen in the three following heat maps, when the algae threshold is increased our model predicts fewer days with algae in the water column. Thus, by raising the threshold our model is able to be more precise and pick out days with larger algae events. The drastic decrease in predicted algae positive days as the threshold increases indicates that under the original 50 cells/L threshold, most of the algae positive days had a low predicted amount of algae. Figure 6 shows the predicted percent of days with algae presence in the water column above the three different thresholds, averaged over months, for six years and eight different sites. By summing over months, we can look at yearly trends. When the threshold is 50 cells/L this figure contains the same information as Fig. 5. The U-shaped pattern present at 50 cells/L, where the earlier and later years studied had higher predicted algae presence in the water column when compared to the middle years, is still prevalent for both the 250 and 500 cells/L threshold for areas 1-5. Areas 6 has a higher predicted algae presence above the threshold in the water column in the beginning and then slowly decreases from 2015 to 2020 for both the additional thresholds while area 7 and 8 vary by year. Heat map of predicted percent of days with algae presence in the water column above the three different thresholds. The x axis is years studied and the y axis is the different threshold levels. Each of the eight boxes represents a different area. For this figure, months were averaged over Figure 7 shows the predicted percent of days with algae presence in the water column above the three different thresholds, averaged over years, for each month and the eight different sites. This heat map focuses on monthly trends as it sums over years. For the 50 cells/L threshold, the number of predicted days where algae is above the threshold in the water column peaks between April and August, often occurring in May. Areas 1-6 all had the highest percentages during the spring and summer months, while areas 7 and 8 also had a high predicted percentage in January. Looking at the 250 cells/L threshold, the highest percentages occurred between April and October, with most peaks happening in April. Areas 1, 2, 6, and 7 also had at least one winter month with a high predicted percentage, however all areas had the most predicted positive days during the spring and summer months. Finally, for the 500 cells/L threshold areas 6 and 7 had the highest predicted percent of days with algae above the threshold in January, however area 6 also has high predictions for April and May. Area 2 has the highest predicted percent of days in February but is closely followed by June and July. The rest of the areas have the highest predicted percent of days between April and July. Heat map of predicted percent of days with algae presence in the water column above the three different thresholds. The x axis is months and the y axis is the different threshold levels. Each of the eight boxes represents a different area. For this figure, years were averaged over Figure 8 shows the predicted percent of days with algae presence in the water column above the three different thresholds, averaged over areas, for each month and year. This heat map sums over geographic differences and lets us examine overall trends, looking at western Andalucía as a whole. Across all years and thresholds (except for the 500 cells/L threshold in 2020), there is a higher predicted percent of days of algae above the threshold in the water column during the spring and summer, although some years had increased algae during the winter as well. Heat map of predicted percent of days with algae presence in the water column above the three different thresholds. The x axis is months and the y axis is the different threshold levels. Each of the six boxes represents a different year. For this figure, areas were averaged over Figures 9 and 10 are similar to Figs. 2 and 3 and show the interpolated Viterbi path for 2016 when the threshold is 250 and 500 cells/L, respectively. Again, as we increase the threshold, the qualification for what counts as algae presence in the water column is harder to meet. Therefore, fewer days are predicted to have algae in the water column, however the predicted algae events are larger. The decrease in number of days that are estimated to have algae in the water column when the threshold is increased can also be seen in these two plots as the Viterbi path in Fig. 9 passes through the 1 state more often than the Viterbi path in Fig. 10 does.  Decoded water column state (S t ) path using the Viterbi algorithm for 2016 for areas 2-8.
The line is the decoded path while the dots indicate absence and presence of observed algae. When algae counts are above the threshold of 250 cells/L, indicated with a black dot at 1, S t must equal one Decoded water column state (S t ) path using the Viterbi algorithm for 2016 for areas 2-8.
The line is the decoded path while the dots indicate absence and presence of observed algae. When algae counts are above the threshold of 500 cells/L, indicated with a black dot at 1, S t must equal one  Decoded water column state (S t ) path using the Viterbi algorithm for 2016 in area 1. The line is the decoded path while the dots indicate absence and presence of observed algae. When algae counts are above 0 we know that S t must equal one, however when algae counts are 0 we cannot say anything about S t . For the figure, this is why when we observe algae (indicated by a dot at 1) the path (the line) must pass through it, but when we do not observe algae (a dot at 0) the path may or may not pass through it Proportion of predicted days with algae presence in the Viterbi path across all areas and years Decoded water column state (S t ) path using the Viterbi algorithm for 2016 for areas 2-8.
The line is the decoded path while the dots indicate absence and presence of observed algae. When algae counts are above the threshold of 50 cells/L, indicated with a black dot at 1, S t must equal one Table 1 Comparison of truth and estimated parameter from three different simulations with 0%, 33%, and 85% missing data  Summary of basic information about the collected algae count and DST state data. DST continuous measurements were discretized into four states: 0 (≤ 40μg of OA equ/kg), 1 (40 ≤ 100μg of OA equ/kg), 2 (100 ≤ 160μg of OA equ/kg), and 3 (≥ 160μg of OA equ/kg)   Table 3 Markov transition probabilities for algae presence/absence in the water column, mean and size negative binomial parameters for water sample algae count, and ordinal logistic regression coefficients for DST state