Dados DOC - Laura

Packages

library(multcomp)

Warning: package ‘multcomp’ was built under R version 4.3.3Carregando pacotes exigidos: mvtnorm
Carregando pacotes exigidos: survival
Carregando pacotes exigidos: TH.data
Warning: package ‘TH.data’ was built under R version 4.3.3Carregando pacotes exigidos: MASS

Attaching package: ‘TH.data’

The following object is masked from ‘package:MASS’:

    geyser

DATA

1 - DMO

data1.1 = subset(data1, DIETA != "D0") #retirando a dieta controle do data
data1.1

#model
mod1 = lmer(DMO~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod1))

shapiro.test(resid(mod1))

    Shapiro-Wilk normality test

data:  resid(mod1)
W = 0.97392, p-value = 0.285

anova(mod1)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
               Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
DIETA           53107 26553.7     2    43 18.0318 2.056e-06 ***
INCLUSAO        10676  5337.9     2    43  3.6248   0.03509 *  
DIETA:INCLUSAO   5916  1478.9     4    43  1.0043   0.41584    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod1.1 = lmer(DMO~DIETA+(1|INOC), data = data1)
hist(resid(mod1.1))

shapiro.test(resid(mod1.1))

    Shapiro-Wilk normality test

data:  resid(mod1.1)
W = 0.98228, p-value = 0.5322

anova(mod1.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
DIETA  71525   23842     3    54  14.675 4.149e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias1.1=emmeans(mod1.1,~ DIETA)
summary(medias1.1) #media do fator dieta

 DIETA emmean    SE   df lower.CL upper.CL
 D0       464 16.58 47.6      431      498
 D1       507  9.72 15.8      487      528
 D2       494  9.72 15.8      474      515
 D3       566  9.72 15.8      546      587

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

tukey1 = pairs(medias1.1, adjust = "tukey")
print(tukey1)

 contrast estimate   SE df t.ratio p.value
 D0 - D1     -43.1 19.0 54  -2.269  0.1181
 D0 - D2     -30.0 19.0 54  -1.578  0.3994
 D0 - D3    -102.1 19.0 54  -5.373  <.0001
 D1 - D2      13.1 13.4 54   0.977  0.7630
 D1 - D3     -59.0 13.4 54  -4.390  0.0003
 D2 - D3     -72.1 13.4 54  -5.367  <.0001

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

tukey_df = as.data.frame(summary(tukey1))
tukey_matrix = matrix(NA, nrow = length(unique(data1$DIETA)), ncol = length(unique(data1$DIETA)))
rownames(tukey_matrix) = levels(data1$DIETA)
colnames(tukey_matrix) = levels(data1$DIETA)

for (i in 1:nrow(tukey_df)) {
  comp = unlist(strsplit(as.character(tukey_df$contrast[i]), " - "))
  tukey_matrix[comp[1], comp[2]] = tukey_df$p.value[i]
  tukey_matrix[comp[2], comp[1]] = tukey_df$p.value[i]
}

letters = multcompLetters(tukey_matrix, reversed = TRUE)$Letters
print(letters) #letras

 D0  D1  D2  D3 
"b" "b" "b" "a" 

#########desdobramento - FATOR INCLUSÃO#########

mod1.2 = lmer(DMO~INCLUSAO+(1|INOC), data = data1)

boundary (singular) fit: see help('isSingular')

hist(resid(mod1.2))

shapiro.test(resid(mod1.2))

    Shapiro-Wilk normality test

data:  resid(mod1.2)
W = 0.98138, p-value = 0.4896

anova(mod1.2) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
INCLUSAO  29093  9697.6     3    56  4.0553 0.01117 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias1.2=emmeans(mod1.2,~ INCLUSAO)
summary(medias1.2) #media do fator inclusão

 INCLUSAO emmean   SE   df lower.CL upper.CL
 0           464 20.0 50.0      424      504
 25          523 11.5 18.5      499      548
 50          540 11.5 18.5      516      564
 100         505 11.5 18.5      481      529

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg1 = lm(DMO~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg1)

Call:
lm(formula = DMO ~ poly(as.numeric(INCLUSAO), degree = 3), data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-88.023 -28.927  -7.087  35.455 104.956 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              516.933      6.313  81.882   <2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1   40.277     48.901   0.824   0.4136    
poly(as.numeric(INCLUSAO), degree = 3)2 -165.577     48.901  -3.386   0.0013 ** 
poly(as.numeric(INCLUSAO), degree = 3)3   -7.421     48.901  -0.152   0.8799    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 48.9 on 56 degrees of freedom
Multiple R-squared:  0.1785,    Adjusted R-squared:  0.1345 
F-statistic: 4.055 on 3 and 56 DF,  p-value: 0.01117

2 - NETGPMO

#model
mod2 = lmer(NETGPMO~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod2))

shapiro.test(resid(mod2))

    Shapiro-Wilk normality test

data:  resid(mod2)
W = 0.97933, p-value = 0.4722

anova(mod2)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
               Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)   
