Packages


library(multcomp) #Complementar ao emmeans, fornece valor de p e letras para os testes de medias.
## Carregando pacotes exigidos: mvtnorm
## Carregando pacotes exigidos: survival
## Carregando pacotes exigidos: TH.data
## Carregando pacotes exigidos: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(emmeans) #Faz anova para modelos mistos.
library(lmerTest) #complementar ao emmeans para funcionar a função lmer
## Carregando pacotes exigidos: lme4
## Carregando pacotes exigidos: Matrix
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(lme4)
library(multcompView)

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack

Database


library(readxl)
data = read_excel("C:/Users/Samsung/OneDrive/Documentos/2 - Material de trabalho/7 - Projetos/7.1 Pesquisa/Leguminosas IZ/Plan_PG_Luciana_VO210826.xlsx", sheet = "R table")
head(data)
#fatores
data$Planta=as.factor(data$Planta)
data$Inclusao=as.numeric(data$Inclusao)
data$DMO=as.numeric(data$DMO)
## Warning: NAs introduzidos por coerção
str(data)
## tibble [84 × 11] (S3: tbl_df/tbl/data.frame)
##  $ Planta    : Factor w/ 8 levels "Calopogonio",..: 8 8 8 5 5 5 3 3 3 1 ...
##  $ Inclusao  : num [1:84] 100 100 100 0 0 0 100 100 100 100 ...
##  $ tempo     : chr [1:84] "24h" "24h" "24h" "24h" ...
##  $ inoculo   : num [1:84] 1 1 1 1 1 1 1 1 1 1 ...
##  $ PG_Total  : num [1:84] 48.7 55.5 46.5 44.3 92.3 ...
##  $ DMO       : num [1:84] NA NA NA 242 250 ...
##  $ Acetico   : num [1:84] 50.3 78.8 22.6 45.7 42.1 ...
##  $ Propionico: num [1:84] 5.82 4.82 4.94 24.47 20.42 ...
##  $ Butirico  : num [1:84] 10.52 11.29 8.67 20.72 17.67 ...
##  $ CH4.24h   : num [1:84] 2.51 1.89 1.72 5.68 5.26 ...
##  $ AP.ratio  : num [1:84] 8.65 16.37 4.58 1.87 2.06 ...
View(data)

###### 1 - Dados Macrotiloma ######


1.1 - ANOVA + REGRESSÃO

# Filtrar apenas CM, M e M_CM
data_M = data %>% dplyr::filter(Planta %in% c("Marandu", "Macrotiloma", "Macrotiloma_Marandu"))
table(data_M$Planta)
## 
##         Calopogonio Calopogonio_Marandu         Macrotiloma Macrotiloma_Marandu 
##                   0                   0                   6                  18 
##             Marandu        Soja_Marandu         Soja_perene         Tifton_LANA 
##                   6                   0                   0                   0
View(data_M)

###################################
data_M$Inclusao_sc = (data_M$Inclusao - 50) / 25
variables = c("PG_Total","DMO","Acetico","Propionico","Butirico","CH4.24h","AP.ratio")

for (var in variables) {
  cat("\n\n##############################\n")
  cat("Variável:", var, "\n")
  cat("################################\n")
  
  # Modelo linear
  formula_lin = as.formula(paste(var, "~ Inclusao_sc + (1|inoculo)"))
  mod_lin = lmer(formula_lin, data = data_M, REML = FALSE)
  
  # Modelo quadrático
  formula_quad = as.formula(paste(var, "~ Inclusao_sc + I(Inclusao_sc^2) + (1|inoculo)"))
  mod_quad = lmer(formula_quad,data = data_M,REML = FALSE)
  
  cat("\n---------------- MODELO LINEAR ----------------\n")
  print(anova(mod_lin))
  cat("\n---------------- MODELO QUADRÁTICO ----------------\n")
  print(anova(mod_quad))
  cat("\n---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------\n")
  print(anova(mod_lin, mod_quad))}
