1 Packages

library(lme4)
library(lmerTest)
library(car)
library(emmeans)
library(multcompView)
library(multcomp)


2 DATA

Data=read.csv("C:/Users/Samsung/OneDrive/Documentos/1 - PósDoc/2 - FAPESP - Silvipastoril/Responsabilidades/Projeto Laura/invivo.csv")
str(Data)
'data.frame':   16 obs. of  8 variables:
 $ Treatment    : chr  "D0" "D0" "D0" "D0" ...
 $ Period       : int  1 2 3 4 1 2 3 4 1 2 ...
 $ Animal       : chr  "1903" "1717" "1417" "IZ1528" ...
 $ CH4_L_d      : num  12.9 13.1 13.3 20.1 26.6 ...
 $ CH4mw_L_Kg   : num  0.541 0.525 0.537 0.746 1.148 ...
 $ DMI          : num  1.36 1.36 1.05 1.11 1.28 ...
 $ Digestibility: num  72.9 71.1 78.1 72.9 55.1 ...
 $ CH4_DMI      : num  9.45 9.58 12.71 18.04 20.72 ...
Data$Treatment=as.factor(Data$Treatment)
Data$Period=as.factor(Data$Period)
Data$Animal=as.factor(Data$Animal)
print(Data)


3 ANOVA + TUKEY

names(Data)
[1] "Treatment"     "Period"        "Animal"        "CH4_L_d"       "CH4mw_L_Kg"    "DMI"          
[7] "Digestibility" "CH4_DMI"      
# Variáveis a serem analisadas
variables <- c("CH4_L_d", "CH4mw_L_Kg", "DMI", "Digestibility", "CH4_DMI")

# Testes de normalidade e homocedasticidade (com lm para simplificação)
testes <- lapply(variables, function(var) {
  formula <- as.formula(paste(var, "~ Treatment + (1|Animal) + (1|Period)"))
  modelo <- lmer(formula, data = Data)  
  res <- residuals(modelo)
  shapiro_p <- shapiro.test(res)$p.value
  levene_p <- leveneTest(res ~ Data$Treatment)[1, "Pr(>F)"]
  
  data.frame(Variavel = var,
             p_Shapiro = round(shapiro_p, 4),
             p_Levene = round(levene_p, 4))
})
boundary (singular) fit: see help('isSingular')
resultado_testes <- do.call(rbind, testes)
print(resultado_testes)

# Análise do modelo misto (quadro latino) e comparação de médias
for (var in variables) {
  cat("\n\n########################### Variável:", var, " ##########################\n")
  formula_mista <- as.formula(paste(var, "~ Treatment + (1|Animal) + (1|Period)"))
  mod_misto <- lmer(formula_mista, data = Data)
  # ANOVA
  cat("ANOVA (modelo misto):\n")
  print(anova(mod_misto))
  # Médias
  cat("\nMédias (emmeans):\n")
  medias <- emmeans(mod_misto, ~Treatment)
  print(summary(medias))
  # Comparações múltiplas 
  cat("\nTukey (glht):\n")
  tukey_glht <- glht(mod_misto, linfct = mcp(Treatment = "Tukey"))
  print(summary(tukey_glht))
  # Letras compactas
  cat("\nLetras do Tukey (glht):\n")
  letras <- cld(tukey_glht)
  print(letras)
}


########################### Variável: CH4_L_d  ##########################
ANOVA (modelo misto):
Type III Analysis of Variance Table with Satterthwaite's method
          Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
Treatment 143.88  47.959     3     6  3.9128 0.07307 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Médias (emmeans):
 Treatment emmean   SE   df lower.CL upper.CL
 D0          14.8 3.83 7.24     5.83     23.8
 D1          19.2 3.83 7.24    10.19     28.2
 D2          22.9 3.83 7.24    13.87     31.9
 D3          16.8 3.83 7.24     7.77     25.8

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