DIETA          657.03  328.51     2    43  7.5055 0.00160 **
INCLUSAO       396.96  198.48     2    43  4.5347 0.01633 * 
DIETA:INCLUSAO  52.91   13.23     4    43  0.3022 0.87489   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod2.1 = lmer(NETGPMO~DIETA+(1|INOC), data = data1)
hist(resid(mod2.1))

shapiro.test(resid(mod2.1))

    Shapiro-Wilk normality test

data:  resid(mod2.1)
W = 0.98965, p-value = 0.8925

anova(mod2.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
DIETA 883.55  294.52     3    54  5.9751 0.001358 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias2.1=emmeans(mod2.1,~ DIETA)
summary(medias2.1) #media do fator dieta

 DIETA emmean   SE   df lower.CL upper.CL
 D0      44.5 4.86 4.21     31.3     57.7
 D1      48.2 4.26 2.50     33.0     63.4
 D2      48.9 4.26 2.50     33.6     64.1
 D3      55.9 4.26 2.50     40.7     71.1

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

tukey2 = pairs(medias2.1, adjust = "tukey") #mudar aqui
print(tukey2) #mudar aqui

 contrast estimate   SE df t.ratio p.value
 D0 - D1    -3.700 3.31 54  -1.118  0.6801
 D0 - D2    -4.334 3.31 54  -1.309  0.5610
 D0 - D3   -11.396 3.31 54  -3.443  0.0060
 D1 - D2    -0.633 2.34 54  -0.271  0.9930
 D1 - D3    -7.696 2.34 54  -3.288  0.0093
 D2 - D3    -7.062 2.34 54  -3.018  0.0197

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

tukey_df = as.data.frame(summary(tukey2)) #mudar aqui
tukey_matrix = matrix(NA, nrow = length(unique(data1$DIETA)), ncol = length(unique(data1$DIETA))) #não mudar
rownames(tukey_matrix) = levels(data1$DIETA) #não mudar
colnames(tukey_matrix) = levels(data1$DIETA) #não mudar

for (i in 1:nrow(tukey_df)) { 
  comp = unlist(strsplit(as.character(tukey_df$contrast[i]), " - "))
  tukey_matrix[comp[1], comp[2]] = tukey_df$p.value[i]
  tukey_matrix[comp[2], comp[1]] = tukey_df$p.value[i]
} #não mudar

letters = multcompLetters(tukey_matrix, reversed = TRUE)$Letters #não mudar
print(letters) #letras

 D0  D1  D2  D3 
"b" "b" "b" "a" 

#########desdobramento - FATOR INCLUSÃO#########

mod2.2 = lmer(NETGPMO~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod2.2))

shapiro.test(resid(mod2.2))

    Shapiro-Wilk normality test

data:  resid(mod2.2)
W = 0.9777, p-value = 0.339

anova(mod2.2) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
INCLUSAO 623.49  207.83     3    54  3.8411 0.0145 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias2.2=emmeans(mod2.2,~ INCLUSAO)
summary(medias2.2) #media do fator inclusão

 INCLUSAO emmean   SE   df lower.CL upper.CL
 0          44.5 4.93 4.46     31.4     57.7
 25         53.2 4.28 2.55     38.1     68.2
 50         52.7 4.28 2.55     37.6     67.7
 100        47.2 4.28 2.55     32.1     62.3

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg2 = lm(NETGPMO~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg2)

Call:
lm(formula = NETGPMO ~ poly(as.numeric(INCLUSAO), degree = 3), 
    data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-22.155  -6.957  -2.375   8.802  20.242 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                               50.346      1.204  41.832   <2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1   -4.988      9.323  -0.535   0.5947    
poly(as.numeric(INCLUSAO), degree = 3)2  -24.172      9.323  -2.593   0.0121 *  
poly(as.numeric(INCLUSAO), degree = 3)3    3.783      9.323   0.406   0.6865    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.323 on 56 degrees of freedom
Multiple R-squared:  0.1136,    Adjusted R-squared:  0.06607 
F-statistic: 2.391 on 3 and 56 DF,  p-value: 0.07822

3 - NETCH4MO

#model
mod3 = lmer(NETCH4MO~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod3))