## 
## 
## ##############################
## Variável: PG_Total 
## ################################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
## Inclusao_sc 1641.8  1641.8     1    28  7.5012 0.0106 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Inclusao_sc      1641.8  1641.8     1    28  12.811 0.0012821 ** 
## I(Inclusao_sc^2) 2540.0  2540.0     1    28  19.820 0.0001238 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_M
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod_lin     4 261.24 266.84 -126.62   253.24                         
## mod_quad    5 248.25 255.26 -119.12   238.25 14.987  1  0.0001083 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ##############################
## Variável: DMO 
## ################################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Inclusao_sc  38950   38950     1 19.017  9.0537 0.007212 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Inclusao_sc       23100   23100     1 19.016  8.5049 0.008850 **
## I(Inclusao_sc^2)  30600   30600     1 19.009 11.2662 0.003313 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_M
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)   
## mod_lin     4 248.24 252.42 -120.12   240.24                       
## mod_quad    5 241.34 246.56 -115.67   231.34 8.903  1   0.002847 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ##############################
## Variável: Acetico 
## ################################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 438.97  438.97     1    28   2.726 0.1099
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc      438.97  438.97     1    28  3.2880 0.08052 .
## I(Inclusao_sc^2) 770.67  770.67     1    28  5.7725 0.02315 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_M
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod_lin     4 252.23 257.83 -122.11   244.23                       
## mod_quad    5 248.98 255.98 -119.49   238.98 5.2484  1    0.02197 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ##############################
## Variável: Propionico 
## ################################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 0.1371  0.1371     1    28  0.0021 0.9634
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Inclusao_sc        0.14    0.14     1    28  0.0028 0.958437   
## I(Inclusao_sc^2) 399.34  399.34     1    28  8.0538 0.008354 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_M
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## mod_lin     4 226.03 231.64 -109.02   218.03                        
## mod_quad    5 220.95 227.96 -105.48   210.95 7.0786  1   0.007801 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ##############################
## Variável: Butirico 
## ################################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc 80.002  80.002     1    28  3.2116 0.08393 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc      80.002  80.002     1    28  3.2304 0.08308 .
## I(Inclusao_sc^2)  4.051   4.051     1    28  0.1636 0.68896  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_M
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod_lin     4 197.78 203.38 -94.890   189.78                     
## mod_quad    5 199.62 206.62 -94.808   189.62 0.1631  1     0.6863
## 
## 
## ##############################
## Variável: CH4.24h 
## ################################
## boundary (singular) fit: see help('isSingular')
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##              Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 0.39718 0.39718     1    30  0.0973 0.7572
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Inclusao_sc       0.397   0.397     1    28  0.1692    0.6839    
## I(Inclusao_sc^2) 49.691  49.691     1    28 21.1744 8.238e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_M
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod_lin     4 135.32 140.93 -63.662   127.32                         
## mod_quad    5 121.53 128.54 -55.766   111.53 15.792  1  7.071e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ##############################
## Variável: AP.ratio 
## ################################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Inclusao_sc 0.3687  0.3687     1    28  12.412 0.001485 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Inclusao_sc      0.36870 0.36870     1    28 12.4884 0.001443 **
## I(Inclusao_sc^2) 0.00512 0.00512     1    28  0.1734 0.680249   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_M
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## mod_lin     4 -5.7998 -0.1950 6.8999  -13.800                     
## mod_quad    5 -3.9727  3.0333 6.9863  -13.973 0.1729  1     0.6775

1.2 - PLOTS - REGRESSÃO QUADRÁTICA

library(lme4)
library(ggplot2)

variables = c("PG_Total","DMO","CH4.24h","AP.ratio")

for (var in variables) {
  # Modelo quadrático
  formula_quad = as.formula(paste(var, "~ Inclusao_sc + I(Inclusao_sc^2) + (1|inoculo)"))
  mod_quad = lmer(formula_quad,data = data_M,REML = FALSE)
  # Dados para construir a curva
  pred = data.frame(Inclusao = seq(
      min(data_M$Inclusao, na.rm = TRUE),
      max(data_M$Inclusao, na.rm = TRUE),length.out = 200))
  # Mesma transformação utilizada no modelo
  pred$Inclusao_sc = (pred$Inclusao - 50) / 25
  # Predição apenas dos efeitos fixos
  pred$y_pred = predict(mod_quad,newdata = pred,re.form = NA)
  
  # Plot
  p = ggplot(data_M,aes(x = Inclusao,y = .data[[var]])) +
    geom_point(size = 2.5,alpha = 0.7) +
    geom_line(data = pred,aes(x = Inclusao,y = y_pred),linewidth = 1) +
    scale_x_continuous(breaks = c(0, 25, 50, 75, 100)) +
    labs(x = "Inclusão (%)",y = var) +
    theme_classic(base_size = 14)
  print(p)}

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

