Honey Bee Colony Population Daily Loss Rate Forecasting and an Early Warning Method Using Temporal Convolutional Networks

The population loss rate of a honey bee colony is a critical index to verify its health condition. Forecasting models for the population loss rate of a honey bee colony can be an essential tool in honey bee health management and pave a way to early warning methods in the understanding of potential abnormalities affecting a honey bee colony. This work presents a forecasting and early warning algorithm for the population daily loss rate of honey bee colonies and determining warning levels based on the predictions. Honey bee colony population daily loss rate data were obtained through embedded image systems to automatically monitor in real-time the in-and-out activity of honey bees at hive entrances. A forecasting model was trained based on temporal convolutional neural networks (TCN) to predict the following day’s population loss rate. The forecasting model was optimized by conducting feature importance analysis, feature selection, and hyperparameter optimization. A warning level determination method using an isolation forest algorithm was applied to classify the population daily loss rate as normal or abnormal. The integrated algorithm was tested on two population loss rate datasets collected from multiple honey bee colonies in a honey bee farm. The test results show that the forecasting model can achieve a weighted mean average percentage error (WMAPE) of 17.1 ± 1.6%, while the warning level determination method reached 90.0 ± 8.5% accuracy. The forecasting model developed through this study can be used to facilitate efficient management of honey bee colonies and prevent colony collapse.


Introduction
Honey bees (Apis mellifera) are essential for global food production [1]. In the last decades, some regions of the world have suffered from a significant loss in the number of active honey bee colonies [2]. The decline of honey bee colonies is partly attributed to a phenomenon known as Colony Collapse Disorder (CCD) [1,2]. To maintain healthy and thriving honey bee colonies, it is essential for beekeepers to regularly monitor the status of the bee hives. Currently, most beekeepers manually estimate the health of a honey bee colony in a relatively time-consuming and imprecise approach by beehive inspection and honey bee flight activity observation at the hive entrance. Automated monitoring system are highly desirable to have a rapid and quantitative access to the health status information of honey bee hives.
Various systems for automatic in-hive and hive-entrance monitoring of honey bee activities have been developed. Reported systems use infrared sensors [3], radio frequency identification (RFID) tags and readers [4,5], and imaging sensors [6]. In our recent work, we described an embedded image monitoring system for automatic and highly reliable its ideal hyperparameters and determining its best input features; (3) implementing an outlier detection method, such as an isolation forest algorithm, to classify the forecasted daily population loss rate as normal or abnormal; and (4) evaluating the developed TCN model and isolation forest algorithm based on test datasets. The developed method can be used by beekeepers, entomologists, and agronomists to obtain fast, essential information about the status and abnormality in honey bee colonies.

Data Collection
An automated system was used for recording the traffic at the beehive entrance and monitoring local environmental conditions. Each set of the system included an NVIDIA Jetson TX2 (NVIDIA Corporation, Santa Clara, CA, USA) for image monitoring in the beehive entrance, and a Raspberry Pi 3 (Raspberry Pi Foundation, Cambridge, UK) embedded system for acquiring sensor data including temperature, humidity, light intensity, and rain level. A black acrylic observation box, which included a web camera and red LED lighting board, was attached to the beehive entrance. Technical details are discussed in our previous work [7].
Data collected by the image monitoring system included daily incoming and outgoing counts, daily differences in incoming and outgoing counts, daily population loss rate, temperature, and relative humidity. The daily population loss rate (LR) of a honey bee colony was calculated using Equation (1): where C in (t) and C out (t), as functions of time (t), are the daily incoming and outgoing counts, respectively. This equation was used in our previous work to assess the effects of pesticides on honey bee colonies [7]. Local weather data were wind speed, wind speed gust, wind direction, wind direction gust, precipitation, and the duration of precipitation and collected from the open-source data platform of the Taiwan Central Weather Bureau. The local weather data were cross verified with the rain level, temperature, and humidity sensor data collected by the monitoring system; as they showed similar values, they were used as additional input features of the forecasting model. All collected data were used for feature selection of the forecasting model and analysis. The nomenclature of the data is shown in Table 1.