shapiro.test(resid(mod3))

    Shapiro-Wilk normality test

data:  resid(mod3)
W = 0.97609, p-value = 0.3512

anova(mod3)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
DIETA          21.7986 10.8993     2    43 12.2274 6.249e-05 ***
INCLUSAO        3.0649  1.5324     2    43  1.7192    0.1913    
DIETA:INCLUSAO  4.9724  1.2431     4    43  1.3946    0.2519    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod3.1 = lmer(NETCH4MO~DIETA+(1|INOC), data = data1)
hist(resid(mod3.1))

shapiro.test(resid(mod3.1)) #essa merda não está normal vou ter que transformar

    Shapiro-Wilk normality test

data:  resid(mod3.1)
W = 0.94452, p-value = 0.008677

mod3.1.1 = lmer(NETCH4MO^0.9~DIETA+(1|INOC), data = data1) #transformação
hist(resid(mod3.1.1))

shapiro.test(resid(mod3.1.1))

    Shapiro-Wilk normality test

data:  resid(mod3.1.1)
W = 0.99171, p-value = 0.9594

anova(mod3.1.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
DIETA 16.676  5.5585     3 53.001   13.25 1.432e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias3.1=emmeans(mod3.1,~ DIETA)
summary(medias3.1) #media do fator dieta, usar o modelo não transformado

 DIETA emmean    SE   df lower.CL upper.CL
 D0      2.54 0.829 3.13  -0.0373     5.11
 D1      3.30 0.764 2.27   0.3647     6.24
 D2      3.24 0.764 2.27   0.2971     6.17
 D3      4.62 0.764 2.27   1.6774     7.55

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

tukey3 = pairs(medias3.1, adjust = "tukey") #mudar aqui
print(tukey3) #mudar aqui

 contrast estimate    SE df t.ratio p.value
 D0 - D1   -0.7665 0.454 54  -1.690  0.3388
 D0 - D2   -0.6989 0.454 54  -1.541  0.4206
 D0 - D3   -2.0792 0.454 54  -4.585  0.0002
 D1 - D2    0.0676 0.321 54   0.211  0.9966
 D1 - D3   -1.3127 0.321 54  -4.094  0.0008
 D2 - D3   -1.3803 0.321 54  -4.304  0.0004

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

tukey_df = as.data.frame(summary(tukey3)) #mudar aqui
tukey_matrix = matrix(NA, nrow = length(unique(data1$DIETA)), ncol = length(unique(data1$DIETA))) #não mudar
rownames(tukey_matrix) = levels(data1$DIETA) #não mudar
colnames(tukey_matrix) = levels(data1$DIETA) #não mudar

for (i in 1:nrow(tukey_df)) { 
  comp = unlist(strsplit(as.character(tukey_df$contrast[i]), " - "))
  tukey_matrix[comp[1], comp[2]] = tukey_df$p.value[i]
  tukey_matrix[comp[2], comp[1]] = tukey_df$p.value[i]
} #não mudar

letters = multcompLetters(tukey_matrix, reversed = TRUE)$Letters #não mudar
print(letters) #letras

 D0  D1  D2  D3 
"b" "b" "b" "a" 

#########desdobramento - FATOR INCLUSÃO#########

mod3.2 = lmer(NETCH4MO~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod3.2))

shapiro.test(resid(mod3.2))

    Shapiro-Wilk normality test

data:  resid(mod3.2)
W = 0.97077, p-value = 0.1592

anova(mod3.2) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
INCLUSAO 10.604  3.5346     3    54  2.7777 0.04989 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias3.2=emmeans(mod3.2,~ INCLUSAO)
summary(medias3.2) #media do fator inclusão

 INCLUSAO emmean    SE   df lower.CL upper.CL
 0          2.54 0.859 3.62   0.0485     5.02
 25         3.98 0.773 2.38   1.1150     6.85
 50         3.77 0.773 2.38   0.8995     6.63
 100        3.41 0.773 2.38   0.5376     6.27

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg3 = lm(NETCH4MO~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg3)

Call:
lm(formula = NETCH4MO ~ poly(as.numeric(INCLUSAO), degree = 3), 
    data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8142 -0.9822 -0.0952  0.7424  3.1865 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                                3.600      0.200  17.997   <2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1    0.312      1.549   0.201   0.8412    
poly(as.numeric(INCLUSAO), degree = 3)2   -2.938      1.549  -1.896   0.0631 .  
poly(as.numeric(INCLUSAO), degree = 3)3    1.370      1.549   0.884   0.3803    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.549 on 56 degrees of freedom
Multiple R-squared:  0.0731,    Adjusted R-squared:  0.02345 
F-statistic: 1.472 on 3 and 56 DF,  p-value: 0.2319

4 - PF

#model
mod4 = lmer(PF~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod4))

