Artigo

Ferdinando Yoshio Agapito Urasaki; Luiz Gustavo Antônio de Souza

DOI: 10.22167/2675-6528-202602215

Artigo derivado de Trabalho de Conclusão de Curso (TCC), com conteúdo baseado no trabalho original do aluno e adaptado ao formato editorial da Revista E&S com apoio da ferramenta ResumeAI, solução de inteligência artificial desenvolvida pelo Instituto Pecege para síntese e organização textual.

Resumo

A estimativa da produtividade é essencial para o planejamento agrícola, mas a estrutura de agrupamento dos dados agronômicos frequentemente viola o pressuposto de independência dos modelos clássicos de regressão. Avaliou-se a acurácia de modelos multiníveis com classificação cruzada na predição do rendimento agrícola da cana-de-açúcar na região sudeste de Mato Grosso do Sul, abrangendo as safras de 2019/2020 a 2023/2024. O estudo analisou 29.674 observações, organizadas sobre 10.336 talhões e 30 variedades. Integraram-se preditores biofísicos e agrometeorológicos, aplicando-se a estratégia “Step-Up” para a construção da estrutura multinível e a decomposição da variância. Os resultados evidenciaram que a abordagem particionou adequadamente a variância espacial, genética e temporal, quantificando com êxito a interação genótipo-ambiente. O ajuste final com inclinações aleatórias (M7) explicou 81,1% da variância total e alcançou uma Raiz do Erro Quadrático Médio (RMSE) de 10,24 t ha⁻¹. A diferença de 9,8% aferida entre o Pseudo-R² marginal e o condicional demonstrou o relevante ganho explicativo gerado pela classificação cruzada. O desempenho preditivo equiparou-se à exatidão reportada por algoritmos de Machine Learning, mantendo o rigor paramétrico e a interpretabilidade dos preditores agronômicos. Concluiu-se que a modelagem multinível consistiu em uma ferramenta analítica robusta e transparente para a estimativa de produtividade e o planejamento do setor sucroenergético.

Palavras-chave: Agrometeorologia; Cana-de-açúcar; Classificação cruzada; Modelagem multinível; Rendimento agrícola.

1. Introdução

A cana-de-açúcar (Saccharum spp.) destaca-se como uma das principais commodities agrícolas do Brasil, ao lado da soja e do milho. Juntas, essas culturas ocupam aproximadamente 6,9% do território nacional e mais de 70% da área agrícola cultivada (EPE, 2024b). Além da produção de açúcar, a cultura é fundamental para a geração de etanol e para o aproveitamento de coprodutos, como a vinhaça na fertirrigação e o bagaço na cogeração de energia elétrica (Avilez et al., 2018; EPE, 2024b). Esse aproveitamento reflete-se na matriz energética brasileira, onde a biomassa da cana-de-açúcar e seus derivados representam 16,8% da oferta interna, sendo o segundo maior recurso energético do país, superados apenas pelo petróleo (EPE, 2024a).

Nesse contexto, a região sudeste de Mato Grosso do Sul configura-se como um polo estratégico para a produção de cana-de-açúcar, impulsionada por sua elevada densidade bioenergética (BIOSUL, 2024) e por condições edafoclimáticas favoráveis (Manzatto et al., 2009). Para garantir a eficiência operacional das unidades produtoras, a acurácia na estimativa da produtividade, mensurada em toneladas de cana-de-açúcar por hectare (TCH), é crucial. Tal precisão é determinante para o planejamento da colheita e para a otimização logística, assegurando o suprimento contínuo à planta industrial (Silva A.S. et al., 2022).

O potencial produtivo da cana-de-açúcar pode atingir 400 t ha⁻¹ ao ano, devido à sua alta eficiência fotossintética. No entanto, o rendimento é condicionado por uma complexa interação de fatores ambientais e operacionais, com destaque para a disponibilidade hídrica e nutricional (Silva et al., 2012). Como cultura semiperene, com um ciclo médio de cinco cortes, a cana-de-açúcar apresenta quatro fases fenológicas distintas (brotação, perfilhamento, crescimento vegetativo e maturação), cada uma com demandas hídricas específicas. As fases iniciais são particularmente sensíveis ao estresse hídrico, e o crescimento vegetativo é o período mais crítico, enquanto uma restrição hídrica moderada na maturação favorece o acúmulo de sacarose (Avilez et al., 2018). A demanda hídrica total varia entre 1.500 e 2.500 mm, e a temperatura ótima situa-se entre 25 e 33 °C (Almeida et al., 2008). Eventos extremos, como geadas, também podem comprometer significativamente a produção, com severidade proporcional à intensidade e duração do fenômeno (Fiorentin et al., 2016).

Com o objetivo de quantificar o impacto desses fatores agrometeorológicos e viabilizar a predição de safras, a modelagem da produtividade da cana-de-açúcar tem sido amplamente explorada. Métodos como a regressão linear múltipla foram utilizados para estimar a produtividade anual com base em variáveis como irrigação, precipitação e temperatura do ar (Silva A.S. et al., 2022), ou para correlacionar a produção de biomassa com índices de vegetação (Silva M.G.S. et al., 2022). Scarpari (2002) também propôs a predição da produtividade agrícola integrando dados de precipitação e graus-dias.

Apesar dos avanços representados por esses modelos, a regressão linear múltipla frequentemente viola o pressuposto de independência das observações, o que pode levar a estimativas de erro incorretas e intervalos de confiança otimistas (Gelman e Hill, 2006). No contexto agrícola, a produtividade não é um evento isolado, mas o resultado de interdependências estruturais complexas. Por exemplo, variações climáticas em uma safra afetam múltiplas áreas simultaneamente, a variedade estabelece potenciais biológicos comuns entre talhões, e a sucessão de ciclos em um mesmo talhão induz autocorrelação temporal.