1.3 - Ranking

cada variavel teve pesos iguais, as médias foram ajustadas via modelo simples Direções: - DMO: quanto maior melhor; - PG_Total, CH4, AP.ratio: quanto menor melhor (invertido);

EXPLICAÇÃO PARA TEXTOS CIENTIFICOS: Foi construído um índice multicritério para ranquear as inclusões de leguminosas com base na degradabilidade da matéria orgânica (DMO), na produção de gas (PG_Total), produção de metano (CH4.24h) e razão acetato:propionato (AP.ratio). As médias ajustadas foram obtidas por modelo linear considerando inclusão como efeito fixo, seguidas de estimativas de médias marginais (emmeans).Para permitir a comparação entre variáveis em escalas distintas, os valores foram padronizados por transformação min–max (0–1). Para DMO, valores maiores foram considerados mais favoráveis. Para PG_Total, CH4.24h e AP.ratio, a escala foi invertida, de modo que menores produções resultassem em maiores escores padronizados.O índice final foi calculado como a média aritmética simples dos escores padronizados (pesos iguais), e os tratamentos foram ranqueados em ordem decrescente do escore composto.

data_M_ranking = data %>% dplyr::filter(Planta %in% c("Macrotiloma_Marandu"))
View(data_M_ranking)

########################### FUNÇÃO MODELO ########################### 
fun_emm = function(data_M_ranking, response){
  f <- as.formula(paste0(response, " ~ Inclusao + factor(inoculo)"))
  mod <- lm(f, data = data_M_ranking)
  emmeans(mod,~ Inclusao,at = list(Inclusao = sort(unique(data_M_ranking$Inclusao)))) %>%
    as.data.frame() %>%
    mutate(variavel = response)}

########################### MÉDIAS AJUSTADAS (EMMEANS) ########################### 
variaveis_modelo = c("PG_Total","DMO","CH4.24h","AP.ratio")
medias1 = bind_rows(lapply(variaveis_modelo, \(v) fun_emm(data_M_ranking, v)))
View(medias1)

########################### ORGANIZAR EM FORMATO LARGO ###########################
emm_wide = medias1 %>% select(Inclusao, variavel, emmean, SE) %>%
  pivot_wider(names_from  = variavel,values_from = c(emmean, SE))
emm_wide
########################### PADRONIZAÇÃO EM ESCALA ###########################
ranking = emm_wide %>% ungroup() %>% mutate(
    # DMO (benefício: maior é melhor)
    DMO_score = (emmean_DMO - min(emmean_DMO, na.rm = TRUE)) /
                (max(emmean_DMO, na.rm = TRUE) - min(emmean_DMO, na.rm = TRUE)),
    # PG_Total (invertido: menor é melhor)
    PG_Total_score = (max(emmean_PG_Total, na.rm = TRUE) - emmean_PG_Total) /
                     (max(emmean_PG_Total, na.rm = TRUE) - min(emmean_PG_Total, na.rm = TRUE)),
    # CH4.24h (invertido: menor é melhor)
    CH4.24h_score = (max(emmean_CH4.24h, na.rm = TRUE) - emmean_CH4.24h) /
                    (max(emmean_CH4.24h, na.rm = TRUE) - min(emmean_CH4.24h, na.rm = TRUE)),
    # AP.ratio (invertido: menor é melhor)
    AP.ratio_score = (max(emmean_AP.ratio, na.rm = TRUE) - emmean_AP.ratio) /
                     (max(emmean_AP.ratio, na.rm = TRUE) - min(emmean_AP.ratio, na.rm = TRUE)))
ranking
########################### ÍNDICE MULTICRITÉRIO (PESOS IGUAIS) ###########################
ranking = ranking %>%
  mutate(Score = (DMO_score + PG_Total_score + CH4.24h_score + AP.ratio_score)/4) %>%
  arrange(desc(Score)) %>% mutate(Rank = row_number())

ranking_final = ranking %>% select(Inclusao,
         emmean_DMO, DMO_score,
         emmean_PG_Total, PG_Total_score,
         emmean_CH4.24h, CH4.24h_score,
         emmean_AP.ratio, AP.ratio_score,
         Score, Rank)
ranking_final

CONCLUSÃO

Para macrotiloma, a melhor inclusão é 25%, com maior DMO e menor PG e CH4.


###### 2 - Dados Calopogonio ######


2.1 - ANOVA + REGRESSÃO