shapiro.test(resid(mod4))

    Shapiro-Wilk normality test

data:  resid(mod4)
W = 0.97604, p-value = 0.3496

anova(mod4)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
DIETA          1.70880 0.85440     2    43 21.8294 2.862e-07 ***
INCLUSAO       0.43731 0.21866     2    43  5.5865   0.00697 ** 
DIETA:INCLUSAO 0.30384 0.07596     4    43  1.9408   0.12093    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod4.1 = lmer(PF~DIETA+(1|INOC), data = data1)
hist(resid(mod4.1))

shapiro.test(resid(mod4.1)) 

    Shapiro-Wilk normality test

data:  resid(mod4.1)
W = 0.99153, p-value = 0.9528

anova(mod4.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
DIETA 2.2922 0.76407     3    54  16.624 8.956e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias4.1=emmeans(mod4.1,~ DIETA)
summary(medias4.1) #media do fator dieta

 DIETA emmean    SE   df lower.CL upper.CL
 D0      2.81 0.144 4.45     2.42     3.19
 D1      3.10 0.125 2.55     2.66     3.54
 D2      2.94 0.125 2.55     2.50     3.38
 D3      3.37 0.125 2.55     2.93     3.81

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

tukey4 = pairs(medias4.1, adjust = "tukey") #mudar aqui
print(tukey4) #mudar aqui

 contrast estimate     SE df t.ratio p.value
 D0 - D1    -0.292 0.1011 54  -2.886  0.0278
 D0 - D2    -0.132 0.1011 54  -1.303  0.5650
 D0 - D3    -0.563 0.1011 54  -5.568  <.0001
 D1 - D2     0.160 0.0715 54   2.239  0.1259
 D1 - D3    -0.271 0.0715 54  -3.792  0.0021
 D2 - D3    -0.431 0.0715 54  -6.031  <.0001

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

tukey_df = as.data.frame(summary(tukey4)) #mudar aqui
tukey_matrix = matrix(NA, nrow = length(unique(data1$DIETA)), ncol = length(unique(data1$DIETA))) #não mudar
rownames(tukey_matrix) = levels(data1$DIETA) #não mudar
colnames(tukey_matrix) = levels(data1$DIETA) #não mudar

for (i in 1:nrow(tukey_df)) { 
  comp = unlist(strsplit(as.character(tukey_df$contrast[i]), " - "))
  tukey_matrix[comp[1], comp[2]] = tukey_df$p.value[i]
  tukey_matrix[comp[2], comp[1]] = tukey_df$p.value[i]
} #não mudar

letters = multcompLetters(tukey_matrix, reversed = TRUE)$Letters #não mudar
print(letters) #letras

  D0   D1   D2   D3 
 "c"  "b" "bc"  "a" 

#########desdobramento - FATOR INCLUSÃO#########

mod4.2 = lmer(PF~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod4.2))

shapiro.test(resid(mod4.2))

    Shapiro-Wilk normality test

data:  resid(mod4.2)
W = 0.98376, p-value = 0.6058

anova(mod4.2) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
INCLUSAO 1.0207 0.34024     3    54   4.895 0.004412 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias4.2=emmeans(mod4.2,~ INCLUSAO)
summary(medias4.2) #media do fator inclusão

 INCLUSAO emmean    SE   df lower.CL upper.CL
 0          2.81 0.156 6.03     2.43     3.19
 25         3.02 0.129 2.85     2.60     3.44
 50         3.24 0.129 2.85     2.82     3.66
 100        3.15 0.129 2.85     2.73     3.57

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg4 = lm(PF~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg4)

Call:
lm(formula = PF ~ poly(as.numeric(INCLUSAO), degree = 3), data = data1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53915 -0.18647  0.01635  0.19400  0.73086 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              3.10367    0.04014  77.326   <2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1  0.77942    0.31090   2.507   0.0151 *  
poly(as.numeric(INCLUSAO), degree = 3)2 -0.57684    0.31090  -1.855   0.0688 .  
poly(as.numeric(INCLUSAO), degree = 3)3 -0.28367    0.31090  -0.912   0.3655    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3109 on 56 degrees of freedom
Multiple R-squared:  0.1586,    Adjusted R-squared:  0.1136 
F-statistic:  3.52 on 3 and 56 DF,  p-value: 0.02074

5 - PH

#model
mod5 = lmer(PH~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod5))