Para abordar adequadamente essas interdependências, a modelagem multinível surge como uma alternativa analítica mais robusta em comparação à regressão clássica. Em dados agrícolas, a organização das observações raramente segue uma estrutura hierárquica estrita, pois fatores de agrupamento como variedade e safra frequentemente se cruzam ao longo do tempo sobre os talhões, configurando uma estrutura de classificação cruzada. Essa abordagem permite isolar a variância intrínseca a cada contexto (espacial, genético e temporal), corrigindo a subestimação dos erros-padrão, problema comum em modelos que ignoram essa complexidade (Gelman e Hill, 2006). Assim, este estudo justifica-se pela necessidade de uma ferramenta analítica que não apenas preveja a produtividade com acurácia, mas também preserve o rigor paramétrico e a interpretabilidade dos preditores agronômicos.

Sob essa premissa, objetivou-se avaliar a acurácia de modelos multiníveis com classificação cruzada na predição da produtividade da cana-de-açúcar em unidades produtoras na região sudeste de Mato Grosso do Sul, abrangendo a série histórica das safras de 2019/2020 a 2023/2024. Especificamente, buscou-se decompor a variância total para quantificar a influência dos componentes de talhão, variedade e safra na explicação da variabilidade produtiva, bem como determinar a contribuição relativa dos preditores biofísicos e agrometeorológicos sobre a variância do TCH.

2. Material e Métodos

O estudo foi conduzido na região sudeste de Mato Grosso do Sul, abrangendo os municípios de Angélica, Deodápolis, Glória de Dourados, Ivinhema, Jateí, Naviraí, Nova Andradina e Novo Horizonte do Sul. A escolha dessa região justificou-se pela sua relevância estratégica e pela robustez das séries históricas de dados agronômicos, fundamentais para o ajuste de modelos de predição de produtividade em nível local.

A base de dados foi estruturada a partir da integração de registros operacionais, dados de sensoriamento remoto multiespectral e variáveis ambientais exógenas, categorizadas como:

• Agronômicos: Registros de unidades de produção sucroenergética, incluindo caracterização edáfica, datas de plantio e colheita, material genético (variedade), número de corte, produtividade agrícola (TCH) e geometria dos talhões.
• Variáveis Ambientais e Espectrais:
  • Sensoriamento Remoto: Imagens do satélite Sentinel-2 (MSI), processadas no Nível 2A, com resolução espacial de 10 m e 20 m e revisita de cinco dias, para extração de 11 índices espectrais.
  • Precipitação: Estimativas da coleção Climate Hazards Group InfraRed Precipitation with Station Data (CHIRPS), com resolução temporal diária e espacial de 0,05° (~5,5 km).
  • Topografia: Modelo Digital de Elevação (MDE) derivado da Shuttle Radar Topography Mission (SRTM), com resolução espacial de 1 arco-segundo (~30 m).
  • Meteorologia: Dados de reanálise da coleção ERA5-Land, com resolução espacial de 0,1° (~11 km) e temporal diária, para extração de temperatura do ar (mínima, média e máxima), temperatura do ponto de orvalho, componentes vetoriais do vento e radiação solar de ondas curtas.

O fluxo de trabalho foi delineado em onze etapas: (i) limpeza e engenharia de atributos agronômicos; (ii) processamento geométrico para estatística zonal; (iii) coleção mensal de índices espectrais; (iv) coleção mensal de variáveis agrometeorológicas; (v) estruturação de séries temporais espectrais e agrometeorológicas; (vi) tratamento de lacunas dos índices espectrais; (vii) sincronização dos ciclos e ajuste das datas de manejo; (viii) consolidação de métricas e integração temporal; (ix) seleção de variáveis e redução de dimensionalidade; (x) especificação do modelo multinível e (xi) validação e avaliação de desempenho dos modelos.

Todo o processamento foi implementado em linguagem Python (versão 3.12), utilizando as bibliotecas Pandas, GeoPandas e NumPy. A integração e o processamento em nuvem dos índices espectrais e variáveis agrometeorológicas foram operacionalizados via API do Google Earth Engine (biblioteca ee) e Geemap. Para a modelagem multinível, utilizou-se a biblioteca Pymer4, que provê interface com os algoritmos do pacote lme4 da linguagem R.

O tratamento dos dados agronômicos brutos (225.992 registros) iniciou-se com a remoção de observações com dados ausentes de produtividade, número de corte e ambiente de produção, restringindo a amostragem às modalidades de produção própria e de parceria. Assegurou-se a consistência temporal e a paridade de área, removendo talhões com plantios multivarietais. O recorte analítico abrangeu as safras 2019/2020 a 2023/2024, do primeiro ao sétimo corte, e as 30 variedades de maior representatividade. Valores atípicos (outliers) de produtividade foram identificados pelo método do Intervalo Interquartil (IQR), com multiplicador de 1,5, estratificado por número de corte e ambiente de produção. A consistência entre atributos alfanuméricos e vetoriais foi validada, descartando-se registros com divergência absoluta superior a 0,025 ha. A base de dados resultante totalizou 32.033 registros.

Para o refinamento da amostragem, selecionaram-se talhões com área superior a 1,0 ha. As geometrias originais foram redefinidas por uma retração espacial de 20 m para mitigar o efeito de borda, assegurando que apenas os pixels internos fossem processados. Os 30.953 registros resultantes serviram como base vetorial para estatísticas zonais na plataforma Google Earth Engine (GEE).

Para a coleção mensal de índices espectrais, foram pré-selecionadas imagens do Sentinel-2 com cobertura de nuvens inferior a 20%, aplicando-se um filtro booleano baseado na camada de classificação de cena (SCL) gerada pelo algoritmo Sen2Cor. Os dados foram organizados em composições mensais, utilizando a mediana de cada pixel para eliminar ruídos atmosféricos. Foram calculados 11 índices espectrais derivados de diferentes bandas, conforme detalhado na Tabela 1.

Tabela 1. Relação de índices espectrais e respectivas formulações.

