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=