shapiro.test(resid(mod5))

    Shapiro-Wilk normality test

data:  resid(mod5)
W = 0.97121, p-value = 0.2178

anova(mod5)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
                 Sum Sq   Mean Sq NumDF DenDF F value    Pr(>F)    
DIETA          0.015648 0.0078241     2    43  4.7078 0.0141632 *  
INCLUSAO       0.037470 0.0187352     2    43 11.2732 0.0001158 ***
DIETA:INCLUSAO 0.001541 0.0003852     4    43  0.2318 0.9190062    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod5.1 = lmer(PH~DIETA+(1|INOC), data = data1)
hist(resid(mod5.1))

shapiro.test(resid(mod5.1)) 

    Shapiro-Wilk normality test

data:  resid(mod5.1)
W = 0.98882, p-value = 0.8584

anova(mod5.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
        Sum Sq   Mean Sq NumDF DenDF F value  Pr(>F)  
DIETA 0.021782 0.0072607     3    54  2.9337 0.04154 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias5.1=emmeans(mod5.1,~ DIETA)
summary(medias5.1) #media do fator dieta

 DIETA emmean      SE    df lower.CL upper.CL
 D0     6.670 0.02056 46.13    6.629    6.711
 D1     6.692 0.01215 14.62    6.666    6.718
 D2     6.692 0.01215 14.62    6.666    6.718
 D3     6.728 0.01215 14.62    6.702    6.754

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

tukey5 = pairs(medias5.1, adjust = "tukey") #mudar aqui
print(tukey5) #mudar aqui

 contrast estimate     SE df t.ratio p.value
 D0 - D1   -0.0217 0.0235 54  -0.924  0.7922
 D0 - D2   -0.0217 0.0235 54  -0.924  0.7922
 D0 - D3   -0.0578 0.0235 54  -2.464  0.0774
 D1 - D2    0.0000 0.0166 54   0.000  1.0000
 D1 - D3   -0.0361 0.0166 54  -2.178  0.1425
 D2 - D3   -0.0361 0.0166 54  -2.178  0.1425

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

tukey_df = as.data.frame(summary(tukey5)) #mudar aqui
tukey_matrix = matrix(NA, nrow = length(unique(data1$DIETA)), ncol = length(unique(data1$DIETA))) #não mudar
rownames(tukey_matrix) = levels(data1$DIETA) #não mudar
colnames(tukey_matrix) = levels(data1$DIETA) #não mudar

for (i in 1:nrow(tukey_df)) { 
  comp = unlist(strsplit(as.character(tukey_df$contrast[i]), " - "))
  tukey_matrix[comp[1], comp[2]] = tukey_df$p.value[i]
  tukey_matrix[comp[2], comp[1]] = tukey_df$p.value[i]
} #não mudar

letters = multcompLetters(tukey_matrix, reversed = TRUE)$Letters #não mudar
print(letters) #letras

 D0  D1  D2  D3 
"a" "a" "a" "a" 

#########desdobramento - FATOR INCLUSÃO#########

mod5.2 = lmer(PH~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod5.2))

shapiro.test(resid(mod5.2))

    Shapiro-Wilk normality test

data:  resid(mod5.2)
W = 0.96041, p-value = 0.04928

anova(mod5.2) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq  Mean Sq NumDF DenDF F value    Pr(>F)    
INCLUSAO 0.043604 0.014535     3    54  7.0188 0.0004529 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias5.2=emmeans(mod5.2,~ INCLUSAO)
summary(medias5.2) #media do fator inclusão

 INCLUSAO emmean      SE    df lower.CL upper.CL
 0         6.670 0.01903 42.19    6.632    6.708
 25        6.681 0.01148 12.09    6.656    6.706
 50        6.690 0.01148 12.09    6.665    6.715
 100       6.741 0.01148 12.09    6.716    6.766

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg5 = lm(PH~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg5)

Call:
lm(formula = PH ~ poly(as.numeric(INCLUSAO), degree = 3), data = data1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.100556 -0.022917  0.009444  0.020000  0.090000 

Coefficients:
                                        Estimate Std. Error  t value Pr(>|t|)    
(Intercept)                             6.700333   0.005926 1130.695  < 2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1 0.190264   0.045901    4.145 0.000116 ***
poly(as.numeric(INCLUSAO), degree = 3)2 0.077108   0.045901    1.680 0.098558 .  
poly(as.numeric(INCLUSAO), degree = 3)3 0.038191   0.045901    0.832 0.408926    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0459 on 56 degrees of freedom
Multiple R-squared:  0.2698,    Adjusted R-squared:  0.2307 
F-statistic: 6.899 on 3 and 56 DF,  p-value: 0.0004921

6 - AGVTotal

#model
mod6 = lmer(AGVTotal~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod6))