Tukey (glht):

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = formula_mista, data = Data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)   
D1 - D0 == 0    4.359      2.476   1.761  0.29248   
D2 - D0 == 0    8.034      2.476   3.245  0.00656 **
D3 - D0 == 0    1.934      2.476   0.781  0.86296   
D2 - D1 == 0    3.675      2.476   1.485  0.44676   
D3 - D1 == 0   -2.424      2.476  -0.979  0.76129   
D3 - D2 == 0   -6.100      2.476  -2.464  0.06576 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


Letras do Tukey (glht):
  D0   D1   D2   D3 
 "a" "ab"  "b" "ab" 


########################### Variável: CH4mw_L_Kg  ##########################
ANOVA (modelo misto):
Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
Treatment 0.25905 0.08635     3     6  4.7772 0.04958 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Médias (emmeans):
 Treatment emmean    SE   df lower.CL upper.CL
 D0         0.587 0.155 6.97    0.220    0.955
 D1         0.783 0.155 6.97    0.416    1.150
 D2         0.932 0.155 6.97    0.565    1.299
 D3         0.686 0.155 6.97    0.319    1.053

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

Tukey (glht):

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = formula_mista, data = Data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)   
D1 - D0 == 0  0.19580    0.09507   2.060  0.16651   
D2 - D0 == 0  0.34460    0.09507   3.625  0.00162 **
D3 - D0 == 0  0.09832    0.09507   1.034  0.72929   
D2 - D1 == 0  0.14880    0.09507   1.565  0.39866   
D3 - D1 == 0 -0.09748    0.09507  -1.025  0.73458   
D3 - D2 == 0 -0.24627    0.09507  -2.591  0.04695 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


Letras do Tukey (glht):
  D0   D1   D2   D3 
 "a" "ab"  "b"  "a" 


########################### Variável: DMI  ##########################
ANOVA (modelo misto):
Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
Treatment 0.033625 0.011208     3     6  2.1524 0.1949

Médias (emmeans):
 Treatment emmean     SE   df lower.CL upper.CL
 D0          1.22 0.0846 4.71     1.00     1.44
 D1          1.24 0.0846 4.71     1.02     1.46
 D2          1.26 0.0846 4.71     1.04     1.48
 D3          1.34 0.0846 4.71     1.12     1.56

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

Tukey (glht):

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = formula_mista, data = Data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)  
D1 - D0 == 0  0.01875    0.05103   0.367   0.9831  
D2 - D0 == 0  0.03950    0.05103   0.774   0.8662  
D3 - D0 == 0  0.12025    0.05103   2.357   0.0856 .
D2 - D1 == 0  0.02075    0.05103   0.407   0.9773  
D3 - D1 == 0  0.10150    0.05103   1.989   0.1921  
D3 - D2 == 0  0.08075    0.05103   1.583   0.3886  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


Letras do Tukey (glht):
 D0  D1  D2  D3 
"a" "a" "a" "a" 


########################### Variável: Digestibility  ##########################
boundary (singular) fit: see help('isSingular')
ANOVA (modelo misto):
Type III Analysis of Variance Table with Satterthwaite's method
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Treatment 36.426  12.142     3     9  0.3402  0.797

Médias (emmeans):
 Treatment emmean   SE   df lower.CL upper.CL
 D0          73.7 3.51 9.75     65.9     81.6
 D1          70.1 3.51 9.75     62.2     78.0
 D2          70.0 3.51 9.75     62.2     77.9
 D3          71.6 3.51 9.75     63.8     79.5

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

Tukey (glht):

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = formula_mista, data = Data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)
D1 - D0 == 0 -3.64275    4.22461  -0.862    0.824
D2 - D0 == 0 -3.71350    4.22461  -0.879    0.816
D3 - D0 == 0 -2.11575    4.22461  -0.501    0.959
D2 - D1 == 0 -0.07075    4.22461  -0.017    1.000
D3 - D1 == 0  1.52700    4.22461   0.361    0.984
D3 - D2 == 0  1.59775    4.22461   0.378    0.982
(Adjusted p values reported -- single-step method)


Letras do Tukey (glht):
 D0  D1  D2  D3 
