Definições da sessão
#-----------------------------------------------------------------------
# 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á
# 2020-out-30 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
library(nlme)
library(multcomp)
library(tidyverse)
library(readxl)
source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/pairwise.R")
#-----------------------------------------------------------------------
# Importação.
tb <- read_xlsx("Dieta3.xlsx", sheet = 1) %>%
mutate(Dieta = factor(Dieta))
str(tb)
## tibble [100 × 3] (S3: tbl_df/tbl/data.frame)
## $ Dias : num [1:100] 35 40 45 50 55 60 65 70 75 35 ...
## $ Dieta: Factor w/ 11 levels "REF","SFA","SFAE",..: 1 1 1 1 1 1 1 1 1 2 ...
## $ Peso : num [1:100] 736 847 1027 1220 1418 ...
#-----------------------------------------------------------------------
# Análise exploratória.
ggplot(data = tb,
mapping = aes(x = Dias, y = Peso)) +
facet_wrap(facets = ~Dieta) +
geom_point() +
geom_line()
ggplot(data = tb,
mapping = aes(x = Dias, y = Peso, color = Dieta)) +
geom_point() +
geom_line()
ggplot(data = filter(tb, Dias >= 40),
mapping = aes(x = Dias, y = Peso)) +
facet_wrap(facets = ~Dieta) +
geom_point() +
geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'
# NOTE: se deixar de fora a primeira data, o comportamento de
# crescimento é linear.
Ajuste do modelo linear
#-----------------------------------------------------------------------
# Ajuste de modelo linear.
m0 <- lm(Peso ~ Dieta * Dias, data = filter(tb, Dias >= 40))
# Inspeciona os resíduos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# NOTE: parece apresentar falta de ajuste.
ggplot(data = cbind(m0$model, r = residuals(m0), f = fitted(m0)),
mapping = aes(x = Dias, y = r)) +
facet_wrap(facets = ~Dieta) +
geom_point() +
geom_smooth(se = FALSE, col = "red", size = 0.25, span = 1) +
geom_hline(yintercept = 0, col = "gray", lty = 2)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# NOTE: não há confirmação visual clara de falta de ajuste.
# Quadro de ANOVA.
anova(m0)
## Analysis of Variance Table
##
## Response: Peso
## Df Sum Sq Mean Sq F value Pr(>F)
## Dieta 10 978022 97802 616.7 < 2.2e-16 ***
## Dias 1 11870631 11870631 74851.6 < 2.2e-16 ***
## Dieta:Dias 10 151294 15129 95.4 < 2.2e-16 ***
## Residuals 67 10625 159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Existe interação entre Dieta e Dias. Ou seja, a taxa de crescimento
# não é a mesma, as retas não são paralelas.
L <- cbind(1, contrasts(tb$Dieta))
a <- m0$assign
# Interceptos.
L %*% coef(m0)[a <= 1]
## [,1]
## REF -664.6315
## SFA -544.9421
## SFAE -420.4936
## SSA -486.2871
## SSF -247.9463
## SSFA -641.6093
## SSFAE -512.3814
## SSFE -414.3286
## SSM -653.9408
## SSMA -593.6279
## SSME -647.9874
# Taxas.
L %*% coef(m0)[a > 1]
## [,1]
## REF 37.80831
## SFA 30.58843
## SFAE 28.03771
## SSA 32.23543
## SSF 24.68976
## SSFA 34.42286
## SSFAE 32.17116
## SSFE 27.07721
## SSM 34.24467
## SSMA 33.36857
## SSME 33.79248
# Teste de hipótese para as taxas.
Ltx <- cbind(0 * L, L)
Ltx %*% coef(m0)
## [,1]
## REF 37.80831
## SFA 30.58843
## SFAE 28.03771
## SSA 32.23543
## SSF 24.68976
## SSFA 34.42286
## SSFAE 32.17116
## SSFE 27.07721
## SSM 34.24467
## SSMA 33.36857
## SSME 33.79248
# Comparações múltiplas para as taxas.
test <- apmc(X = Ltx, model = m0, focus = "Dieta", test = "fdr") %>%
arrange(fit) %>%
mutate(cld2 = ordered_cld(let = cld, means = fit))
test
## Dieta fit lwr upr cld cld2
## 1 SSF 24.68976 23.91404 25.46548 f g
## 2 SSFE 27.07721 26.30150 27.85293 g f
## 3 SFAE 28.03771 27.26200 28.81343 g f
## 4 SFA 30.58843 29.81271 31.36415 e e
## 5 SSFAE 32.17116 31.47990 32.86243 b d
## 6 SSA 32.23543 31.45971 33.01115 bc cd
## 7 SSMA 33.36857 32.59285 34.14429 cd bc
## 8 SSME 33.79248 33.01676 34.56819 d b
## 9 SSM 34.24467 33.46895 35.02039 d b
## 10 SSFA 34.42286 33.64714 35.19858 d b
## 11 REF 37.80831 37.03259 38.58403 a a
# Estimativas intervalares com resultado do teste de hipótese.
ggplot(data = test,
mapping = aes(y = reorder(Dieta, fit), x = fit)) +
geom_point() +
geom_errorbarh(mapping = aes(xmin = lwr, xmax = upr),
height = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.2f%s", fit, cld2)),
vjust = 0, nudge_y = 0.1) +
labs(x = expression("Taxa de crescimento" ~ (g ~ dia^{-1})),
y = "Dieta")
#
Ajuste do modelo logístico
#-----------------------------------------------------------------------
# Ajuste do modelo não linear logístico.
# Estrutura de dados para a `gnls()`.
tbg <- groupedData(Peso ~ Dias | Dieta,
data = filter(tb, Dias > 38),
order.groups = FALSE)
# Ajustes individuais primeiro.
n00 <- nlsList(Peso ~ SSlogis(Dias, A, C, S) | Dieta,
data = tbg)
coef(n00)
## A C S
## REF 2871.175 55.47320 17.89316
## SFA 2641.496 61.08485 20.55304
## SFAE 2529.890 60.17943 21.61182
## SSA 2480.912 53.43466 18.01786
## SSF 3530.413 79.94650 31.07577
## SSFA 2658.068 57.23572 18.29408
## SSFAE 2968.875 62.18390 22.03741
## SSFE 2899.738 68.72097 24.78058
## SSM 2769.104 59.59333 19.20956
## SSMA 2541.641 55.81140 17.96512
## SSME 2660.631 58.58033 18.65765
# plot(intervals(n00))
# Ajuste conjunto.
n0 <- gnls(Peso ~ A/(1 + exp(-4 * R * (Dias - C)/A)),
params = A + C + R ~ Dieta,
start = c(c(3000, rep(0, 10)),
c(60, rep(0, 10)),
c(20, rep(0, 10))),
data = tbg)
# Teste hipótese sobre os parâmetros.
anova(n0, type = "marginal")
## Denom. DF: 56
## numDF F-value p-value
## A.(Intercept) 1 2291.841 <.0001
## A.Dieta 10 4.312 2e-04
## C.(Intercept) 1 4558.302 <.0001
## C.Dieta 10 6.638 <.0001
## R.(Intercept) 1 13061.911 <.0001
## R.Dieta 10 83.285 <.0001
# Valores preditos para verificar qualidade do ajuste.
pred <- with(tbg,
crossing(Dieta = levels(Dieta),
Dias = seq(0, 100, length.out = 81)))
pred$Peso <- predict(n0, newdata = pred)
# Observa as curvas ajustadas.
ggplot(data = tbg,
mapping = aes(x = Dias, y = Peso)) +
facet_wrap(facets = ~Dieta) +
geom_point() +
geom_line(data = pred)
# Teste de hipótese para a taxa R.
Lvel <- cbind(0 * L, 0 * L, L)
Lvel %*% coef(n0)
## [,1]
## REF 40.11554
## SFA 32.13024
## SFAE 29.26512
## SSA 34.42297
## SSF 28.40166
## SSFA 36.32414
## SSFAE 33.67995
## SSFE 29.25415
## SSM 36.03810
## SSMA 35.36910
## SSME 35.65068
# Intervalo de confiança para as taxas.
# ci <- (stats::confint(glht(n0, linfct = Lvel),
# calpha = multcomp::univariate_calpha())$confint) %>%
# as.data.frame() %>%
# rownames_to_column(var = "Dieta") %>%
# rename(fit = Estimate)
# ci
# ATTENTION: esse é o truque porque objetos de classe `gnls` não tem
# `terms()`. Esse componente é necessário para a função
# `multcomp:::extr()` que é chamada dentro da `multcomp::cld()` que é
# chamada pela `apmc()`. Foi necessário abrir as funções para entender
# isso.
#
# str(m0$model)
# str(lm(Peso ~ Dieta * Dias, data = tbg)$model)
n0$model <- m0$model
# Comparações múltiplas para as taxas (R).
test <- apmc(X = Lvel, model = n0, focus = "Dieta", test = "fdr") %>%
arrange(fit) %>%
mutate(cld2 = ordered_cld(let = cld, means = fit))
test
## Dieta fit lwr upr cld cld2
## 1 SSF 28.40166 24.15302 32.65030 ef ef
## 2 SSFE 29.25415 27.59584 30.91246 f f
## 3 SFAE 29.26512 28.57039 29.95985 f f
## 4 SFA 32.13024 31.40275 32.85773 e e
## 5 SSFAE 33.67995 33.03281 34.32709 b d
## 6 SSA 34.42297 33.43340 35.41254 bc cd
## 7 SSMA 35.36910 34.52168 36.21651 cd bc
## 8 SSME 35.65068 34.92530 36.37606 cd bc
## 9 SSM 36.03810 35.33032 36.74587 d b
## 10 SSFA 36.32414 35.55069 37.09759 d b
## 11 REF 40.11554 39.24826 40.98282 a a
# Gráfico com as estimativas intervalares e resultado dos testes de hipótese.
ggplot(data = test, mapping = aes(y = reorder(Dieta, fit), x = fit)) +
geom_point() +
geom_errorbarh(mapping = aes(xmin = lwr, xmax = upr),
height = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.2f%s", fit, cld2)),
vjust = 0, nudge_y = 0.1) +
labs(x = expression("Taxa de crescimento" ~ (g ~ dia^{-1})),
y = "Dieta")
#