shapiro.test(resid(mod6))

    Shapiro-Wilk normality test

data:  resid(mod6)
W = 0.93516, p-value = 0.00589

mod6.0 = lmer(AGVTotal^3.3~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1)#Transformação

Warning: Model may not have converged with 1 eigenvalue close to zero: 1.2e-10

hist(resid(mod6.0))

shapiro.test(resid(mod6.0))

    Shapiro-Wilk normality test

data:  resid(mod6.0)
W = 0.95721, p-value = 0.05172

anova(mod6.0)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
                   Sum Sq    Mean Sq NumDF DenDF F value Pr(>F)
DIETA          3.4277e+12 1.7138e+12     2   Inf  1.1499 0.3167
INCLUSAO       5.2322e+12 2.6161e+12     2   Inf  1.7552 0.1729
DIETA:INCLUSAO 4.1721e+12 1.0430e+12     4   Inf  0.6998 0.5920

#########desdobramento - FATOR DIETAS#########

mod6.1 = lmer(AGVTotal~DIETA+(1|INOC), data = data1)
hist(resid(mod6.1))

shapiro.test(resid(mod6.1)) 

    Shapiro-Wilk normality test

data:  resid(mod6.1)
W = 0.95738, p-value = 0.03509

mod6.1.1 = lmer(AGVTotal^1.4~DIETA+(1|INOC), data = data1)
hist(resid(mod6.1.1))

shapiro.test(resid(mod6.1.1))

    Shapiro-Wilk normality test

data:  resid(mod6.1.1)
W = 0.96197, p-value = 0.05872

anova(mod6.1.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
DIETA 8951.9    2984     3    54   1.068 0.3704

medias6.1=emmeans(mod6.1,~ DIETA)
summary(medias6.1) #media do fator dieta

 DIETA emmean   SE   df lower.CL upper.CL
 D0       124 4.17 3.64      112      136
 D1       124 3.75 2.38      111      138
 D2       127 3.75 2.38      113      141
 D3       127 3.75 2.38      113      141

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

#########desdobramento - FATOR INCLUSÃO#########

mod6.2 = lmer(AGVTotal~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod6.2))

shapiro.test(resid(mod6.2))

    Shapiro-Wilk normality test

data:  resid(mod6.2)
W = 0.91353, p-value = 0.0004252

anova(mod6.2) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
INCLUSAO 126.52  42.174     3    54    1.42  0.247

medias6.2=emmeans(mod6.2,~ INCLUSAO)
summary(medias6.2) #media do fator inclusão

 INCLUSAO emmean   SE   df lower.CL upper.CL
 0           124 4.16 3.61      112      136
 25          127 3.75 2.37      113      140
 50          128 3.75 2.37      114      142
 100         124 3.75 2.37      110      138

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg6 = lm(AGVTotal~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg6)

Call:
lm(formula = AGVTotal ~ poly(as.numeric(INCLUSAO), degree = 3), 
    data = data1)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.5696  -4.7641  -0.1848   5.7795  17.4441 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                             125.9199     0.9681 130.072   <2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1  -2.3730     7.4987  -0.316    0.753    
poly(as.numeric(INCLUSAO), degree = 3)2 -10.6493     7.4987  -1.420    0.161    
poly(as.numeric(INCLUSAO), degree = 3)3  -2.7354     7.4987  -0.365    0.717    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.499 on 56 degrees of freedom
Multiple R-squared:  0.03863,   Adjusted R-squared:  -0.01287 
F-statistic:  0.75 on 3 and 56 DF,  p-value: 0.5269

7 - AcAcetico

#model
mod7 = lmer(AcAcetico~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod7))

shapiro.test(resid(mod7))

    Shapiro-Wilk normality test

data:  resid(mod7)
W = 0.98409, p-value = 0.6885

anova(mod7)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
               Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
DIETA           1.095  0.5477     2    43  0.1050 0.9005  
INCLUSAO       42.123 21.0615     2    43  4.0388 0.0247 *
DIETA:INCLUSAO 28.434  7.1084     4    43  1.3631 0.2626  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod7.1 = lmer(AcAcetico~DIETA+(1|INOC), data = data1)
hist(resid(mod7.1))

shapiro.test(resid(mod7.1)) 

    Shapiro-Wilk normality test

data:  resid(mod7.1)
W = 0.97058, p-value = 0.1557