Experimental Setup
Experimental data were collected from two experiments in a honey bee farm located in Hsinchu City, Taiwan. The details of each experiment are summarized in Table 2. The honey bees studied in this research were Apis mellifera. In each experiment, four healthy colonies with four honey bee combs and one queen each were prepared. The four replicated beehives were used to test the adaptability of the proposed forecasting algorithm. The hives were checked every two weeks to ensure the healthy condition of the colony. The data obtained from the two experiments were referred to as dataset E An and E Bn , respectively, where n is the beehive number.

Honey Bee Colony Population Loss Rate Forecasting and Early Warning Algorithm
The flowchart of the forecasting and early warning algorithm is presented in Figure 1. The algorithm was implemented using Python 3.5 with the support of Keras v2.2.6 deep learning library [23] and SciKit-Learn machine learning library [24].

Experimental Setup
Experimental data were collected from two experiments in a honey bee farm located in Hsinchu City, Taiwan. The details of each experiment are summarized in Table 2. The honey bees studied in this research were Apis mellifera. In each experiment, four healthy colonies with four honey bee combs and one queen each were prepared. The four replicated beehives were used to test the adaptability of the proposed forecasting algorithm. The hives were checked every two weeks to ensure the healthy condition of the colony. The data obtained from the two experiments were referred to as dataset EAn and EBn, respectively, where n is the beehive number.

Honey Bee Colony Population Loss Rate Forecasting and Early Warning Algorithm
The flowchart of the forecasting and early warning algorithm is presented in Figure  1. The algorithm was implemented using Python 3.5 with the support of Keras v2.2.6 deep learning library [23] and SciKit-Learn machine learning library [24].

Data Pre-Processing
LR data collected by the image monitoring system were pre-processed by interpolation and data normalization before feeding into the TCN forecasting model. Sample results are shown in Figure 1. Interpolation is a mathematical method that fits a function to a dataset and uses the function to fill in missing data based on the nearest past and future

Data Pre-Processing
LR data collected by the image monitoring system were pre-processed by interpolation and data normalization before feeding into the TCN forecasting model. Sample results are shown in Figure 1. Interpolation is a mathematical method that fits a function to a dataset and uses the function to fill in missing data based on the nearest past and future values. Each set of missing data was filled in with the mean computed from its four nearest values.
Afterwards, the completed data were normalized by scaling the input feature into values from 0 to 1 based on the minimum and maximum values of each input feature, respectively. Dataset E A had five days of missing data while dataset E B had eight days of missing data, both due to power outage and maintenance of the image monitoring system.

Temporal Convolutional Network for Daily Population Loss Rate Forecasting
A temporal convolutional network (TCN) [20] was used to forecast the future honey bee daily population loss rate and predict potential colony collapse. The key characteristic of TCNs is its usage of convolutions. Convolutions are causal and does not depend on any future timestep data. Unlike other deep learning models for forecasting like GRU and LSTM. TCN has longer memory and can process time series data of any length and generate similar long predictions [25]. For these reasons, TCN was selected as the forecasting model in this work. The TCN forecasting model comprised a series of blocks, which individually contained a sequence of convolutional layers. Each layer was composed of dilated convolutions, associated with a dilation factor d, with rectified linear unit (ReLU) as non-linear activation function. Dilation introduces a fixed step between every adjacent filter. Larger dilations and larger filters of size k effectively expand the receptive field [19,20]. In these convolutions, exponential increments in the value of d increases the depth of the network. This guarantees the presence of a filter that hits each input within the effective history [19]. A residual connection was added for each dilated convolution to integrate the convolutional result with the input layer. In this work, the input of the TCN model was defined as x t and the output was represented by y t , where x t contains n-dimensional parameters. The output of the TCN model was the LR of the next day.