Indice	Sigla	Equação	Referência
índice de vegetação por diferença normalizada	NDVI	(PNIR – PR) / (PNIR + PR)	Rouse et al. (1974)
índice de vegetação por diferença normalizada verde	GNDVI	(PNIR – PG) / (PNIR + PG)	Gitelson et al. (1996)
índice de vegetação ajustado ao solo	SAVI	(PNIR – PR / (PNIR + PR + L)) * (1 + L)	Huete (1988)
índice de vegetação ajustado ao solo modificado 2	MSAVI2	(2PNIR + 1 – sqrt((2PNIR + 1)^2 – 8(PNIR – PR))) / 2	Qi et al. (1994)
índice de vegetação realçado	EVI	2.5 * (PNIR – PR / (PNIR + 6PR – 7.5PB + 1))	Huete et al. (2002)
índice de água por diferença normalizada	NDWI	(PNIR – PSWIR1) / (PNIR + PSWIR1)	Gao (1996)
índice de borda do vermelho por diferença normalizada	NDRE	(PRE2 – PRE1) / (PRE2 + PRE1)	Barnes et al. (2000)
índice de clorofila da borda do vermelho	CIre	(PRE3 – 1) / PRE1	Gitelson et al. (2005)
índice de clorofila da borda do vermelho invertido	IRECI	(PRE3 – PR) / (PRE1/PRE2)	Frampton et al. (2013)
índice de área foliar do sentinel-2	SeLI	(PNIR – PRE1) / (PNIR + PRE1)	Pasqualotto et al. (2019)
índice de vegetação de ampla faixa dinâmica	WDRVI	(PNIR – PR) / (aPNIR + PR)	Gitelson (2004)

Fonte: Dados originais da pesquisa. Nota: α = 0,1; L = 0,5. Bandas Sentinel-2: B2 (B), B3 (G), B4 (R), B5 (RE1), B6 (RE2), B7 (RE3), B8 (NIR), B11 (SWIR1).

A coleção mensal de variáveis agrometeorológicas utilizou dados de reanálise e estimativas de satélite (CHIRPS e ERA5-Land) para garantir continuidade temporal e cobertura espacial. Indicadores foram estruturados em cinco eixos: oferta hídrica (precipitação – P), regime térmico (temperaturas mínima, média e máxima – Tmin, Tmed, Tmax), disponibilidade energética (radiação solar incidente de ondas curtas – Rs), demanda evaporativa (evapotranspiração de referência – ET0) e estresse atmosférico (déficit de pressão de vapor – VPD). O acúmulo térmico foi quantificado por graus-dia diários (GDD) com Tbase de 18 °C e T corte de 35 °C. A ET0 foi estimada pelo método de Penman-Monteith. Após a computação diária, as séries matriciais foram agregadas em composições mensais, utilizando o redutor de soma para variáveis de fluxo (P, ET0, GDD e Rs) e a média para variáveis de estado (Tmin, Tmed, Tmax e VPD). O conjunto de indicadores e métodos de agregação é apresentado na Tabela 2.

Tabela 2. Indicadores agrometeorológicos e métodos de agregação.

Indicador	Sigla	Unidade	Representação	Agregação	Referência
precipitação	P	mm	oferta hídrica pluvial	soma	CHIRPS
temperatura mínima	Tmin	°C	limite térmico inferior	média	ERA5-Land
temperatura média	Tmed	°C	estado térmico médio	média	ERA5-Land
temperatura máxima	Tmax	°C	limite térmico superior	média	ERA5-Land
radiação solar de ondas curtas	Rs	MJ · m⁻²	oferta energética solar	soma	ERA5-Land
evapotranspiração de referência	ET₀	mm	demanda hídrica potencial	soma	Allen et al. (1998)
graus-dias acumulados	GDA	°C	acúmulo térmico fenológico	soma	Villa Nova et al. (1972)
déficit de pressão de vapor	VPD	kPa	estresse higrométrico	média	ERA5-Land

Fonte: Dados originais da pesquisa.

Para a sistematização das séries temporais, os dados matriciais mensais das coleções de índices espectrais e de variáveis agrometeorológicas foram integrados às geometrias dos talhões via estatística zonal na plataforma GEE. Os valores médios no talhão foram computados com um redutor aritmético, com escala espacial de saída em 20 m, convertendo os dados matriciais em atributos tabulares.

As lacunas internas dos índices espectrais foram tratadas via interpolação linear. Para as extremidades das séries, onde a interpolação não era viável, aplicou-se a imputação por persistência (“forward fill” e “backward fill”), propagando o último valor válido observado. As séries foram estendidas com uma janela temporal de 100 dias anteriores e posteriores aos limites de cada ciclo para prevenir distorções matemáticas.

As séries temporais dos índices espectrais foram suavizadas com o filtro polinomial Savitzky-Golay (polinômio de segunda ordem, janela móvel de cinco observações). Os limites do ciclo foram ajustados com base na resposta espectral do NDVI suavizado, utilizando a detecção de mínimos locais dentro de uma janela de tolerância de 90 dias ao redor das datas operacionais registradas.

Com os limites dos ciclos ajustados, as séries temporais mensais foram agregadas em janelas cronológicas fixas: F1 (30 a 120 dias após o início do ciclo) e F2 (120 a 270 dias). Para os 11 índices espectrais, calculou-se a integral temporal (AUC) pela regra trapezoidal. As variáveis agrometeorológicas de fluxo (P, ET0, GDD e Rs) foram integradas pela soma, e as variáveis de estado (Tmin, Tmed, Tmax e VPD) pela média aritmética. Cada atributo foi identificado por uma nomenclatura específica. Essa abordagem resultou em um conjunto de 38 características biofísicas e agrometeorológicas em 29.674 observações.

Para a seleção de variáveis e redução de dimensionalidade, utilizou-se o método de Máxima Relevância e Mínima Redundância (mRMR), que prioriza variáveis com maior dependência estatística em relação à variável-alvo e menor redundância. Após ranqueamento, selecionaram-se as 10 variáveis de maior relevância. Em seguida, para mitigar a multicolinearidade, conduziu-se o diagnóstico via Fator de Inflação da Variância (VIF), excluindo iterativamente variáveis com VIF superior a cinco, resultando em um conjunto de seis variáveis. As variáveis quantitativas foram centralizadas na grande média.