"a" "a" "a" "a" 


########################### Variável: CH4_DMI  ##########################
ANOVA (modelo misto):
Type III Analysis of Variance Table with Satterthwaite's method
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Treatment 93.165  31.055     3     6  3.1556 0.1074

Médias (emmeans):
 Treatment emmean   SE   df lower.CL upper.CL
 D0          12.4 3.47 7.25     4.30     20.6
 D1          15.7 3.47 7.25     7.57     23.9
 D2          18.5 3.47 7.25    10.38     26.7
 D3          13.0 3.47 7.25     4.90     21.2

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

Tukey (glht):

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = formula_mista, data = Data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)  
D1 - D0 == 0   3.2749     2.2183   1.476   0.4518  
D2 - D0 == 0   6.0812     2.2183   2.741   0.0315 *
D3 - D0 == 0   0.5983     2.2183   0.270   0.9931  
D2 - D1 == 0   2.8062     2.2183   1.265   0.5853  
D3 - D1 == 0  -2.6767     2.2183  -1.207   0.6226  
D3 - D2 == 0  -5.4829     2.2183  -2.472   0.0643 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


Letras do Tukey (glht):
  D0   D1   D2   D3 
 "a" "ab"  "b" "ab" 
LS0tDQp0aXRsZTogIkRhZG9zIExhdXJhIEluIHZpdm8iDQphdXRob3I6ICJWYWduZXIgT3ZhbmkiDQpkYXRlOiAiMzAvMDUvMjAyNSINCm91dHB1dDoNCiAgaHRtbF9ub3RlYm9vazoNCiAgICB0b2M6IFRSVUUNCiAgICB0b2NfZGVwdGg6IDINCiAgICB0aGVtZTogdW5pdGVkDQotLS0NCg0KKioqDQoqKioNCg0KIyAqKjEgUGFja2FnZXMqKg0KDQpgYGB7cn0NCmxpYnJhcnkobG1lNCkNCmxpYnJhcnkobG1lclRlc3QpDQpsaWJyYXJ5KGNhcikNCmxpYnJhcnkoZW1tZWFucykNCmxpYnJhcnkobXVsdGNvbXBWaWV3KQ0KbGlicmFyeShtdWx0Y29tcCkNCg0KYGBgDQoNCioqKg0KKioqDQoNCiMgKioyIERBVEEqKg0KDQpgYGB7cn0NCkRhdGE9cmVhZC5jc3YoIkM6L1VzZXJzL1NhbXN1bmcvT25lRHJpdmUvRG9jdW1lbnRvcy8xIC0gUMOzc0RvYy8yIC0gRkFQRVNQIC0gU2lsdmlwYXN0b3JpbC9SZXNwb25zYWJpbGlkYWRlcy9Qcm9qZXRvIExhdXJhL2ludml2by5jc3YiKQ0Kc3RyKERhdGEpDQpEYXRhJFRyZWF0bWVudD1hcy5mYWN0b3IoRGF0YSRUcmVhdG1lbnQpDQpEYXRhJFBlcmlvZD1hcy5mYWN0b3IoRGF0YSRQZXJpb2QpDQpEYXRhJEFuaW1hbD1hcy5mYWN0b3IoRGF0YSRBbmltYWwpDQpwcmludChEYXRhKQ0KYGBgDQoNCioqKg0KKioqDQoNCiMgKiozIEFOT1ZBICsgVFVLRVkqKg0KDQpgYGB7cn0NCm5hbWVzKERhdGEpDQojIFZhcmnDoXZlaXMgYSBzZXJlbSBhbmFsaXNhZGFzDQp2YXJpYWJsZXMgPC0gYygiQ0g0X0xfZCIsICJDSDRtd19MX0tnIiwgIkRNSSIsICJEaWdlc3RpYmlsaXR5IiwgIkNINF9ETUkiKQ0KDQojIFRlc3RlcyBkZSBub3JtYWxpZGFkZSBlIGhvbW9jZWRhc3RpY2lkYWRlIChjb20gbG0gcGFyYSBzaW1wbGlmaWNhw6fDo28pDQp0ZXN0ZXMgPC0gbGFwcGx5KHZhcmlhYmxlcywgZnVuY3Rpb24odmFyKSB7DQogIGZvcm11bGEgPC0gYXMuZm9ybXVsYShwYXN0ZSh2YXIsICJ+IFRyZWF0bWVudCArICgxfEFuaW1hbCkgKyAoMXxQZXJpb2QpIikpDQogIG1vZGVsbyA8LSBsbWVyKGZvcm11bGEsIGRhdGEgPSBEYXRhKSAgDQogIHJlcyA8LSByZXNpZHVhbHMobW9kZWxvKQ0KICBzaGFwaXJvX3AgPC0gc2hhcGlyby50ZXN0KHJlcykkcC52YWx1ZQ0KICBsZXZlbmVfcCA8LSBsZXZlbmVUZXN0KHJlcyB+IERhdGEkVHJlYXRtZW50KVsxLCAiUHIoPkYpIl0NCiAgDQogIGRhdGEuZnJhbWUoVmFyaWF2ZWwgPSB2YXIsDQogICAgICAgICAgICAgcF9TaGFwaXJvID0gcm91bmQoc2hhcGlyb19wLCA0KSwNCiAgICAgICAgICAgICBwX0xldmVuZSA9IHJvdW5kKGxldmVuZV9wLCA0KSkNCn0pDQoNCnJlc3VsdGFkb190ZXN0ZXMgPC0gZG8uY2FsbChyYmluZCwgdGVzdGVzKQ0KcHJpbnQocmVzdWx0YWRvX3Rlc3RlcykNCg0KIyBBbsOhbGlzZSBkbyBtb2RlbG8gbWlzdG8gKHF1YWRybyBsYXRpbm8pIGUgY29tcGFyYcOnw6NvIGRlIG3DqWRpYXMNCmZvciAodmFyIGluIHZhcmlhYmxlcykgew0KICBjYXQoIlxuXG4jIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMgVmFyacOhdmVsOiIsIHZhciwgIiAjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjI1xuIikNCiAgZm9ybXVsYV9taXN0YSA8LSBhcy5mb3JtdWxhKHBhc3RlKHZhciwgIn4gVHJlYXRtZW50ICsgKDF8QW5pbWFsKSArICgxfFBlcmlvZCkiKSkNCiAgbW9kX21pc3RvIDwtIGxtZXIoZm9ybXVsYV9taXN0YSwgZGF0YSA9IERhdGEpDQogICMgQU5PVkENCiAgY2F0KCJBTk9WQSAobW9kZWxvIG1pc3RvKTpcbiIpDQogIHByaW50KGFub3ZhKG1vZF9taXN0bykpDQogICMgTcOpZGlhcw0KICBjYXQoIlxuTcOpZGlhcyAoZW1tZWFucyk6XG4iKQ0KICBtZWRpYXMgPC0gZW1tZWFucyhtb2RfbWlzdG8sIH5UcmVhdG1lbnQpDQogIHByaW50KHN1bW1hcnkobWVkaWFzKSkNCiAgIyBDb21wYXJhw6fDtWVzIG3Dumx0aXBsYXMgDQogIGNhdCgiXG5UdWtleSAoZ2xodCk6XG4iKQ0KICB0dWtleV9nbGh0IDwtIGdsaHQobW9kX21pc3RvLCBsaW5mY3QgPSBtY3AoVHJlYXRtZW50ID0gIlR1a2V5IikpDQogIHByaW50KHN1bW1hcnkodHVrZXlfZ2xodCkpDQogICMgTGV0cmFzIGNvbXBhY3Rhcw0KICBjYXQoIlxuTGV0cmFzIGRvIFR1a2V5IChnbGh0KTpcbiIpDQogIGxldHJhcyA8LSBjbGQodHVrZXlfZ2xodCkNCiAgcHJpbnQobGV0cmFzKQ0KfQ0KDQpgYGANCg0KDQo=