# Filtrar apenas CM, M e M_CM
data_C = data %>% dplyr::filter(Planta %in% c("Marandu", "Calopogonio", "Calopogonio_Marandu"))
table(data_C$Planta)
## 
##         Calopogonio Calopogonio_Marandu         Macrotiloma Macrotiloma_Marandu 
##                   6                  18                   0                   0 
##             Marandu        Soja_Marandu         Soja_perene         Tifton_LANA 
##                   6                   0                   0                   0
View(data_C)

###################################
data_C$Inclusao_sc = (data_C$Inclusao - 50) / 25
variables = c("PG_Total","DMO","Acetico","Propionico","Butirico","CH4.24h","AP.ratio")

for (var in variables) {
  cat("\n\n###########################\n")
  cat("Variável:", var, "\n")
  cat("##############################\n")
  
  # Modelo linear
  formula_lin = as.formula(paste(var, "~ Inclusao_sc + (1|inoculo)"))
  mod_lin = lmer(formula_lin, data = data_C, REML = FALSE)
  
  # Modelo quadrático
  formula_quad = as.formula(paste(var, "~ Inclusao_sc + I(Inclusao_sc^2) + (1|inoculo)"))
  mod_quad = lmer(formula_quad,data = data_C,REML = FALSE)
  
  cat("\n---------------- MODELO LINEAR ----------------\n")
  print(anova(mod_lin))
  cat("\n---------------- MODELO QUADRÁTICO ----------------\n")
  print(anova(mod_quad))
  cat("\n---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------\n")
  print(anova(mod_lin, mod_quad))}
## 
## 
## ###########################
## Variável: PG_Total 
## ##############################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc 1458.8  1458.8     1    28  5.1057 0.03182 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc      1458.8  1458.8     1    28  5.9789 0.02103 *
## I(Inclusao_sc^2) 1168.6  1168.6     1    28  4.7892 0.03714 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_C
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod_lin     4 267.99 273.59 -130.00   259.99                       
## mod_quad    5 265.57 272.57 -127.78   255.57 4.4211  1     0.0355 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ###########################
## Variável: DMO 
## ##############################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Inclusao_sc  84207   84207     1    19  12.367 0.002307 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)   
## Inclusao_sc       76213   76213     1 18.998 11.8834 0.0027 **
## I(Inclusao_sc^2)   7928    7928     1 19.014  1.2361 0.2801   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_C
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar   AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)
## mod_lin     4 257.6 261.78 -124.8    249.6                     
## mod_quad    5 258.4 263.62 -124.2    248.4 1.2004  1     0.2732
## 
## 
## ###########################
## Variável: Acetico 
## ##############################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 57.753  57.753     1    28  0.4737  0.497
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc      57.753  57.753     1    28  0.4737 0.4969
## I(Inclusao_sc^2)  0.124   0.124     1    28  0.0010 0.9748
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_C
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod_lin     4 244.63 250.24 -118.32   236.63                    
## mod_quad    5 246.63 253.64 -118.32   236.63 0.001  1     0.9746
## 
## 
## ###########################
## Variável: Propionico 
## ##############################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 81.571  81.571     1    28  1.0696 0.3099
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc      81.571  81.571     1    28  1.0968 0.3039
## I(Inclusao_sc^2) 52.924  52.924     1    28  0.7116 0.4061
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_C
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod_lin     4 230.85 236.46 -111.42   222.85                     
## mod_quad    5 232.15 239.15 -111.07   222.15 0.7027  1     0.4019
## 
## 
## ###########################
## Variável: Butirico 
## ##############################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 29.348  29.348     1    28  1.1839 0.2858
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc      29.3477 29.3477     1    28  1.1854 0.2855
## I(Inclusao_sc^2)  0.8858  0.8858     1    28  0.0358 0.8513
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_C
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod_lin     4 197.96 203.56 -94.980   189.96                     
## mod_quad    5 199.92 206.93 -94.962   189.92 0.0358  1       0.85
## 
## 
## ###########################
## Variável: CH4.24h 
## ##############################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 2.3857  2.3857     1    28  0.5707 0.4563
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc       2.3857  2.3857     1    28  0.6815 0.41604  
## I(Inclusao_sc^2) 19.0289 19.0289     1    28  5.4360 0.02715 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_C
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
## mod_lin     4 139.15 144.76 -65.577   131.15                      
## mod_quad    5 136.19 143.19 -63.093   126.19 4.968  1    0.02582 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ###########################
## Variável: AP.ratio 
## ##############################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 0.012339 0.012339     1    28  0.5355 0.4704
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc      0.012339 0.012339     1    28  0.6571 0.42442  
## I(Inclusao_sc^2) 0.119383 0.119383     1    28  6.3576 0.01766 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_C
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar     AIC      BIC logLik deviance  Chisq Df Pr(>Chisq)  
## mod_lin     4 -13.843  -8.2381 10.921  -21.843                       
## mod_quad    5 -17.572 -10.5662 13.786  -27.572 5.7293  1    0.01668 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2 - PLOTS - REGRESSÃO QUADRÁTICA