Warning Level Determination
LR data were fitted to an isolation forest model [26] to categorize the honey bee daily loss rate according to different warning levels. Isolation forest is an unsupervised learning classification method based on decision trees. It recursively creates partitions by randomly selecting a feature and then picks a random split value between the minimum and maximum value of the selected feature. It produces smaller paths for the outlier values. Unlike other outlier detection methods, isolation forest explicitly identifies anomalies instead of profiling normal data points. Thus, the outlier values can be easily distinguished from non-outlier data. Isolation forest performs better than most anomaly detection algorithms across different datasets, based on receiver operating characteristic (ROC) performance and precision [26] and therefore chosen for the given problematic. Anomalies from LR data were detected by isolation forest algorithm and categorized as abnormal, while the rest of the data were classified as normal. Normal level indicates a honey bee population within normal conditions at a natural loss rate without the need of manual interception. An abnormal level suggests detailed monitoring to prevent potential colony damage.

Data Characterization
Forecasting accuracy depends strongly on the data characteristic. Determining the data characteristic can help in developing an optimized forecasting strategy and selecting a suitable model [27]. Two indices were computed to identify the demand pattern of the data: average inter-demand interval (ADI) and square of the coefficient of variation (CV 2 ). ADI is a measure of the demand regularity in time based on the average interval between demands while CV 2 is an index that measures the variation in quantities. The two indices were used to classify LR data according to four different categories: smooth, intermittent, erratic, and lumpy [28].

Feature Importance Analysis and Selection
The objective of feature selection is to minimize the input features to optimize forecasting performance. LR data were fit to a random forest classifier to compute the importance score of each feature using functions provided in Scikit-Learn machine learning library [24]. Afterwards, correlation analysis was carried out to generate a correlation heatmap that describes relationships between features. The correlation coefficient computed between features ranges from -1 to 1, where values higher than 0.5 indicate significant linear correlation. Feature groups (FG) were formed based on the feature analysis results for comparison.

Model Training and Optimization
TCN models were trained based on the different feature groups obtained from the feature selection step for comparison. Grid search was used to find the best values of the hyperparameters, including learning rate, number of epochs, number of TCN layers, batch size, window size, number of filters, kernel size, dilation, and dropout rate. Each model was optimized using the Adam optimizer, while minimizing the mean squared error (MSE), with He normal as the kernel initializer [29].

Algorithm Evaluation and Statistical Analysis
To evaluate the performance of the forecasting model, the weighted mean average percentage error (WMAPE) was computed. WMAPE measures the absolute percentage error in prediction and was calculated by Equation (2): where y obs i denotes the real value, y pred i denotes the predicted value, and N is the total amount of data.
There is no standard in determining warnings of the daily population loss rate of honey bee colonies. Thus, the warning level determination method was validated based on the daily population loss rate values found in related literature. In the work of Rumkee et al. [30], it was found that a daily mortality of forager honey bees of more than 15% led to a rapid loss in population and a colony survival rate of only 50%. Meanwhile, Dukas [31] discovered a natural a daily mortality rate of forager honey bees of about 13.4% due to aging and natural deaths. Therefore, here, colonies with LR values higher than 13% were considered to be abnormal, for verification purposes. The accuracy of the warning level determination method was computed using Equation (3): where A predicted and A true are the number of predicted abnormal LR by isolation forest algorithm and number of true abnormal LR greater than 13%, respectively. Student paired t-tests (p < 0.05) were used to confirm the significant differences between the results and data acquired.

Data Characterization Results
LR data obtained from the beehives of each dataset were characterized by computing ADI and CV 2 , as summarized in Table 3. The results showed a smooth demand pattern for most LR data, except for dataset E B3 . On average, the ADI and CV 2 of all the datasets indicated a smooth characteristic and the data were regular in time and quantity. High regularity indicates high forecastability [27]. After each experiment, the honey bee colonies were also manually inspected to determine if their condition was normal or collapsed. A colony with thousands of in-and-out activities at the hive entrance was classified to be normal, otherwise, as collapsed. In experiment E A , all four beehives were normal while there were two normal beehives (E B2 and E B4 ) and two collapsed beehives (E B1 and E B3 ) in experiment E B . The higher values of CV 2 obtained from datasets E B1 and E B3 compared to the other data might be indicative. The two collapsed beehives in experiment E B were further analyzed and are discussed in later sections.