A construção do modelo multinível adotou a estratégia “Step-Up”, que consiste na adição de complexidade paramétrica crescente. O procedimento estatístico iniciou-se com a estimação de modelos incondicionais para particionar a variância total da produtividade (TCH) associada aos contextos espacial (talhão), genético (variedade) e temporal (safra).

O modelo M0 estabeleceu a linha de base da variância, incluindo apenas o intercepto aleatório para o talhão:

Y_i = beta_0 + u_t[i] + e_i (1)

Em que: Y_i: produtividade (TCH) na observação i; beta_0: intercepto fixo global; u_t[i]: efeito aleatório do talhão; e_i: erro residual.

Os modelos M1 e M2 testaram, independentemente, a inclusão dos componentes de variância associados à variedade e à safra, respectivamente:

M1: Intercepto Aleatório para Talhão e Variedade

Y_i = beta_0 + u_t[i] + v_v[i] + e_i (2)

Em que: v_v[i]: efeito aleatório da variedade.

M2: Intercepto Aleatório para Talhão e Safra

Y_i = beta_0 + u_t[i] + w_s[i] + e_i (3)

Em que: w_s[i]: efeito aleatório da safra.

No modelo M3, os interceptos aleatórios de talhão, variedade e safra foram incorporados simultaneamente, consolidando a estrutura de classificação cruzada:

M3: de Classificação Cruzada

Y_i = beta_0 + u_t[i] + v_v[i] + w_s[i] + e_i (4)

A estimação dos parâmetros e componentes de variância ocorreu pelo método de Máxima Verossimilhança Restrita (REML). Para avaliação do ganho de ajuste estatístico, utilizaram-se o Critério de Informação Bayesiano (BIC) e o Teste de Razão de Verossimilhança (LRT), com reajustes por Máxima Verossimilhança (ML) para comparação de modelos com diferentes estruturas de efeitos fixos.

Consolidados os componentes de variância no M3, adicionaram-se as variáveis categóricas de efeito fixo. O modelo M4 incorporou o número de corte e o ambiente de produção:

M4: Condicional com Controle Agronômico

Y_i = beta_0 + beta_1 corte_i + beta_2 ambiente_t[i] + u_t[i] + v_v[i] + w_s[i] + e_i (5)

Em que: corte_i: número de corte; ambiente_t[i]: ambiente de produção do talhão t; beta_1, beta_2: coeficientes de efeitos fixos.

Na sequência, o modelo M5 incorporou a idade da cultura (em dias) como covariável, centralizada na grande média:

M5: Condicional com Controle de Idade

Y_i = beta_0 + beta_1 corte_i + beta_2 ambiente_t[i] + beta_3 idade_cent_i + u_t[i] + v_v[i] + w_s[i] + e_i (6)

Em que: idade_cent_i: idade da cana-de-açúcar centralizada; beta_3: coeficiente de inclinação para a idade.

A integração das variáveis biofísicas e agrometeorológicas previamente selecionadas ocorreu via procedimento “Forward Stepwise” no modelo M6:

M6: Condicional com Efeitos Fixos

Y_i = beta_0 + beta_1 corte_i + beta_2 ambiente_t[i] + beta_3 idade_cent_i + Sum(lambda_q Z_q,i) from q=1 to Q + u_t[i] + v_v[i] + w_s[i] + e_i (7)

Em que: Z_q,i: valor do q-ésimo preditor biofísico ou agrometeorológico; lambda_q: coeficiente fixo parcial.

Por fim, para avaliar variações na sensibilidade de um contexto em resposta a variáveis específicas, testaram-se efeitos aleatórios de inclinação no M6, priorizando variáveis de estresse higrométrico (VPD), acúmulo térmico (GDA) e estado nutricional (NDRE e Clre), nas janelas F1 e F2. A modelagem seguiu a estratégia “Forward Stepwise” para cada variável candidata, especificando um modelo basal e incorporando-a como efeito aleatório de inclinação para comparação do critério de informação. A validação exigiu redução do BIC, ausência de falhas na convergência e VPC da inclinação de no mínimo 1% da variância total. O modelo M7 incluiu inclinações aleatórias:

M7: Inclinações Aleatórias (Interação Genótipo × Ambiente)

Y_i = beta_0 + beta_1 corte_i + beta_2 ambiente_t[i] + beta_3 idade_cent_i + Sum(lambda_q Z_q,i) from q=1 to Q + u_t[i] + v_0,v[i] + Sum(v_p,v[i] X_p,i) from p=1 to P + w_s[i] + e_i (8)

Em que: v_0,v[i]: intercepto aleatório da variedade v; X_p,i: valor do p-ésimo preditor submetido ao efeito aleatório; v_p,v[i]: efeito aleatório de inclinação da variedade sobre o preditor X_p.

Para a validação e avaliação de desempenho dos modelos, adotaram-se duas métricas de exatidão preditiva: Erro Médio Absoluto (MAE) e Raiz do Erro Quadrático Médio (RMSE). O ganho explicativo foi mensurado pela proporção de variância retida, utilizando o Pseudo-R² Marginal (R²m), que representa a variância explicada pelos efeitos fixos, e o Pseudo-R² Condicional (R²c), que representa a variância total explicada pelo modelo completo (efeitos fixos combinados aos componentes aleatórios).

3. Resultados e Discussão

Decomposição da Variância: M0 ao M3

A análise inicial, utilizando o modelo nulo (M0), revelou que 21,83% da variabilidade da produtividade da cana-de-açúcar (TCH) estava associada à localização espacial (talhão), validando a necessidade da abordagem multinível, conforme detalhado na Tabela 3.

Tabela 3. Evolução do ajuste estatístico dos modelos na etapa incondicional (M0 a M3).

Modelo	σ² (talhão)	σ² (variedade)	σ² (safra)	σ² (resíduos)	BIC	ΔBIC
M0	131,00			468,99	272.658
M1	95,20	166,90		441,85	269.967	2.691
M2	169,87		58,87	395,48	269.755	2.903
M3	182,73	303,13	94,86	318,25	265.139	7.519

