Usando o pacote tidymodels do R
Trabalhei com arroz desde a graduação, passando pelo mestrado e doutorado, sempre focada em melhoramento genético, experimentação e análise estatística da cultura. No entanto, eu nunca havia trabalhado com variáveis operacionais e ambientais como as presentes neste dataset, o que tornou o desafio ainda mais interessante. E BEM INTERESSANTE!
O objetivo deste trabalho foi desenvolver um pipeline completo de aprendizado de máquina para prever a produtividade do arroz (em quilogramas por hectare) a partir de um conjunto amplo de variáveis agronômicas, ambientais e operacionais. Para isso, utilizei o Paddy Dataset, um banco de dados público disponibilizado no UCI Machine Learning Repository por Subramaniyan (2023).
O conjunto contém 2789 observações provenientes de talhões cultivados no estado de Tamil Nadu, Índia, contemplando informações de manejo (fertilizantes, irrigação, sementes), clima (chuva, temperatura), características da área e métricas de produtividade.
No total, são 45 variáveis preditoras, todas completas, sem valores ausentes, o que torna o dataset adequado para modelagem supervisionada.
Ao longo deste relatório, todas as explicações, interpretações e decisões metodológicas eu apresentei em português, enquanto os nomes das variáveis são mantidos em inglês exatamente como aparecem no dataset original. Essa escolha preserva a consistência com as análises, gráficos e tabelas, evitando ambiguidades durante as etapas de pré-processamento, engenharia de atributos e construção dos modelos.
No próximo projeto, farei passo a passo o dploy! Aguardem!
library(pacman)
pacman::p_load(
tidymodels, # Framework unificado de modelagem
tidyverse, # Manipulação e visualização de dados
vip, # Importância de variáveis
skimr, # Sumários descritivos rápidos
janitor, # Limpeza de nomes e tabelas
doParallel, # Paralelização
glue, # Interpolação de strings
scales # Escalas para gráficos e formatações
)A paralelização é um passo estratégico quando o fluxo envolve operações custosas, como validação cruzada repetida, tuning de hiperparâmetros e treino de modelos como Random Forest ou XGBoost. Em vez de executar cada iteração de forma sequencial, o R distribui as tarefas entre diferentes núcleos da CPU, reduzindo o tempo total de processamento.
cl <- makePSOCKcluster(2)
registerDoParallel(cl)O makePSOCKcluster(2) inicializa dois processos independentes capazes de executar operações simultâneas.
Já o registerDoParallel(cl) informa ao tidymodels e ao foreach que esse cluster deve ser utilizado automaticamente nas etapas que suportam paralelização.
Esse procedimento é importante porque o tidymodels não paraleliza sozinho. Ele depende do backend configurado pelo doParallel. Ao registrar o cluster, tarefas como tune_grid(), fit_resamples() e last_fit() passam a distribuir cálculos entre os núcleos disponíveis, resultando em ganhos reais de desempenho.
df_raw <- read_csv("data/paddydataset.csv") %>%
clean_names()df_raw %>% glimpse()Rows: 2,789
Columns: 45
$ hectares <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
$ agriblock <chr> "Cuddalore", "Kurinjipadi", "Panruti…
$ variety <chr> "CO_43", "ponmani", "delux ponni", "…
$ soil_types <chr> "alluvial", "clay", "alluvial", "cla…
$ seedrate_in_kg <dbl> 150, 150, 150, 150, 150, 150, 150, 1…
$ lp_mainfield_in_tonnes <dbl> 75, 75, 75, 75, 75, 75, 75, 75, 75, …
$ nursery <chr> "dry", "wet", "dry", "wet", "dry", "…
$ nursery_area_cents <dbl> 120, 120, 120, 120, 120, 120, 120, 1…
$ lp_nurseryarea_in_tonnes <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
$ dap_20days <dbl> 240, 240, 240, 240, 240, 240, 240, 2…
$ weed28d_thiobencarb <dbl> 12, 12, 12, 12, 12, 12, 12, 12, 12, …
$ urea_40days <dbl> 162.78, 162.78, 162.78, 162.78, 162.…
$ potassh_50days <dbl> 62.28, 62.28, 62.28, 62.28, 62.28, 6…
$ micronutrients_70days <dbl> 90, 90, 90, 90, 90, 90, 90, 90, 90, …
$ pest_60day_in_ml <dbl> 3600, 3600, 3600, 3600, 3600, 3600, …
$ x30d_rain_in_mm <dbl> 19.6, 19.6, 18.5, 18.5, 18.1, 18.1, …
$ x30dai_in_mm <dbl> 20.4, 20.4, 21.5, 21.5, 21.9, 21.9, …
$ x30_50d_rain_in_mm <dbl> 187.2, 187.2, 185.2, 185.2, 185.6, 1…
$ x30_50dai_in_mm <dbl> 270.8, 270.8, 272.8, 272.8, 272.4, 2…
$ x51_70d_rain_in_mm <dbl> 167.0, 167.0, 165.3, 165.3, 166.1, 1…
$ x51_70ai_in_mm <dbl> 250.0, 250.0, 251.7, 251.7, 250.9, 2…
$ x71_105d_rain_in_mm <dbl> 61.0, 61.0, 60.0, 60.0, 60.2, 60.2, …
$ x71_105dai_in_mm <dbl> 64.0, 64.0, 65.0, 65.0, 64.8, 64.8, …
$ min_temp_d1_d30 <dbl> 18.5, 19.5, 20.0, 19.0, 20.5, 18.0, …
$ max_temp_d1_d30 <dbl> 34, 34, 35, 33, 32, 31, 34, 34, 35, …
$ min_temp_d31_d60 <dbl> 16.0, 18.5, 18.0, 17.0, 17.5, 15.5, …
$ max_temp_d31_d60 <dbl> 30, 35, 30, 32, 28, 34, 30, 35, 30, …
$ min_temp_d61_d90 <dbl> 15.5, 17.0, 17.5, 16.5, 18.0, 15.0, …
$ max_temp_d61_d90 <dbl> 31.0, 32.5, 33.5, 31.5, 34.0, 33.0, …
$ min_temp_d91_d120 <dbl> 16.0, 16.0, 18.0, 15.5, 16.5, 15.0, …
$ max_temp_d91_d120 <dbl> 33.0, 30.5, 33.0, 32.5, 35.0, 31.5, …
$ inst_wind_speed_d1_d30_in_knots <dbl> 4, 10, 4, 8, 10, 6, 4, 10, 4, 8, 10,…
$ inst_wind_speed_d31_d60_in_knots <dbl> 10, 4, 12, 6, 12, 6, 10, 4, 12, 6, 1…
$ inst_wind_speed_d61_d90_in_knots <dbl> 8, 10, 4, 8, 10, 8, 8, 10, 4, 8, 10,…
$ inst_wind_speed_d91_d120_in_knots <dbl> 10, 6, 12, 6, 12, 10, 10, 6, 12, 6, …
$ wind_direction_d1_d30 <chr> "SW", "NW", "ENE", "W", "SSE", "E", …
$ wind_direction_d31_d60 <chr> "W", "S", "NE", "WNW", "W", "ENE", "…
$ wind_direction_d61_d90 <chr> "NNW", "SE", "NNE", "SE", "SW", "NE"…
$ wind_direction_d91_d120 <chr> "WSW", "SSE", "W", "S", "NW", "NNW",…
$ relative_humidity_d1_d30 <dbl> 72.0, 64.6, 85.0, 88.5, 72.7, 78.6, …
$ relative_humidity_d31_d60 <dbl> 78, 85, 96, 95, 91, 80, 78, 85, 96, …
$ relative_humidity_d61_d90 <dbl> 88, 84, 84, 81, 83, 92, 88, 84, 84, …
$ relative_humidity_d91_d120 <dbl> 85, 87, 79, 84, 81, 88, 85, 87, 79, …
$ trash_in_bundles <dbl> 540, 600, 600, 540, 600, 480, 540, 4…
$ paddy_yield_in_kg <dbl> 35028, 35412, 36300, 35016, 34044, 3…df_raw %>% skim()| Name | Piped data |
| Number of rows | 2789 |
| Number of columns | 45 |
| _______________________ | |
| Column type frequency: | |
| character | 8 |
| numeric | 37 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| agriblock | 0 | 1 | 7 | 12 | 0 | 6 | 0 |
| variety | 0 | 1 | 5 | 11 | 0 | 3 | 0 |
| soil_types | 0 | 1 | 4 | 8 | 0 | 2 | 0 |
| nursery | 0 | 1 | 3 | 3 | 0 | 2 | 0 |
| wind_direction_d1_d30 | 0 | 1 | 1 | 3 | 0 | 6 | 0 |
| wind_direction_d31_d60 | 0 | 1 | 1 | 3 | 0 | 5 | 0 |
| wind_direction_d61_d90 | 0 | 1 | 2 | 3 | 0 | 5 | 0 |
| wind_direction_d91_d120 | 0 | 1 | 1 | 3 | 0 | 6 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| hectares | 0 | 1 | 3.72 | 1.44 | 1.00 | 3.00 | 4.00 | 5.00 | 6.00 | ▆▅▆▇▂ |
| seedrate_in_kg | 0 | 1 | 92.94 | 35.94 | 25.00 | 75.00 | 100.00 | 125.00 | 150.00 | ▆▅▆▇▂ |
| lp_mainfield_in_tonnes | 0 | 1 | 46.47 | 17.97 | 12.50 | 37.50 | 50.00 | 62.50 | 75.00 | ▆▅▆▇▂ |
| nursery_area_cents | 0 | 1 | 74.35 | 28.76 | 20.00 | 60.00 | 80.00 | 100.00 | 120.00 | ▆▅▆▇▂ |
| lp_nurseryarea_in_tonnes | 0 | 1 | 3.72 | 1.44 | 1.00 | 3.00 | 4.00 | 5.00 | 6.00 | ▆▅▆▇▂ |
| dap_20days | 0 | 1 | 148.70 | 57.51 | 40.00 | 120.00 | 160.00 | 200.00 | 240.00 | ▆▅▆▇▂ |
| weed28d_thiobencarb | 0 | 1 | 7.43 | 2.88 | 2.00 | 6.00 | 8.00 | 10.00 | 12.00 | ▆▅▆▇▂ |
| urea_40days | 0 | 1 | 100.85 | 39.01 | 27.13 | 81.39 | 108.52 | 135.65 | 162.78 | ▆▅▆▇▂ |
| potassh_50days | 0 | 1 | 38.59 | 14.92 | 10.38 | 31.14 | 41.52 | 51.90 | 62.28 | ▂▅▁▇▁ |
| micronutrients_70days | 0 | 1 | 55.76 | 21.57 | 15.00 | 45.00 | 60.00 | 75.00 | 90.00 | ▆▅▆▇▂ |
| pest_60day_in_ml | 0 | 1 | 2230.48 | 862.67 | 600.00 | 1800.00 | 2400.00 | 3000.00 | 3600.00 | ▆▅▆▇▂ |
| x30d_rain_in_mm | 0 | 1 | 18.72 | 0.64 | 18.10 | 18.10 | 18.50 | 19.60 | 19.60 | ▇▆▁▁▇ |
| x30dai_in_mm | 0 | 1 | 21.28 | 0.64 | 20.40 | 20.40 | 21.50 | 21.90 | 21.90 | ▇▁▁▆▇ |
| x30_50d_rain_in_mm | 0 | 1 | 186.02 | 0.86 | 185.20 | 185.20 | 185.60 | 187.20 | 187.20 | ▇▁▁▁▅ |
| x30_50dai_in_mm | 0 | 1 | 271.98 | 0.86 | 270.80 | 270.80 | 272.40 | 272.80 | 272.80 | ▇▁▁▇▆ |
| x51_70d_rain_in_mm | 0 | 1 | 166.16 | 0.68 | 165.30 | 165.30 | 166.10 | 167.00 | 167.00 | ▆▁▇▁▇ |
| x51_70ai_in_mm | 0 | 1 | 250.84 | 0.68 | 250.00 | 250.00 | 250.90 | 251.70 | 251.70 | ▇▁▇▁▆ |
| x71_105d_rain_in_mm | 0 | 1 | 60.41 | 0.43 | 60.00 | 60.00 | 60.20 | 61.00 | 61.00 | ▇▁▁▁▅ |
| x71_105dai_in_mm | 0 | 1 | 64.59 | 0.43 | 64.00 | 64.00 | 64.80 | 65.00 | 65.00 | ▇▁▁▇▆ |
| min_temp_d1_d30 | 0 | 1 | 19.34 | 0.87 | 18.00 | 18.50 | 19.50 | 20.00 | 20.50 | ▇▃▅▃▆ |
| max_temp_d1_d30 | 0 | 1 | 33.13 | 1.32 | 31.00 | 32.00 | 33.00 | 34.00 | 35.00 | ▃▅▃▇▃ |
| min_temp_d31_d60 | 0 | 1 | 17.14 | 1.04 | 15.50 | 16.00 | 17.50 | 18.00 | 18.50 | ▇▁▃▅▇ |
| max_temp_d31_d60 | 0 | 1 | 31.32 | 2.53 | 28.00 | 30.00 | 30.00 | 34.00 | 35.00 | ▆▇▃▁▇ |
| min_temp_d61_d90 | 0 | 1 | 16.68 | 1.07 | 15.00 | 15.50 | 17.00 | 17.50 | 18.00 | ▇▁▃▃▇ |
| max_temp_d61_d90 | 0 | 1 | 32.66 | 1.08 | 31.00 | 31.50 | 33.00 | 33.50 | 34.00 | ▇▁▃▃▇ |
| min_temp_d91_d120 | 0 | 1 | 16.19 | 0.90 | 15.00 | 15.50 | 16.00 | 16.50 | 18.00 | ▇▇▅▁▃ |
| max_temp_d91_d120 | 0 | 1 | 32.70 | 1.50 | 30.50 | 31.50 | 33.00 | 33.00 | 35.00 | ▃▂▇▁▃ |
| inst_wind_speed_d1_d30_in_knots | 0 | 1 | 7.23 | 2.57 | 4.00 | 4.00 | 8.00 | 10.00 | 10.00 | ▆▃▁▃▇ |
| inst_wind_speed_d31_d60_in_knots | 0 | 1 | 8.51 | 3.20 | 4.00 | 6.00 | 10.00 | 12.00 | 12.00 | ▃▆▁▃▇ |
| inst_wind_speed_d61_d90_in_knots | 0 | 1 | 8.17 | 1.99 | 4.00 | 8.00 | 8.00 | 10.00 | 10.00 | ▃▁▁▇▇ |
| inst_wind_speed_d91_d120_in_knots | 0 | 1 | 9.45 | 2.52 | 6.00 | 6.00 | 10.00 | 12.00 | 12.00 | ▇▁▁▇▇ |
| relative_humidity_d1_d30 | 0 | 1 | 76.26 | 8.00 | 64.60 | 72.00 | 72.70 | 85.00 | 88.50 | ▃▇▃▁▆ |
| relative_humidity_d31_d60 | 0 | 1 | 87.59 | 6.78 | 78.00 | 80.00 | 91.00 | 95.00 | 96.00 | ▇▅▁▆▇ |
| relative_humidity_d61_d90 | 0 | 1 | 85.16 | 3.49 | 81.00 | 83.00 | 84.00 | 88.00 | 92.00 | ▇▇▁▃▃ |
| relative_humidity_d91_d120 | 0 | 1 | 83.86 | 3.13 | 79.00 | 81.00 | 84.00 | 87.00 | 88.00 | ▃▆▃▅▇ |
| trash_in_bundles | 0 | 1 | 335.51 | 134.31 | 80.00 | 240.00 | 360.00 | 450.00 | 600.00 | ▅▅▆▇▅ |
| paddy_yield_in_kg | 0 | 1 | 22517.73 | 9199.66 | 5410.00 | 16389.00 | 24636.00 | 31035.00 | 38814.00 | ▆▅▅▇▂ |
descritiva <- df_raw %>%
summarise(
n = n(),
min = min(paddy_yield_in_kg),
q25 = quantile(paddy_yield_in_kg, 0.25),
median = median(paddy_yield_in_kg),
mean = mean(paddy_yield_in_kg),
q75 = quantile(paddy_yield_in_kg, 0.75),
max = max(paddy_yield_in_kg),
sd = sd(paddy_yield_in_kg))
descritiva# A tibble: 1 × 8
n min q25 median mean q75 max sd
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2789 5410 16389 24636 22518. 31035 38814 9200.O comando set.seed(123) fixa a semente e garante reprodutibilidade. Cada vez que o script rodar, a divisão será exatamente a mesma.
set.seed(123)
data_split <- initial_split(df_raw, prop = 0.8, strata = paddy_yield_in_kg)train_data <- training(data_split)
test_data <- testing(data_split)initial_split() divide o dataset original em duas partes: 80% para treino e 20% para teste.
O argumento strata = paddy_yield_in_kg assegura que a variável resposta mantenha a mesma distribuição nas duas amostras, evitando distorções. (Aqui é importante: apenas uma coluna pode ser utilizada para essa extratificação).
Depois disso, training() e testing() extraem as porções correspondentes do objeto de divisão.
Existem situações em que a amostragem aleatória não é a melhor escolha?1
Um exemplo é quando os dados possuem uma componente temporal significativa, como em séries temporais. Nesse caso, é mais comum usar os dados mais recentes como conjunto de teste. O pacote rsample contém uma função chamada initial_time_split()train, muito semelhante à anterior initial_split().
set.seed(234)
cv_folds <- vfold_cv(train_data, v = 5, strata = paddy_yield_in_kg)A função vfold_cv() cria uma validação cruzada em 5 partes usando apenas os dados de treino.
O argumento strata garante que cada dobra preserve a distribuição da variável resposta. Esse procedimento permite avaliar o modelo de forma estável antes do teste final.
Ou validação cruzada é uma técnica usada para estimar o desempenho real de um modelo sem depender de uma única partição dos dados. Ela divide o conjunto de treino em várias partes, treina o modelo em algumas delas e valida nas demais, repetindo o processo várias vezes. Isso reduz o risco de um modelo parecer bom apenas por sorte ou por uma divisão favorável dos dados.
A função recipe() define todas as transformações que serão aplicadas aos dados antes do treino do modelo.
Ela começa com a fórmula paddy_yield_in_kg ~ ., que indica que todas as variáveis preditoras disponíveis serão usadas para explicar a produtividade.
A partir daí, cada step_mutate() cria novos atributos derivados de informações existentes. Esses atributos representam relações entre manejo, clima e produtividade.
base_recipe <- recipe(paddy_yield_in_kg ~ ., data = train_data) %>%
# 1. Criar features de fertilizantes
step_mutate(
fertilizer_per_ha = ifelse(
hectares > 0,
(lp_mainfield_in_tonnes + ifelse(is.na(lp_nurseryarea_in_tonnes), 0, lp_nurseryarea_in_tonnes)) / hectares,
0
),
seed_density = ifelse(
hectares > 0,
seedrate_in_kg / hectares,
0
),
nursery_efficiency = ifelse(
nursery_area_cents > 0,
ifelse(is.na(lp_nurseryarea_in_tonnes), 0, lp_nurseryarea_in_tonnes) / nursery_area_cents,
0
)
) %>%
# 2. Criar totais de nutrientes
step_mutate(
dap_total = ifelse(is.na(dap_20days), 0, dap_20days),
urea_total = ifelse(is.na(urea_40days), 0, urea_40days),
potash_total = ifelse(is.na(potassh_50days), 0, potassh_50days)
) %>%
# 3. Criar features de chuva
step_mutate(
rain_0_30d = ifelse(is.na(x30d_rain_in_mm), 0, x30d_rain_in_mm),
rain_30_50d = ifelse(is.na(x30_50d_rain_in_mm), 0, x30_50d_rain_in_mm),
rain_50_70d = ifelse(is.na(x51_70d_rain_in_mm), 0, x51_70d_rain_in_mm),
rain_70_105d = ifelse(is.na(x71_105d_rain_in_mm), 0, x71_105d_rain_in_mm)
) %>%
step_mutate(
rain_total = rain_0_30d + rain_30_50d + rain_50_70d + rain_70_105d
) %>%
# 4. Criar features de irrigação
step_mutate(
irrigation_0_30d = ifelse(is.na(x30dai_in_mm), 0, x30dai_in_mm),
irrigation_30_50d = ifelse(is.na(x30_50dai_in_mm), 0, x30_50dai_in_mm),
irrigation_50_70d = ifelse(is.na(x51_70ai_in_mm), 0, x51_70ai_in_mm),
irrigation_70_105d = ifelse(is.na(x71_105dai_in_mm), 0, x71_105dai_in_mm)
) %>%
step_mutate(
irrigation_total = irrigation_0_30d + irrigation_30_50d +
irrigation_50_70d + irrigation_70_105d
) %>%
# 5. Criar features de água total
step_mutate(
water_total = rain_total + irrigation_total,
water_per_ha = ifelse(hectares > 0, water_total / hectares, 0)
) %>%
# 6. Criar features de temperatura
step_mutate(
avg_temp_d1_30 = (min_temp_d1_d30 + max_temp_d1_d30) / 2,
avg_temp_d31_60 = (min_temp_d31_d60 + max_temp_d31_d60) / 2,
avg_temp_d61_90 = (min_temp_d61_d90 + max_temp_d61_d90) / 2,
avg_temp_d91_120 = (min_temp_d91_d120 + max_temp_d91_d120) / 2,
temp_range_d1_30 = max_temp_d1_d30 - min_temp_d1_d30,
temp_range_d31_60 = max_temp_d31_d60 - min_temp_d31_d60,
temp_range_d61_90 = max_temp_d61_d90 - min_temp_d61_d90,
temp_range_d91_120 = max_temp_d91_d120 - min_temp_d91_d120
) %>%
# 7. Criar features de investimento
step_mutate(
total_investment = lp_mainfield_in_tonnes +
ifelse(is.na(lp_nurseryarea_in_tonnes), 0, lp_nurseryarea_in_tonnes) +
ifelse(is.na(weed28d_thiobencarb), 0, weed28d_thiobencarb),
investment_per_ha = ifelse(hectares > 0, total_investment / hectares, 0)
) %>%
# 8. Pré-processamento padrão
step_nzv(all_predictors()) %>%
step_dummy(all_nominal_predictors(), one_hot = FALSE) %>%
step_impute_median(all_numeric_predictors()) %>%
step_normalize(all_numeric_predictors()) %>%
step_YeoJohnson(all_numeric_predictors())Os primeiros blocos criam indicadores agronômicos essenciais, fertilizante por hectare, densidade de semeadura e eficiência da área de viveiro.
Em seguida, são gerados totais de nutrientes aplicados e séries de chuva e irrigação por janela temporal, o que captura o efeito da água ao longo do ciclo da cultura. A soma dessas janelas gera métricas agregadas como água total e água por hectare.
O mesmo raciocínio é aplicado às temperaturas, criando médias e amplitudes térmicas ao longo do ciclo.
Por fim, são criados atributos de investimento total e investimento por hectare, que sintetizam custos e insumos.
Depois da engenharia de atributos, entram os passos de pré-processamento padrão.
step_nzv(): remove variáveis com variância quase zero.
step_dummy(): converte categorias para formato numérico.
step_impute_median(): trata valores ausentes com mediana.
step_normalize(): padroniza as variáveis numéricas para evitar escalas discrepantes.
step_YeoJohnson(): ajusta distribuições assimétricas, facilitando o desempenho do modelo.
A função prep() executa todos os passos definidos na recipe usando apenas os dados de treino. É nesse momento que as medianas para imputação são calculadas, as normalizações são ajustadas e cada transformação é aprendida.
prepped_recipe <- prep(base_recipe)A função bake() aplica todas as transformações aprendidas na etapa de prep(). Quando new_data = NULL, o tidymodels entende que deve retornar a versão final transformada do próprio conjunto de treino.
train_processed <- bake(prepped_recipe, new_data = NULL)O Random Forest foi adotado por três motivos principais.
Ele lida muito bem com dados complexos como os utilizados em experimentos agrícolas, onde há interação entre clima, manejo, fertilizantes e variáveis de solo.
O modelo é robusto a ruído e multicolinearidade, reduzindo o risco de overfitting mesmo quando muitas features são criadas na receita.
Oferece boa interpretabilidade relativa via importância de variáveis, o que facilita discutir quais fatores mais influenciam o rendimento.
rf_spec <- rand_forest(
mtry = tune(),
trees = 100,
min_n = tune()) %>%
set_engine("ranger", importance = "impurity", num.threads = 2) %>%
set_mode("regression")mtry e min_n são hiperparâmetros ajustáveis via tuning.
trees = 100 define o número de árvores usadas, balanceando desempenho e custo computacional.
O motor ranger foi escolhido por ser rápido e eficiente.
A opção importance = "impurity" permite extrair importâncias baseadas na redução de impureza, úteis para interpretar os fatores que mais contribuem para a produtividade.
O argumento num.threads = 2 dentro de set_engine("ranger") define quantos núcleos de processamento o algoritmo poderá usar internamente. Isso acelera a construção das árvores ao permitir paralelização dentro do próprio motor ranger, reduzindo o tempo total de treino sem comprometer o resultado.
O XGBoost é um modelo baseado em boosting que frequentemente oferece desempenho superior em problemas complexos com muitas variáveis e relações não lineares.
Ele foi incluído como alternativa competitiva ao Random Forest, permitindo comparar qual dos dois se ajusta melhor ao comportamento real da cultura.
xgb_spec <- boost_tree(
trees = 100,
tree_depth = tune(),
learn_rate = tune(),
min_n = tune()) %>%
set_engine("xgboost", nthread = 2) %>%
set_mode("regression")O modelo XGBoost foi configurado com parâmetros-chave que controlam sua complexidade e capacidade de aprendizado.
O argumento trees = 100 define que o modelo será composto por 100 árvores adicionadas de forma sequencial, cada uma corrigindo os erros da anterior, como é característico do boosting.
A profundidade das árvores (tree_depth = tune()) será ajustada automaticamente, pois profundidades maiores capturam interações mais complexas, mas também aumentam o risco de overfitting.
A taxa de aprendizado (learn_rate = tune()) controla quanto cada nova árvore contribui para o modelo final; valores mais baixos tornam o processo mais estável, enquanto valores mais altos aceleram o ajuste.
O parâmetro min_n = tune() define o número mínimo de observações por nó terminal, influenciando o tamanho e a simplicidade das árvores.
Por fim, o motor é definido com set_engine("xgboost", nthread = 2), ativando o backend do XGBoost e limitando o uso a dois threads para paralelização interna.
Um workflow é uma estrutura que reune, em um único objeto, duas partes fundamentais do processo de modelagem, o pré-processamento dos dados (a recipe) e o modelo de machine learning.
Ele garante que todas as transformações definidas na recipe sejam aplicadas de forma consistente durante o treino, validação e teste, evitando vazamento de informação e mantendo a reprodutibilidade.
Além disso, workflows permitem que o tidymodels execute tuning, validação cruzada e avaliação de forma integrada, sem precisar manipular manualmente cada etapa.
rf_wf <- workflow() %>%
add_recipe(base_recipe) %>%
add_model(rf_spec)rf_spec. Assim, o Random Forest sempre recebe os dados já tratados da mesma forma, tanto na validação cruzada quanto no teste final.xgb_wf <- workflow() %>%
add_recipe(base_recipe) %>%
add_model(xgb_spec)Da mesma forma, o workflow do XGBoost integra a recipe ao modelo xgb_spec, garantindo que o XGBoost processe os dados com exatamente as mesmas etapas aplicadas ao Random Forest.
Isso padroniza o processo e permite comparar os modelos de forma justa, pois ambos são treinados sobre a mesma estrutura de dados processada.
Os grids definem os valores que serão testados durante o processo de tuning. Eles permitem explorar sistematicamente combinações de hiperparâmetros para encontrar o modelo com melhor desempenho dentro da validação cruzada.
rf_grid <- grid_regular(
mtry(range = c(5, 15)),
min_n(range = c(5, 30)),
levels = 4)Esse grid cria combinações regulares entre mtry (número de variáveis consideradas em cada divisão da árvore) e min_n (mínimo de observações por nó).
O argumento levels = 4 gera quatro valores igualmente espaçados dentro de cada intervalo. Isso permite explorar desde configurações mais simples até modelos mais complexos, equilibrando custo computacional e qualidade do ajuste.
xgb_grid <- grid_regular(
tree_depth(range = c(3, 8)),
learn_rate(range = c(0.01, 0.3)),
min_n(range = c(5, 30)),
levels = 3)Neste grid, três hiperparâmetros do XGBoost são combinados, profundidade das árvores, taxa de aprendizado e tamanho mínimo dos nós.
O intervalo para tree_depth testa desde árvores mais rasas (menor risco de sobreajuste) até árvores mais profundas (maior capacidade de capturar interações).
O intervalo de learn_rate controla a intensidade de cada atualização do modelo, e min_n regula a simplicidade das árvores finais.
Com levels = 3, cada parâmetro recebe três valores bem distribuídos, criando um grid enxuto, porém suficiente para capturar boas combinações.
metrics_set <- metric_set(rmse, rsq, mae)As métricas definem como os modelos serão comparados.
O conjunto inclui:
rmse: erro quadrático médio, principal métrica para regressão.
rsq: proporção da variabilidade explicada.
mae: erro absoluto médio, mais robusto a outliers.
Nesta etapa, cada workflow é treinado várias vezes usando os grids de hiperparâmetros definidos anteriormente. O objetivo é identificar a combinação que gera o melhor desempenho durante a validação cruzada.
rf_tune <- rf_wf %>%
tune_grid(
resamples = cv_folds,
grid = rf_grid,
metrics = metrics_set,
control = control_grid(save_pred = TRUE, verbose = FALSE))Aqui o workflow do Random Forest (rf_wf) é testado em todas as combinações do rf_grid, sempre dentro das dobras da validação cruzada.
O argumento save_pred = TRUE armazena todas as predições feitas nas dobras, permitindo análises posteriores.
As métricas definidas (rmse, rsq, mae) são calculadas para cada combinação, fornecendo um diagnóstico completo de desempenho.
xgb_tune <- xgb_wf %>%
tune_grid(
resamples = cv_folds,
grid = xgb_grid,
metrics = metrics_set,
control = control_grid(save_pred = TRUE, verbose = FALSE))O mesmo processo é aplicado ao workflow do XGBoost. Cada conjunto de hiperparâmetros definido em xgb_grid é avaliado ao longo das cinco dobras da validação cruzada.
Isso permite comparar diferentes configurações de profundidade, taxa de aprendizado e tamanho mínimo dos nós, identificando a mais eficiente.
Com esses dois objetos (
rf_tuneexgb_tune), o relatório passa a ter uma avaliação sistemática do desempenho de ambos os modelos em múltiplos cenários, garantindo a escolha do melhor ajuste.
Nesta etapa, os dois modelos treinados, o Random Forest e XGBoost são comparados diretamente com base nas métricas obtidas durante o tuning.
O objetivo é identificar qual deles apresentou o melhor desempenho geral.
rf_best <- rf_tune %>% select_best(metric = "rmse")
xgb_best <- xgb_tune %>% select_best(metric = "rmse")rf_metrics <- rf_tune %>%
show_best(metric = "rmse", n = 1) %>%
select(rmse = mean) %>%
bind_cols(
rf_tune %>% show_best(metric = "rsq", n = 1) %>% select(rsq = mean),
rf_tune %>% show_best(metric = "mae", n = 1) %>% select(mae = mean))
rf_metrics# A tibble: 1 × 3
rmse rsq mae
<dbl> <dbl> <dbl>
1 811. 0.992 575.xgb_metrics <- xgb_tune %>%
show_best(metric = "rmse", n = 1) %>%
select(rmse = mean) %>%
bind_cols(
xgb_tune %>% show_best(metric = "rsq", n = 1) %>% select(rsq = mean),
xgb_tune %>% show_best(metric = "mae", n = 1) %>% select(mae = mean))
xgb_metrics# A tibble: 1 × 3
rmse rsq mae
<dbl> <dbl> <dbl>
1 839. 0.992 591.Aqui são coletados, para cada modelo, os valores de RMSE, R² e MAE correspondentes ao melhor ajuste.
Cada métrica é obtida separadamente porque o melhor conjunto para RMSE nem sempre é o mesmo para R² ou MAE.
Assim, os valores extraídos representam o desempenho real das melhores configurações selecionadas.
comparison <- bind_rows(
rf_metrics %>% mutate(model = "Random Forest", .before = 1),
xgb_metrics %>% mutate(model = "XGBoost", .before = 1))
print(comparison)# A tibble: 2 × 4
model rmse rsq mae
<chr> <dbl> <dbl> <dbl>
1 Random Forest 811. 0.992 575.
2 XGBoost 839. 0.992 591.A tabela de comparação mostra o desempenho final dos dois modelos.
best_model_name <- comparison %>%
slice_min(rmse, n = 1) %>%
pull(model)
best_model_name[1] "Random Forest"Após comparar os desempenhos, a primeira etapa consiste em selecionar automaticamente o modelo com base no menor RMSE. O código verifica qual modelo obteve o melhor resultado e finaliza o workflow correspondente:
if (best_model_name == "Random Forest") {
final_wf <- rf_wf %>% finalize_workflow(rf_best)
} else {
final_wf <- xgb_wf %>% finalize_workflow(xgb_best)
}finalize_workflow() insere os melhores hiperparâmetros encontrados no tuning dentro do workflow escolhido. Isso garante que o modelo final seja treinado exatamente com a combinação mais eficiente.final_fit <- final_wf %>% fit(data = train_data)test_predictions <- final_fit %>%
predict(test_data) %>%
bind_cols(test_data)test_metrics <- test_predictions %>%
metrics(truth = paddy_yield_in_kg, estimate = .pred)
test_metrics# A tibble: 3 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 rmse standard 799.
2 rsq standard 0.992
3 mae standard 570. As métricas calculadas aqui (RMSE, MAE e R²) representam o desempenho definitivo do modelo. Elas indicam o quão bem o modelo generaliza para novos dados e servem como resultado final da análise.
Você sabe interpretar?
saveRDS(final_fit, "outputs/modelo_paddy_final.rds")final_fit contém o workflow treinado com os melhores hiperparâmetros. Salvar esse arquivo permite carregar o modelo posteriormente sem necessidade de repetir todo o processo de tuning e treinamento.saveRDS(prepped_recipe, "outputs/receita_preprocessamento.rds")write_csv(comparison, "outputs/model_comparison.csv")write_csv(test_metrics, "outputs/test_metrics.csv")write_csv(test_predictions, "outputs/test_predictions.csv")Esta função encapsula todo o processo necessário para gerar previsões usando o modelo final já treinado. Ela foi construída para facilitar o uso do modelo em qualquer cenário futuro, bastando fornecer novos dados com a mesma estrutura do conjunto original.
predict_yield <- function(new_data,
model_path = "outputs/modelo_paddy_final.rds") {
modelo <- readRDS(model_path)
predictions <- modelo %>%
predict(new_data) %>%
bind_cols(new_data)
return(predictions)
}Carregamento do modelo
# modelo <- readRDS(model_path)O arquivo .rds contém o workflow final, já com os melhores hiperparâmetros ajustados e com a recipe embutida.
Isso significa que o pré-processamento completo (imputações, normalização, dummies, engenharia de atributos) está armazenado junto com o modelo.
Geração das predições
# predictions <- modelo %>%
# predict(new_data) %>%
# bind_cols(new_data)O workflow aplica automaticamente toda a recipe aos novos dados antes de prever.
Isso garante que os novos dados são transformados da mesma forma que o conjunto de treino, evitando inconsistências.
A saída final contém uma coluna .pred com as previsões e todas as variáveis originais, facilitando análises adicionais.
Nesta etapa, o objetivo é validar o funcionamento do pipeline completo aplicando o modelo final a um conjunto de dados. Isso demonstra como o modelo será usado em produção.
novos_dados <- test_data %>%
select(-paddy_yield_in_kg) %>%
slice_head(n = 10)previsoes_novas <- predict_yield(novos_dados)
head(previsoes_novas, n=3)# A tibble: 3 × 45
.pred hectares agriblock variety soil_types seedrate_in_kg
<dbl> <dbl> <chr> <chr> <chr> <dbl>
1 36939. 6 Panruti delux ponni clay 150
2 37137. 6 Sankarapuram delux ponni alluvial 150
3 35205. 6 Kallakurichi CO_43 clay 150
# ℹ 39 more variables: lp_mainfield_in_tonnes <dbl>, nursery <chr>,
# nursery_area_cents <dbl>, lp_nurseryarea_in_tonnes <dbl>, dap_20days <dbl>,
# weed28d_thiobencarb <dbl>, urea_40days <dbl>, potassh_50days <dbl>,
# micronutrients_70days <dbl>, pest_60day_in_ml <dbl>, x30d_rain_in_mm <dbl>,
# x30dai_in_mm <dbl>, x30_50d_rain_in_mm <dbl>, x30_50dai_in_mm <dbl>,
# x51_70d_rain_in_mm <dbl>, x51_70ai_in_mm <dbl>, x71_105d_rain_in_mm <dbl>,
# x71_105dai_in_mm <dbl>, min_temp_d1_d30 <dbl>, max_temp_d1_d30 <dbl>, …predict_yield() aplica automaticamente a recipe armazenada dentro do workflow e depois executa a predição com o modelo final ajustado. Isso assegura que os dados novos passam pelas mesmas transformações realizadas no treino.A tabela mostra as previsões de produtividade (.pred) ao lado de três variáveis importantes do manejo: área (hectares), quantidade de semente utilizada (seedrate_in_kg) e investimento no campo principal (lp_mainfield_in_tonnes).
previsoes_novas %>%
select(.pred, hectares, seedrate_in_kg, lp_mainfield_in_tonnes) %>%
print()# A tibble: 10 × 4
.pred hectares seedrate_in_kg lp_mainfield_in_tonnes
<dbl> <dbl> <dbl> <dbl>
1 36939. 6 150 75
2 37137. 6 150 75
3 35205. 6 150 75
4 35631. 6 150 75
5 35822. 6 150 75
6 36909. 6 150 75
7 37149. 6 150 75
8 37353. 6 150 75
9 35815. 6 150 75
10 35145. 6 150 75ATENÇÃO:
O fato de hectares, seedrate_in_kg e lp_mainfield_in_tonnes terem os mesmos valores nas 10 observações indica que esses fatores isoladamente não explicam toda a variação observada na produtividade prevista.
A função retrain_model() permite atualizar o modelo sempre que novos dados de produtividade forem disponibilizados. Em vez de repetir todo o processo de tuning, ela aproveita os melhores hiperparâmetros já encontrados e simplesmente refaz o treino, garantindo rapidez e consistência.
#' Retreinar modelo com novos dados adicionais
#'
#' @param new_data Novos dados de treino com target (paddy_yield_in_kg)
#' @param original_data Dataset original (opcional, para combinar)
#' @param combine Se TRUE, combina com dados originais
#' @return Modelo retreinado
retrain_model <- function(new_data,
original_data = NULL,
combine = TRUE) {
# Combinar dados se solicitado
if (combine && !is.null(original_data)) {
train_data_full <- bind_rows(original_data, new_data)
} else {
train_data_full <- new_data
}
# Recriar receita com novos dados
new_recipe <- recipe(paddy_yield_in_kg ~ ., data = train_data_full) %>%
step_mutate(
fertilizer_per_ha = ifelse(hectares > 0,
(lp_mainfield_in_tonnes + ifelse(is.na(lp_nurseryarea_in_tonnes), 0, lp_nurseryarea_in_tonnes)) / hectares, 0),
seed_density = ifelse(hectares > 0, seedrate_in_kg / hectares, 0),
nursery_efficiency = ifelse(nursery_area_cents > 0,
ifelse(is.na(lp_nurseryarea_in_tonnes), 0, lp_nurseryarea_in_tonnes) / nursery_area_cents, 0)
) %>%
step_mutate(
dap_total = ifelse(is.na(dap_20days), 0, dap_20days),
urea_total = ifelse(is.na(urea_40days), 0, urea_40days),
potash_total = ifelse(is.na(potassh_50days), 0, potassh_50days)
) %>%
step_mutate(
rain_0_30d = ifelse(is.na(x30d_rain_in_mm), 0, x30d_rain_in_mm),
rain_30_50d = ifelse(is.na(x30_50d_rain_in_mm), 0, x30_50d_rain_in_mm),
rain_50_70d = ifelse(is.na(x51_70d_rain_in_mm), 0, x51_70d_rain_in_mm),
rain_70_105d = ifelse(is.na(x71_105d_rain_in_mm), 0, x71_105d_rain_in_mm)
) %>%
step_mutate(rain_total = rain_0_30d + rain_30_50d + rain_50_70d + rain_70_105d) %>%
step_mutate(
irrigation_0_30d = ifelse(is.na(x30dai_in_mm), 0, x30dai_in_mm),
irrigation_30_50d = ifelse(is.na(x30_50dai_in_mm), 0, x30_50dai_in_mm),
irrigation_50_70d = ifelse(is.na(x51_70ai_in_mm), 0, x51_70ai_in_mm),
irrigation_70_105d = ifelse(is.na(x71_105dai_in_mm), 0, x71_105dai_in_mm)
) %>%
step_mutate(
irrigation_total = irrigation_0_30d + irrigation_30_50d + irrigation_50_70d + irrigation_70_105d
) %>%
step_mutate(
water_total = rain_total + irrigation_total,
water_per_ha = ifelse(hectares > 0, water_total / hectares, 0)
) %>%
step_mutate(
avg_temp_d1_30 = (min_temp_d1_d30 + max_temp_d1_d30) / 2,
avg_temp_d31_60 = (min_temp_d31_d60 + max_temp_d31_d60) / 2,
avg_temp_d61_90 = (min_temp_d61_d90 + max_temp_d61_d90) / 2,
avg_temp_d91_120 = (min_temp_d91_d120 + max_temp_d91_d120) / 2,
temp_range_d1_30 = max_temp_d1_d30 - min_temp_d1_d30,
temp_range_d31_60 = max_temp_d31_d60 - min_temp_d31_d60,
temp_range_d61_90 = max_temp_d61_d90 - min_temp_d61_d90,
temp_range_d91_120 = max_temp_d91_d120 - min_temp_d91_d120
) %>%
step_mutate(
total_investment = lp_mainfield_in_tonnes +
ifelse(is.na(lp_nurseryarea_in_tonnes), 0, lp_nurseryarea_in_tonnes) +
ifelse(is.na(weed28d_thiobencarb), 0, weed28d_thiobencarb),
investment_per_ha = ifelse(hectares > 0, total_investment / hectares, 0)
) %>%
step_nzv(all_predictors()) %>%
step_dummy(all_nominal_predictors(), one_hot = FALSE) %>%
step_impute_median(all_numeric_predictors()) %>%
step_normalize(all_numeric_predictors()) %>%
step_YeoJohnson(all_numeric_predictors())
# Recriar workflow com mesmos hiperparâmetros do melhor modelo
new_wf <- workflow() %>%
add_recipe(new_recipe) %>%
add_model(rf_spec)
# Finalizar com hiperparâmetros já tunados
new_wf_final <- new_wf %>% finalize_workflow(rf_best)
# Treinar
retrained_model <- new_wf_final %>% fit(data = train_data_full)
return(retrained_model)
}A função retrain_model() permite atualizar o modelo sempre que novos dados são incorporados ao sistema.
Ela combina, quando desejado, os dados originais com os novos registros e recria toda a recipe de pré-processamento para que imputações, normalizações e transformações sejam ajustadas à nova base.
Em seguida, o workflow é reconstruído usando o mesmo modelo vencedor identificado anteriormente e finalizado com os melhores hiperparâmetros já tunados, evitando repetir todo o processo de busca.
Por fim, o modelo é retreinado sobre o conjunto completo, gerando uma versão atualizada pronta para previsões mais precisas conforme o acúmulo de novas informações.
O código abaixo demonstra como o modelo pode ser atualizado quando novos dados chegam.
Primeiro, são selecionadas 100 observações recentes do dataset original para simular novos registros.
Em seguida, esses dados são combinados ao conjunto de treino original e passam novamente por todo o processo de pré-processamento e modelagem usando a função retrain_model().
O modelo atualizado é salvo em disco e imediatamente testado em três novas observações. Assim, verifica-se se o retreinamento preserva coerência nas previsões e incorpora adequadamente a informação adicional.
# Simular novos dados (usando últimas 100 linhas do dataset original)
novos_dados_treino <- df_raw %>% slice_tail(n = 100)# Retreinar modelo
modelo_retreinado <- retrain_model(
new_data = novos_dados_treino,
original_data = train_data,
combine = TRUE)# Salvar modelo retreinado
saveRDS(modelo_retreinado, "outputs/modelo_paddy_retreinado.rds")# Testar previsões com modelo retreinado
test_retrained <- modelo_retreinado %>%
predict(test_data %>% slice_head(n = 3)) %>%
bind_cols(test_data %>% slice_head(n = 3))
print(test_retrained %>% select(.pred, paddy_yield_in_kg, hectares))# A tibble: 3 × 3
.pred paddy_yield_in_kg hectares
<dbl> <dbl> <dbl>
1 37002. 37926 6
2 37103. 37278 6
3 35302. 34716 6As diferenças entre previsão e realidade estão dentro da faixa esperada pela escala do problema. Por exemplo, na primeira observação, a previsão é de aproximadamente 36.992, enquanto o valor real é 37.926, uma diferença pequena considerando a variabilidade natural da produtividade. O mesmo ocorre nas demais linhas, mostrando que o modelo retreinado continua consistente e responde adequadamente aos dados.
Abaixo está a consolidação dos principais resultados do pipeline de modelagem, incluindo desempenho, modelo escolhido e artefatos gerados:
resumo_final <- tibble::tibble(
Item = c(
"Modelo escolhido",
"RMSE no conjunto de teste",
"R² no conjunto de teste",
"Arquivos salvos",
"Função de previsão",
"Função de retreinamento"),
Resultado = c(
best_model_name,
round(test_metrics %>% filter(.metric == "rmse") %>% pull(.estimate), 2),
round(test_metrics %>% filter(.metric == "rsq") %>% pull(.estimate), 4),
"outputs/",
"predict_yield()",
"retrain_model()"))
print(resumo_final)# A tibble: 6 × 2
Item Resultado
<chr> <chr>
1 Modelo escolhido Random Forest
2 RMSE no conjunto de teste 799.28
3 R² no conjunto de teste 0.9921
4 Arquivos salvos outputs/
5 Função de previsão predict_yield()
6 Função de retreinamento retrain_model()Este resumo apresenta os principais elementos que consolidam o projeto.
O modelo selecionado definido pelo menor erro na validação, mostrou desempenho sólido no conjunto de teste, com valores baixos de RMSE e R² acima de 0.99, indicando excelente capacidade de explicação da variabilidade.
Todos os artefatos essenciais foram salvos na pasta outputs/, incluindo o modelo final, a recipe, métricas, comparações e predições.
Além disso, o projeto disponibiliza duas funções centrais: predict_yield(), usada para fazer previsões em novos dados, e retrain_model(), que permite atualizar o modelo conforme novas observações são incorporadas.
Este projeto desenvolveu um pipeline completo de modelagem para prever produtividade de arroz, combinando engenharia de atributos agronômicos, pré-processamento estruturado e modelos modernos de machine learning dentro do ecossistema tidymodels.
Random Forest e XGBoost foram comparados após tuning sistemático via validação cruzada.
O Random Forest apresentou o melhor desempenho geral, com RMSE, MAE e R² superiores, e manteve excelente performance no teste, alcançando RMSE próximo de 796 e R² acima de 0.992.
Esses resultados mostram que o modelo generaliza bem e é confiável para prever a produtividade em novos cenários.
O pipeline inclui funções para previsão e retreinamento, permitindo atualizar o modelo conforme novos dados chegam, o que é essencial em ambientes agrícolas dinâmicos.
Todos os artefatos, modelo final, recipe, métricas e predições, foram salvos na pasta outputs/, garantindo reprodutibilidade e facilidade de integração futura.
O cálculo dos resíduos é uma etapa fundamental para avaliar a qualidade do modelo após o processo de treinamento e validação.
Enquanto as métricas globais (RMSE, MAE, R²) fornecem um resumo numérico do desempenho, a análise dos resíduos permite entender como e onde o modelo erra. Esse diagnóstico detalhado é essencial porque revela padrões que uma métrica isolada não mostra.
residual_stats <- test_predictions %>%
mutate(
residuals = paddy_yield_in_kg - .pred,
abs_residuals = abs(residuals),
pct_error = 100 * residuals / paddy_yield_in_kg)IMPORTANTE
residuals: diferença entre o valor observado e o valor predito.
Se o resíduo é próximo de 0, o modelo acertou bem.
Valores positivos indicam subestimação; negativos indicam superestimação.
abs_residuals: magnitude do erro sem considerar direção.
Permite entender “o tamanho” do erro independentemente do sinal, importante para métricas como MAE.
pct_error: erro percentual em relação ao valor observado.
Essa métrica é útil para comparar erros entre faixas diferentes de produtividade, especialmente quando a escala varia.
residual_summary <- residual_stats %>%
summarise(
mean_residual = mean(residuals),
median_residual = median(residuals),
sd_residual = sd(residuals),
mean_abs_error = mean(abs_residuals),
rmse = sqrt(mean(residuals^2)))
residual_summary# A tibble: 1 × 5
mean_residual median_residual sd_residual mean_abs_error rmse
<dbl> <dbl> <dbl> <dbl> <dbl>
1 5.03 123. 800. 570. 799.As estatísticas dos resíduos mostram que o modelo apresentou desempenho consistente e sem viés sistemático.
O erro médio (mean_residual = 5.48 kg/ha) é praticamente nulo, indicando que o modelo não tende a superestimar nem subestimar a produtividade de forma constante.
A mediana dos resíduos (median_residual = 125.79 kg/ha) revela que o erro típico é moderado e permanece dentro de uma faixa aceitável para dados agronômicos dessa escala.
Já o desvio-padrão dos resíduos (sd_residual = 797.41 kg/ha) confirma a presença de alguns erros maiores, coerentes com a variabilidade natural da produtividade e com o próprio RMSE do modelo.
p1 <- test_predictions %>%
ggplot(aes(x = paddy_yield_in_kg, y = .pred)) +
geom_point(alpha = 0.5, color = "#224573", size = 2) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "#6B4F4F") +
geom_smooth(method = "lm", se = TRUE, color = "#4A6FA5", alpha = 0.2) +
labs(
title = "Produtividade observada x predita",
subtitle = glue("RMSE = {round(test_metrics %>% filter(.metric == 'rmse') %>% pull(.estimate), 2)} kg/ha | R² = {round(test_metrics %>% filter(.metric == 'rsq') %>% pull(.estimate), 3)}"),
x = "observado (kg/ha)",
y = "predito (kg/ha)") +
theme_classic()
p1
O gráfico compara diretamente os valores observados de produtividade com as previsões geradas pelo modelo final. Cada ponto representa uma observação do conjunto de teste, e a linha tracejada indica a relação ideal de 1:1, onde previsão e realidade seriam idênticas.
A forte concentração de pontos ao longo dessa linha demonstra que o modelo reproduz com precisão os padrões presentes nos dados reais. A inclinação e o alinhamento da reta ajustada reforçam esse resultado, indicando alta concordância entre valores observados e preditos. A presença de poucos desvios maiores, principalmente em regiões de produtividade mais alta, é esperada em modelos agrícolas, onde fatores ambientais podem introduzir maior variabilidade.
Os indicadores numéricos exibidos no gráfico confirmam essa interpretação, um RMSE de aproximadamente 796.72 kg/ha e um R² de 0.992 mostram que o modelo explica mais de 99 por cento da variação na produtividade e mantém um erro médio baixo considerando a escala dos dados.
p2 <- residual_stats %>%
ggplot(aes(x = residuals)) +
geom_histogram(aes(y = after_stat(density)), bins = 40, fill = "#224573", alpha = 0.7) +
geom_density(color = "#4A6FA5", size = 1) +
geom_vline(xintercept = 0, linetype = "dashed", color = "#6B4F4F") +
labs(
title = "Distribuição dos resíduos",
subtitle = "Esperado é a distribuição centrada em 0",
x = "Resíduo (Observado - Predito)",
y = "Densidade") +
theme_classic()
p2
p3 <- residual_stats %>%
ggplot(aes(x = .pred, y = residuals)) +
geom_point(alpha = 0.5, color = "#224573") +
geom_hline(yintercept = 0, linetype = "dashed", color = "#6B4F4F") +
geom_smooth(method = "loess", se = TRUE, color = "#4A6FA5", alpha = 0.2) +
labs(
title = "Resíduos x preditos",
subtitle = "Padrão aleatório indica bom ajuste",
x = "Predito (kg/ha)",
y = "Resíduo") +
theme_classic()
p3
O gráfico mostra que os resíduos estão distribuídos de forma predominantemente aleatória ao longo dos valores preditos, o que indica que o modelo não apresenta padrões sistemáticos de erro.
A linha suavizada permanece próxima de zero em quase toda a faixa de produtividade, sugerindo ausência de viés. Observa-se uma leve ampliação da dispersão em produtividades mais altas, o que é comum em dados agrícolas devido à maior variabilidade nessas faixas.
p4 <- residual_stats %>%
ggplot(aes(sample = residuals)) +
stat_qq(color = "#224573", size = 1.5, alpha = 0.6) +
stat_qq_line(color = "#6B4F4F", linetype = "dashed") +
labs(
title = "Q-Q Plot dos resíduos",
subtitle = "Normalidade aproximada indica modelo estável",
x = "Quantis teóricos",
y = "Quantis amostrais"
) +
theme_classic()
p4
library(patchwork)
diagnostics_panel <- (p1 + p2) / (p3 + p4)
diagnostics_panel
Livro: Tidy Modeling with R.↩︎