variables = c("PG_Total","DMO","CH4.24h","AP.ratio")

for (var in variables) {
  # Modelo quadrático
  formula_quad = as.formula(paste(var, "~ Inclusao_sc + I(Inclusao_sc^2) + (1|inoculo)"))
  mod_quad = lmer(formula_quad,data = data_C,REML = FALSE)
  # Dados para construir a curva
  pred = data.frame(Inclusao = seq(
      min(data_C$Inclusao, na.rm = TRUE),
      max(data_C$Inclusao, na.rm = TRUE),length.out = 200))
  # Mesma transformação utilizada no modelo
  pred$Inclusao_sc = (pred$Inclusao - 50) / 25
  # Predição apenas dos efeitos fixos
  pred$y_pred = predict(mod_quad,newdata = pred,re.form = NA)
  
  # Plot
  p = ggplot(data_C,aes(x = Inclusao,y = .data[[var]])) +
    geom_point(size = 2.5,alpha = 0.7) +
    geom_line(data = pred,aes(x = Inclusao,y = y_pred),linewidth = 1) +
    scale_x_continuous(breaks = c(0, 25, 50, 75, 100)) +
    labs(x = "Inclusão (%)",y = var) +
    theme_classic(base_size = 14)
  print(p)}

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

2.3 - Ranking

cada variavel teve pesos iguais, as médias foram ajustadas via modelo simples Direções: - DMO: quanto maior melhor; - PG_Total, CH4, AP.ratio: quanto menor melhor (invertido);

EXPLICAÇÃO PARA TEXTOS CIENTIFICOS: Foi construído um índice multicritério para ranquear as inclusões de leguminosas com base na degradabilidade da matéria orgânica (DMO), na produção de gas (PG_Total), produção de metano (CH4.24h) e razão acetato:propionato (AP.ratio). As médias ajustadas foram obtidas por modelo linear considerando inclusão como efeito fixo, seguidas de estimativas de médias marginais (emmeans).Para permitir a comparação entre variáveis em escalas distintas, os valores foram padronizados por transformação min–max (0–1). Para DMO, valores maiores foram considerados mais favoráveis. Para PG_Total, CH4.24h e AP.ratio, a escala foi invertida, de modo que menores produções resultassem em maiores escores padronizados.O índice final foi calculado como a média aritmética simples dos escores padronizados (pesos iguais), e os tratamentos foram ranqueados em ordem decrescente do escore composto.

data_C_ranking = data %>% dplyr::filter(Planta %in% c("Calopogonio_Marandu"))
View(data_C_ranking)

########################### FUNÇÃO MODELO ########################### 
fun_emm = function(data_C_ranking, response){
  f <- as.formula(paste0(response, " ~ Inclusao + factor(inoculo)"))
  mod <- lm(f, data = data_C_ranking)
  emmeans(mod,~ Inclusao,at = list(Inclusao = sort(unique(data_C_ranking$Inclusao)))) %>%
    as.data.frame() %>%
    mutate(variavel = response)}

########################### MÉDIAS AJUSTADAS (EMMEANS) ########################### 
variaveis_modelo = c("PG_Total","DMO","CH4.24h","AP.ratio")
medias1 = bind_rows(lapply(variaveis_modelo, \(v) fun_emm(data_C_ranking, v)))
View(medias1)

########################### ORGANIZAR EM FORMATO LARGO ###########################
emm_wide = medias1 %>% select(Inclusao, variavel, emmean, SE) %>%
  pivot_wider(names_from  = variavel,values_from = c(emmean, SE))