anova(mod7.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
DIETA 12.549  4.1831     3    54  0.6277 0.6003

medias7.1=emmeans(mod7.1,~ DIETA)
summary(medias7.1) #media do fator dieta

 DIETA emmean   SE   df lower.CL upper.CL
 D0      70.6 1.93 3.73     65.1     76.1
 D1      72.1 1.72 2.40     65.7     78.4
 D2      71.9 1.72 2.40     65.6     78.3
 D3      72.3 1.72 2.40     65.9     78.6

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

#########desdobramento - FATOR INCLUSÃO#########

mod7.2 = lmer(AcAcetico~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod7.2))

shapiro.test(resid(mod7.2))

    Shapiro-Wilk normality test

data:  resid(mod7.2)
W = 0.97475, p-value = 0.2474

anova(mod7.2) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
INCLUSAO 53.577  17.859     3    54  3.0246 0.03734 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias7.2=emmeans(mod7.2,~ INCLUSAO)
summary(medias7.2) #media do fator inclusão

 INCLUSAO emmean   SE   df lower.CL upper.CL
 0          70.6 1.90 3.51     65.1     76.2
 25         70.9 1.71 2.35     64.5     77.3
 50         73.0 1.71 2.35     66.5     79.4
 100        72.4 1.71 2.35     66.0     78.9

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg7 = lm(AcAcetico~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg7)

Call:
lm(formula = AcAcetico ~ poly(as.numeric(INCLUSAO), degree = 3), 
    data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.4559 -2.2357  0.1723  2.1300  7.4831 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              71.9434     0.4379 164.308   <2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1   5.7850     3.3916   1.706   0.0936 .  
poly(as.numeric(INCLUSAO), degree = 3)2  -2.0316     3.3916  -0.599   0.5516    
poly(as.numeric(INCLUSAO), degree = 3)3  -3.9979     3.3916  -1.179   0.2435    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.392 on 56 degrees of freedom
Multiple R-squared:  0.07679,   Adjusted R-squared:  0.02733 
F-statistic: 1.553 on 3 and 56 DF,  p-value: 0.2111

8 - AcPropionico

#model
mod8 = lmer(AcPropionico~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod8))

shapiro.test(resid(mod8))

    Shapiro-Wilk normality test

data:  resid(mod8)
W = 0.95287, p-value = 0.03322

mod8.0 = lmer(AcPropionico^1.5~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1)
hist(resid(mod8.0))

shapiro.test(resid(mod8.0))

    Shapiro-Wilk normality test

data:  resid(mod8.0)
W = 0.95853, p-value = 0.05925

anova(mod8.0)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
               Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
DIETA           848.4  424.18     2    43  1.6993 0.1948624    
INCLUSAO       5264.2 2632.11     2    43 10.5443 0.0001878 ***
DIETA:INCLUSAO 1954.1  488.54     4    43  1.9571 0.1182769    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod8.1 = lmer(AcPropionico~DIETA+(1|INOC), data = data1)
hist(resid(mod8.1))

shapiro.test(resid(mod8.1)) 

    Shapiro-Wilk normality test

data:  resid(mod8.1)
W = 0.98189, p-value = 0.5134

anova(mod8.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
      Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
DIETA  23.46  7.8202     3    54  1.6261 0.1941

medias8.1=emmeans(mod8.1,~ DIETA)
summary(medias8.1) #media do fator dieta

 DIETA emmean    SE    df lower.CL upper.CL
 D0      30.0 1.127 10.07     27.5     32.5
 D1      31.6 0.858  3.58     29.1     34.1
 D2      32.0 0.858  3.58     29.5     34.5
 D3      30.9 0.858  3.58     28.4     33.4

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

#########desdobramento - FATOR INCLUSÃO#########

mod8.2 = lmer(AcPropionico~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod8.2))

shapiro.test(resid(mod8.2))

    Shapiro-Wilk normality test

data:  resid(mod8.2)
W = 0.9604, p-value = 0.04919

mod8.2.1 = lmer(AcPropionico^1.1~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod8.2.1))

shapiro.test(resid(mod8.2.1))

    Shapiro-Wilk normality test

data:  resid(mod8.2.1)
W = 0.96118, p-value = 0.05371

anova(mod8.2.1) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
INCLUSAO 208.53  69.512     3    54   7.921 0.0001808 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias8.2=emmeans(mod8.2,~ INCLUSAO)
summary(medias8.2) #media do fator inclusão

 INCLUSAO emmean    SE   df lower.CL upper.CL
 0          30.0 1.047 7.68     27.6     32.5
 25         32.1 0.831 3.16     29.5     34.7
 50         32.6 0.831 3.16     30.0     35.2
 100        29.9 0.831 3.16     27.3     32.5

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg8 = lm(AcPropionico~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg8)