Fonte: Dados originais da pesquisa. Nota: N=29.674 observações. Estrutura de agrupamento: talhão (n = 10.336), variedade (n = 30) e safra (n = 5). ΔBIC calculado em relação ao modelo nulo (M0).

A inclusão da variedade no modelo M1 resultou em uma redução de 27,3% na variância do talhão, indicando um viés de posicionamento varietal. Por outro lado, a inserção independente da safra no modelo M2 aumentou a variância do talhão em 29,7%, sugerindo que a flutuação temporal pode mascarar a heterogeneidade espacial intrínseca. No modelo M3, com classificação cruzada, a variedade emergiu como o principal componente de variância da produtividade (33,72%), superando o talhão (20,33%) e a safra (10,55%). Essa estrutura foi ratificada pela redução do BIC, de 272.659 no M0 para 265.134, com uma variância residual de 318,25 (t ha⁻¹)².

Variáveis de Controle Agronômico: M4 e M5

A inserção dos controles agronômicos de número de corte e ambiente de produção no modelo M4 reduziu a variância da variedade de 303,13 para 51,14 (t ha⁻¹)², corrigindo um viés de variável omitida e demonstrando que a variabilidade produtiva atribuída à variedade no modelo incondicional era, em parte, decorrente desses fatores (Tabela 4).

Tabela 4. Evolução do ajuste estatístico dos modelos na etapa condicional (M4 a M5).

Modelo	σ²t (talhão)	σ²v (variedade)	σ²s (safra)	σ²e (resíduos)	BIC	ΔBIC
M3	182,73	303,13	94,86	318,25	265.139	7.519
M4	80,83	51,14	41,25	212,81	250.926	14.213
M5	92,53	48,85	31,89	177,27	247.267	3.659

Fonte: Dados originais da pesquisa. Nota: ΔBIC calculado em relação ao modelo anterior.

A adição da idade da cultura no modelo M5 reduziu a variância residual para 177,27 (t ha⁻¹)² e revelou um efeito de supressão sobre a variância intrínseca do talhão, que aumentou para 92,53 (t ha⁻¹)². Ambas as adições (M4 e M5) foram validadas por sucessivas reduções no BIC.

Seleção de Variáveis e Redução de Dimensionalidade

O ranqueamento das 38 características biofísicas e agrometeorológicas via mRMR (Figura 1) indicou a radiação solar na janela F2 como a de maior relevância informacional. Contudo, a avaliação de VIF revelou multicolinearidade severa entre as 10 variáveis pré-selecionadas, com a radiação solar na F2 apresentando um VIF de 84,36. Após quatro rodadas de exclusão, seis variáveis remanescentes de baixa colinearidade foram selecionadas: IRECI na F2 (VIF = 1,28), EVI na F1 (VIF = 1,65), radiação solar na F1 (VIF = 1,87), precipitação na F1 (VIF = 2,10) e na F2 (VIF = 2,88), e evapotranspiração de referência na F2 (VIF = 2,92).

Figura 1. Ranking das dez variáveis pré-selecionadas via mRMR

Fonte: Dados originais da pesquisa.

Efeitos Fixos: M6

A inserção das seis variáveis selecionadas no modelo M5, centralizadas na grande média, foi realizada pelo procedimento Forward Stepwise (Tabela 5). O índice biofísico IRECI (F2) foi o que mais contribuiu, promovendo uma redução de 9.829 pontos no BIC. As adições subsequentes resultaram em reduções marginais decrescentes, culminando com a precipitação (F1). Ao todo, as inclusões reduziram o BIC de 247.267 para 235.774, resultando no modelo M6.

Tabela 5. Inclusão iterativa de variáveis biofísicas e agrometeorológicas via procedimento Forward Stepwise (M6).

Passo	Variável Adicionada	Log-Likelihood	LRT	BIC	ΔBIC
M5		-123.551	< 0,001	247.267
1	ireci_auc_f2_cent	– 118.632	< 0,001	237.438	9.829
2	rs_acc_f1_cent	– 118.354	< 0,001	236.893	545
3	evi_auc_f1_cent	– 117.927	< 0,001	236.050	843
4	prec_mm_acc_f2_cent	– 117.853	< 0,001	235.911	139
5	eto_mm_acc_f2_cent	– 117.799	< 0,001	235.815	96
6	prec_mm_acc_f1_cent	– 117.774	< 0,001	235.774	41

Fonte: Dados originais da pesquisa. Nota: ΔBIC calculado em relação ao passo iterativo imediatamente anterior. O sufixo _cent indica que a variável foi centralizada na grande média.

Essa estruturação resultou no maior ganho explicativo do estudo (Tabela 6). A inclusão simultânea dos preditores biofísicos e agrometeorológicos reduziu a variância residual em 26,1% e a variância da safra em 90,7%. Esse decréscimo indicou que a flutuação interanual da produtividade foi quase integralmente explicada pelos preditores de efeito fixo, restringindo o efeito aleatório da safra a 1,43% da variância total do modelo. Esse comportamento corrobora a premissa de que a variabilidade temporal da produtividade é predominantemente governada pelas condições agrometeorológicas e pela resposta biofísica da cultura em cada ciclo (Marin et al., 2008; Silva, M.G.S. et al., 2022).

Tabela 6. Evolução do ajuste estatístico do modelo na adição de efeitos fixos (M6).

Modelo	σ²t (talhão)	σ²v (variedade)	σ²s (safra)	σ²e (resíduos)	BIC	ΔBIC
M5	92,53	48,85	31,89	177,27	247.267	3.659
M6	42,40	29,66	2,95	130,91	235.774	11.493

Fonte: Dados originais da pesquisa. Nota: ΔBIC calculado em relação ao modelo anterior.

Inclinações Aleatórias: M7

A etapa final de estruturação avaliou a sensibilidade da cultura a variáveis específicas por meio da inclusão de inclinações aleatórias para o contexto da variedade. A estratégia Forward Stepwise identificou o índice NDRE na F1 e o déficit de pressão de vapor na F2 como as duas variáveis que atenderam aos critérios de validação (Tabela 7).