emm_wide
########################### PADRONIZAÇÃO EM ESCALA ###########################
ranking = emm_wide %>% ungroup() %>% mutate(
    # DMO (benefício: maior é melhor)
    DMO_score = (emmean_DMO - min(emmean_DMO, na.rm = TRUE)) /
                (max(emmean_DMO, na.rm = TRUE) - min(emmean_DMO, na.rm = TRUE)),
    # PG_Total (invertido: menor é melhor)
    PG_Total_score = (max(emmean_PG_Total, na.rm = TRUE) - emmean_PG_Total) /
                     (max(emmean_PG_Total, na.rm = TRUE) - min(emmean_PG_Total, na.rm = TRUE)),
    # CH4.24h (invertido: menor é melhor)
    CH4.24h_score = (max(emmean_CH4.24h, na.rm = TRUE) - emmean_CH4.24h) /
                    (max(emmean_CH4.24h, na.rm = TRUE) - min(emmean_CH4.24h, na.rm = TRUE)),
    # AP.ratio (invertido: menor é melhor)
    AP.ratio_score = (max(emmean_AP.ratio, na.rm = TRUE) - emmean_AP.ratio) /
                     (max(emmean_AP.ratio, na.rm = TRUE) - min(emmean_AP.ratio, na.rm = TRUE)))
ranking
########################### ÍNDICE MULTICRITÉRIO (PESOS IGUAIS) ###########################
ranking = ranking %>%
  mutate(Score = (DMO_score + PG_Total_score + CH4.24h_score + AP.ratio_score)/4) %>%
  arrange(desc(Score)) %>% mutate(Rank = row_number())

ranking_final = ranking %>% select(Inclusao,
         emmean_DMO, DMO_score,
         emmean_PG_Total, PG_Total_score,
         emmean_CH4.24h, CH4.24h_score,
         emmean_AP.ratio, AP.ratio_score,
         Score, Rank)
ranking_final

CONCLUSÃO

Para Calopognio, não foi possivel encontrar uma diferença de score, observamos diferenças no score individual de cada variavel, mas o score médio (as 4 variaveis dividido por 4) ficou 0.5 para os tres niveis. Então a melhor inclusão terá que ser escolhida por um fator de maior peso, eu escolheria 25% de inclusão, por ser o nivel com maior DMO e menor A:P ratio.


###### 3 - Dados Sojá Perene ######


3.1 - ANOVA + REGRESSÃO

# Filtrar apenas CM, M e M_CM
data_SP = data %>% dplyr::filter(Planta %in% c("Marandu", "Soja_perene", "Soja_Marandu"))
table(data_SP$Planta)
## 
##         Calopogonio Calopogonio_Marandu         Macrotiloma Macrotiloma_Marandu 
##                   0                   0                   0                   0 
##             Marandu        Soja_Marandu         Soja_perene         Tifton_LANA 
##                   6                  18                   6                   0
View(data_SP)

###################################
data_SP$Inclusao_sc = (data_SP$Inclusao - 50) / 25
variables = c("PG_Total","DMO","Acetico","Propionico","Butirico","CH4.24h","AP.ratio")

for (var in variables) {
  cat("\n\n###########################\n")
  cat("Variável:", var, "\n")
  cat("#########################\n")
  
  # Modelo linear
  formula_lin = as.formula(paste(var, "~ Inclusao_sc + (1|inoculo)"))
  mod_lin = lmer(formula_lin, data = data_SP, REML = FALSE)
  
  # Modelo quadrático
  formula_quad = as.formula(paste(var, "~ Inclusao_sc + I(Inclusao_sc^2) + (1|inoculo)"))
  mod_quad = lmer(formula_quad,data = data_SP,REML = FALSE)
  
  cat("\n---------------- MODELO LINEAR ----------------\n")
  print(anova(mod_lin))
  cat("\n---------------- MODELO QUADRÁTICO ----------------\n")
  print(anova(mod_quad))
  cat("\n---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------\n")
  print(anova(mod_lin, mod_quad))}