Feature Correlation and Importance Analysis
The feature importance analysis results obtained from the data on each beehive are presented in Figure 2. The feature importance scores computed from the data of each beehive varied. In general, it was found that precipitation, hourly precipitation, incoming counts, outgoing counts, and daily count difference had the highest scores. However, it was also found that the collapsed beehives (Figure 2e,g) had different results from the normal beehives (Figure 2a-d,f,h). The feature importance score of normal hives were high in terms of precipitation, precipitation duration, incoming counts, outgoing counts, difference in count (Figure 2a-d,f,h). Meanwhile, the incoming count, difference in count had high scores for collapsed hives (Figure 2e,g), suggesting that the incoming counts and daily count difference were also important for forecasting.
The correlation heatmap of the features from dataset E A2 is shown in Figure 3a and is representative for a normal case. The correlation heatmap is highly correlated with the computed feature importance scores (greater than 0.7) for groups of features, such as temperature, humidity, and precipitation, most likely due to their similarities. Highly correlated features are often considered redundant because they do not add useful information to the model. Hence, one feature was selected from each set of related features, forming a feature group called FG2, as shown in Table 4. The computed correlation values of each feature in relation to LR were sorted (Figure 3b). It was found that some features had low correlation values of −0.25 to 0.25 such as WS gust , WS, WD, T min , T max , and WD gust signifying that these features had neither a positive nor negative correlation to LR. Therefore, the features were removed to form another feature group called FG3.

Feature Group Selection
For comparison, a feature group named FG1, which included all the available input features, was prepared. A univariate model, using LR data as the only feature, was also trained. Individual forecasting models were trained with a training dataset of 80 days, with the last 20% of the training data as validation set. Other training parameters were not yet considered and were optimized later. The results are shown in Figure 4.
Based on the testing results, the WMAPE of the univariate model was about 42.5 ± 3.4%. This extensively high error demonstrates that using only the daily loss rate as input feature is inadequate for forecasting. The multivariate model using FG1 performs better with a WMAPE of 20.1 ± 2.8% as supported by T-test results with a significant difference between both approaches. There is no significant difference between FG1, FG2, and FG3, this indicates that reducing the number of features does not affect the performance as expected from apparent redundancies. Interestingly, adding the forecasted weather data as an additional set of features improves the performance and reduces the error to 17.1 ± 1.6%. This reflects that forecasted weather-data as supplemental information contribute to an enhanced forecasting model. Based on these findings, FG4 was selected as the input feature group of the forecasting model throughout this work.  The weather forecast data from the weather station were also utilized to examine their potential effect on the forecasting results. This formed another feature group called FG4, which included the following day's forecasted values of precipitation, relative humidity, and temperature. Based on the related literature, temperature also has an effect to the behavior of honey bees [32]. Yet, it can be seen from the feature importance score and correlation analysis results that its effect was not considerable for the current datasets. Even so, temperature was still added in all feature groups to ensure that the forecasting model would also be trained to adapt to drastic changes in the environmental condition.

Feature Group Selection
For comparison, a feature group named FG1, which included all the available input features, was prepared. A univariate model, using LR data as the only feature, was also trained. Individual forecasting models were trained with a training dataset of 80 days, with the last 20% of the training data as validation set. Other training parameters were not yet considered and were optimized later. The results are shown in Figure 4.
Based on the testing results, the WMAPE of the univariate model was about 42.5 ± 3.4%. This extensively high error demonstrates that using only the daily loss rate as input feature is inadequate for forecasting. The multivariate model using FG1 performs better with a WMAPE of 20.1 ± 2.8% as supported by T-test results with a significant difference between both approaches. There is no significant difference between FG1, FG2, and FG3, this indicates that reducing the number of features does not affect the performance as expected from apparent redundancies. Interestingly, adding the forecasted weather data as an additional set of features improves the performance and reduces the error to 17.1 ± 1.6%. This reflects that forecasted weather-data as supplemental information contribute to an enhanced forecasting model. Based on these findings, FG4 was selected as the input feature group of the forecasting model throughout this work.