Tabela 7. Resumo das iterações válidas para inclusão de inclinações aleatórias.

Rodada	Preditor sob Efeito Aleatório	Efeito Fixo Substituído	ΔBIC	VPC (%)
1	vpd_kpa_avg_f1_cent	eto_mm_acc_f2_cent	82	2,37
1	vpd_kpa_avg_f2_cent	eto_mm_acc_f2_cent	282	4,46
1	ndre_auc_f1_cent	evi_auc_f1_cent	378	3,22
1	ndre_auc_f2_cent	ireci_auc_f2_cent	157	3,74
2	vpd_kpa_avg_f1_cent	eto_mm_acc_f2_cent	121	1,99
2	vpd_kpa_avg_f2_cent	eto_mm_acc_f2_cent	223	4,06

Fonte: Dados originais da pesquisa. Nota: ΔBIC calculado em relação ao respectivo modelo basal. Tabela somente com os preditores avaliados que superaram os critérios de validação.

No escopo dos efeitos fixos, a inserção dessas novas inclinações exigiu a remoção prévia do EVI (F1), da evapotranspiração de referência (F2) e da precipitação (F1) devido à multicolinearidade e perda de significância estatística (p = 0,062). A análise dos efeitos aleatórios do modelo M7 evidenciou a interação genótipo-ambiente. As inclinações do NDRE na F1 e do déficit de pressão de vapor na F2 retiveram, respectivamente, 2,41% e 4,06% da variância total. Esses preditores atuaram como indicadores quantitativos da heterogeneidade genética, com a dinâmica do NDRE na F1 refletindo o fechamento do dossel e a eficiência inicial de absorção de nitrogênio (Shendryk et al., 2020), e a sensibilidade ao déficit de pressão de vapor na F2 sinalizando a tolerância ao estresse higrométrico (Inman-Bamber e Smith, 2005). A modelagem dessa interação reduziu a variância residual para 125,79 (t ha⁻¹)² e atenuou a variância da safra para 1,85 (t ha⁻¹)², representando 0,86% da variância total (Tabela 8).

Tabela 8. Evolução do ajuste estatístico do modelo na adição de efeitos aleatórios de inclinação (M7).

Modelo	σ²t (talhão)	σ²v (variedade)	σ²v0	σ²v1	σ²v2	σ²w (safra)	σ²e (resíduos)	BIC	ΔBIC
M6	42,40	29,66				2,95	130,91	235.774	11.493
M7	45,94	31,00	5,20	8,74		1,85	125,79	235.032	742

Fonte: Dados originais da pesquisa. Nota: ΔBIC calculado em relação ao modelo anterior. σ²v0: variância do intercepto aleatório da variedade; σ²v1: variância da inclinação aleatória do NDRE (F1); σ²v2: variância da inclinação aleatória do VPD (F2).

Validação e Avaliação de Desempenho dos Modelos

Para quantificar a evolução do poder preditivo, o modelo M3 serviu como linha de base comparativa (Tabela 9). Seu Pseudo-R² Condicional (R²c) indicou que os contextos de talhão, variedade e safra explicaram 45,7% da variância total da produtividade, com MAE de 12,63 t ha⁻¹ e RMSE de 15,88 t ha⁻¹.

Tabela 9. Métricas de avaliação de desempenho preditivo e critérios de informação dos modelos multinível ajustados para estimativa de TCH.

Modelo	BIC	LRT	R²m	R²c	MAE (t ha⁻¹)	RMSE (t ha⁻¹)
M3: de Classificação Cruzada	265.139	< 0,001	0,000	0,457	12,63	15,88
M4: com Controle Agronômico	250.926	< 0,001	0,475	0,664	10,33	13,26
M5: com Controle de Idade	247.267	< 0,001	0,521	0,730	9,25	11,91
M6: com Efeitos Fixos	235.774	< 0,001	0,709	0,802	7,89	10,48
M7: com Inclinação Aleatória	235.033	< 0,001	0,713	0,811	7,72	10,24

Fonte: Dados originais da pesquisa

A adição dos controles agronômicos e da idade da cultura (M4 e M5) incrementou a capacidade explicativa marginal. O M5 capturou 52,1% da variância por meio de seus efeitos fixos (R²m), elevando a explicação global (R²c) para 73,0% e reduzindo o RMSE para 11,91 t ha⁻¹. O maior ganho marginal ocorreu no M6, com a inclusão das variáveis biofísicas e agrometeorológicas, que elevou a variância explicada pelos efeitos fixos (R²m) para 70,9% e reduziu o MAE para 7,89 t ha⁻¹.

A inclusão de inclinações aleatórias no modelo M7 consolidou a interação genótipo-ambiente, validada pelo LRT (p < 0,001) e com o menor BIC (235.033). O modelo definitivo alcançou MAE de 7,72 t ha⁻¹ e RMSE de 10,24 t ha⁻¹, com uma explicação global (R²c) de 81,1%. A diferença de 9,8% entre o R²c e o R²m (0,713) representou a variabilidade atribuída aos contextos de talhão, variedade e safra, justificando a abordagem multinível (Gelman e Hill, 2006).

A verificação dos pressupostos de homocedasticidade e linearidade do modelo M7, por meio de análise gráfica, mostrou que a dispersão dos resíduos padronizados não evidenciou padrões severos de afunilamento (Figura 2), validando a constância da variância. A linha de tendência não paramétrica manteve-se próxima a zero, atestando a ausência de viés sistemático nas predições (Zuur et al., 2009).

Figura 2. Dispersão dos resíduos padronizados em função dos valores ajustados para diagnóstico de homocedasticidade do modelo M7.

Fonte: Dados originais da pesquisa.

A distribuição dos resíduos apresentou aderência aos quantis teóricos na porção central (Figura 3), embora com desvios nas caudas (leptocurtose). Tais desvios não inviabilizam a inferência estatística em modelos multiníveis, dada a robustez assintótica garantida pelo Teorema do Limite Central frente ao tamanho expressivo da amostra (Zuur et al., 2009).