## 
## 
## ###########################
## Variável: PG_Total 
## #########################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc 969.48  969.48     1    28  5.1001 0.03191 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Inclusao_sc       969.48  969.48     1    28  7.6451 0.0099592 ** 
## I(Inclusao_sc^2) 1771.86 1771.86     1    28 13.9724 0.0008445 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_SP
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod_lin     4 257.36 262.96 -124.68   249.36                         
## mod_quad    5 248.02 255.03 -119.01   238.02 11.335  1  0.0007608 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ###########################
## Variável: DMO 
## #########################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Inclusao_sc  80111   80111     1    20  31.754 1.626e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Inclusao_sc       67732   67732     1    20  49.764 7.686e-07 ***
## I(Inclusao_sc^2)  23235   23235     1    20  17.071  0.000517 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_SP
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod_lin     4 247.62 251.98 -119.81   239.62                         
## mod_quad    5 237.27 242.73 -113.64   227.27 12.342  1  0.0004428 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ###########################
## Variável: Acetico 
## #########################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 115.66  115.66     1    28  0.3935 0.5355
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Inclusao_sc       115.66  115.66     1    28  0.5307 0.472374   
## I(Inclusao_sc^2) 2126.42 2126.42     1    28  9.7566 0.004129 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_SP
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## mod_lin     4 267.35 272.96 -129.68   259.35                        
## mod_quad    5 260.98 267.99 -125.49   250.98 8.3707  1   0.003813 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ###########################
## Variável: Propionico 
## #########################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 0.032648 0.032648     1    28   3e-04 0.9852
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Inclusao_sc         0.03    0.03     1    28  0.0006 0.9809934    
## I(Inclusao_sc^2) 1034.31 1034.31     1    28 18.3050 0.0001984 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_SP
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod_lin     4 235.44 241.04 -113.72   227.44                         
## mod_quad    5 223.35 230.36 -106.67   213.35 14.085  1  0.0001747 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ###########################
## Variável: Butirico 
## #########################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc 99.311  99.311     1    28  2.9097 0.09912 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Inclusao_sc       99.311  99.311     1    28  3.3651 0.07724 .
## I(Inclusao_sc^2) 129.324 129.324     1    28  4.3821 0.04550 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_SP
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod_lin     4 205.82 211.42 -98.908   197.82                       
## mod_quad    5 203.75 210.75 -96.873   193.75 4.0712  1    0.04362 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## ###########################
## Variável: CH4.24h 
## #########################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 0.000663 0.000663     1    28   1e-04 0.9906
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc      0.0007  0.0007     1    28  0.0002 0.9902
## I(Inclusao_sc^2) 8.8480  8.8480     1    28  2.0379 0.1645
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_SP
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod_lin     4 142.43 148.04 -67.216   134.43                     
## mod_quad    5 142.47 149.47 -66.233   132.47 1.9671  1     0.1608
## 
## 
## ###########################
## Variável: AP.ratio 
## #########################
## 
## ---------------- MODELO LINEAR ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##              Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc 0.03929 0.03929     1    28  0.4863 0.4913
## 
## ---------------- MODELO QUADRÁTICO ----------------
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Inclusao_sc      0.039290 0.039290     1    28  0.4951 0.4875
## I(Inclusao_sc^2) 0.040199 0.040199     1    28  0.5065 0.4825
## 
## ---------------- COMPARAÇÃO LINEAR vs QUADRÁTICO ----------------
## Data: data_SP
## Models:
## mod_lin: formula_lin
## mod_quad: formula_quad
##          npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod_lin     4 22.204 27.808 -7.1018   14.204                    
## mod_quad    5 23.702 30.708 -6.8508   13.702 0.502  1     0.4786

3.2 - PLOTS - REGRESSÃO QUADRÁTICA

variables = c("PG_Total","DMO","CH4.24h","AP.ratio")

for (var in variables) {
  # Modelo quadrático
  formula_quad = as.formula(paste(var, "~ Inclusao_sc + I(Inclusao_sc^2) + (1|inoculo)"))
  mod_quad = lmer(formula_quad,data = data_SP,REML = FALSE)
  # Dados para construir a curva
  pred = data.frame(Inclusao = seq(
      min(data_SP$Inclusao, na.rm = TRUE),
      max(data_SP$Inclusao, na.rm = TRUE),length.out = 200))
  # Mesma transformação utilizada no modelo
  pred$Inclusao_sc = (pred$Inclusao - 50) / 25
  # Predição apenas dos efeitos fixos
  pred$y_pred = predict(mod_quad,newdata = pred,re.form = NA)
  
  # Plot
  p = ggplot(data_SP,aes(x = Inclusao,y = .data[[var]])) +
    geom_point(size = 2.5,alpha = 0.7) +
    geom_line(data = pred,aes(x = Inclusao,y = y_pred),linewidth = 1) +
    scale_x_continuous(breaks = c(0, 25, 50, 75, 100)) +
    labs(x = "Inclusão (%)",y = var) +
    theme_classic(base_size = 14)
  print(p)}

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