Model Training and Hyperparameter Optimization
The dataset was split into training, validation, and testing sets according to time. The training set was used to train the model based on the validation set. The testing set was used for assessing the model performance in forecasting data outside the training and validation set. For normal colonies (EA1~EA4; EB2, EB4), the training set (Ytrain) in each time series was the first N-day observations, where N is the number of training days. The validation set (Yval) included the last 20% of the training set. The rest of the dataset was used for testing (Ytest). For the collapsed colonies (EB1, EB3), the daily LR during the first N days were normal and started to become abnormal at day 105 and 90 for EB1 and EB3, respectively. From prior testing, it was not possible to use normal LR data to forecast the abnormal LR data. Therefore, a model was trained using the data of EB1 for forecasting the data of EB3 and vice versa. Similar to the normal colonies, the last 20% of the Ytrain were used for validation.
The effect of different number of days N in training the forecasting model was evaluated to know how much data were needed to attain reliable forecasting performance. TCN models, using FG4 as input feature group, were trained with values of N, ranging from 20 to 100, as presented in Table 5. The results show that the model trained with fewer than 80 training days had poor performances with WMAPE as high as 50.2%. Using 80 training days and above, it led to better forecasting performance with WMAPE as low as 21.6%. Interestingly, it also shows that increasing the value of N from 80 to 100 did not considerably improve the forecasting performance. Based on these findings, the number of training days was set to 80 for later analysis.
The results of finding the best model hyperparameters via grid search, using N = 80, is shown in Table 6. It was found that the model can perform best by setting the input window size to 10; this means that the model needs at most a length of 10 days to achieve satisfactory forecasting performance. It was also observed that 2 TCN layers were sufficient to yield good forecasting performance; this can be mainly attributed to the smooth characteristic of the data. By applying the tuned hyperparameters, the WMAPE of the TCN model was reduced to 16.2%, from 21.6% of the model using the default values.

Model Training and Hyperparameter Optimization
The dataset was split into training, validation, and testing sets according to time. The training set was used to train the model based on the validation set. The testing set was used for assessing the model performance in forecasting data outside the training and validation set. For normal colonies (E A1~EA4 ; E B2 , E B4 ), the training set (Y train ) in each time series was the first N-day observations, where N is the number of training days. The validation set (Y val ) included the last 20% of the training set. The rest of the dataset was used for testing (Y test ). For the collapsed colonies (E B1 , E B3 ), the daily LR during the first N days were normal and started to become abnormal at day 105 and 90 for E B1 and E B3 , respectively. From prior testing, it was not possible to use normal LR data to forecast the abnormal LR data. Therefore, a model was trained using the data of E B1 for forecasting the data of E B3 and vice versa. Similar to the normal colonies, the last 20% of the Y train were used for validation.
The effect of different number of days N in training the forecasting model was evaluated to know how much data were needed to attain reliable forecasting performance. TCN models, using FG4 as input feature group, were trained with values of N, ranging from 20 to 100, as presented in Table 5. The results show that the model trained with fewer than 80 training days had poor performances with WMAPE as high as 50.2%. Using 80 training days and above, it led to better forecasting performance with WMAPE as low as 21.6%. Interestingly, it also shows that increasing the value of N from 80 to 100 did not considerably improve the forecasting performance. Based on these findings, the number of training days was set to 80 for later analysis. The results of finding the best model hyperparameters via grid search, using N = 80, is shown in Table 6. It was found that the model can perform best by setting the input window size to 10; this means that the model needs at most a length of 10 days to achieve satisfactory forecasting performance. It was also observed that 2 TCN layers were sufficient to yield good forecasting performance; this can be mainly attributed to the smooth characteristic of the data. By applying the tuned hyperparameters, the WMAPE of the TCN model was reduced to 16.2%, from 21.6% of the model using the default values.

