Análise das variáveis de ingestão e ganho de peso
Prof. Dr. Walmes Zeviani
2024/Jan/11
# Prof. Dr. Walmes M. Zeviani
# leg.ufpr.br/~walmes · github.com/walmes
# walmes@ufpr.br · @walmeszeviani
# Laboratory of Statistics and Geoinformation (LEG)
# Department of Statistics · Federal University of Paraná
# 2024-jan-11 · Curitiba/PR/Brazil1 Informações da sessão
# Pacotes.library(emmeans)
library(rstatix)
library(segmented)
library(nlme)
library(tidyverse)# Funções.
sig_stars <- function(x) {
stars <- symnum(x,
corr = FALSE,
na = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", " "))
attr(stars, "legend") <- NULL
stars
}
# Para não exibir a correlação entre estimativas.
assignInNamespace("print.correlation",
function(x, title) {
return()
},
ns = "nlme")2 Consumo de ração e água
2.1 Importação e análise exploratória dos dados
# Consumo de ração e água.
# Caminho para o arquivo.
url <- "PESO corporal e _CONSUMO de agua e ração_Estatística(Recuperado Automaticamente).xlsx"
# Lê as duas abas.
tb <-
imap(c("ingestao_racao" = 1, "ingestao_agua" = 2),
.f = function(.x, .y) {
tb <- readxl::read_excel(path = url,
sheet = .x,
range = "A2:E16")
names(tb) <- c("dias",
"macho_2000",
"femea_2000",
"macho_0",
"femea_0")
tb <- tb |>
pivot_longer(cols = -dias,
names_to = c("sexo", "trat"),
values_to = .y,
names_sep = "_") |>
mutate(sexo = factor(sexo, levels = c("macho", "femea")),
trat = factor(trat))
return(tb)
})
# Coloca uma ao lado da outra e melhora a disposição de linhas e
# colunas.
tb <- full_join(tb$ingestao_racao,
tb$ingestao_agua,
by = c("dias", "sexo", "trat")) |>
arrange(sexo, trat, dias) |>
relocate(sexo, trat, dias)
# # Exibe os dados.
# tb |>
# reactable::reactable(pagination = FALSE,
# sortable = TRUE,
# width = 800)
gg1 <-
ggplot(data = tb,
mapping = aes(x = dias,
y = ingestao_racao,
color = factor(trat))) +
facet_wrap(facets = ~sexo) +
geom_point() +
# stat_summary(geom = "line", fun = "mean") +
geom_smooth(se = FALSE, span = 0.95, linewidth = 0.5) +
# geom_smooth(method = "lm",
# se = FALSE,
# formula = y ~ poly(x, degree = 2)) +
labs(x = "Dias",
y = "Ingestão de ração (g/dia)",
color = "Tratamento")
gg2 <-
ggplot(data = tb,
mapping = aes(x = dias,
y = ingestao_agua,
color = factor(trat))) +
facet_wrap(facets = ~sexo) +
geom_point() +
# stat_summary(geom = "line", fun = "mean") +
geom_smooth(se = FALSE, span = 0.95, linewidth = 0.5) +
# geom_smooth(method = "lm",
# se = FALSE,
# formula = y ~ poly(x, degree = 2)) +
labs(x = "Dias",
y = "Ingestão de água (ml/dia)",
color = "Tratamento")gridExtra::grid.arrange(gg1, gg2, nrow = 2) |>
suppressMessages()
Figure 2.1: Ingestão de ração e água para os animais de cada sexo em função dos dias para cada tratamento.
#=======================================================================
# Ajuste com modelos polinomiais.2.2 Análise da ingestão de ração
# Ingestão de ração.
# Especificação e ajuste do modelo.
k <- 2
m0 <- lm(ingestao_racao ~ sexo * trat * poly(dias, degree = k),
data = tb)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: ingestao_racao
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 1 5401.8 5401.8 61.7511 6.453e-10 ***
## trat 1 848.6 848.6 9.7014 0.003236 **
## poly(dias, degree = k) 2 640.5 320.3 3.6610 0.033827 *
## sexo:trat 1 52.1 52.1 0.5953 0.444516
## sexo:poly(dias, degree = k) 2 460.3 230.2 2.6313 0.083290 .
## trat:poly(dias, degree = k) 2 2423.2 1211.6 13.8505 2.160e-05 ***
## sexo:trat:poly(dias, degree = k) 2 74.7 37.3 0.4269 0.655188
## Residuals 44 3849.0 87.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Abandona termos não significativos.
m1 <- update(m0, ~ sexo + trat * poly(dias, degree = k))
anova(m1, m0)## Analysis of Variance Table
##
## Model 1: ingestao_racao ~ sexo + trat + poly(dias, degree = k) + trat:poly(dias,
## degree = k)
## Model 2: ingestao_racao ~ sexo * trat * poly(dias, degree = k)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 49 4436.1
## 2 44 3849.0 5 587.11 1.3423 0.2645# Modelo após o abandono de termos.
anova(m1)## Analysis of Variance Table
##
## Response: ingestao_racao
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 1 5401.8 5401.8 59.6669 5.052e-10 ***
## trat 1 848.6 848.6 9.3739 0.003568 **
## poly(dias, degree = k) 2 640.5 320.3 3.5375 0.036725 *
## trat:poly(dias, degree = k) 2 2423.2 1211.6 13.3830 2.305e-05 ***
## Residuals 49 4436.1 90.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary(m0)$r.squared
# summary(m1)$r.squared
# Cria as combinações de sexo x tratamento e dias com espaçamento fino
# para fazer os gráficos.
tb_pred <- tb |>
expand(nesting(sexo, trat),
dias = seq(min(dias), max(dias), by = 0.25))
tb_pred <-
predict(m1, newdata = tb_pred, interval = "confidence") |>
as.data.frame() |>
bind_cols(tb_pred)
gg <-
ggplot(data = tb,
mapping = aes(x = dias,
y = ingestao_racao,
color = trat)) +
facet_wrap(facets = ~sexo) +
geom_point() +
geom_ribbon(data = tb_pred,
inherit.aes = FALSE,
show.legend = FALSE,
mapping = aes(x = dias,
ymin = lwr,
ymax = upr,
group = trat,
fill = trat),
alpha = 0.2) +
geom_line(data = tb_pred,
inherit.aes = FALSE,
show.legend = FALSE,
mapping = aes(x = dias,
y = fit,
group = trat,
color = trat)) +
labs(x = "Dias",
y = "Ingestão de ração (g/dia)",
color = "Tratamento")
gg
# Reparametriza para facilitar a interpretação.
# m2 <- update(m0, ~ sexo + trat/poly(dias, degree = k, raw = TRUE))
m2 <- update(m0, ~ sexo + trat/(dias + I(dias^2)))
summary(m2)##
## Call:
## lm(formula = ingestao_racao ~ sexo + trat + trat:dias + trat:I(dias^2),
## data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.9990 -4.7953 0.4098 5.8502 24.6470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.2692 6.3964 19.584 < 2e-16 ***
## sexofemea -19.6429 2.5430 -7.724 5.05e-10 ***
## trat2000 31.5082 8.8653 3.554 0.000851 ***
## trat0:dias 10.1073 1.9226 5.257 3.19e-06 ***
## trat2000:dias -3.8790 1.9226 -2.018 0.049134 *
## trat0:I(dias^2) -0.6212 0.1247 -4.983 8.21e-06 ***
## trat2000:I(dias^2) 0.2837 0.1247 2.275 0.027313 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.515 on 49 degrees of freedom
## Multiple R-squared: 0.6774, Adjusted R-squared: 0.6379
## F-statistic: 17.15 on 6 and 49 DF, p-value: 1.478e-10# Coeficientes da equação de segundo grau.
tb_coef <-
coef(m2) |>
as.list() |>
with({
rbind("macho:0" = c(b0 = `(Intercept)`,
b1 = `trat0:dias`,
b2 = `trat0:I(dias^2)`),
"femea:0" = c(b0 = `(Intercept)` + sexofemea,
b1 = `trat0:dias`,
b2 = `trat0:I(dias^2)`),
"macho:2000" = c(b0 = `(Intercept)` + trat2000,
b1 = `trat2000:dias`,
b2 = `trat2000:I(dias^2)`),
"femea:2000" = c(b0 = `(Intercept)` + sexofemea + trat2000,
b1 = `trat2000:dias`,
b2 = `trat2000:I(dias^2)`))
}) |>
as.data.frame() |>
rownames_to_column(var = "sexo_trat") |>
separate(sexo_trat, into = c("sexo", "trat"), sep = ":")
# tb_coef
# Adiciona o R2 em cada combinação de sexo e tratamento.
tb_coef <- tb |>
add_column(fitted = fitted(m1)) |>
group_by(sexo, trat) |>
summarise(R2 = cor(ingestao_racao, fitted)^2,
.groups = "drop") |>
ungroup() |>
left_join(tb_coef, y = _, by = c("sexo", "trat"))
tb_coef## sexo trat b0 b1 b2 R2
## 1 macho 0 125.2692 10.107349 -0.6212225 0.6975827
## 2 femea 0 105.6264 10.107349 -0.6212225 0.5678814
## 3 macho 2000 156.7775 -3.878984 0.2836538 0.1885979
## 4 femea 2000 137.1346 -3.878984 0.2836538 0.1393804# # Para conferir os coeficientes.
# gg <-
# ggplot(data = tb,
# mapping = aes(x = dias,
# y = ingestao_racao,
# color = trat)) +
# facet_wrap(facets = ~sexo) +
# geom_point() +
# geom_ribbon(data = tb_pred,
# inherit.aes = FALSE,
# show.legend = FALSE,
# mapping = aes(x = dias,
# ymin = lwr,
# ymax = upr,
# group = trat,
# fill = trat),
# alpha = 0.2) +
# labs(x = "Dias",
# y = "Ingestão de ração (g/dia)",
# color = "Tratamento")
# fun <- function(x, b0, b1, b2) {
# b0 + b1 * x + b2 * x^2
# }
# for (i in seq(nrow(tb_coef))) {
# gg <- gg +
# geom_function(data = tb_coef[i, c("sexo", "trat")],
# inherit.aes = FALSE,
# fun = fun,
# args = as.list(tb_coef[i, c("b0", "b1", "b2")]))
# }
# gg2.2.1 Resultados principais
# Conteúdo que vai no relatório.
# Quadro de análise de variância completo.
anova(m0) |>
as.data.frame(check.names = FALSE) |>
rownames_to_column(var = "Effect") |>
mutate(Effect = gsub("poly\\(dias, degree = k\\)", "dias", Effect)) |>
mutate(" " = sig_stars(`Pr(>F)`))## Effect Df Sum Sq Mean Sq F value Pr(>F)
## 1 sexo 1 5401.78571 5401.78571 61.7511429 6.452950e-10 ***
## 2 trat 1 848.64286 848.64286 9.7013597 3.235723e-03 **
## 3 dias 2 640.51133 320.25567 3.6610400 3.382715e-02 *
## 4 sexo:trat 1 52.07143 52.07143 0.5952606 4.445160e-01
## 5 sexo:dias 2 460.34815 230.17407 2.6312617 8.328993e-02 .
## 6 trat:dias 2 2423.18880 1211.59440 13.8504826 2.159777e-05 ***
## 7 sexo:trat:dias 2 74.69155 37.34578 0.4269226 6.551880e-01
## 8 Residuals 44 3848.97445 87.47669 NA NA# Resumo do quadro de análise de variância.
rstatix::anova_test(m0) |>
map(~.x) |>
# as_tibble() |>
as.data.frame(check.names = FALSE) |>
dplyr::select(-ges) |>
mutate(Effect = gsub("poly\\(dias, degree = k\\)", "dias", Effect))## Effect DFn DFd F p p<.05
## 1 sexo 1 44 61.751 6.45e-10 *
## 2 trat 1 44 9.701 3.00e-03 *
## 3 dias 2 44 3.661 3.40e-02 *
## 4 sexo:trat 1 44 0.595 4.45e-01
## 5 sexo:dias 2 44 2.631 8.30e-02
## 6 trat:dias 2 44 13.850 2.16e-05 *
## 7 sexo:trat:dias 2 44 0.427 6.55e-01# Resultados do modelo ajustado.
summary(m2)##
## Call:
## lm(formula = ingestao_racao ~ sexo + trat + trat:dias + trat:I(dias^2),
## data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.9990 -4.7953 0.4098 5.8502 24.6470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.2692 6.3964 19.584 < 2e-16 ***
## sexofemea -19.6429 2.5430 -7.724 5.05e-10 ***
## trat2000 31.5082 8.8653 3.554 0.000851 ***
## trat0:dias 10.1073 1.9226 5.257 3.19e-06 ***
## trat2000:dias -3.8790 1.9226 -2.018 0.049134 *
## trat0:I(dias^2) -0.6212 0.1247 -4.983 8.21e-06 ***
## trat2000:I(dias^2) 0.2837 0.1247 2.275 0.027313 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.515 on 49 degrees of freedom
## Multiple R-squared: 0.6774, Adjusted R-squared: 0.6379
## F-statistic: 17.15 on 6 and 49 DF, p-value: 1.478e-10# Coeficientes.
tb_coef## sexo trat b0 b1 b2 R2
## 1 macho 0 125.2692 10.107349 -0.6212225 0.6975827
## 2 femea 0 105.6264 10.107349 -0.6212225 0.5678814
## 3 macho 2000 156.7775 -3.878984 0.2836538 0.1885979
## 4 femea 2000 137.1346 -3.878984 0.2836538 0.1393804# Gráfico final.ggFigure 2.2: Ajuste do modelo para ingestão de ração para os animais de cada sexo em função dos dias para cada tratamento.
2.3 Análise da ingestão de água
# Ingestão de água.
# Especificação e ajuste do modelo.
k <- 2
m0 <- lm(ingestao_agua ~ sexo * trat * poly(dias, degree = k),
data = tb)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: ingestao_agua
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 1 11429 11428.6 8.3233 0.00604 **
## trat 1 350 350.0 0.2549 0.61616
## poly(dias, degree = k) 2 9672 4836.0 3.5220 0.03812 *
## sexo:trat 1 7 7.1 0.0052 0.94283
## sexo:poly(dias, degree = k) 2 5991 2995.4 2.1816 0.12492
## trat:poly(dias, degree = k) 2 296 148.1 0.1079 0.89799
## sexo:trat:poly(dias, degree = k) 2 4011 2005.7 1.4607 0.24310
## Residuals 44 60415 1373.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Abandona termos não significativos.
m1 <- update(m0, ~ sexo + trat + poly(dias, degree = k))
anova(m1, m0)## Analysis of Variance Table
##
## Model 1: ingestao_agua ~ sexo + trat + poly(dias, degree = k)
## Model 2: ingestao_agua ~ sexo * trat * poly(dias, degree = k)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 51 70721
## 2 44 60415 7 10306 1.0722 0.3971# Modelo após o abandono de termos.
anova(m1)## Analysis of Variance Table
##
## Response: ingestao_agua
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 1 11429 11428.6 8.2417 0.005947 **
## trat 1 350 350.0 0.2524 0.617552
## poly(dias, degree = k) 2 9672 4836.0 3.4875 0.038055 *
## Residuals 51 70721 1386.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary(m0)$r.squared
# summary(m1)$r.squared
# Cria as combinações de sexo x tratamento e dias com espaçamento fino
# para fazer os gráficos.
tb_pred <- tb |>
expand(nesting(sexo, trat),
dias = seq(min(dias), max(dias), by = 0.25))
tb_pred <-
predict(m1, newdata = tb_pred, interval = "confidence") |>
as.data.frame() |>
bind_cols(tb_pred)
gg <-
ggplot(data = tb,
mapping = aes(x = dias,
y = ingestao_agua,
color = trat)) +
facet_wrap(facets = ~sexo) +
geom_point() +
geom_ribbon(data = tb_pred,
inherit.aes = FALSE,
show.legend = FALSE,
mapping = aes(x = dias,
ymin = lwr,
ymax = upr,
group = trat,
fill = trat),
alpha = 0.2) +
geom_line(data = tb_pred,
inherit.aes = FALSE,
show.legend = FALSE,
mapping = aes(x = dias,
y = fit,
group = trat,
color = trat)) +
labs(x = "Dias",
y = "Ingestão de água (ml/dia)",
color = "Tratamento")
# Reparametriza para facilitar a interpretação.
# m2 <- update(m0, ~ sexo + trat/poly(dias, degree = k, raw = TRUE))
m2 <- update(m0, ~ sexo + trat + (dias + I(dias^2)))
summary(m2)##
## Call:
## lm(formula = ingestao_agua ~ sexo + trat + dias + I(dias^2),
## data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -70.519 -27.535 -1.173 22.570 89.955
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 247.7610 18.7211 13.234 < 2e-16 ***
## sexofemea -28.5714 9.9523 -2.871 0.00595 **
## trat2000 5.0000 9.9523 0.502 0.61755
## dias 13.1885 5.3207 2.479 0.01653 *
## I(dias^2) -0.9049 0.3450 -2.623 0.01148 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.24 on 51 degrees of freedom
## Multiple R-squared: 0.2327, Adjusted R-squared: 0.1725
## F-statistic: 3.867 on 4 and 51 DF, p-value: 0.008077# Coeficientes da equação de segundo grau.
tb_coef <-
coef(m2) |>
as.list() |>
with({
rbind("macho:0" = c(b0 = `(Intercept)`,
b1 = `dias`,
b2 = `I(dias^2)`),
"femea:0" = c(b0 = `(Intercept)` + sexofemea,
b1 = `dias`,
b2 = `I(dias^2)`),
"macho:2000" = c(b0 = `(Intercept)` + trat2000,
b1 = `dias`,
b2 = `I(dias^2)`),
"femea:2000" = c(b0 = `(Intercept)` + sexofemea + trat2000,
b1 = `dias`,
b2 = `I(dias^2)`))
}) |>
as.data.frame() |>
rownames_to_column(var = "sexo_trat") |>
separate(sexo_trat, into = c("sexo", "trat"), sep = ":")
# tb_coef
# Adiciona o R2 em cada combinação de sexo e tratamento.
tb_coef <- tb |>
add_column(fitted = fitted(m1)) |>
group_by(sexo, trat) |>
summarise(R2 = cor(ingestao_racao, fitted)^2,
.groups = "drop") |>
ungroup() |>
left_join(tb_coef, y = _, by = c("sexo", "trat"))
tb_coef## sexo trat b0 b1 b2 R2
## 1 macho 0 247.7610 13.18853 -0.9048764 0.4731468
## 2 femea 0 219.1896 13.18853 -0.9048764 0.6086095
## 3 macho 2000 252.7610 13.18853 -0.9048764 0.1438854
## 4 femea 2000 224.1896 13.18853 -0.9048764 0.1777327# # Para conferir os coeficientes.
# gg <-
# ggplot(data = tb,
# mapping = aes(x = dias,
# y = ingestao_agua,
# color = trat)) +
# facet_wrap(facets = ~sexo) +
# geom_point() +
# geom_ribbon(data = tb_pred,
# inherit.aes = FALSE,
# show.legend = FALSE,
# mapping = aes(x = dias,
# ymin = lwr,
# ymax = upr,
# group = trat,
# fill = trat),
# alpha = 0.2) +
# labs(x = "Dias",
# y = "Ingestão de água (ml/dia)",
# color = "Tratamento")
# fun <- function(x, b0, b1, b2) {
# b0 + b1 * x + b2 * x^2
# }
# for (i in seq(nrow(tb_coef))) {
# gg <- gg +
# geom_function(data = tb_coef[i, c("sexo", "trat")],
# inherit.aes = FALSE,
# fun = fun,
# args = as.list(tb_coef[i, c("b0", "b1", "b2")]))
# }
# gg2.3.1 Resultados principais
# Conteúdo que vai no relatório.
# Quadro de análise de variância completo.
anova(m0) |>
as.data.frame(check.names = FALSE) |>
rownames_to_column(var = "Effect") |>
mutate(Effect = gsub("poly\\(dias, degree = k\\)", "dias", Effect)) |>
mutate(" " = sig_stars(`Pr(>F)`))## Effect Df Sum Sq Mean Sq F value Pr(>F)
## 1 sexo 1 11428.571429 11428.571429 8.323348256 0.006040248 **
## 2 trat 1 350.000000 350.000000 0.254902540 0.616163475
## 3 dias 2 9672.012363 4836.006181 3.522029316 0.038121020 *
## 4 sexo:trat 1 7.142857 7.142857 0.005202093 0.942828869
## 5 sexo:dias 2 5990.858516 2995.429258 2.181550078 0.124924603
## 6 trat:dias 2 296.188187 148.094093 0.107855887 0.897993576
## 7 sexo:trat:dias 2 4011.407967 2005.703984 1.460740116 0.243098001
## 8 Residuals 44 60415.247253 1373.073801 NA NA# Resumo do quadro de análise de variância.
rstatix::anova_test(m0) |>
map(~.x) |>
# as_tibble() |>
as.data.frame(check.names = FALSE) |>
dplyr::select(-ges) |>
mutate(Effect = gsub("poly\\(dias, degree = k\\)", "dias", Effect))## Effect DFn DFd F p p<.05
## 1 sexo 1 44 8.323 0.006 *
## 2 trat 1 44 0.255 0.616
## 3 dias 2 44 3.522 0.038 *
## 4 sexo:trat 1 44 0.005 0.943
## 5 sexo:dias 2 44 2.182 0.125
## 6 trat:dias 2 44 0.108 0.898
## 7 sexo:trat:dias 2 44 1.461 0.243# Resultados do modelo ajustado.
summary(m2)##
## Call:
## lm(formula = ingestao_agua ~ sexo + trat + dias + I(dias^2),
## data = tb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -70.519 -27.535 -1.173 22.570 89.955
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 247.7610 18.7211 13.234 < 2e-16 ***
## sexofemea -28.5714 9.9523 -2.871 0.00595 **
## trat2000 5.0000 9.9523 0.502 0.61755
## dias 13.1885 5.3207 2.479 0.01653 *
## I(dias^2) -0.9049 0.3450 -2.623 0.01148 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.24 on 51 degrees of freedom
## Multiple R-squared: 0.2327, Adjusted R-squared: 0.1725
## F-statistic: 3.867 on 4 and 51 DF, p-value: 0.008077# Coeficientes.
tb_coef## sexo trat b0 b1 b2 R2
## 1 macho 0 247.7610 13.18853 -0.9048764 0.4731468
## 2 femea 0 219.1896 13.18853 -0.9048764 0.6086095
## 3 macho 2000 252.7610 13.18853 -0.9048764 0.1438854
## 4 femea 2000 224.1896 13.18853 -0.9048764 0.1777327# Gráfico final.gg
Figure 2.3: Ajuste do modelo para ingestão de água para os animais de cada sexo em função dos dias para cada tratamento.
3 Curva de crescimento
3.1 Importação e análise exploratória dos dados
#=======================================================================
# Curva de crescimento.
# Lê os dados de crescimento.
url <- "PESO corporal e _CONSUMO de agua e ração_Estatística(Recuperado Automaticamente).xlsx"
tb <- readxl::read_excel(path = url,
sheet = 3,
range = "A4:U16",
col_names = FALSE) |>
suppressMessages()
# str(tb)
# Cria nome para as colunas.
col_names <-
expand.grid(animal = 1:5,
sexo = c("macho", "femea"),
trat = c(2000, 0),
KEEP.OUT.ATTRS = FALSE) |>
unite("grupo", c("animal", "sexo", "trat")) |>
add_row(grupo = "dias", .before = 1) |>
pull(grupo)
names(tb) <- col_names
# str(tb)
# Empilha os dados.
tb <- tb |>
pivot_longer(cols = -dias,
names_to = c("animal", "sexo", "trat"),
values_to = "peso", names_sep = "_") |>
mutate(sexo = factor(sexo, levels = c("macho", "femea")),
trat = factor(trat),
animal = str_c(sexo, trat, animal, sep = "_")) |>
arrange(trat, sexo, animal)
str(tb)## tibble [260 × 5] (S3: tbl_df/tbl/data.frame)
## $ dias : num [1:260] 1 2 3 4 5 6 7 8 9 10 ...
## $ animal: chr [1:260] "macho_0_1" "macho_0_1" "macho_0_1" "macho_0_1" ...
## $ sexo : Factor w/ 2 levels "macho","femea": 1 1 1 1 1 1 1 1 1 1 ...
## $ trat : Factor w/ 2 levels "0","2000": 1 1 1 1 1 1 1 1 1 1 ...
## $ peso : num [1:260] 263 264 266 269 272 275 278 281 286 289 ...# Subsitui o valor de um animal em um dia específico pois parece ser um
# erro de digitação.
tb <- tb |>
mutate(peso = replace(peso,
animal == "macho_2000_1" & dias == 7,
295))
# # Confere.
# tb |>
# filter(animal == "macho_2000_1")
# tb |>
# reactable::reactable(pagination = FALSE,
# sortable = TRUE,
# width = 600)
ggplot(data = tb,
mapping = aes(x = dias,
y = peso,
color = factor(trat),
group = animal)) +
facet_wrap(facets = ~sexo) +
geom_point() +
geom_line() +
labs(x = "Dias",
y = "Peso (g)",
color = "Tratamento")
3.2 Ajuste do modelo
# Modelo de efeitos aleatórios.
# Define uma estrutura de dados agrupados para explorar os gráficos
# automáticos.
tb_grp <- groupedData(peso ~ dias | animal, data = tb)
# Ajusta o modelo aos dados.
m0 <- lme(peso ~ sexo * trat * dias,
random = ~1 + dias | animal,
correlation = corAR1(form = ~dias | animal),
data = tb_grp,
method = "ML")
summary(m0)## Linear mixed-effects model fit by maximum likelihood
## Data: tb_grp
## AIC BIC logLik
## 1078.722 1125.011 -526.3612
##
## Random effects:
## Formula: ~1 + dias | animal
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 5.1046572 (Intr)
## dias 0.2219867 0.284
## Residual 3.6483261
##
## Correlation Structure: AR(1)
## Formula: ~dias | animal
## Parameter estimate(s):
## Phi
## 0.8987075
## Fixed effects: peso ~ sexo * trat * dias
## Value Std.Error DF t-value p-value
## (Intercept) 266.16282 2.898131 236 91.83946 0.0000
## sexofemea -93.44734 4.098577 16 -22.79995 0.0000
## trat2000 4.84760 4.098577 16 1.18275 0.2542
## dias 3.22819 0.194092 236 16.63224 0.0000
## sexofemea:trat2000 13.50321 5.796263 16 2.32964 0.0332
## sexofemea:dias -0.23063 0.274488 236 -0.84021 0.4016
## trat2000:dias 0.26104 0.274488 236 0.95101 0.3426
## sexofemea:trat2000:dias -0.70534 0.388184 236 -1.81703 0.0705
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.87304690 -0.29261197 0.06511862 0.46843246 2.39306000
##
## Number of Observations: 260
## Number of Groups: 20# Quadro de testes de Wald.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 236 27814.273 <.0001
## sexo 1 16 985.226 <.0001
## trat 1 16 15.986 0.0010
## dias 1 236 998.821 <.0001
## sexo:trat 1 16 3.890 0.0661
## sexo:dias 1 236 9.032 0.0029
## trat:dias 1 236 0.223 0.6373
## sexo:trat:dias 1 236 3.302 0.0705# plot(augPred(m0, level = 0:1))
# plot(ACF(m0))
# Simplifica o modelo abadonando termos não significativos.
m1 <- update(m0, ~ trat + sexo * dias)
anova(m0, m1)## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 13 1078.722 1125.011 -526.3612
## m1 2 10 1079.222 1114.829 -529.6111 1 vs 2 6.499927 0.0897# Quadro de anova.
anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 238 25535.347 <.0001
## trat 1 17 13.371 0.0020
## sexo 1 17 877.946 <.0001
## dias 1 238 854.433 <.0001
## sexo:dias 1 238 7.718 0.0059# Cria as combinações de sexo x tratamento e dias com espaçamento fino
# para fazer os gráficos.
dias_seq <- with(tb, seq(min(dias), max(dias), by = 0.25))
# tb_pred <- tb |>
# expand(nesting(sexo, trat),
# dias = dias_seq)
# tb_pred$fit <- predict(m1, newdata = tb_pred, level = 0)
tb_pred <-
emmeans(m1, ~ sexo + trat + dias,
at = list(dias = dias_seq)) |>
as.data.frame() |>
rename(fit = emmean,
lwr = lower.CL,
upr = upper.CL)
gg <-
ggplot(data = tb,
mapping = aes(x = dias,
y = peso)) +
facet_wrap(facets = ~sexo) +
geom_line(mapping = aes(group = animal),
linewidth = 0.3,
alpha = 0.5,
color = "gray30") +
geom_point(mapping = aes(color = trat),
alpha = 0.5,
cex = 0.75,
pch = 19) +
geom_ribbon(data = tb_pred,
inherit.aes = FALSE,
show.legend = FALSE,
mapping = aes(x = dias,
ymin = lwr,
ymax = upr,
group = trat,
fill = trat),
alpha = 0.2) +
geom_line(data = tb_pred,
inherit.aes = FALSE,
show.legend = FALSE,
linewidth = 1,
mapping = aes(x = dias,
y = fit,
group = trat,
color = trat)) +
labs(x = "Dias",
y = "Peso (g)",
color = "Tratamento")
gg
# Reparametriza para facilitar a interpretação.
m2 <- update(m1, ~ trat + sexo/dias)
summary(m2)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 263.029537 2.762136 238 95.226847 1.769224e-191
## trat2000 11.105813 3.037209 17 3.656585 1.953576e-03
## sexofemea -86.680315 3.262894 17 -26.565474 2.770129e-15
## sexomacho:dias 3.360280 0.148464 238 22.633640 2.858764e-61
## sexofemea:dias 2.776989 0.148464 238 18.704799 1.238473e-48summary(m1)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 263.0295372 2.7621364 238 95.226847 1.769224e-191
## trat2000 11.1058132 3.0372092 17 3.656585 1.953576e-03
## sexofemea -86.6803153 3.2628936 17 -26.565474 2.770129e-15
## dias 3.3602798 0.1484640 238 22.633640 2.858764e-61
## sexofemea:dias -0.5832913 0.2099597 238 -2.778110 5.903578e-03# Coeficientes da equação de segundo grau.
tb_coef <-
fixef(m2) |>
as.list() |>
with({
rbind("macho:0" = c(b0 = `(Intercept)`,
b1 = `sexomacho:dias`),
"femea:0" = c(b0 = `(Intercept)` + sexofemea,
b1 = `sexofemea:dias`),
"macho:2000" = c(b0 = `(Intercept)` + trat2000,
b1 = `sexomacho:dias`),
"femea:2000" = c(b0 = `(Intercept)` + sexofemea + trat2000,
b1 = `sexofemea:dias`))
}) |>
as.data.frame() |>
rownames_to_column(var = "sexo_trat") |>
separate(sexo_trat, into = c("sexo", "trat"), sep = ":")
tb_coef## sexo trat b0 b1
## 1 macho 0 263.0295 3.360280
## 2 femea 0 176.3492 2.776989
## 3 macho 2000 274.1354 3.360280
## 4 femea 2000 187.4550 2.776989# Adiciona o R2 em cada combinação de sexo e tratamento.
tb_coef <- tb |>
add_column(fitted = fitted(m1)) |>
group_by(sexo, trat) |>
summarise(R2 = cor(peso, fitted)^2) |>
ungroup() |>
left_join(tb_coef, y = _, by = c("sexo", "trat"))## `summarise()` has grouped output by 'sexo'. You can override using the
## `.groups` argument.tb_coef## sexo trat b0 b1 R2
## 1 macho 0 263.0295 3.360280 0.9879994
## 2 femea 0 176.3492 2.776989 0.9852488
## 3 macho 2000 274.1354 3.360280 0.9551186
## 4 femea 2000 187.4550 2.776989 0.98208053.2.1 Resultados principais
# Conteúdo que vai no relatório.
# Quadro de testes de Wald para os termos de efeito fixo.
anova(m0) |>
as.data.frame() |>
rownames_to_column(var = "Effect") |>
as_tibble() |>
mutate(stars = sig_stars(`p-value`)) |>
# mutate(across(c("F-value", "p-value"),
# ~formatC(.x, digits = 4))) |>
filter(Effect != "(Intercept)")## # A tibble: 7 × 6
## Effect numDF denDF `F-value` `p-value` stars
## <chr> <int> <dbl> <dbl> <dbl> <noquote>
## 1 sexo 1 16 985. 8.88e-16 ***
## 2 trat 1 16 16.0 1.04e- 3 **
## 3 dias 1 236 999. 0 ***
## 4 sexo:trat 1 16 3.89 6.61e- 2 .
## 5 sexo:dias 1 236 9.03 2.94e- 3 **
## 6 trat:dias 1 236 0.223 6.37e- 1
## 7 sexo:trat:dias 1 236 3.30 7.05e- 2 .# Estimativas dos coeficientes do modelo final.
tb_coef## sexo trat b0 b1 R2
## 1 macho 0 263.0295 3.360280 0.9879994
## 2 femea 0 176.3492 2.776989 0.9852488
## 3 macho 2000 274.1354 3.360280 0.9551186
## 4 femea 2000 187.4550 2.776989 0.9820805# (opcional) Estimativas das variâncias dos termos de efeito aleatório.
VarCorr(m0)## animal = pdLogChol(1 + dias)
## Variance StdDev Corr
## (Intercept) 26.05752543 5.1046572 (Intr)
## dias 0.04927811 0.2219867 0.284
## Residual 13.31028335 3.6483261# Gráfico final.ggFigure 3.1: Ajuste do modelo para o peso dos animais de cada sexo em função dos dias para cada tratamento.