2.3 - Ranking

cada variavel teve pesos iguais, as médias foram ajustadas via modelo simples Direções: - DMO: quanto maior melhor; - PG_Total, CH4, AP.ratio: quanto menor melhor (invertido);

EXPLICAÇÃO PARA TEXTOS CIENTIFICOS: Foi construído um índice multicritério para ranquear as inclusões de leguminosas com base na degradabilidade da matéria orgânica (DMO), na produção de gas (PG_Total), produção de metano (CH4.24h) e razão acetato:propionato (AP.ratio). As médias ajustadas foram obtidas por modelo linear considerando inclusão como efeito fixo, seguidas de estimativas de médias marginais (emmeans).Para permitir a comparação entre variáveis em escalas distintas, os valores foram padronizados por transformação min–max (0–1). Para DMO, valores maiores foram considerados mais favoráveis. Para PG_Total, CH4.24h e AP.ratio, a escala foi invertida, de modo que menores produções resultassem em maiores escores padronizados.O índice final foi calculado como a média aritmética simples dos escores padronizados (pesos iguais), e os tratamentos foram ranqueados em ordem decrescente do escore composto.

data_SP_ranking = data %>% dplyr::filter(Planta %in% c("Soja_Marandu"))
View(data_SP_ranking)

########################### FUNÇÃO MODELO ########################### 
fun_emm = function(data_SP_ranking, response){
  f <- as.formula(paste0(response, " ~ Inclusao + factor(inoculo)"))
  mod <- lm(f, data = data_SP_ranking)
  emmeans(mod,~ Inclusao,at = list(Inclusao = sort(unique(data_SP_ranking$Inclusao)))) %>%
    as.data.frame() %>%
    mutate(variavel = response)}

########################### MÉDIAS AJUSTADAS (EMMEANS) ########################### 
variaveis_modelo = c("PG_Total","DMO","CH4.24h","AP.ratio")
medias1 = bind_rows(lapply(variaveis_modelo, \(v) fun_emm(data_SP_ranking, v)))
View(medias1)

########################### ORGANIZAR EM FORMATO LARGO ###########################
emm_wide = medias1 %>% select(Inclusao, variavel, emmean, SE) %>%
  pivot_wider(names_from  = variavel,values_from = c(emmean, SE))
emm_wide
########################### PADRONIZAÇÃO EM ESCALA ###########################
ranking = emm_wide %>% ungroup() %>% mutate(
    # DMO (benefício: maior é melhor)
    DMO_score = (emmean_DMO - min(emmean_DMO, na.rm = TRUE)) /
                (max(emmean_DMO, na.rm = TRUE) - min(emmean_DMO, na.rm = TRUE)),
    # PG_Total (invertido: menor é melhor)
    PG_Total_score = (max(emmean_PG_Total, na.rm = TRUE) - emmean_PG_Total) /
                     (max(emmean_PG_Total, na.rm = TRUE) - min(emmean_PG_Total, na.rm = TRUE)),
    # CH4.24h (invertido: menor é melhor)
    CH4.24h_score = (max(emmean_CH4.24h, na.rm = TRUE) - emmean_CH4.24h) /
                    (max(emmean_CH4.24h, na.rm = TRUE) - min(emmean_CH4.24h, na.rm = TRUE)),
    # AP.ratio (invertido: menor é melhor)
    AP.ratio_score = (max(emmean_AP.ratio, na.rm = TRUE) - emmean_AP.ratio) /
                     (max(emmean_AP.ratio, na.rm = TRUE) - min(emmean_AP.ratio, na.rm = TRUE)))
ranking
########################### ÍNDICE MULTICRITÉRIO (PESOS IGUAIS) ###########################
ranking = ranking %>%
  mutate(Score = (DMO_score + PG_Total_score + CH4.24h_score + AP.ratio_score)/4) %>%
  arrange(desc(Score)) %>% mutate(Rank = row_number())

ranking_final = ranking %>% select(Inclusao,
         emmean_DMO, DMO_score,
         emmean_PG_Total, PG_Total_score,
         emmean_CH4.24h, CH4.24h_score,
         emmean_AP.ratio, AP.ratio_score,
         Score, Rank)
ranking_final

CONCLUSÃO

Para sojá perene, o melhor nivel foi 75%, com maior DMO e menor CH4 e A:P ratio.