Daily Population Loss Rate Forecasting Results
By exploiting the optimized model, the long-term forecasting results of the honey bee daily population loss rate for normal, and collapsed beehives are presented in Figures 5 and 6, respectively, while the WMAPE boxplots of the forecasting model is shown in Figure 7.

Daily Population Loss Rate Forecasting Results
By exploiting the optimized model, the long-term forecasting results of the honey bee daily population loss rate for normal, and collapsed beehives are presented in Figures 5  and 6, respectively, while the WMAPE boxplots of the forecasting model is shown in Figure 7.   The trained TCN model was able to successfully predict the testing LR data of normal beehives ( Figure 5). Most importantly, it was able to predict sudden changes in LR, as particularly seen in Figure 5a-d. The TCN model performed well on the collapsed beehives by forecasting the sudden changes in LR on both datasets, as particularly seen in Figure 6 from day 260 and so on. In general, the TCN model performed well in all datasets with WMAPE from 17% to 19% (Figure 7), with paired t-tests showing no significant difference between the original and predicted values per normal LR dataset. The results also indicate that the selected features in FG4 were sufficient for forecasting, even though there were differences in the condition of the beehives. Paired t-test showed a significant difference between normal and collapsed beehives (EB1, EB3) which most likely was due to the difference in data characteristics as mentioned in Table 3.  The trained TCN model was able to successfully predict the testing LR data of normal beehives ( Figure 5). Most importantly, it was able to predict sudden changes in LR, as particularly seen in Figure 5a-d. The TCN model performed well on the collapsed beehives by forecasting the sudden changes in LR on both datasets, as particularly seen in Figure 6 from day 260 and so on. In general, the TCN model performed well in all datasets with WMAPE from 17% to 19% (Figure 7), with paired t-tests showing no significant difference between the original and predicted values per normal LR dataset. The results also indicate that the selected features in FG4 were sufficient for forecasting, even though there were differences in the condition of the beehives. Paired t-test showed a significant difference between normal and collapsed beehives (EB1, EB3) which most likely was due to the difference in data characteristics as mentioned in Table 3. The trained TCN model was able to successfully predict the testing LR data of normal beehives ( Figure 5). Most importantly, it was able to predict sudden changes in LR, as particularly seen in Figure 5a-d. The TCN model performed well on the collapsed beehives by forecasting the sudden changes in LR on both datasets, as particularly seen in Figure 6 from day 260 and so on. In general, the TCN model performed well in all datasets with WMAPE from 17% to 19% (Figure 7), with paired t-tests showing no significant difference between the original and predicted values per normal LR dataset. The results also indicate that the selected features in FG4 were sufficient for forecasting, even though there were differences in the condition of the beehives. Paired t-test showed a significant difference between normal and collapsed beehives (E B1 , E B3 ) which most likely was due to the difference in data characteristics as mentioned in Table 3.

Daily Loss Rate Warning Threshold Results
The results of applying isolation forest algorithm on normal beehives to determine abnormal LR values are presented in Figures 8 and 9. Figure 8 applies for dataset E A1~EA4 and Figure 9 applies for datasets E B1~EB4 . It can be easily distinguished that the collapsed beehives (Figure 9e,g) had the most anomalies detected. Particularly, most of the results showed that LR values less than 13% were classified as normal otherwise as abnormal. Meanwhile, the results of the normal colonies showed only a few days with abnormal LR values. Upon validation, the accuracy of the fitted isolation forest model was about 90.0 ± 8.5%, as shown in Table 7.