Figura 3. Distribuição de densidade de probabilidade e gráfico de quantis-quantis para verificação da normalidade dos resíduos do modelo M7.

Fonte: Dados originais da pesquisa.

Na avaliação de desempenho preditivo, a dispersão das estimativas (Figura 4) validou a capacidade explicativa do modelo, apresentando alta aderência à reta de identidade. As métricas globais de erro (MAE = 7,72 t ha⁻¹ e RMSE = 10,24 t ha⁻¹) e o coeficiente de determinação (R²preditivo = 0,824) atestaram a adequação prática do modelo.

Figura 4. Relação entre a produtividade de cana-de-açúcar (TCH) observada e os valores preditos pelo modelo M7.

Fonte: Dados originais da pesquisa.

Em comparação a técnicas de Machine Learning, os resultados obtidos equipararam-se aos de algoritmos de Random Forest, que registraram RMSE entre 11,70 e 15,56 t ha⁻¹ e MAE de 4,57 a 7,53 t ha⁻¹ (Amaro, 2023; Bocca e Rodrigues, 2016). A abordagem multinível apresentou precisão preditiva equivalente, com a vantagem inferencial de permitir a interpretação direta dos parâmetros agronômicos, quantificando a partição de variância e explicitando a interação genótipo-ambiente.

A extração dos coeficientes de efeitos fixos revelou a magnitude estimada dos impactos de manejo, ambiente, clima e estado biofísico (Tabela 10). O intercepto global situou-se em 98,06 ± 1,35 t ha⁻¹, representando a produtividade média basal teórica de um canavial em primeiro corte (cana-planta), no ambiente de maior potencial (ambiente A), sob preditores contínuos centralizados nas suas respectivas médias globais.

Tabela 10. Estimativas de parâmetros de efeitos fixos do modelo M7 para predição da produtividade (TCH).

Parâmetro	Média Global	Estimativa ± EP	IC (95%)	t	p-valor
produtividade média basal		98,06 ± 1,35	95,41; 100,70	72,73	< 0,001
Controles de Manejo
2º corte		-8,99 ± 0,29	-9,55; -8,42	-31,11	< 0,001
3º corte		-18,82 ± 0,30	-19,40; -18,25	-63,87	< 0,001
4º corte		-23,60 ± 0,31	-24,21; -22,99	-76,19	< 0,001
5º corte		-25,00 ± 0,34	-25,67; -24,32	-72,62	< 0,001
6º corte		-26,91 ± 0,42	-27,73; -26,09	-64,44	< 0,001
7º corte		-28,94 ± 0,65	-30,22; -27,66	-44,42	< 0,001
ambiente B		-4,22 ± 0,57	-5,35; -3,10	-7,35	<0,001
ambiente C		-9,69 ± 0,58	-10,83; -8,55	-16,61	< 0,001
ambiente D		-12,27 ± 0,59	-13,43; -11,12	-20,81	< 0,001
ambiente E		-15,32 ± 0,63	-16,56; -14,08	-24,18	< 0,001
Preditores Contínuos
idade (dias)	387	0,079 ± 0,001	0,076; 0,081	62,69	< 0,001
NDRE (F1)	25,73	0,464 ± 0,053	0,360; 0,568	8,74	< 0,001
IRECI (F2)	94,66	0,292 ± 0,004	0,286; 0,299	82,67	< 0,001
déficit de pressão de vapor (F2)	1,32	-11,40 ± 2,75	-16,79; -6,00	-4,14	< 0,001
precipitação (F2)	527,79	0,007 ± 0,001	0,006; 0,008	10,08	< 0,001
radiação solar (F1)	1.899,53	-0,015 ± 0,001	-0,015; -0,014	-35,35	< 0,001

Fonte: Dados originais da pesquisa. Nota: EP = Erro Padrão; IC = Intervalo de Confiança; t = Estatística t; Variáveis contínuas centralizadas na grande média.

A sucessão dos ciclos evidenciou um decréscimo não linear na produtividade, com o segundo corte resultando em uma redução estimada de 8,99 ± 0,29 t ha⁻¹, e o sétimo corte em uma penalização de 28,94 ± 0,65 t ha⁻¹. A transição do ambiente A para o E (restritivo) subtraiu 15,32 ± 0,63 t ha⁻¹ da produtividade média basal. Quanto à idade da cultura, o modelo estimou um ganho de 0,079 ± 0,001 t ha⁻¹ para cada dia adicional, indicando que um incremento de 1 t ha⁻¹ demanda aproximadamente 13 dias no ciclo.

No escopo biofísico, o NDRE (F1) apresentou coeficiente de 0,464 ± 0,053 t ha⁻¹, validando sua relação com a eficiência inicial de absorção de nitrogênio (Shendryk et al., 2020) e a velocidade de fechamento do dossel (Dong et al., 2019). O IRECI (F2) corroborou a predição do crescimento vegetativo com coeficiente de 0,292 ± 0,004 t ha⁻¹, ratificando a capacidade da banda Red-Edge em quantificar a biomassa com resistência à saturação ótica em dosséis densos (Dong et al., 2019). Para esse preditor, o ganho unitário de produtividade está associado a um aumento médio de 0,023 no índice durante a janela F2.

Entre as variáveis agroclimatológicas, o déficit de pressão de vapor (F2) destacou-se como fator restritivo, com coeficiente de -11,40 ± 2,75 t ha⁻¹, indicando que o estresse higrométrico restringe a assimilação de carbono (Inman-Bamber e Smith, 2005). Em oposição, a precipitação (F2) apresentou coeficiente de 0,007 ± 0,001 t ha⁻¹, sugerindo que um volume suplementar de 143 mm (28,6 mm mês⁻¹) é necessário para o mesmo incremento produtivo. A radiação solar na F1 apresentou coeficiente negativo de -0,015 ± 0,001 t ha⁻¹, confirmando sua função como proxy para o estresse termohídrico de dias excessivamente limpos, que aceleram a desidratação (Inman-Bamber e Smith, 2005).