Call:
lm(formula = AcPropionico ~ poly(as.numeric(INCLUSAO), degree = 3), 
    data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.2583 -1.2119 -0.2708  1.5801  4.7142 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              31.3897     0.2796 112.264  < 2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1  -3.0231     2.1658  -1.396 0.168273    
poly(as.numeric(INCLUSAO), degree = 3)2  -8.6484     2.1658  -3.993 0.000192 ***
poly(as.numeric(INCLUSAO), degree = 3)3  -1.5565     2.1658  -0.719 0.475343    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.166 on 56 degrees of freedom
Multiple R-squared:  0.2474,    Adjusted R-squared:  0.2071 
F-statistic: 6.137 on 3 and 56 DF,  p-value: 0.001105

9 - AcButirico

#model
mod9 = lmer(AcButirico~DIETA+INCLUSAO+DIETA*INCLUSAO+(1|INOC), data = data1.1) #modelo misto, dois fatores fixos e um aleatorio
hist(resid(mod9))

shapiro.test(resid(mod9))

    Shapiro-Wilk normality test

data:  resid(mod9)
W = 0.97356, p-value = 0.2974

anova(mod9)#P VALOR INTERAÇÃO

Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
DIETA           0.0331  0.0165     2 41.000  0.0196 0.980552   
INCLUSAO        2.3392  1.1696     2 41.001  1.3898 0.260619   
DIETA:INCLUSAO 18.8716  4.7179     4 41.001  5.6059 0.001077 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#########desdobramento - FATOR DIETAS#########

mod9.1 = lmer(AcButirico~DIETA+(1|INOC), data = data1)
hist(resid(mod9.1))

shapiro.test(resid(mod9.1)) 

    Shapiro-Wilk normality test

data:  resid(mod9.1)
W = 0.96413, p-value = 0.08409

anova(mod9.1) #p value do fator dieta

Type III Analysis of Variance Table with Satterthwaite's method
       Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
DIETA 0.48485 0.16162     3    52  0.1464 0.9315

medias9.1=emmeans(mod9.1,~ DIETA)
summary(medias9.1) #media do fator dieta

 DIETA emmean   SE   df lower.CL upper.CL
 D0      21.6 1.39 2.39     16.5     26.7
 D1      21.3 1.34 2.12     15.8     26.8
 D2      21.3 1.34 2.10     15.7     26.8
 D3      21.3 1.34 2.10     15.8     26.8

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

#########desdobramento - FATOR INCLUSÃO#########

mod9.2 = lmer(AcButirico~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod9.2))

shapiro.test(resid(mod9.2))

    Shapiro-Wilk normality test

data:  resid(mod9.2)
W = 0.9521, p-value = 0.02273

mod9.2.1 = lmer(AcButirico^0.4~INCLUSAO+(1|INOC), data = data1)
hist(resid(mod9.2.1))

shapiro.test(resid(mod9.2.1))

    Shapiro-Wilk normality test

data:  resid(mod9.2.1)
W = 0.96054, p-value = 0.05669

anova(mod9.2.1) #p value do fator inclusão

Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
INCLUSAO 0.012833 0.0042777     3 52.001  1.0217 0.3906

medias9.2=emmeans(mod9.2,~ INCLUSAO)
summary(medias9.2) #media do fator inclusão

 INCLUSAO emmean   SE   df lower.CL upper.CL
 0          21.6 1.39 2.37     16.4     26.7
 25         21.6 1.35 2.09     16.0     27.1
 50         21.3 1.35 2.10     15.8     26.8
 100        21.0 1.35 2.10     15.5     26.5

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

reg9 = lm(AcButirico~poly(as.numeric(INCLUSAO), degree = 3), data = data1)
summary(reg9)

Call:
lm(formula = AcButirico ~ poly(as.numeric(INCLUSAO), degree = 3), 
    data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.6024 -2.0162  0.1866  1.7187  4.8000 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              21.3234     0.2861  74.527   <2e-16 ***
poly(as.numeric(INCLUSAO), degree = 3)1  -1.3607     2.2076  -0.616    0.540    
poly(as.numeric(INCLUSAO), degree = 3)2   0.1137     2.2089   0.051    0.959    
poly(as.numeric(INCLUSAO), degree = 3)3   0.6362     2.2071   0.288    0.774    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.178 on 54 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.008579,  Adjusted R-squared:  -0.0465 
F-statistic: 0.1558 on 3 and 54 DF,  p-value: 0.9255