Daily Loss Rate Warning Threshold Results
The results of applying isolation forest algorithm on normal beehives to determine abnormal LR values are presented in Figures 8 and 9. Figure 8 applies for dataset EA1~EA4 and Figure 9 applies for datasets EB1~EB4. It can be easily distinguished that the collapsed beehives (Figure 9e,g) had the most anomalies detected. Particularly, most of the results showed that LR values less than 13% were classified as normal otherwise as abnormal. Meanwhile, the results of the normal colonies showed only a few days with abnormal LR values. Upon validation, the accuracy of the fitted isolation forest model was about 90.0 ± 8.5%, as shown in Table 7.

Discussions
The abnormal daily LR of the honey bee colonies in this study was about 13%, in agreement with related literature [30,31]. In the work of Rumkee et al. [30], it was found that a colony will collapse if the daily mortality of forager honey bees exceeds 15% and was equally observed for both collapsed beehives (E B1 , E B3 ) in this study (Figure 9e,g).
Our experiments show that honey bee colony collapses can be predicted. Previous studies identified the collapse of honey bee colonies was most likely caused by in-hive colony behavior and health of the honey bees [31]. A typical in-hive anomaly is caused by the loss of a queen bee. Indicators of a queen-less colony include missing eggs and brood, increased number of drones, a significant drop in population, and stored pollen and honey in the brood cells. These indicators were observed from the two collapsed beehives in the present study. Based on our recorded data, the loss of the queen occurred at about day 105 of E B1 and day 95 of dataset E B3 . According to Lopes et al. [33], workers of queen-less colonies can live up to 80 days. This was similarly observed in our study which shows that the colony collapsed about 60-80 days after the queen disappeared (Figure 6a,b).
There are no strict conventions for defining the LR threshold of a collapsed honey bee colony. Based on the experimental results obtained, the daily loss rate of the beehives became abnormal or several consecutive days, which was about 50 and 60 days for dataset E B1 and E B3 , respectively. The presence of an abnormal status as automatically classified is a strong indication for the required interference of beekeepers at early stage to prevent further damage to the honey bee colony.

Conclusions
A reliable method for forecasting and early detection of abnormal honey bee colony population loss rate is presented. The proposed method was optimized by appropriately selecting informative features and tuning the hyperparameters of a TCN forecasting model. The forecasting model performed best by using selected features, such as daily population loss rate, incoming counts, difference in counts, temperature, humidity, precipitation, and the forecasted following day's temperature, humidity, and precipitation obtained from the local weather station. Upon comparison with other selected feature groups, it was discovered that using the forecasted weather data as supplemental input feature improved the model performance. It was also found that the forecasting model performed well by training it with at least 80 days of historical population loss rate data. The forecasting model was able to accurately forecast the following day's honey bee colony population loss rate with a WMAPE of 17% to 19%. The forecasting model optimization strategy can be used as a reference for researchers to effectively improve their model.
The population loss rate data output of the forecasting model was further utilized to define the warning levels using an isolation forest algorithm, thus enabling the proposed method to determine whether the beehive colony status to be normal or abnormal. By validating the obtained values with the results of related works, it was found that the warning level determination method was able to yield an accuracy of 90.0 ± 8.6%.
Although the exact reason why a honey bee colony collapsed remained unclear, the data of the image monitoring system were successfully utilized to train a forecasting model that can potentially explain the phenomenon. If the daily loss rates were abnormal for a certain number of consecutive days, it indicated that honey bee collapse may potentially occur. The proposed method, together with the image monitoring system, was proven to assist in the early detection of collapsed honey bee colonies. This work can be used to deliver fast and reliable data-driven information to honey bee colony managers so that they can ensure the health of their honey bee colonies. To actualize the system, an early warning can be determined by detecting abnormal daily loss rates while a notification can be sent to the beekeeper via mobile phone. Upon receiving the message, beekeepers can check the beehive and prevent potential harm to the honey bee colony. Beekeepers can install the monitoring system in selected beehives and monitor the beehive health index remotely.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, T.-T.L., upon reasonable request.