As inclinações aleatórias evidenciaram a heterogeneidade varietal. A variedade IAC SP01 5503 apresentou o maior desvio positivo de intercepto (11,31 t ha⁻¹), elevando seu patamar de referência para 109,37 t ha⁻¹ sob condições de referência. Em contrapartida, a CTC 7515BT registrou um desvio negativo severo (-14,69 t ha⁻¹), indicando menor teto produtivo basal. A Figura 5 ilustra a relação entre o potencial produtivo basal e a resposta varietal ao NDRE na F1.

Figura 5. Relação entre o potencial produtivo basal e a resposta varietal ao NDRE na F1.

Fonte: Dados originais da pesquisa.

Na janela de grande crescimento vegetativo (F2), a avaliação do déficit de pressão de vapor mapeou o perfil de vulnerabilidade ao estresse higrométrico. A IAC SP01 5503 demonstrou extrema sensibilidade, com BLUP de -36,74 t ha⁻¹, resultando em uma penalização de 48,14 t ha⁻¹ para cada 1 kPa de incremento nesse déficit (Figura 6), o que sugere sua inadequação para áreas com deficiência hídrica. Em contraste, genótipos como a CTC 20 exibiram alta rusticidade fisiológica, com inclinação aleatória fortemente positiva (24,49 t ha⁻¹), anulando a penalidade global imposta pelo estresse termohídrico e justificando sua alocação estratégica em ambientes mais restritivos.

Figura 6. Relação entre o potencial produtivo basal e a sensibilidade varietal ao déficit de pressão de vapor na F2.

Fonte: Dados originais da pesquisa.

Em síntese, a modelagem multinível com classificação cruzada demonstrou ser uma abordagem robusta para a predição da produtividade da cana-de-açúcar, permitindo a decomposição da variância em componentes espaciais, genéticos e temporais. A integração de preditores biofísicos e agrometeorológicos capturou a dinâmica fenológica da cultura, e o modelo final atingiu exatidão preditiva comparável a algoritmos de Machine Learning, mantendo a interpretabilidade dos coeficientes agronômicos. Os resultados indicaram a influência significativa do número de corte, ambiente de produção, idade da cultura e variáveis biofísicas e agrometeorológicas na produtividade, além de evidenciar a heterogeneidade varietal na resposta a fatores como o NDRE e o déficit de pressão de vapor.

4. Conclusão

A modelagem multinível com classificação cruzada demonstrou ser uma abordagem analítica robusta e metodologicamente adequada para a predição da produtividade da cana-de-açúcar na região sudeste de Mato Grosso do Sul, abrangendo as safras de 2019/2020 a 2023/2024. Verificou-se que a abordagem particionou adequadamente a variância espacial, genética e temporal, quantificando com êxito a interação genótipo-ambiente. O ajuste final com inclinações aleatórias (M7) explicou 81,1% da variância total e alcançou uma Raiz do Erro Quadrático Médio de 10,24 t ha⁻¹. A diferença de 9,8% entre o Pseudo-R² marginal e o condicional evidenciou o relevante ganho explicativo gerado pela classificação cruzada, equiparando seu desempenho preditivo à exatidão reportada por algoritmos de Machine Learning. Essa abordagem, contudo, manteve o rigor paramétrico e a interpretabilidade dos preditores agronômicos, como o número de corte, ambiente de produção, idade da cultura, índices biofísicos (NDRE e IRECI) e variáveis agrometeorológicas (déficit de pressão de vapor, precipitação e radiação solar). A principal contribuição reside na oferta de uma ferramenta analítica transparente e robusta para a estimativa de produtividade e o planejamento do setor sucroenergético, permitindo a compreensão detalhada da influência de fatores contextuais e a heterogeneidade varietal.

Apesar da robustez atestada, o estudo apresenta limitações inerentes ao escopo dos modelos empíricos. A extrapolação das inferências restringe-se ao portfólio varietal e às condições edafoclimáticas específicas da área de estudo. Observaram-se desvios na distribuição dos resíduos nas caudas extremas, o que sinaliza a ocorrência de anomalias fenológicas, variáveis agrometeorológicas omitidas ou interações complexas que a estrutura linear paramétrica não conseguiu absorver em sua totalidade. Para desenvolvimentos futuros, sugere-se a expansão do domínio amostral, a incorporação de estruturas de autocorrelação espacial para capturar a dependência entre talhões próximos e a transição para Modelos Aditivos Generalizados Mistos (GAMM). Essas iniciativas permitirão refinar o mapeamento de extremos climáticos e consolidar a técnica como uma valiosa ferramenta de inteligência analítica no setor sucroenergético.

Referências Bibliográficas

Associação de Produtores de Bioenergia de Mato Grosso do Sul [BIOSUL]. 2024. Mapa da Bioenergia de Mato Grosso do Sul. Disponível em: . Acesso em: 20 mar. 2025.

Avilez, A.Μ.Α.; Hernandez, F.B.T.; Bispo, R.C.; Giovanelli, L.B. 2018. Necessidade hídrica da cana-de-açúcar no Noroeste Paulista. Irriga 1(1): 171-188.

Empresa de Pesquisa Energética [EPE]. 2024a. Balanço Energético Nacional 2024: Ano base 2023. EPE, Rio de Janeiro, RJ, Brasil.

Empresa de Pesquisa Energética [EPE]. 2024b. Nota técnica: Análise de conjuntura dos biocombustíveis – ano 2023. EPE, Rio de Janeiro, RJ, Brasil.

Manzatto, C.V.; Assad, E.D.; Bacca, J.F.M.; Zaroni, M.J.; Pereira, S.E.M. 2009. Zoneamento agroecológico da cana-de-açúcar: expandir a produção, preservar a vida, garantir o futuro.

Artigo oriundo de Trabalho de Conclusão de Curso da Especialização em Data Science e Analytics do MBA USP/Esalq

Para saber mais sobre o curso, clique aqui e acesse a plataforma MBX Academy