AnĂ¡lise de experimento fatorial com tratamentos adicionais
Andre Felipe Hermann & Walmes Marques Zeviani
2020-10-13
1 ImportaĂ§Ă£o e preparo
#-----------------------------------------------------------------------
# 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-08 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(lme4)
library(lmerTest)
library(emmeans)
# library(multcomp)
# library(multcompView)
library(tidyverse)
library(readxl)
as_letters <- function(x) {
xs <- strsplit(trimws(as.character(x)), "")
xs <- lapply(xs, as.integer)
m <- max(unlist(xs))
sapply(xs, function(i) paste(rev(letters[m - i + 1]), collapse = ""))
}# LĂª 3 abas.
tbs <- map(1:3, read_excel, path = "Dados para Walmes.xlsx")
# str(tbs)
# Faz a juntaĂ§Ă£o de todas as tabelas.
tb <- Reduce(inner_join, tbs)## Joining, by = c("Parcela", "Cultivar", "PopulaĂ§Ă£o")
## Joining, by = c("Parcela", "Cultivar", "PopulaĂ§Ă£o")# Encurta os nomes.
# dput(names(tb))
names(tb) <- c("ue", "cult", "pop", "rendimento",
"altura", "vagens")
# Modifica as variĂ¡veis.
tb <- tb %>%
mutate(cult = str_remove(cult, "Cultivar "),
pop = pop/1000,
reg = ifelse(is.element(pop, c(120, 160, 200)), "fat", "adic"))
# Delineamento fatorial irregular.
xtabs(~cult + pop, data = tb)## pop
## cult 120 160 200 240 300 340
## A 4 4 4 0 4 0
## B 4 4 4 0 4 0
## C 4 4 4 4 0 0
## D 4 4 4 0 4 0
## E 4 4 4 0 0 4# VisualizaĂ§Ă£o dos dados.
tb %>%
gather(key = "var", value = "val", rendimento:vagens) %>%
ggplot(mapping = aes(x = pop, y = val, color = reg)) +
facet_grid(facets = var ~ cult, scales = "free_y") +
geom_point() +
geom_smooth(mapping = aes(group = 1),
color = "black",
size = 0.25,
se = FALSE)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
2 AnĂ¡lise 1: conversĂ£o para fatorial 5x4
Nessa abordagem de anĂ¡lise, os nĂveis de populaĂ§Ă£o adicionais de cada cultivar receberĂ£o o mesmo rĂ³tulo. Ou seja, os nĂveis de populaĂ§Ă£o serĂ£o \(\{120, 160, 200, \text{Recomendada}\}\). Com isso, a estrutura de tratamentos do experimento passa a ser um fatorial regular \(5 \times 4\) para o qual a anĂ¡lise Ă© bem simples.
No caso de nĂ£o haver interaĂ§Ă£o, faz-se o estudo dos efeitos principais. Quando houver interaĂ§Ă£o, pode-se estudar o efeito de populaĂ§Ă£o dentro de cada cultivar por meio de comparaĂ§Ă£o mĂºltipla de mĂ©dias ou contrastes planejados contra o nĂvel adicional.
Uma questĂ£o que esse modelo deixa em aberto Ă© que, havendo interaĂ§Ă£o, nĂ£o Ă© possĂvel discrinar o quanto dela se deve a parte regular, ou seja, proveniente do cruzamento com os nĂveis comuns de populaĂ§Ă£o (120, 160 e 200) com cultivar (\((3 - 1)(5 - 1) = 8\) graus de liberdade), e o quanto se deve a parte adicional. Para separar isso, deve-se representar os fatores experimentais de outra forma para o modelo.
# Converte variĂ¡veis para fator.
tb <- tb %>%
mutate(Cult = factor(cult),
fat = is.element(pop, c(120, 160, 200)),
Pop = factor(ifelse(fat, pop, "Recomendada")),
Trat = interaction(Cult, pop, drop = TRUE, sep = ":"),
Adic = factor(ifelse(!fat, as.character(Trat), "fat"))) %>%
group_by(Cult, pop) %>%
mutate(rept = 1:n()) %>%
ungroup()
str(tb)## tibble [80 Ă— 13] (S3: tbl_df/tbl/data.frame)
## $ ue : num [1:80] 1 2 3 4 5 6 7 8 9 10 ...
## $ cult : chr [1:80] "A" "A" "A" "A" ...
## $ pop : num [1:80] 120 120 120 120 160 160 160 160 200 200 ...
## $ rendimento: num [1:80] 4889 5305 5131 5148 4962 ...
## $ altura : num [1:80] 91.8 100.4 92.4 97.2 100.2 ...
## $ vagens : num [1:80] 56.2 53 64.2 55.4 32.8 34.6 30.4 30.2 31.6 23.2 ...
## $ reg : chr [1:80] "fat" "fat" "fat" "fat" ...
## $ Cult : Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ fat : logi [1:80] TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ Pop : Factor w/ 4 levels "120","160","200",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ Trat : Factor w/ 20 levels "A:120","B:120",..: 1 1 1 1 6 6 6 6 11 11 ...
## $ Adic : Factor w/ 6 levels "A:300","B:300",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ rept : int [1:80] 1 2 3 4 1 2 3 4 1 2 ...# ATTENTION: para criar a variĂ¡vel `rept` estĂ¡ sendo assumido que a
# as repetições estĂ£o em ordem.
# Combinações experimentais.
xtabs(~Cult + Pop, data = tb)## Pop
## Cult 120 160 200 Recomendada
## A 4 4 4 4
## B 4 4 4 4
## C 4 4 4 4
## D 4 4 4 4
## E 4 4 4 4xtabs(~Adic, data = tb)## Adic
## A:300 B:300 C:240 D:300 E:340 fat
## 4 4 4 4 4 60xtabs(~Trat, data = tb)## Trat
## A:120 B:120 C:120 D:120 E:120 A:160 B:160 C:160 D:160 E:160 A:200 B:200
## 4 4 4 4 4 4 4 4 4 4 4 4
## C:200 D:200 E:200 C:240 A:300 B:300 D:300 E:340
## 4 4 4 4 4 4 4 42.1 Rendimento
# AnĂ¡lise para rendimento.
mm_r <- lmer(rendimento ~ Cult * Pop + (1 | Cult:rept),
data = tb)
anova(mm_r)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 1152296 288074 4 15 1.6672 0.2097
## Pop 785964 261988 3 45 1.5162 0.2232
## Cult:Pop 2931371 244281 12 45 1.4137 0.1952# Faz a inspeĂ§Ă£o dos pressupostos.
mm <- mm_r
par(mfrow = c(2, 2))
plot(residuals(mm) ~ fitted(mm))
plot(residuals(mm) ~ jitter(as.integer(tb$Cult)))
plot(residuals(mm) ~ jitter(as.integer(tb$Pop)))
qqnorm(residuals(mm))
qqline(residuals(mm), lty = 2)
layout(1)
# ComparaĂ§Ă£o mĂºltipla das mĂ©dias marginais de Cultivar.
tb_cult <- emmeans(mm, specs = ~Cult) %>%
multcomp::cld() %>%
as.data.frame() %>%
mutate(cld = as_letters(.group))## NOTE: Results may be misleading due to involvement in interactionstb_cult## Cult emmean SE df lower.CL upper.CL .group cld
## 1 D 4947.051 108.2826 15 4716.253 5177.850 1 a
## 2 A 5156.678 108.2826 15 4925.879 5387.477 1 a
## 3 E 5234.099 108.2826 15 5003.301 5464.898 1 a
## 4 C 5250.771 108.2826 15 5019.972 5481.570 1 a
## 5 B 5304.324 108.2826 15 5073.525 5535.122 1 agg1 <-
ggplot(data = tb_cult,
mapping = aes(x = reorder(Cult, emmean), y = emmean)) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f%s", emmean, cld)),
hjust = 0, nudge_x = 0.05) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL),
width = 0.05) +
labs(x = "Cultivares",
y = "Rendimento")
# ComparaĂ§Ă£o mĂºltipla das mĂ©dias marginais de PopulaĂ§Ă£o.
tb_pop <- emmeans(mm, specs = ~Pop) %>%
multcomp::cld() %>%
as.data.frame() %>%
mutate(cld = as_letters(.group))## NOTE: Results may be misleading due to involvement in interactionstb_pop## Pop emmean SE df lower.CL upper.CL .group cld
## 1 160 5061.646 93.94028 59.92091 4873.732 5249.560 1 a
## 2 120 5103.188 93.94028 59.92091 4915.274 5291.101 1 a
## 3 Recomendada 5247.896 93.94028 59.92091 5059.983 5435.810 1 a
## 4 200 5301.608 93.94028 59.92091 5113.694 5489.522 1 agg2 <-
ggplot(data = tb_pop,
mapping = aes(x = Pop, y = emmean)) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f%s", emmean, cld)),
hjust = 0, nudge_x = 0.05) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL),
width = 0.05) +
labs(x = "PopulaĂ§Ă£o (/1000)",
y = "Rendimento")
gridExtra::grid.arrange(gg1, gg2, ncol = 1)
2.2 Altura de plantas
# AnĂ¡lise para altura.
mm_a <- lmer(altura ~ Cult * Pop + (1 | Cult:rept),
data = tb)## boundary (singular) fit: see ?isSingularanova(mm_a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 10799.5 2699.87 4 60 267.6393 < 2.2e-16 ***
## Pop 1113.1 371.04 3 60 36.7813 1.29e-13 ***
## Cult:Pop 104.7 8.73 12 60 0.8652 0.5854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Faz a inspeĂ§Ă£o dos pressupostos.
mm <- mm_a
par(mfrow = c(2, 2))
plot(residuals(mm) ~ fitted(mm))
plot(residuals(mm) ~ jitter(as.integer(tb$Cult)))
plot(residuals(mm) ~ jitter(as.integer(tb$Pop)))
qqnorm(residuals(mm))
qqline(residuals(mm), lty = 2)
layout(1)
# ComparaĂ§Ă£o mĂºltipla das mĂ©dias marginais de Cultivar.
tb_cult <- emmeans(mm, specs = ~Cult) %>%
multcomp::cld() %>%
as.data.frame() %>%
mutate(cld = as_letters(.group))## NOTE: Results may be misleading due to involvement in interactionstb_cult## Cult emmean SE df lower.CL upper.CL .group cld
## 1 E 88.13725 0.7940288 15 86.44482 89.82968 1 d
## 2 B 100.98750 0.7940288 15 99.29507 102.67993 2 c
## 3 A 103.00000 0.7940288 15 101.30757 104.69243 2 c
## 4 C 108.83750 0.7940288 15 107.14507 110.52993 3 b
## 5 D 123.86250 0.7940288 15 122.17007 125.55493 4 agg1 <-
ggplot(data = tb_cult,
mapping = aes(x = reorder(Cult, emmean), y = emmean)) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f%s", emmean, cld)),
hjust = 0, nudge_x = 0.05) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL),
width = 0.05) +
labs(x = "Cultivares",
y = "Altura de plantas")
# ComparaĂ§Ă£o mĂºltipla das mĂ©dias marginais de PopulaĂ§Ă£o.
tb_pop <- emmeans(mm, specs = ~Pop) %>%
multcomp::cld() %>%
as.data.frame() %>%
mutate(cld = as_letters(.group))## NOTE: Results may be misleading due to involvement in interactionstb_pop## Pop emmean SE df lower.CL upper.CL .group cld
## 1 120 100.1600 0.7102009 60 98.73939 101.5806 1 c
## 2 160 103.6100 0.7102009 60 102.18939 105.0306 2 b
## 3 200 105.6200 0.7102009 60 104.19939 107.0406 2 b
## 4 Recomendada 110.4698 0.7102009 60 109.04919 111.8904 3 agg2 <-
ggplot(data = tb_pop,
mapping = aes(x = Pop, y = emmean)) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f%s", emmean, cld)),
hjust = 0, nudge_x = 0.05) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL),
width = 0.05) +
labs(x = "PopulaĂ§Ă£o (/1000)",
y = "Altura de plantas")
gridExtra::grid.arrange(gg1, gg2, ncol = 1)
2.3 NĂºmero de vagens
A anĂ¡lise foi feita com o logarĂtmo para que nĂ£o ocorresse afastamento dos pressupostos.
# AnĂ¡lise para nĂºmero de vagens.
mm_v <- lmer(log(vagens) ~ Cult * Pop + (1 | Cult:rept),
data = tb)
anova(mm_v)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 7.313 1.8282 4 15 17.0358 1.888e-05 ***
## Pop 44.112 14.7040 3 45 137.0143 < 2.2e-16 ***
## Cult:Pop 8.954 0.7462 12 45 6.9529 6.558e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Faz a inspeĂ§Ă£o dos pressupostos.
mm <- mm_v
par(mfrow = c(2, 2))
plot(residuals(mm) ~ fitted(mm))
plot(residuals(mm) ~ jitter(as.integer(tb$Cult)))
plot(residuals(mm) ~ jitter(as.integer(tb$Pop)))
qqnorm(residuals(mm))
qqline(residuals(mm), lty = 2)
layout(1)
# ComparaĂ§Ă£o mĂºltipla das mĂ©dias marginais de Cultivar.
tb_popcult <- emmeans(mm, specs = ~Pop | Cult) %>%
multcomp::cld() %>%
as.data.frame() %>%
group_by(Cult) %>%
mutate(cld = as_letters(.group)) %>%
ungroup() %>%
mutate_at(c("emmean", "lower.CL", "upper.CL"), "exp")
tb_popcult## # A tibble: 20 x 9
## Pop Cult emmean SE df lower.CL upper.CL .group cld
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 Recomendada A 10.3 0.188 51.3 7.04 15.0 " 1 " c
## 2 200 A 24.7 0.188 51.3 17.0 36.0 " 2 " b
## 3 160 A 31.9 0.188 51.3 21.9 46.6 " 23" ab
## 4 120 A 57.1 0.188 51.3 39.1 83.2 " 3" a
## 5 Recomendada B 1.26 0.188 51.3 0.866 1.84 " 1 " d
## 6 200 B 7.87 0.188 51.3 5.40 11.5 " 2 " c
## 7 160 B 19.0 0.188 51.3 13.0 27.6 " 3 " b
## 8 120 B 57.5 0.188 51.3 39.5 83.8 " 4" a
## 9 Recomendada C 17.5 0.188 51.3 12.0 25.6 " 1 " c
## 10 200 C 31.0 0.188 51.3 21.3 45.2 " 12 " bc
## 11 160 C 49.2 0.188 51.3 33.8 71.7 " 23" ab
## 12 120 C 67.6 0.188 51.3 46.4 98.6 " 3" a
## 13 Recomendada D 12.8 0.188 51.3 8.79 18.7 " 1 " c
## 14 200 D 30.6 0.188 51.3 21.0 44.6 " 2 " b
## 15 160 D 39.7 0.188 51.3 27.2 57.9 " 23" ab
## 16 120 D 59.2 0.188 51.3 40.6 86.2 " 3" a
## 17 Recomendada E 11.7 0.188 51.3 8.04 17.1 " 1 " c
## 18 200 E 24.2 0.188 51.3 16.6 35.3 " 2 " b
## 19 160 E 44.6 0.188 51.3 30.6 64.9 " 23" ab
## 20 120 E 63.1 0.188 51.3 43.3 91.9 " 3" agg3 <-
ggplot(data = tb_popcult,
mapping = aes(x = Pop, y = emmean)) +
facet_wrap(facets = ~Cult, nrow = 3) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f%s", emmean, cld)),
hjust = 0, nudge_x = 0.05) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL),
width = 0.05) +
expand_limits(x = c(NA, 5)) +
labs(x = "PopulaĂ§Ă£o",
y = "Altura de plantas")
gg3
3 AnĂ¡lise 2: conversĂ£o para fatorial 5x3 + 5
3.1 Rendimento
# AnĂ¡lise para rendimento.
mm_r <- lmer(rendimento ~ Adic + Cult * Pop + (1 | Cult:rept),
data = tb)## fixed-effect model matrix is rank deficient so dropping 5 columns / coefficientsanova(mm_r)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Adic 872764 174553 5 45.000 1.0102 0.42285
## Cult 1530288 382572 4 24.844 2.2140 0.09662 .
## Pop 657854 328927 2 45.000 1.9036 0.16084
## Cult:Pop 1875415 234427 8 45.000 1.3567 0.24142
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2 Altura de plantas
# AnĂ¡lise para altura.
mm_a <- lmer(altura ~ Adic + Cult * Pop + (1 | Cult:rept),
data = tb)## fixed-effect model matrix is rank deficient so dropping 5 columns / coefficients## boundary (singular) fit: see ?isSingularanova(mm_a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Adic 1133.3 226.65 5 60 22.468 1.270e-12 ***
## Cult 7972.9 1993.24 4 60 197.591 < 2.2e-16 ***
## Pop 305.0 152.51 2 60 15.119 4.819e-06 ***
## Cult:Pop 59.9 7.48 8 60 0.742 0.6543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.3 NĂºmero de vagens
# AnĂ¡lise para nĂºmero de vagens.
mm_v <- lmer(log(vagens) ~ Adic + Cult * Pop + (1 | Cult:rept),
data = tb)## fixed-effect model matrix is rank deficient so dropping 5 columns / coefficientsanova(mm_v)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Adic 49.049 9.8098 5 45.000 91.4088 < 2.2e-16 ***
## Cult 2.631 0.6578 4 19.627 6.1295 0.002264 **
## Pop 10.907 5.4534 2 45.000 50.8161 2.865e-12 ***
## Cult:Pop 2.490 0.3112 8 45.000 2.9002 0.010745 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14 AnĂ¡lise 3: uso de regressĂ£o polinomial
Nessa abordagem de anĂ¡lise, o efeito de populaĂ§Ă£o Ă© representado com termos de funções polinomiais. Pela anĂ¡lise exploratĂ³ria, o uso de polinĂ´mio de grau de 2 Ă© suficiente para acomodar o efeito de populaĂ§Ă£o. Se os termos de maior ordem nĂ£o forem significativos, eles serĂ£o abandonados e um modelo mais simples serĂ¡ ajustado. As curvas podem ser comparadas por meio de testes de hipĂ³tese sobre os parĂ¢metros.
4.1 Rendimento
# Reduz a escala.
tb <- tb %>%
mutate(popul = pop/100)
# AnĂ¡lise para rendimento.
mm_r <- lmer(rendimento ~ Cult * (popul + I(popul^2)) + (1 | Cult:rept),
data = tb)## boundary (singular) fit: see ?isSingularanova(mm_r)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 223504 55876 4 65 0.2939 0.8809
## popul 71669 71669 1 65 0.3769 0.5414
## I(popul^2) 30215 30215 1 65 0.1589 0.6915
## Cult:popul 366003 91501 4 65 0.4812 0.7494
## Cult:I(popul^2) 469427 117357 4 65 0.6172 0.6518mm_r <- update(mm_r, . ~ Cult * popul + (1 | Cult:rept))
anova(mm_r)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 1380447 345112 4 65.082 1.8847 0.12369
## popul 532417 532417 1 55.000 2.9076 0.09381 .
## Cult:popul 1174441 293610 4 55.000 1.6034 0.18646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Coeficientes estimados.
summary(mm_r)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rendimento ~ Cult + popul + (1 | Cult:rept) + Cult:popul
## Data: tb
##
## REML criterion at convergence: 1070.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.54885 -0.36213 -0.02188 0.43638 3.12344
##
## Random effects:
## Groups Name Variance Std.Dev.
## Cult:rept (Intercept) 1123 33.5
## Residual 183112 427.9
## Number of obs: 80, groups: Cult:rept, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4907.303 330.108 65.567 14.866 <2e-16 ***
## CultB 767.396 466.843 65.567 1.644 0.105
## CultC -593.641 553.260 62.917 -1.073 0.287
## CultD -206.469 466.843 65.567 -0.442 0.660
## CultE 169.192 436.639 66.699 0.387 0.700
## popul 127.885 159.919 55.000 0.800 0.427
## CultB:popul -317.821 226.160 55.000 -1.405 0.166
## CultC:popul 392.731 287.744 55.000 1.365 0.178
## CultD:popul -1.619 226.160 55.000 -0.007 0.994
## CultE:popul -51.004 205.477 55.000 -0.248 0.805
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CultB CultC CultD CultE popul CltB:p CltC:p
## CultB -0.707
## CultC -0.597 0.422
## CultD -0.707 0.500 0.422
## CultE -0.756 0.535 0.451 0.535
## popul -0.945 0.668 0.564 0.668 0.714
## CultB:popul 0.668 -0.945 -0.399 -0.472 -0.505 -0.707
## CultC:popul 0.525 -0.371 -0.960 -0.371 -0.397 -0.556 0.393
## CultD:popul 0.668 -0.472 -0.399 -0.945 -0.505 -0.707 0.500 0.393
## CultE:popul 0.735 -0.520 -0.439 -0.520 -0.936 -0.778 0.550 0.433
## CltD:p
## CultB
## CultC
## CultD
## CultE
## popul
## CultB:popul
## CultC:popul
## CultD:popul
## CultE:popul 0.550# Valores preditos com intervalos de confiança.
popul_seq <- extendrange(tb$popul, f = 0.1)
popul_seq <- seq(popul_seq[1], popul_seq[2], length.out = 31)
tb_curves <- emmeans(mm_r,
specs = ~Cult * popul,
at = list(popul = popul_seq)) %>%
as.data.frame()
# Filtra para valores compatĂveis com a regiĂ£o experimental explorada.
tb_curves <- tb_curves %>%
group_by(Cult) %>%
mutate(popul = {
m <- 1.1 * max(tb$popul[tb$Cult == Cult[1]])
ifelse(popul > m, NA, popul)
}) %>%
filter(!is.na(popul)) %>%
mutate(pop = popul * 100)
ggplot(data = tb_curves,
mapping = aes(x = pop, y = emmean)) +
facet_wrap(facets = ~Cult) +
geom_ribbon(mapping = aes(ymin = lower.CL, ymax = upper.CL),
fill = "gray75") +
geom_line() +
geom_point(data = tb,
mapping = aes(x = pop, y = rendimento)) +
labs(x = "PopulaĂ§Ă£o (/1000)",
y = "Rendimento")
4.2 Altura de plantas
# AnĂ¡lise para altura de plantas.
mm_a <- lmer(altura ~ Cult * (popul + I(popul^2)) + (1 | Cult:rept),
data = tb)
anova(mm_a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 105.064 26.266 4 50.643 2.8329 0.03388 *
## popul 51.170 51.170 1 50.000 5.5189 0.02280 *
## I(popul^2) 13.055 13.055 1 50.000 1.4080 0.24100
## Cult:popul 32.800 8.200 4 50.000 0.8844 0.48011
## Cult:I(popul^2) 23.611 5.903 4 50.000 0.6366 0.63878
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mm_a <- update(mm_a, . ~ Cult * popul + (1 | Cult:rept))
anova(mm_a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 1021.06 255.26 4 65.515 28.3703 1.094e-13 ***
## popul 1061.75 1061.75 1 55.000 118.0033 2.648e-15 ***
## Cult:popul 103.67 25.92 4 55.000 2.8806 0.03083 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Coeficientes estimados.
summary(mm_a)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: altura ~ Cult + popul + (1 | Cult:rept) + Cult:popul
## Data: tb
##
## REML criterion at convergence: 377
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2835 -0.3894 0.0572 0.4828 2.2248
##
## Random effects:
## Groups Name Variance Std.Dev.
## Cult:rept (Intercept) 0.1999 0.4471
## Residual 8.9976 2.9996
## Number of obs: 80, groups: Cult:rept, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 87.3891 2.3218 66.0238 37.639 < 2e-16 ***
## CultB 5.0821 3.2835 66.0238 1.548 0.12646
## CultC 7.9259 3.8876 63.3220 2.039 0.04565 *
## CultD 24.8687 3.2835 66.0238 7.574 1.53e-10 ***
## CultE -6.8638 3.0725 67.1482 -2.234 0.02882 *
## popul 8.0056 1.1210 55.0000 7.141 2.21e-09 ***
## CultB:popul -3.6383 1.5853 55.0000 -2.295 0.02558 *
## CultC:popul -0.4931 2.0170 55.0000 -0.244 0.80778
## CultD:popul -2.0545 1.5853 55.0000 -1.296 0.20042
## CultE:popul -4.2924 1.4404 55.0000 -2.980 0.00428 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CultB CultC CultD CultE popul CltB:p CltC:p
## CultB -0.707
## CultC -0.597 0.422
## CultD -0.707 0.500 0.422
## CultE -0.756 0.534 0.451 0.534
## popul -0.941 0.666 0.562 0.666 0.711
## CultB:popul 0.666 -0.941 -0.398 -0.471 -0.503 -0.707
## CultC:popul 0.523 -0.370 -0.958 -0.370 -0.395 -0.556 0.393
## CultD:popul 0.666 -0.471 -0.398 -0.941 -0.503 -0.707 0.500 0.393
## CultE:popul 0.733 -0.518 -0.438 -0.518 -0.933 -0.778 0.550 0.433
## CltD:p
## CultB
## CultC
## CultD
## CultE
## popul
## CultB:popul
## CultC:popul
## CultD:popul
## CultE:popul 0.550# Valores preditos com intervalos de confiança.
popul_seq <- extendrange(tb$popul, f = 0.1)
popul_seq <- seq(popul_seq[1], popul_seq[2], length.out = 31)
tb_curves <- emmeans(mm_a,
specs = ~Cult * popul,
at = list(popul = popul_seq)) %>%
as.data.frame()
# Filtra para valores compatĂveis com a regiĂ£o experimental explorada.
tb_curves <- tb_curves %>%
group_by(Cult) %>%
mutate(popul = {
m <- 1.1 * max(tb$popul[tb$Cult == Cult[1]])
ifelse(popul > m, NA, popul)
}) %>%
filter(!is.na(popul)) %>%
mutate(pop = popul * 100)
ggplot(data = tb_curves,
mapping = aes(x = pop, y = emmean)) +
facet_wrap(facets = ~Cult) +
geom_ribbon(mapping = aes(ymin = lower.CL, ymax = upper.CL),
fill = "gray75") +
geom_line() +
geom_point(data = tb,
mapping = aes(x = pop, y = altura)) +
labs(x = "PopulaĂ§Ă£o (/1000)",
y = "Altura de plantas")
4.3 NĂºmero de vagens
# AnĂ¡lise para nĂºmero de vagens.
mm_v <- lmer(log(vagens) ~ Cult * (popul + I(popul^2)) + (1 | Cult:rept),
data = tb)
anova(mm_v)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 0.61388 0.15347 4 51.457 1.5392 0.204723
## popul 0.91443 0.91443 1 50.000 9.1710 0.003883 **
## I(popul^2) 0.03951 0.03951 1 50.000 0.3963 0.531884
## Cult:popul 0.62279 0.15570 4 50.000 1.5615 0.199067
## Cult:I(popul^2) 0.30452 0.07613 4 50.000 0.7635 0.554009
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mm_v <- update(mm_v, . ~ Cult * popul + (1 | Cult:rept))
anova(mm_v)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cult 2.192 0.548 4 69.388 5.355 0.0008117 ***
## popul 39.675 39.675 1 55.000 387.755 < 2.2e-16 ***
## Cult:popul 8.932 2.233 4 55.000 21.824 7.929e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Coeficientes estimados.
summary(mm_v)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(vagens) ~ Cult + popul + (1 | Cult:rept) + Cult:popul
## Data: tb
##
## REML criterion at convergence: 75.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.95452 -0.43076 0.09476 0.55495 1.97380
##
## Random effects:
## Groups Name Variance Std.Dev.
## Cult:rept (Intercept) 0.03482 0.1866
## Residual 0.10232 0.3199
## Number of obs: 80, groups: Cult:rept, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.04736 0.26351 69.99582 19.154 < 2e-16 ***
## CultB 1.32924 0.37267 69.99582 3.567 0.000657 ***
## CultC 0.58460 0.43375 68.81843 1.348 0.182147
## CultD 0.01940 0.37267 69.99582 0.052 0.958636
## CultE -0.10836 0.35161 69.62804 -0.308 0.758869
## popul -0.91598 0.11954 55.00000 -7.662 3.10e-10 ***
## CultB:popul -1.16304 0.16906 55.00000 -6.879 5.94e-09 ***
## CultC:popul -0.21160 0.21509 55.00000 -0.984 0.329542
## CultD:popul 0.07841 0.16906 55.00000 0.464 0.644617
## CultE:popul 0.16393 0.15360 55.00000 1.067 0.290508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CultB CultC CultD CultE popul CltB:p CltC:p
## CultB -0.707
## CultC -0.608 0.430
## CultD -0.707 0.500 0.430
## CultE -0.749 0.530 0.455 0.530
## popul -0.885 0.626 0.537 0.626 0.663
## CultB:popul 0.626 -0.885 -0.380 -0.442 -0.469 -0.707
## CultC:popul 0.492 -0.348 -0.916 -0.348 -0.368 -0.556 0.393
## CultD:popul 0.626 -0.442 -0.380 -0.885 -0.469 -0.707 0.500 0.393
## CultE:popul 0.688 -0.487 -0.418 -0.487 -0.869 -0.778 0.550 0.433
## CltD:p
## CultB
## CultC
## CultD
## CultE
## popul
## CultB:popul
## CultC:popul
## CultD:popul
## CultE:popul 0.550# Valores preditos com intervalos de confiança.
popul_seq <- extendrange(tb$popul, f = 0.1)
popul_seq <- seq(popul_seq[1], popul_seq[2], length.out = 31)
tb_curves <- emmeans(mm_v,
specs = ~Cult * popul,
at = list(popul = popul_seq)) %>%
as.data.frame()
# Filtra para valores compatĂveis com a regiĂ£o experimental explorada.
tb_curves <- tb_curves %>%
group_by(Cult) %>%
mutate(popul = {
m <- 1.1 * max(tb$popul[tb$Cult == Cult[1]])
ifelse(popul > m, NA, popul)
}) %>%
filter(!is.na(popul)) %>%
mutate(pop = popul * 100)
ggplot(data = tb_curves,
mapping = aes(x = pop, y = emmean)) +
facet_wrap(facets = ~Cult) +
geom_ribbon(mapping = aes(ymin = lower.CL, ymax = upper.CL),
fill = "gray75") +
geom_line() +
geom_point(data = tb,
mapping = aes(x = pop, y = log(vagens))) +
labs(x = "PopulaĂ§Ă£o (/1000)",
y = "NĂºmero de vagens (log)")