- Introduction
Aquaculture, understood as the cultivation of aquatic organisms in controlled environments, has become one of the pillars of global food security throughout the 21st century. The continuous growth in demand for high–biological-value protein, combined with the stagnation of capture fisheries at sustainable production levels, has driven aquaculture to assume a central role in the supply of aquatic foods. According to the Food and Agriculture Organization of the United Nations (FAO, 2024), aquaculture production has already surpassed capture fisheries as the main source of fish for human consumption, standing out for its high feed conversion efficiency, potential for large-scale production, and adaptability to diverse socio-environmental contexts.
In Brazil, this expansion scenario is particularly evident. The country benefits from favorable edaphoclimatic conditions, abundant water resources, an extensive coastline, and large hydrographic basins, all of which support the growth of national aquaculture. According to data from the Brazilian Fish Farming Association (Peixe BR, 2025), tilapia farming leads Brazilian aquaculture production, consolidating itself as the sector’s main segment. At the same time, shrimp farming and the cultivation of native species, such as tambaqui (Colossoma macropomum), play a strategic role, especially in specific regions of the country. Despite the progress achieved, the FAO (2024) emphasizes that Brazil, although ranked among the largest aquaculture producers in Latin America, still has significant underexploited productive potential.
In this context, the use of data analysis tools has become increasingly relevant to the sustainable development of aquaculture. The growing availability of information on production, markets, and environmental variables enables the application of statistical and machine-learning methods capable of identifying patterns, modeling complex relationships, and supporting decision-making processes (Presenza et al., 2025). Time-series models, such as ARIMA, are widely employed to capture trends, seasonality, and autocorrelation structures in longitudinal data (Box et al., 2015; Panagiotelis et al., 2021).
Additionally, linear regression methods and tree-based algorithms, such as Random Forest, stand out for their robustness in modeling nonlinear relationships and handling high-dimensional datasets (Biau & Scornet, 2016). Multivariate approaches, such as Principal Component Analysis (PCA), allow for dimensionality reduction and the identification of dominant patterns of variability (Jolliffe & Cadima, 2016), while clustering techniques, such as the Partitioning Around Medoids (PAM) method, enable the segmentation of regions or production systems based on similar characteristics (Kaufman & Rousseeuw, 2009).
Thus, the integrated application of these statistical and computational tools goes beyond mere scenario forecasting, contributing to strategic planning, sustainable aquaculture management, and the formulation of evidence-based public policies.
In light of this framework, the present study aims to analyze the evolution of Brazilian aquaculture production between 2013 and 2023, with emphasis on the main farmed species, tilapia, tambaqui, and shrimp, identify regional patterns of productive specialization, and project prospective scenarios through 2030.
We hypothesized that Brazilian aquaculture exhibits distinct regional production patterns and that time-series models, particularly ARIMA, provide more accurate forecasts than traditional regression and machine-learning approaches.
2. Methodology
The analytical procedure was conducted using the R programming language, employing a set of statistical and data-manipulation packages to ensure the robustness and integrity of the analysis. The study adopted a quantitative approach, with a methodological focus on exploratory analysis of secondary data. The database was obtained from the Brazilian Institute of Geography and Statistics (IBGE) and contains annual aquaculture production data (kg) covering the period from 2013 to 2023 (https://sidra.ibge.gov.br/tabela/3940).
After exporting the original data in CSV format, a tabulation and filtering stage was performed, in which the 20 most representative and relevant species in national production were selected. These data were organized into a new file, named Aquicultura.csv, which served as the structured basis for subsequent analyses.
During the preprocessing stage, the absence of production records for a given species in a specific region was not treated as random missing data; but rather interpreted as an indication of non-production. This pattern falls under the category of Missing Not at Random (MNAR), according to the typology proposed by Little and Rubin (2002), thereby justifying the assumption that “not reported” values correspond to zero production. Accordingly, to ensure that the analysis reflected only actual production, null records were excluded, retaining only observations with Producao_kg > 0.
Data standardization was carried out by harmonizing variable nomenclature to Region, Year, Species, and Production_kg, as well as by converting numerical values to the appropriate format.
Following data tabulation and preprocessing, initial descriptive analyses were conducted, followed by the application of Principal Component Analysis (PCA) to reduce data dimensionality and facilitate the interpretation of underlying patterns. In the segmentation stage, clustering was performed using the Partitioning Around Medoids (PAM) method, with the elbow and silhouette criteria applied to determine the optimal number of clusters. Finally, predictive models, including ARIMA, decision trees, and linear regression, were implemented to model and forecast the behavior of the analyzed variables (Kaufman & Rousseeuw, 2009; Jolliffe & Cadima, 2016; Box et al., 2015; Biau & Scornet, 2016).
The analyses were supported by statistical packages such as tidyverse, readr, reshape2, and lubridate for data acquisition and preprocessing; forecast, e1071, caret, and randomForest for predictive modeling; and cluster, factoextra, and corrplot for clustering and principal component analyses. All code developed for this study is publicly available at: https://github.com/mtank6691-coder/C-digo-TCC-Murilo-/blob/main/DATA%20SCIENCE%20MURILO%20TANK.
2.1 Exploratory Data Analysis (EDA)
Exploratory Data Analysis (EDA) was conducted with the purpose of statistically describing the dataset and guiding the formulation of hypotheses to be tested in subsequent stages of predictive modeling. For this purpose, the ggplot2 library was employed.
Among the visualizations produced, bar charts representing total production by region and by species stand out. Categories were ordered in descending order based on production volume using the reorder() function. The choice of this type of chart is grounded in its effectiveness in representing categorical variables with numerical aggregations, as discussed by Cleveland (1985) in his work on best practices in data visualization.
2.2 Principal Component Analysis (PCA)
Principal Component Analysis (PCA) is a statisti-cal technique widely used for dimensionality reduction, whose objective is to transform a set of correlated variables into a new set of orthogonal (uncorrelated) variables, known as principal com-ponents (Jolliffe & Cadima, 2016). Its application in this study is justified by the need to extract the main variance structures from the data and to identify associations between regions and aqua-culture species based on their production profiles.
Mathematically, the first principal component (PC1) is defined as a linear combination of the original variablesnX_1,X_2,…,X_p , expressed as:
PC1=a11X1+a12X2 +⋯+a1p Xp x'
where the coefficients a_1jcorrespond to the elements of the eigenvector associated with the largest eigenvalue of the covariance (or correlation) matrix and represent the contribution of the variable X_jto the formation of the principal component. These eigenvectors indicate the directions of maximum variance in the multivariate data space (Jolliffe & Cadima, 2016). The second principal component (PC2) is also a linear combination of the original variables, given by:
PC2=a21X1+a22X2+⋯+a2p Xp
being defined based on the eigenvector asso-ciated with the second largest eigenvalue of the covariance or correlation matrix. PC2 is orthogonal to PC1 and explains the largest possible proportion of the remaining variance not captured by the first component. Thus, while PC1 represents the primary gradient of variation in the data, PC2 reveals secondary patterns of differentiation among observations, contributing to a more in-depth understanding of the associations between regions and aquaculture species.
To prevent variables with larger numerical magnitudes from disproportionately influencing the PCA results, the data were previously standardized using the function prcomp (..., scale. = TRUE), ensuring that all variables had a mean of zero and a standard deviation of one (Rencher & Christensen, 2012). Standardization is particu-larly important when variables are measured on different scales, which is common in studies involving production data.
The interpretation of the results was enhanced using the factoextra package, which provides effective visualization tools for principal compo-nent analysis. The plots generated using the fviz_pca_ind() function enabled the visualization of similarity patterns among regions (Kassambara, 2017).
2.3 Clustering analysis
To group regions based on similar production profiles, a clustering approach centered on the Partitioning Around Medoids (PAM) method was employed. This technique was chosen due to its greater robustness to the presence of outliers when compared to the traditional K-means algorithm. While K-means relies on centroids, hypothetical mean values that may not correspond to actual observations, PAM operates with medoids, which are real data points from the dataset and therefore more accurately represent the central elements of each group (Kaufman & Rousseeuw, 2009).
The PAM algorithm functions by minimizing the sum of dissimilarities between each observation and the medoid of the cluster to which it belongs. To compute these dissimilarities, the Euclidean distance was used, defined by the following mathematical expression:
where x and y represent vectors of observations and n denotes the number of variables conside-red. This metric allows quantifying the degree of similarity (or dissimilarity) between two regions based on their production profiles.
The determination of the optimal number of clusters (k) was performed and validated using internal statistical validation methods, notably the Elbow Method and Silhouette Analysis. The Elbow Method consists of evaluating the variation in the sum of intra-cluster dissimilarities as a function of the number of clusters, with the optimal number of groups identified at the point where a less pronounced marginal reduction in this measure occurs, characterizing an “elbow” in the curve.
Complementarily, Silhouette Analysis was emplo-yed to assess clustering quality based on intra-cluster cohesion and inter-cluster separation. The silhouette coefficient for an observation i is given by:
where a(i)represents the average distance from observation i to all other observations belonging to the same cluster, and b(i)corresponds to the smallest average distance from i to observations in any other cluster. The value of s(i)ranges from −1 to 1, with values close to 1 indicating an appropriate assignment of the observation to its cluster and good separation from the remaining clusters, thus evidencing a satisfactory clustering configuration.
Complementary to the PAM partitioning technique, a Hierarchical Clustering Analysis was conducted in order to explore the hierarchical structure of the groupings in a visual and intuitive manner. The Ward.D2 linkage method was applied, which seeks to minimize total within-cluster variance at each step of cluster fusion. Euclidean distance was likewise adopted as the dissimilarity metric.
The results of this analysis were represented through a dendrogram, allowing the visualization of similarity relationships among fish-producing regions at different levels of aggregation and serving as a complementary tool for the inter-pretation and validation of the clusters obtained.
2.4 Predictive modeling methodology
The final phase of this study focused on the application and comparison of three predictive modeling approaches aimed at forecasting the future production of the main aquaculture species. The use of multiple models, combined with the application of formal performance evaluation crite-ria, ensured greater robustness of the generated estimates and enabled a comprehensive compa-rative analysis of the adopted methodologies.
2.5 Performance metric and validation strategy
To assess model accuracy, the Mean Absolute Error (MAE) was adopted as the performance metric. This measure was selected due to its interpretability, as its values are expressed in the same unit as the original data (kg), thereby facilitating direct comparison between observed and predicted values (Hyndman & Koehler, 2006).
The MAE is formally defined as:
where Y_trepresents the actual production value at observation t, Y ̂_tcorresponds to the value predicted by the model for the same observation, and ndenotes the total number of observations considered.
2.6 Training and testing strategy
The evaluation and forecasting methodology was based on predictive time-series analysis. The historical dataset covering the period from 2013 to 2022 was used as the training set, while data for the year 2023 were reserved as the test set. This strategy ensures an unbiased assessment of the models, as the validation data were not used during the training process (Hyndman & Athanasopoulos, 2018).
Because the time series consisted of only eleven annual observations, more complex resampling strategies, such as repeated cross-validation or rolling-origin evaluation, were not considered appropriate. Under these circumstances, the hold-out approach represents a robust and widely accepted validation strategy for short annual time series. The same train-test split approach, as well as performance comparison, was consistently applied to both the total production series and the species-specific production time series.
2.7 Evaluated predictive models
Three distinct predictive models were evaluated for each analyzed time series (total production and species-specific production): Linear Regression (LM), Random Forest (RF), and ARIMA.
2.8 Linear Regression (LM)
The first model evaluated was Linear Regression (LM), a classical statistical model that assumes a linear relationship between the response variable (production) and the explanatory variable time (year). Its main advantage lies in its ease of implementation and interpretability. However, its fundamental limitation is the inability to capture more complex patterns, such as nonlinear behaviors or structural changes, which are frequently observed in economic and production data (Box et al., 2015).
2.9 Random Forest (RF)
The second model considered was Random Forest (RF), a machine-learning technique based on an ensemble of decision trees, where the aggregation of individual predictions contributes to increased robustness and predictive power (Biau & Scornet, 2016). This model was included to capture nonlinear relationships present in the data. To mitigate the risk of overfitting, five-fold cross-validation was employed using the trainControl function. This strategy ensured that the model was evaluated across different subsets of the training data, promoting greater generalization of predictions.
2.10 ARIMA (AutoRegressive Integrated Moving Average)
The third model evaluated was ARIMA, widely used in time-series analysis due to its ability to capture trends, seasonal patterns, and autoco-rrelations, that is, the dependence between present and past values of the series (Shumway & Stoffer, 2017).
An ARIMA model is denoted as ARIMA(p, d, q), combining three fundamental components: autoregressive (AR), integration (I), and moving average (MA). Its general formulation is expressed as:
(1-ϕ_1 B-…-ϕ_p B^p)(1-B)^d Y_t=c+(1+θ_1 B+⋯+θ_q B^q)ε_t (5)
The term Y_trepresents the value of the time series at time t; ϕ_iare the autoregressive parameters; θ_iare the moving average parameters; dcorresponds to the number of differences applied to the series to achieve stationarity; Bis the backshift (lag) operator; and ε_tis the random error term, assumed to be white noise.
The auto.arima() function from the R forecast package was used to automatically select the optimal values of pand q, based on information criteria such as the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). This approach allowed the fitting of models tailored to the dynamic characteristics of each individual production time series, following the methodology described Hyndman & Athanasopoulos (2018).
2.11 Forecasting Strategy and Scope of Analyses
Among the models evaluated for the total production series, ARIMA exhibited the best performance according to the MAE and was therefore selected as the reference model for supporting the final conclusions of the study.
However, for the production series of the key species, forecasts were generated using all three models (LM, RF, and ARIMA), with the objective of enabling a comparative analysis between linear, nonlinear, and time-series–specific approaches. Projections were made for a seven-year horizon (2024–2030), covering the following categories:
Tilapia: the species with the highest production volume;
Tambaqui: one of the main native species;
Shrimp: representative of shrimp farming;
Oysters, scallops, and mussels: a group characterized by high added value.
3. Results and discussion
3.1 Exploratory analysis
The exploratory analysis of the data revealed a significant growth in Brazilian aquaculture over the last decade, with annual production increasing from approximately 484,000 tons in 2013 to over 887,000 tons in 2023. However, this production (Figure 1A) is not evenly distributed across the country, being concentrated primarily in key regions such as the South, followed by the Northeast and Southeast. This regional predomi-nance is explained by factors such as favorable climatic conditions, developed infrastructure, and historical investment in the sector. Additionally, production is dominated by a few species, Figure 1B, with Tilapia, Tambaqui, and Shrimp accoun-ting for the majority of total output.
Principal Component Analysis (PCA) was applied to identify latent patterns and group regions with similar production characteristics, Figure 2.A. Native species such as Tambaqui, Pirarucu (Arapaima gigas), and Matrinxã (Brycon amazonicus) are concentrated in the North and Central-West regions, reflecting the relationship between local biodiversity and regional production potential. These species are adapted to the ecological conditions of the Amazonian and Pantanal rivers, favoring their adoption as production focuses. The Northeast shows a predominance of Shrimp, consolidating shrimp farming as a strategic activity due to the limitation of inland water resources. Meanwhile, the South and Southeast cluster around Tilapia, highlighting historical investments in intensive systems, infrastructure, and specific incentive policies.
The analysis of species’ contributions to the principal components shows that Dim1 explains 48.4% of the data variability, while Dim2 explains 24.7%, totaling approximately 73%. Species such as Tilapia, Traíra, and Trairão are strongly associated with the positive Dim1 axis, indicating intensive and consolidated production systems (Figure 2B).
