Preparação e análise exploratória

library(nlme)      # Modelos de efeitos mistos.

# Para não exibir o bloco de correlação entre estimativas dos efeitos
# fixos na saída do `summary()`. ATTENTION: tem que vir antes da
# `lmerTest`.
assignInNamespace("print.correlation",
                  function(x, title) return(),
                  ns = "nlme")

library(lmerTest)  # Testes de hipótese para classe lme e similares.
library(emmeans)   # Médias ajustadas.
library(tidyverse) # Manipulação e visualização de dados.

# Para apresentar letras ordenadas com o "a" vínculado a maior média.
numbers_to_letters <- function(.group) {
    list_num <- trimws(.group) %>%
        strsplit(split = "")
    n <- max(as.integer(unlist(list_num)))
    list_num %>%
        lapply(FUN = function(x) sort(letters[n:1][as.integer(x)])) %>%
        sapply(FUN = paste, collapse = "")
}
# Leitura da tabela de dados.
url <- "Dados níveis de erva-mate - PARA ENVIAR ESTATISTICA.xlsx"
tb <- gdata::read.xls(xls = url,
                      sheet = 1,
                      na.strings = ".",
                      check.names = FALSE)

# Adequação dos nomes das variáveis.
names(tb) <- names(tb) %>%
    str_to_lower() %>%
    str_replace_all("[[:punct:]]", "")
str(tb)
## 'data.frame':    120 obs. of  8 variables:
##  $ mate     : num  0 0 0 0 0 0 1.5 1.5 1.5 1.5 ...
##  $ rep      : int  1 2 3 4 5 6 1 2 3 4 ...
##  $ per      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ consumo  : num  98.8 98.3 104.1 97.5 97.5 ...
##  $ postura  : num  94 94 94 93.5 NA ...
##  $ massaovos: num  55.6 50.2 53.7 55.6 43.5 ...
##  $ caduzia  : num  7.5 7.46 7.91 7.45 9.52 ...
##  $ camovos  : num  1.78 1.96 1.94 1.75 2.24 ...
# Faz ordenação dos registros.
tb <- tb %>%
    arrange(mate, rep, per)

# Versão empilhada ou verticalizada das respostas.
tb_long <- tb %>%
    gather(key = "variable", value = "value", consumo:camovos)
str(tb_long)
## 'data.frame':    600 obs. of  5 variables:
##  $ mate    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ rep     : int  1 1 1 1 2 2 2 2 3 3 ...
##  $ per     : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ variable: chr  "consumo" "consumo" "consumo" "consumo" ...
##  $ value   : num  98.8 81.8 83.7 83.3 98.3 ...
ggplot(data = tb_long,
       mapping = aes(x = mate, y = value)) +
    facet_grid(facets = variable ~ per, scale = "free_y") +
    # geom_point() +
    geom_jitter(width = 0.1, height = 0) +
    stat_summary(mapping = aes(group = 1),
                 fun.y = "mean",
                 geom = "line") +
    labs(x = "Concentração de mate",
         y = "Valor de cada resposta")
## Warning: Removed 4 rows containing non-finite values (stat_summary).
## Warning: Removed 4 rows containing missing values (geom_point).

ggplot(data = tb_long,
       mapping = aes(x = per, y = value)) +
    facet_grid(facets = variable ~ mate, scale = "free_y") +
    # geom_point() +
    geom_jitter(width = 0.1, height = 0) +
    stat_summary(mapping = aes(group = 1),
                 fun.y = "mean",
                 geom = "line") +
    labs(x = "Período de postura",
         y = "Valor de cada resposta")
## Warning: Removed 4 rows containing non-finite values (stat_summary).

## Warning: Removed 4 rows containing missing values (geom_point).

# ATTENTION: `massaovos` deve ser o valor total ou médio porque está
# repetindo no período.

# Conta os valores presentes (tamanho de amostra).
tb %>%
    group_by(mate, per) %>%
    summarise_at(vars(consumo:camovos),
                 function(x) sum(is.finite(x)))
## # A tibble: 20 x 7
## # Groups:   mate [5]
##     mate   per consumo postura massaovos caduzia camovos
##    <dbl> <int>   <int>   <int>     <int>   <int>   <int>
##  1   0       1       6       5         6       6       6
##  2   0       2       6       5         6       6       6
##  3   0       3       6       5         6       6       6
##  4   0       4       6       5         6       6       6
##  5   1.5     1       6       6         6       6       6
##  6   1.5     2       6       6         6       6       6
##  7   1.5     3       6       6         6       6       6
##  8   1.5     4       6       6         6       6       6
##  9   3       1       6       6         6       6       6
## 10   3       2       6       6         6       6       6
## 11   3       3       6       6         6       6       6
## 12   3       4       6       6         6       6       6
## 13   4.5     1       6       6         6       6       6
## 14   4.5     2       6       6         6       6       6
## 15   4.5     3       6       6         6       6       6
## 16   4.5     4       6       6         6       6       6
## 17   6       1       6       6         6       6       6
## 18   6       2       6       6         6       6       6
## 19   6       3       6       6         6       6       6
## 20   6       4       6       6         6       6       6
# Conta os valores ausentes.
tb %>%
    group_by(mate, per) %>%
    summarise_at(vars(consumo:camovos),
                 function(x) sum(is.na(x)))
## # A tibble: 20 x 7
## # Groups:   mate [5]
##     mate   per consumo postura massaovos caduzia camovos
##    <dbl> <int>   <int>   <int>     <int>   <int>   <int>
##  1   0       1       0       1         0       0       0
##  2   0       2       0       1         0       0       0
##  3   0       3       0       1         0       0       0
##  4   0       4       0       1         0       0       0
##  5   1.5     1       0       0         0       0       0
##  6   1.5     2       0       0         0       0       0
##  7   1.5     3       0       0         0       0       0
##  8   1.5     4       0       0         0       0       0
##  9   3       1       0       0         0       0       0
## 10   3       2       0       0         0       0       0
## 11   3       3       0       0         0       0       0
## 12   3       4       0       0         0       0       0
## 13   4.5     1       0       0         0       0       0
## 14   4.5     2       0       0         0       0       0
## 15   4.5     3       0       0         0       0       0
## 16   4.5     4       0       0         0       0       0
## 17   6       1       0       0         0       0       0
## 18   6       2       0       0         0       0       0
## 19   6       3       0       0         0       0       0
## 20   6       4       0       0         0       0       0
# Cria variável para representar a parcela e cria cópias como fator.
tb <- tb %>%
    mutate(parcela = interaction(mate, rep, drop = TRUE),
           Mate = factor(mate),
           Per = factor(per))

#

Especificação do modelo

O modelo para análise das variáveis medidas ao nível de animais, ou seja, nas unidades experimentais que são as gaiolas é \[ \begin{aligned} y_{ijk}| ijk &\sim \text{Normal}(\mu_{ij}, \sigma^2) \\ \mu_{ij} &= \mu + \alpha_i + e_{ik} + \beta_j + \gamma_{ij} \\ e_{ik} &\sim \text{Normal}(0, \sigma_{e}^2) \end{aligned} \]

Colocado de outra forma, esse modelo pode ser escrito como soma da parte determinística e aleatória que terá dois estratos: o de parcela e subparcela. Na parcela foram casualizados os níveis de mate. E na subparcela assume-se que foram casualizados os níveis de tempo. Todavia, de fato não foi feita essa casualização porque trata-se de um experimento longitudinal ou de medidas repetidas, também chamado de parcelas subdivididas no tempo.

O modelo escrito como soma de termos determinístico e de erro é \[ \begin{aligned} y_{ijk} &= \mu + \tau_i + \theta_j + \gamma_{ij} + \varepsilon_{ij} + \epsilon_{ijk}\\ \varepsilon_{ij} &\sim \text{Normal}(0, \sigma_\varepsilon^2) \\ \epsilon_{ijk} &\sim \text{Normal}(0, \sigma_\epsilon^2), \\ \end{aligned} \] em que \(\varepsilon_{ij}\) e \(\epsilon_{ijk}\) são os erros ao nível de parcela e subparcela (residual), respectivamente. Estes erros são independentes entre si.

Nas seções a seguir será feito a análise para uma das variáveis resposta. Essa análise será feita usando o método ANOVA, o modelo misto (REML) e por último o modelo misto com estrutura de dependência entre observações ao longo do tempo.

Análise da variável consumo

Método ANOVA

O método ANOVA é muito utilizado. Todavia, para experimentos desbalanceados ocorrem complicações para a obtenção das somas de quadrados. Por isso o método REML é preferível. A título de ilustração, o código abaixo ajusta o modelo e exibe o quadro de análise de variância de dois estratos: com erro ao nível de parcela e erro ao nível de subparcela.

# Modelo de um estrato apenas para conferir pressupostos.
m0 <- lm(terms(consumo ~ Mate + parcela + Per + Mate:Per,
               keep.order = TRUE),
         data = tb)

# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Algumas estatísticas F estão erradas.
anova(m0)
## Analysis of Variance Table
## 
## Response: consumo
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Mate       4  215.76   53.94  2.2021   0.07679 .  
## parcela   25  805.14   32.21  1.3148   0.18241    
## Per        3 2790.02  930.01 37.9682 4.935e-15 ***
## Mate:Per  12 1831.95  152.66  6.2326 1.743e-07 ***
## Residuals 75 1837.08   24.49                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo de 2 estratos.
m0 <- aov(consumo ~ Mate + Error(parcela) + Per + Mate:Per,
          data = tb)

# Quadro de análise de variância de 2 estratos.
summary(m0)
## 
## Error: parcela
##           Df Sum Sq Mean Sq F value Pr(>F)
## Mate       4  215.8   53.94   1.675  0.187
## Residuals 25  805.1   32.21               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Per        3   2790   930.0  37.968 4.94e-15 ***
## Mate:Per  12   1832   152.7   6.233 1.74e-07 ***
## Residuals 75   1837    24.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#

Modelos de efeitos aleatórios REML

O modelo ajustado com as instruções a seguir corresponde ao mesmo modelo usado no método ANOVA. Para dados balanceados, os resultados serão os mesmo, exceto somente que o método REML só produz estimativas dentro do espaço paramétrico, ou seja, as estimativa para os componentes serão positivas.

mm0 <- lme(consumo ~ Mate + Per + Mate:Per,
           random = ~1 | parcela,
           method = "ML",
           data = tb)
anova(mm0)
##             numDF denDF   F-value p-value
## (Intercept)     1    75 26473.649  <.0001
## Mate            4    25     1.675  0.1872
## Per             3    75    37.968  <.0001
## Mate:Per       12    75     6.233  <.0001
#

Para verificar se os resídos apresentam algum padrão de correlação serial, será obtida a correlação entre resíduos de observações vizinhas de primeira e segunda ordem dentro de cada parcela.

# Extrai dados envolvidos no ajuste e adiciona resíduos e ajustados.
tb_fit <- mm0$data
tb_fit <- tb_fit %>%
    mutate(res = residuals(mm0),
           fit = fitted(mm0))

tb_cor <- tb_fit %>%
    select(parcela, per, res) %>%
    spread(key = "per", value = "res") %>%
    select(-1)

lattice::splom(tb_cor, type = c("p", "r"))

tb_cor %>%
    cor()
##             1           2           3           4
## 1  1.00000000 -0.09830009 -0.08129783 -0.24978273
## 2 -0.09830009  1.00000000 -0.18133727 -0.01598283
## 3 -0.08129783 -0.18133727  1.00000000  0.35518427
## 4 -0.24978273 -0.01598283  0.35518427  1.00000000
# Calcula correlação média entre vizinhos dentro de cada parcela.
tb_cor <- tb_fit %>%
    mutate(res_1 = lead(res, n = 1),
           res_2 = lead(res, n = 2)) %>%
    group_by(parcela) %>%
    summarise(viz1 = cor(res, res_1, method = "pearson"),
              viz2 = cor(res, res_2, method = "pearson"))

tb_cor %>%
    summarise_if(is.numeric, mean, na.rm = TRUE)
## # A tibble: 1 x 2
##     viz1   viz2
##    <dbl>  <dbl>
## 1 -0.103 -0.245
tb_cor %>%
    gather(key = "vizinho", value = "pearson", viz1:viz2) %>%
    ggplot(mapping = aes(x = pearson, color = vizinho)) +
    stat_ecdf() +
    geom_vline(xintercept = 0, col = "orange")
## Warning: Removed 2 rows containing non-finite values (stat_ecdf).

#

Pela análise exploratória, parece não haver evidências, pelo menos não gráficas, para contra a hipótese nula de que os erros são independentes entre observações da mesma parcela.

Modelos de efeitos aleatórios com estrutura AR1 nos erros

O modelo a seguir considera uma estrutura de erro autoregressivo de primeira ordem. Ou seja, os erros dentro da mesma parcela apresentam correção decrescente em função do grau de vizinhança. A correlação entre vizinhos de ordem \(k\) é \(\rho^k\). A matriz para 4 medidas no tempo é, \[ \begin{bmatrix} 1 & \rho & \rho^2 & \rho^3 \\ & 1 & \rho & \rho^2 \\ & & 1 & \rho \\ & & & 1 \end{bmatrix}. \]

# Modelo com correlação AR(1) entre observações.
mm1 <- lme(consumo ~ Mate + Per + Mate:Per,
           random = ~1 | parcela,
           correlation = corAR1(value = 0.1, form = ~per | parcela),
           method = "ML",
           data = tb)
anova(mm1)
##             numDF denDF   F-value p-value
## (Intercept)     1    75 26339.562  <.0001
## Mate            4    25     1.852  0.1505
## Per             3    75    39.215  <.0001
## Mate:Per       12    75     6.030  <.0001
summary(mm1)
## Linear mixed-effects model fit by maximum likelihood
##  Data: tb 
##        AIC     BIC    logLik
##   756.3027 820.415 -355.1513
## 
## Random effects:
##  Formula: ~1 | parcela
##         (Intercept) Residual
## StdDev:    0.763366  4.63111
## 
## Correlation Structure: AR(1)
##  Formula: ~per | parcela 
##  Parameter estimate(s):
##        Phi 
## 0.09537041 
## Fixed effects: consumo ~ Mate + Per + Mate:Per 
##                  Value Std.Error DF  t-value p-value
## (Intercept)  100.08135  2.099043 75 47.67951  0.0000
## Mate1.5       -6.56053  2.968495 25 -2.21005  0.0365
## Mate3         -3.09327  2.968495 25 -1.04203  0.3074
## Mate4.5      -12.49008  2.968495 25 -4.20755  0.0003
## Mate6        -21.41568  2.968495 25 -7.21432  0.0000
## Per2         -22.03768  2.785804 75 -7.91071  0.0000
## Per3         -17.80853  2.915621 75 -6.10797  0.0000
## Per4         -21.70833  2.927701 75 -7.41481  0.0000
## Mate1.5:Per2  11.18255  3.939721 75  2.83841  0.0058
## Mate3:Per2     8.07540  3.939721 75  2.04974  0.0439
## Mate4.5:Per2  14.36210  3.939721 75  3.64546  0.0005
## Mate6:Per2    25.38590  3.939721 75  6.44358  0.0000
## Mate1.5:Per3   8.32155  4.123310 75  2.01817  0.0471
## Mate3:Per3     8.29067  4.123310 75  2.01068  0.0480
## Mate4.5:Per3  18.12798  4.123310 75  4.39646  0.0000
## Mate6:Per3    26.87700  4.123310 75  6.51831  0.0000
## Mate1.5:Per4   5.56850  4.140394 75  1.34492  0.1827
## Mate3:Per4     3.09328  4.140394 75  0.74710  0.4573
## Mate4.5:Per4  12.49012  4.140394 75  3.01665  0.0035
## Mate6:Per4    24.43652  4.140394 75  5.90198  0.0000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.94167994 -0.65613082 -0.09842584  0.59321265  3.70606830 
## 
## Number of Observations: 120
## Number of Groups: 30
# Estimativa de \rho.
intervals(mm1)$corStruct
##          lower       est.     upper
## Phi -0.1564151 0.09537041 0.3355168
## attr(,"label")
## [1] "Correlation structure:"
# LRT entre modelos aninhados.
anova(mm1, mm0)
##     Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## mm1     1 23 756.3027 820.4150 -355.1513                         
## mm0     2 22 754.6911 816.0159 -355.3455 1 vs 2 0.3883899  0.5331
#

Pelo teste de razão de verossimilhanças entre modelos encaixados, não rejeitou-se a hipótese nula de que \(\rho = 0\).

Médias ajustadas e comparações múltiplas

# Médias ajustadas para todas as combinações entre níveis.
emm <- emmeans(mm0, specs = ~Per + Mate)
emm
##  Per Mate emmean  SE df lower.CL upper.CL
##  1   0     100.1 2.1 29     95.8    104.4
##  2   0      78.0 2.1 29     73.8     82.3
##  3   0      82.3 2.1 29     78.0     86.6
##  4   0      78.4 2.1 29     74.1     82.7
##  1   1.5    93.5 2.1 25     89.2     97.8
##  2   1.5    82.7 2.1 25     78.3     87.0
##  3   1.5    84.0 2.1 25     79.7     88.4
##  4   1.5    77.4 2.1 25     73.1     81.7
##  1   3      97.0 2.1 25     92.7    101.3
##  2   3      83.0 2.1 25     78.7     87.3
##  3   3      87.5 2.1 25     83.1     91.8
##  4   3      78.4 2.1 25     74.1     82.7
##  1   4.5    87.6 2.1 25     83.3     91.9
##  2   4.5    79.9 2.1 25     75.6     84.2
##  3   4.5    87.9 2.1 25     83.6     92.2
##  4   4.5    78.4 2.1 25     74.1     82.7
##  1   6      78.7 2.1 25     74.3     83.0
##  2   6      82.0 2.1 25     77.7     86.3
##  3   6      87.7 2.1 25     83.4     92.1
##  4   6      81.4 2.1 25     77.1     85.7
## 
## d.f. method: containment 
## Confidence level used: 0.95
# Gráfico com intervalos de confiança.
gg1 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Per,
                     y = emmean)) +
    facet_wrap(facets = ~Mate,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Período de postura",
         y = "Consumo") +
    ggtitle(label = NULL,
            subtitle = "Níveis de mate")

# Gráfico com intervalos de confiança.
gg2 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Níveis de mate",
         y = "Consumo") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

gridExtra::grid.arrange(gg1, gg2, ncol = 1)

# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Per | Mate)
emm
## Mate = 0:
##  Per emmean  SE df lower.CL upper.CL
##  1    100.1 2.1 29     95.8    104.4
##  2     78.0 2.1 29     73.8     82.3
##  3     82.3 2.1 29     78.0     86.6
##  4     78.4 2.1 29     74.1     82.7
## 
## Mate = 1.5:
##  Per emmean  SE df lower.CL upper.CL
##  1     93.5 2.1 25     89.2     97.8
##  2     82.7 2.1 25     78.3     87.0
##  3     84.0 2.1 25     79.7     88.4
##  4     77.4 2.1 25     73.1     81.7
## 
## Mate = 3:
##  Per emmean  SE df lower.CL upper.CL
##  1     97.0 2.1 25     92.7    101.3
##  2     83.0 2.1 25     78.7     87.3
##  3     87.5 2.1 25     83.1     91.8
##  4     78.4 2.1 25     74.1     82.7
## 
## Mate = 4.5:
##  Per emmean  SE df lower.CL upper.CL
##  1     87.6 2.1 25     83.3     91.9
##  2     79.9 2.1 25     75.6     84.2
##  3     87.9 2.1 25     83.6     92.2
##  4     78.4 2.1 25     74.1     82.7
## 
## Mate = 6:
##  Per emmean  SE df lower.CL upper.CL
##  1     78.7 2.1 25     74.3     83.0
##  2     82.0 2.1 25     77.7     86.3
##  3     87.7 2.1 25     83.4     92.1
##  4     81.4 2.1 25     77.1     85.7
## 
## d.f. method: containment 
## Confidence level used: 0.95
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Mate = 0:
##  contrast estimate   SE df t.ratio p.value
##  1 - 2      22.038 2.86 75  7.712  <.0001 
##  1 - 3      17.809 2.86 75  6.232  <.0001 
##  1 - 4      21.708 2.86 75  7.597  <.0001 
##  2 - 3      -4.229 2.86 75 -1.480  0.4545 
##  2 - 4      -0.329 2.86 75 -0.115  0.9994 
##  3 - 4       3.900 2.86 75  1.365  0.5252 
## 
## Mate = 1.5:
##  contrast estimate   SE df t.ratio p.value
##  1 - 2      10.855 2.86 75  3.799  0.0016 
##  1 - 3       9.487 2.86 75  3.320  0.0074 
##  1 - 4      16.140 2.86 75  5.648  <.0001 
##  2 - 3      -1.368 2.86 75 -0.479  0.9636 
##  2 - 4       5.285 2.86 75  1.849  0.2588 
##  3 - 4       6.653 2.86 75  2.328  0.1007 
## 
## Mate = 3:
##  contrast estimate   SE df t.ratio p.value
##  1 - 2      13.962 2.86 75  4.886  <.0001 
##  1 - 3       9.518 2.86 75  3.331  0.0072 
##  1 - 4      18.615 2.86 75  6.515  <.0001 
##  2 - 3      -4.444 2.86 75 -1.555  0.4102 
##  2 - 4       4.653 2.86 75  1.628  0.3692 
##  3 - 4       9.097 2.86 75  3.184  0.0111 
## 
## Mate = 4.5:
##  contrast estimate   SE df t.ratio p.value
##  1 - 2       7.676 2.86 75  2.686  0.0432 
##  1 - 3      -0.319 2.86 75 -0.112  0.9995 
##  1 - 4       9.218 2.86 75  3.226  0.0099 
##  2 - 3      -7.995 2.86 75 -2.798  0.0324 
##  2 - 4       1.543 2.86 75  0.540  0.9490 
##  3 - 4       9.538 2.86 75  3.338  0.0071 
## 
## Mate = 6:
##  contrast estimate   SE df t.ratio p.value
##  1 - 2      -3.348 2.86 75 -1.172  0.6464 
##  1 - 3      -9.068 2.86 75 -3.174  0.0115 
##  1 - 4      -2.728 2.86 75 -0.955  0.7753 
##  2 - 3      -5.720 2.86 75 -2.002  0.1965 
##  2 - 4       0.620 2.86 75  0.217  0.9964 
##  3 - 4       6.340 2.86 75  2.219  0.1275 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
multcomp::cld(emm)
## Mate = 0:
##  Per emmean  SE df lower.CL upper.CL .group
##  2     78.0 2.1 29     73.8     82.3  1    
##  4     78.4 2.1 29     74.1     82.7  1    
##  3     82.3 2.1 29     78.0     86.6  1    
##  1    100.1 2.1 29     95.8    104.4   2   
## 
## Mate = 1.5:
##  Per emmean  SE df lower.CL upper.CL .group
##  4     77.4 2.1 25     73.1     81.7  1    
##  2     82.7 2.1 25     78.3     87.0  1    
##  3     84.0 2.1 25     79.7     88.4  1    
##  1     93.5 2.1 25     89.2     97.8   2   
## 
## Mate = 3:
##  Per emmean  SE df lower.CL upper.CL .group
##  4     78.4 2.1 25     74.1     82.7  1    
##  2     83.0 2.1 25     78.7     87.3  12   
##  3     87.5 2.1 25     83.1     91.8   2   
##  1     97.0 2.1 25     92.7    101.3    3  
## 
## Mate = 4.5:
##  Per emmean  SE df lower.CL upper.CL .group
##  4     78.4 2.1 25     74.1     82.7  1    
##  2     79.9 2.1 25     75.6     84.2  1    
##  1     87.6 2.1 25     83.3     91.9   2   
##  3     87.9 2.1 25     83.6     92.2   2   
## 
## Mate = 6:
##  Per emmean  SE df lower.CL upper.CL .group
##  1     78.7 2.1 25     74.3     83.0  1    
##  4     81.4 2.1 25     77.1     85.7  12   
##  2     82.0 2.1 25     77.7     86.3  12   
##  3     87.7 2.1 25     83.4     92.1   2   
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Mate | Per)
emm
## Per = 1:
##  Mate emmean  SE df lower.CL upper.CL
##  0     100.1 2.1 29     95.8    104.4
##  1.5    93.5 2.1 25     89.2     97.8
##  3      97.0 2.1 25     92.7    101.3
##  4.5    87.6 2.1 25     83.3     91.9
##  6      78.7 2.1 25     74.3     83.0
## 
## Per = 2:
##  Mate emmean  SE df lower.CL upper.CL
##  0      78.0 2.1 29     73.8     82.3
##  1.5    82.7 2.1 25     78.3     87.0
##  3      83.0 2.1 25     78.7     87.3
##  4.5    79.9 2.1 25     75.6     84.2
##  6      82.0 2.1 25     77.7     86.3
## 
## Per = 3:
##  Mate emmean  SE df lower.CL upper.CL
##  0      82.3 2.1 29     78.0     86.6
##  1.5    84.0 2.1 25     79.7     88.4
##  3      87.5 2.1 25     83.1     91.8
##  4.5    87.9 2.1 25     83.6     92.2
##  6      87.7 2.1 25     83.4     92.1
## 
## Per = 4:
##  Mate emmean  SE df lower.CL upper.CL
##  0      78.4 2.1 29     74.1     82.7
##  1.5    77.4 2.1 25     73.1     81.7
##  3      78.4 2.1 25     74.1     82.7
##  4.5    78.4 2.1 25     74.1     82.7
##  6      81.4 2.1 25     77.1     85.7
## 
## d.f. method: containment 
## Confidence level used: 0.95
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Per = 1:
##  contrast   estimate   SE df t.ratio p.value
##  0 - 1.5    6.560533 2.97 25  2.211  0.2086 
##  0 - 3      3.093267 2.97 25  1.042  0.8333 
##  0 - 4.5   12.490083 2.97 25  4.209  0.0025 
##  0 - 6     21.415683 2.97 25  7.216  <.0001 
##  1.5 - 3   -3.467267 2.97 25 -1.168  0.7688 
##  1.5 - 4.5  5.929550 2.97 25  1.998  0.2959 
##  1.5 - 6   14.855150 2.97 25  5.006  0.0003 
##  3 - 4.5    9.396817 2.97 25  3.166  0.0301 
##  3 - 6     18.322417 2.97 25  6.174  <.0001 
##  4.5 - 6    8.925600 2.97 25  3.008  0.0429 
## 
## Per = 2:
##  contrast   estimate   SE df t.ratio p.value
##  0 - 1.5   -4.622017 2.97 25 -1.557  0.5368 
##  0 - 3     -4.982133 2.97 25 -1.679  0.4643 
##  0 - 4.5   -1.872017 2.97 25 -0.631  0.9686 
##  0 - 6     -3.970217 2.97 25 -1.338  0.6710 
##  1.5 - 3   -0.360117 2.97 25 -0.121  0.9999 
##  1.5 - 4.5  2.750000 2.97 25  0.927  0.8839 
##  1.5 - 6    0.651800 2.97 25  0.220  0.9994 
##  3 - 4.5    3.110117 2.97 25  1.048  0.8306 
##  3 - 6      1.011917 2.97 25  0.341  0.9969 
##  4.5 - 6   -2.098200 2.97 25 -0.707  0.9530 
## 
## Per = 3:
##  contrast   estimate   SE df t.ratio p.value
##  0 - 1.5   -1.761017 2.97 25 -0.593  0.9748 
##  0 - 3     -5.197400 2.97 25 -1.751  0.4227 
##  0 - 4.5   -5.637900 2.97 25 -1.900  0.3434 
##  0 - 6     -5.461317 2.97 25 -1.840  0.3741 
##  1.5 - 3   -3.436383 2.97 25 -1.158  0.7745 
##  1.5 - 4.5 -3.876883 2.97 25 -1.306  0.6899 
##  1.5 - 6   -3.700300 2.97 25 -1.247  0.7247 
##  3 - 4.5   -0.440500 2.97 25 -0.148  0.9999 
##  3 - 6     -0.263917 2.97 25 -0.089  1.0000 
##  4.5 - 6    0.176583 2.97 25  0.060  1.0000 
## 
## Per = 4:
##  contrast   estimate   SE df t.ratio p.value
##  0 - 1.5    0.992033 2.97 25  0.334  0.9971 
##  0 - 3     -0.000017 2.97 25  0.000  1.0000 
##  0 - 4.5   -0.000033 2.97 25  0.000  1.0000 
##  0 - 6     -3.020833 2.97 25 -1.018  0.8447 
##  1.5 - 3   -0.992050 2.97 25 -0.334  0.9971 
##  1.5 - 4.5 -0.992067 2.97 25 -0.334  0.9971 
##  1.5 - 6   -4.012867 2.97 25 -1.352  0.6624 
##  3 - 4.5   -0.000017 2.97 25  0.000  1.0000 
##  3 - 6     -3.020817 2.97 25 -1.018  0.8447 
##  4.5 - 6   -3.020800 2.97 25 -1.018  0.8447 
## 
## P value adjustment: tukey method for comparing a family of 5 estimates
multcomp::cld(emm)
## Per = 1:
##  Mate emmean  SE df lower.CL upper.CL .group
##  6      78.7 2.1 25     74.3     83.0  1    
##  4.5    87.6 2.1 25     83.3     91.9   2   
##  1.5    93.5 2.1 25     89.2     97.8   23  
##  3      97.0 2.1 25     92.7    101.3    3  
##  0     100.1 2.1 29     95.8    104.4    3  
## 
## Per = 2:
##  Mate emmean  SE df lower.CL upper.CL .group
##  0      78.0 2.1 29     73.8     82.3  1    
##  4.5    79.9 2.1 25     75.6     84.2  1    
##  6      82.0 2.1 25     77.7     86.3  1    
##  1.5    82.7 2.1 25     78.3     87.0  1    
##  3      83.0 2.1 25     78.7     87.3  1    
## 
## Per = 3:
##  Mate emmean  SE df lower.CL upper.CL .group
##  0      82.3 2.1 29     78.0     86.6  1    
##  1.5    84.0 2.1 25     79.7     88.4  1    
##  3      87.5 2.1 25     83.1     91.8  1    
##  6      87.7 2.1 25     83.4     92.1  1    
##  4.5    87.9 2.1 25     83.6     92.2  1    
## 
## Per = 4:
##  Mate emmean  SE df lower.CL upper.CL .group
##  1.5    77.4 2.1 25     73.1     81.7  1    
##  0      78.4 2.1 29     74.1     82.7  1    
##  3      78.4 2.1 25     74.1     82.7  1    
##  4.5    78.4 2.1 25     74.1     82.7  1    
##  6      81.4 2.1 25     77.1     85.7  1    
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05
grid <- multcomp::cld(emm)
grid <- grid %>%
    group_by(Per) %>%
    mutate(cld = numbers_to_letters(.group)) %>%
    ungroup()

ggplot(data = grid,
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    geom_text(mapping = aes(label = sprintf("%0.2f %s", emmean, cld)),
              angle = 90,
              vjust = 0,
              nudge_x = -0.05) +
    labs(x = "Níveis de mate",
         y = "Consumo") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

#

Na seções a seguir serão apenas considerados os modelos de efeito aleatório, o teste para o para o parâmetro de autocorrelação de primeira ordem, as médias ajustadas e constrastes entre elas.

Análise da variável postura

# ATTENTION: a variável `postura` contém valores ausentes.

# Modelo de efeitos aleatórios.
mm0 <- lme(postura ~ Mate + Per + Mate:Per,
           random = ~1 | parcela,
           method = "ML",
           na.action = na.omit,
           data = tb)

# Modelo com correlação AR(1) entre observações.
mm1 <- update(mm0,
              correlation = corAR1(value = 0.1, form = ~per | parcela))

# Estimativa de \rho.
intervals(mm1)$corStruct
##          lower        est.     upper
## Phi -0.3688275 -0.06999067 0.2419598
## attr(,"label")
## [1] "Correlation structure:"
# LRT entre modelos aninhados.
anova(mm1, mm0)
##     Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## mm1     1 23 366.9966 430.3291 -160.4983                         
## mm0     2 22 365.1943 425.7733 -160.5971 1 vs 2 0.1977131  0.6566
# Estimativas dos parâmetros do modelo escolhido.
summary(mm0)
## Linear mixed-effects model fit by maximum likelihood
##  Data: tb 
##        AIC      BIC    logLik
##   365.1943 425.7733 -160.5971
## 
## Random effects:
##  Formula: ~1 | parcela
##         (Intercept)  Residual
## StdDev:   0.8321831 0.7796921
## 
## Fixed effects: postura ~ Mate + Per + Mate:Per 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  93.92856 0.5606024 72 167.54934  0.0000
## Mate1.5      -1.36906 0.7590592 24  -1.80363  0.0839
## Mate3        -1.17063 0.7590592 24  -1.54221  0.1361
## Mate4.5      -1.07144 0.7590592 24  -1.41154  0.1709
## Mate6        -1.66666 0.7590592 24  -2.19569  0.0380
## Per2         -0.83332 0.5420589 72  -1.53732  0.1286
## Per3          0.00000 0.5420589 72   0.00000  1.0000
## Per4          0.00000 0.5420589 72   0.00000  1.0000
## Mate1.5:Per2  0.33732 0.7339512 72   0.45959  0.6472
## Mate3:Per2    1.82539 0.7339512 72   2.48707  0.0152
## Mate4.5:Per2  1.23017 0.7339512 72   1.67609  0.0981
## Mate6:Per2    1.03172 0.7339512 72   1.40571  0.1641
## Mate1.5:Per3  0.59523 0.7339512 72   0.81100  0.4200
## Mate3:Per3    0.79363 0.7339512 72   1.08132  0.2832
## Mate4.5:Per3  0.00003 0.7339512 72   0.00005  1.0000
## Mate6:Per3    0.09920 0.7339512 72   0.13516  0.8929
## Mate1.5:Per4  0.59525 0.7339512 72   0.81102  0.4200
## Mate3:Per4    0.29762 0.7339512 72   0.40550  0.6863
## Mate4.5:Per4 -0.09920 0.7339512 72  -0.13516  0.8929
## Mate6:Per4    0.29760 0.7339512 72   0.40548  0.6863
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.69310930 -0.38900664 -0.01593018  0.58917028  2.58239450 
## 
## Number of Observations: 116
## Number of Groups: 29
# Quadro de testes para os termos de efeito fixo.
anova(mm0)
##             numDF denDF   F-value p-value
## (Intercept)     1    72 245745.03  <.0001
## Mate            4    24      1.37  0.2735
## Per             3    72      0.76  0.5205
## Mate:Per       12    72      1.23  0.2827
# Médias ajustadas para todas as combinações entre níveis.
emm <- emmeans(mm0, specs = ~Per + Mate)
emm
##  Per Mate emmean    SE df lower.CL upper.CL
##  1   0      93.9 0.561 28     92.8     95.1
##  2   0      93.1 0.561 28     91.9     94.2
##  3   0      93.9 0.561 28     92.8     95.1
##  4   0      93.9 0.561 28     92.8     95.1
##  1   1.5    92.6 0.512 24     91.5     93.6
##  2   1.5    92.1 0.512 24     91.0     93.1
##  3   1.5    93.2 0.512 24     92.1     94.2
##  4   1.5    93.2 0.512 24     92.1     94.2
##  1   3      92.8 0.512 24     91.7     93.8
##  2   3      93.8 0.512 24     92.7     94.8
##  3   3      93.6 0.512 24     92.5     94.6
##  4   3      93.1 0.512 24     92.0     94.1
##  1   4.5    92.9 0.512 24     91.8     93.9
##  2   4.5    93.3 0.512 24     92.2     94.3
##  3   4.5    92.9 0.512 24     91.8     93.9
##  4   4.5    92.8 0.512 24     91.7     93.8
##  1   6      92.3 0.512 24     91.2     93.3
##  2   6      92.5 0.512 24     91.4     93.5
##  3   6      92.4 0.512 24     91.3     93.4
##  4   6      92.6 0.512 24     91.5     93.6
## 
## d.f. method: containment 
## Confidence level used: 0.95
# Gráfico com intervalos de confiança.
gg1 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Per,
                     y = emmean)) +
    facet_wrap(facets = ~Mate,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Período de postura",
         y = "Postura") +
    ggtitle(label = NULL,
            subtitle = "Níveis de mate")

# Gráfico com intervalos de confiança.
gg2 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Níveis de mate",
         y = "Postura") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

gridExtra::grid.arrange(gg1, gg2, ncol = 1)

# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Per | Mate)
multcomp::cld(emm)
## Mate = 0:
##  Per emmean    SE df lower.CL upper.CL .group
##  2     93.1 0.561 28     91.9     94.2  1    
##  1     93.9 0.561 28     92.8     95.1  1    
##  3     93.9 0.561 28     92.8     95.1  1    
##  4     93.9 0.561 28     92.8     95.1  1    
## 
## Mate = 1.5:
##  Per emmean    SE df lower.CL upper.CL .group
##  2     92.1 0.512 24     91.0     93.1  1    
##  1     92.6 0.512 24     91.5     93.6  1    
##  3     93.2 0.512 24     92.1     94.2  1    
##  4     93.2 0.512 24     92.1     94.2  1    
## 
## Mate = 3:
##  Per emmean    SE df lower.CL upper.CL .group
##  1     92.8 0.512 24     91.7     93.8  1    
##  4     93.1 0.512 24     92.0     94.1  1    
##  3     93.6 0.512 24     92.5     94.6  1    
##  2     93.8 0.512 24     92.7     94.8  1    
## 
## Mate = 4.5:
##  Per emmean    SE df lower.CL upper.CL .group
##  4     92.8 0.512 24     91.7     93.8  1    
##  1     92.9 0.512 24     91.8     93.9  1    
##  3     92.9 0.512 24     91.8     93.9  1    
##  2     93.3 0.512 24     92.2     94.3  1    
## 
## Mate = 6:
##  Per emmean    SE df lower.CL upper.CL .group
##  1     92.3 0.512 24     91.2     93.3  1    
##  3     92.4 0.512 24     91.3     93.4  1    
##  2     92.5 0.512 24     91.4     93.5  1    
##  4     92.6 0.512 24     91.5     93.6  1    
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Mate = 0:
##  contrast  estimate    SE df t.ratio p.value
##  1 - 2     8.33e-01 0.542 72  1.537  0.4209 
##  1 - 3     0.00e+00 0.542 72  0.000  1.0000 
##  1 - 4     0.00e+00 0.542 72  0.000  1.0000 
##  2 - 3    -8.33e-01 0.542 72 -1.537  0.4209 
##  2 - 4    -8.33e-01 0.542 72 -1.537  0.4209 
##  3 - 4     0.00e+00 0.542 72  0.000  1.0000 
## 
## Mate = 1.5:
##  contrast  estimate    SE df t.ratio p.value
##  1 - 2     4.96e-01 0.495 72  1.002  0.7485 
##  1 - 3    -5.95e-01 0.495 72 -1.203  0.6270 
##  1 - 4    -5.95e-01 0.495 72 -1.203  0.6270 
##  2 - 3    -1.09e+00 0.495 72 -2.205  0.1316 
##  2 - 4    -1.09e+00 0.495 72 -2.205  0.1316 
##  3 - 4    -1.67e-05 0.495 72  0.000  1.0000 
## 
## Mate = 3:
##  contrast  estimate    SE df t.ratio p.value
##  1 - 2    -9.92e-01 0.495 72 -2.005  0.1958 
##  1 - 3    -7.94e-01 0.495 72 -1.604  0.3830 
##  1 - 4    -2.98e-01 0.495 72 -0.601  0.9313 
##  2 - 3     1.98e-01 0.495 72  0.401  0.9780 
##  2 - 4     6.94e-01 0.495 72  1.403  0.5014 
##  3 - 4     4.96e-01 0.495 72  1.002  0.7485 
## 
## Mate = 4.5:
##  contrast  estimate    SE df t.ratio p.value
##  1 - 2    -3.97e-01 0.495 72 -0.802  0.8533 
##  1 - 3    -3.33e-05 0.495 72  0.000  1.0000 
##  1 - 4     9.92e-02 0.495 72  0.200  0.9971 
##  2 - 3     3.97e-01 0.495 72  0.802  0.8533 
##  2 - 4     4.96e-01 0.495 72  1.002  0.7484 
##  3 - 4     9.92e-02 0.495 72  0.201  0.9971 
## 
## Mate = 6:
##  contrast  estimate    SE df t.ratio p.value
##  1 - 2    -1.98e-01 0.495 72 -0.401  0.9780 
##  1 - 3    -9.92e-02 0.495 72 -0.200  0.9971 
##  1 - 4    -2.98e-01 0.495 72 -0.601  0.9313 
##  2 - 3     9.92e-02 0.495 72  0.200  0.9971 
##  2 - 4    -9.92e-02 0.495 72 -0.200  0.9971 
##  3 - 4    -1.98e-01 0.495 72 -0.401  0.9780 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Mate | Per)
multcomp::cld(emm)
## Per = 1:
##  Mate emmean    SE df lower.CL upper.CL .group
##  6      92.3 0.512 24     91.2     93.3  1    
##  1.5    92.6 0.512 24     91.5     93.6  1    
##  3      92.8 0.512 24     91.7     93.8  1    
##  4.5    92.9 0.512 24     91.8     93.9  1    
##  0      93.9 0.561 28     92.8     95.1  1    
## 
## Per = 2:
##  Mate emmean    SE df lower.CL upper.CL .group
##  1.5    92.1 0.512 24     91.0     93.1  1    
##  6      92.5 0.512 24     91.4     93.5  1    
##  0      93.1 0.561 28     91.9     94.2  1    
##  4.5    93.3 0.512 24     92.2     94.3  1    
##  3      93.8 0.512 24     92.7     94.8  1    
## 
## Per = 3:
##  Mate emmean    SE df lower.CL upper.CL .group
##  6      92.4 0.512 24     91.3     93.4  1    
##  4.5    92.9 0.512 24     91.8     93.9  1    
##  1.5    93.2 0.512 24     92.1     94.2  1    
##  3      93.6 0.512 24     92.5     94.6  1    
##  0      93.9 0.561 28     92.8     95.1  1    
## 
## Per = 4:
##  Mate emmean    SE df lower.CL upper.CL .group
##  6      92.6 0.512 24     91.5     93.6  1    
##  4.5    92.8 0.512 24     91.7     93.8  1    
##  3      93.1 0.512 24     92.0     94.1  1    
##  1.5    93.2 0.512 24     92.1     94.2  1    
##  0      93.9 0.561 28     92.8     95.1  1    
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Per = 1:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5     1.3691 0.759 24  1.804  0.3946 
##  0 - 3       1.1706 0.759 24  1.542  0.5466 
##  0 - 4.5     1.0714 0.759 24  1.412  0.6265 
##  0 - 6       1.6667 0.759 24  2.196  0.2153 
##  1.5 - 3    -0.1984 0.724 24 -0.274  0.9987 
##  1.5 - 4.5  -0.2976 0.724 24 -0.411  0.9936 
##  1.5 - 6     0.2976 0.724 24  0.411  0.9936 
##  3 - 4.5    -0.0992 0.724 24 -0.137  0.9999 
##  3 - 6       0.4960 0.724 24  0.685  0.9577 
##  4.5 - 6     0.5952 0.724 24  0.822  0.9211 
## 
## Per = 2:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5     1.0317 0.759 24  1.359  0.6583 
##  0 - 3      -0.6548 0.759 24 -0.863  0.9077 
##  0 - 4.5    -0.1587 0.759 24 -0.209  0.9995 
##  0 - 6       0.6349 0.759 24  0.836  0.9165 
##  1.5 - 3    -1.6865 0.724 24 -2.330  0.1700 
##  1.5 - 4.5  -1.1905 0.724 24 -1.645  0.4849 
##  1.5 - 6    -0.3968 0.724 24 -0.548  0.9811 
##  3 - 4.5     0.4960 0.724 24  0.685  0.9577 
##  3 - 6       1.2897 0.724 24  1.782  0.4064 
##  4.5 - 6     0.7937 0.724 24  1.097  0.8065 
## 
## Per = 3:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5     0.7738 0.759 24  1.019  0.8439 
##  0 - 3       0.3770 0.759 24  0.497  0.9869 
##  0 - 4.5     1.0714 0.759 24  1.411  0.6265 
##  0 - 6       1.5675 0.759 24  2.065  0.2673 
##  1.5 - 3    -0.3968 0.724 24 -0.548  0.9811 
##  1.5 - 4.5   0.2976 0.724 24  0.411  0.9936 
##  1.5 - 6     0.7936 0.724 24  1.097  0.8065 
##  3 - 4.5     0.6944 0.724 24  0.959  0.8703 
##  3 - 6       1.1905 0.724 24  1.645  0.4849 
##  4.5 - 6     0.4960 0.724 24  0.685  0.9577 
## 
## Per = 4:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5     0.7738 0.759 24  1.019  0.8439 
##  0 - 3       0.8730 0.759 24  1.150  0.7786 
##  0 - 4.5     1.1706 0.759 24  1.542  0.5465 
##  0 - 6       1.3691 0.759 24  1.804  0.3946 
##  1.5 - 3     0.0992 0.724 24  0.137  0.9999 
##  1.5 - 4.5   0.3968 0.724 24  0.548  0.9811 
##  1.5 - 6     0.5952 0.724 24  0.822  0.9211 
##  3 - 4.5     0.2976 0.724 24  0.411  0.9936 
##  3 - 6       0.4960 0.724 24  0.685  0.9577 
##  4.5 - 6     0.1984 0.724 24  0.274  0.9987 
## 
## P value adjustment: tukey method for comparing a family of 5 estimates
grid <- multcomp::cld(emm)
grid <- grid %>%
    group_by(Per) %>%
    mutate(cld = numbers_to_letters(.group)) %>%
    ungroup()

ggplot(data = grid,
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    geom_text(mapping = aes(label = sprintf("%0.2f %s", emmean, cld)),
              angle = 90,
              vjust = 0,
              nudge_x = -0.05) +
    labs(x = "Níveis de mate",
         y = "Postura") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

#

Análise da variável caduzia

# Modelo de efeitos aleatórios.
mm0 <- lme(caduzia ~ Mate + Per + Mate:Per,
           random = ~1 | parcela,
           method = "ML",
           na.action = na.omit,
           data = tb)

# Modelo com correlação AR(1) entre observações.
mm1 <- update(mm0,
              correlation = corAR1(value = 0.1, form = ~per | parcela))

# Estimativa de \rho.
intervals(mm1)$corStruct
##          lower       est.     upper
## Phi -0.2205656 0.08816632 0.3808404
## attr(,"label")
## [1] "Correlation structure:"
# LRT entre modelos aninhados.
anova(mm1, mm0)
##     Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## mm1     1 23 172.7175 236.8298 -63.35876                         
## mm0     2 22 171.0476 232.3724 -63.52381 1 vs 2 0.3301042  0.5656
# Estimativas dos parâmetros do modelo escolhido.
summary(mm0)
## Linear mixed-effects model fit by maximum likelihood
##  Data: tb 
##        AIC      BIC    logLik
##   171.0476 232.3724 -63.52381
## 
## Random effects:
##  Formula: ~1 | parcela
##         (Intercept)  Residual
## StdDev:   0.3057077 0.3437366
## 
## Fixed effects: caduzia ~ Mate + Per + Mate:Per 
##                  Value Std.Error DF  t-value p-value
## (Intercept)   7.960317 0.2057241 75 38.69414  0.0000
## Mate1.5      -0.741700 0.2909378 25 -2.54934  0.0173
## Mate3        -0.491400 0.2909378 25 -1.68902  0.1037
## Mate4.5      -1.222217 0.2909378 25 -4.20095  0.0003
## Mate6        -1.868567 0.2909378 25 -6.42256  0.0000
## Per2         -1.730050 0.2173981 75 -7.95798  0.0000
## Per3         -1.412483 0.2173981 75 -6.49722  0.0000
## Per4         -1.723167 0.2173981 75 -7.92632  0.0000
## Mate1.5:Per2  0.933267 0.3074474 75  3.03553  0.0033
## Mate3:Per2    0.588650 0.3074474 75  1.91464  0.0594
## Mate4.5:Per2  1.111217 0.3074474 75  3.61433  0.0005
## Mate6:Per2    1.973517 0.3074474 75  6.41904  0.0000
## Mate1.5:Per3  0.637850 0.3074474 75  2.07466  0.0414
## Mate3:Per3    0.623067 0.3074474 75  2.02658  0.0463
## Mate4.5:Per3  1.444167 0.3074474 75  4.69728  0.0000
## Mate6:Per3    2.105933 0.3074474 75  6.84974  0.0000
## Mate1.5:Per4  0.444083 0.3074474 75  1.44442  0.1528
## Mate3:Per4    0.271617 0.3074474 75  0.88346  0.3798
## Mate4.5:Per4  1.021150 0.3074474 75  3.32138  0.0014
## Mate6:Per4    1.912717 0.3074474 75  6.22128  0.0000
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.576517355 -0.596609987 -0.004727694  0.488033630  2.861621803 
## 
## Number of Observations: 120
## Number of Groups: 30
# Quadro de testes para os termos de efeito fixo.
anova(mm0)
##             numDF denDF  F-value p-value
## (Intercept)     1    75 8673.713  <.0001
## Mate            4    25    0.941  0.4567
## Per             3    75   40.566  <.0001
## Mate:Per       12    75    6.684  <.0001
# Médias ajustadas para todas as combinações entre níveis.
emm <- emmeans(mm0, specs = ~Per + Mate)
emm
##  Per Mate emmean    SE df lower.CL upper.CL
##  1   0      7.96 0.206 29     7.54     8.38
##  2   0      6.23 0.206 29     5.81     6.65
##  3   0      6.55 0.206 29     6.13     6.97
##  4   0      6.24 0.206 29     5.82     6.66
##  1   1.5    7.22 0.206 25     6.79     7.64
##  2   1.5    6.42 0.206 25     6.00     6.85
##  3   1.5    6.44 0.206 25     6.02     6.87
##  4   1.5    5.94 0.206 25     5.52     6.36
##  1   3      7.47 0.206 25     7.05     7.89
##  2   3      6.33 0.206 25     5.90     6.75
##  3   3      6.68 0.206 25     6.26     7.10
##  4   3      6.02 0.206 25     5.59     6.44
##  1   4.5    6.74 0.206 25     6.31     7.16
##  2   4.5    6.12 0.206 25     5.70     6.54
##  3   4.5    6.77 0.206 25     6.35     7.19
##  4   4.5    6.04 0.206 25     5.61     6.46
##  1   6      6.09 0.206 25     5.67     6.52
##  2   6      6.34 0.206 25     5.91     6.76
##  3   6      6.79 0.206 25     6.36     7.21
##  4   6      6.28 0.206 25     5.86     6.70
## 
## d.f. method: containment 
## Confidence level used: 0.95
# Gráfico com intervalos de confiança.
gg1 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Per,
                     y = emmean)) +
    facet_wrap(facets = ~Mate,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Período de postura",
         y = "???") +
    ggtitle(label = NULL,
            subtitle = "Níveis de mate")

# Gráfico com intervalos de confiança.
gg2 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Níveis de mate",
         y = "???") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

gridExtra::grid.arrange(gg1, gg2, ncol = 1)

# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Per | Mate)
multcomp::cld(emm)
## Mate = 0:
##  Per emmean    SE df lower.CL upper.CL .group
##  2     6.23 0.206 29     5.81     6.65  1    
##  4     6.24 0.206 29     5.82     6.66  1    
##  3     6.55 0.206 29     6.13     6.97  1    
##  1     7.96 0.206 29     7.54     8.38   2   
## 
## Mate = 1.5:
##  Per emmean    SE df lower.CL upper.CL .group
##  4     5.94 0.206 25     5.52     6.36  1    
##  2     6.42 0.206 25     6.00     6.85  1    
##  3     6.44 0.206 25     6.02     6.87  1    
##  1     7.22 0.206 25     6.79     7.64   2   
## 
## Mate = 3:
##  Per emmean    SE df lower.CL upper.CL .group
##  4     6.02 0.206 25     5.59     6.44  1    
##  2     6.33 0.206 25     5.90     6.75  12   
##  3     6.68 0.206 25     6.26     7.10   2   
##  1     7.47 0.206 25     7.05     7.89    3  
## 
## Mate = 4.5:
##  Per emmean    SE df lower.CL upper.CL .group
##  4     6.04 0.206 25     5.61     6.46  1    
##  2     6.12 0.206 25     5.70     6.54  1    
##  1     6.74 0.206 25     6.31     7.16   2   
##  3     6.77 0.206 25     6.35     7.19   2   
## 
## Mate = 6:
##  Per emmean    SE df lower.CL upper.CL .group
##  1     6.09 0.206 25     5.67     6.52  1    
##  4     6.28 0.206 25     5.86     6.70  12   
##  2     6.34 0.206 25     5.91     6.76  12   
##  3     6.79 0.206 25     6.36     7.21   2   
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Mate = 0:
##  contrast estimate    SE df t.ratio p.value
##  1 - 2     1.73005 0.217 75  7.958  <.0001 
##  1 - 3     1.41248 0.217 75  6.497  <.0001 
##  1 - 4     1.72317 0.217 75  7.926  <.0001 
##  2 - 3    -0.31757 0.217 75 -1.461  0.4662 
##  2 - 4    -0.00688 0.217 75 -0.032  1.0000 
##  3 - 4     0.31068 0.217 75  1.429  0.4854 
## 
## Mate = 1.5:
##  contrast estimate    SE df t.ratio p.value
##  1 - 2     0.79678 0.217 75  3.665  0.0025 
##  1 - 3     0.77463 0.217 75  3.563  0.0035 
##  1 - 4     1.27908 0.217 75  5.884  <.0001 
##  2 - 3    -0.02215 0.217 75 -0.102  0.9996 
##  2 - 4     0.48230 0.217 75  2.219  0.1276 
##  3 - 4     0.50445 0.217 75  2.320  0.1025 
## 
## Mate = 3:
##  contrast estimate    SE df t.ratio p.value
##  1 - 2     1.14140 0.217 75  5.250  <.0001 
##  1 - 3     0.78942 0.217 75  3.631  0.0028 
##  1 - 4     1.45155 0.217 75  6.677  <.0001 
##  2 - 3    -0.35198 0.217 75 -1.619  0.3743 
##  2 - 4     0.31015 0.217 75  1.427  0.4869 
##  3 - 4     0.66213 0.217 75  3.046  0.0165 
## 
## Mate = 4.5:
##  contrast estimate    SE df t.ratio p.value
##  1 - 2     0.61883 0.217 75  2.847  0.0285 
##  1 - 3    -0.03168 0.217 75 -0.146  0.9989 
##  1 - 4     0.70202 0.217 75  3.229  0.0098 
##  2 - 3    -0.65052 0.217 75 -2.992  0.0192 
##  2 - 4     0.08318 0.217 75  0.383  0.9808 
##  3 - 4     0.73370 0.217 75  3.375  0.0063 
## 
## Mate = 6:
##  contrast estimate    SE df t.ratio p.value
##  1 - 2    -0.24347 0.217 75 -1.120  0.6785 
##  1 - 3    -0.69345 0.217 75 -3.190  0.0110 
##  1 - 4    -0.18955 0.217 75 -0.872  0.8193 
##  2 - 3    -0.44998 0.217 75 -2.070  0.1725 
##  2 - 4     0.05392 0.217 75  0.248  0.9946 
##  3 - 4     0.50390 0.217 75  2.318  0.1031 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Mate | Per)
multcomp::cld(emm)
## Per = 1:
##  Mate emmean    SE df lower.CL upper.CL .group
##  6      6.09 0.206 25     5.67     6.52  1    
##  4.5    6.74 0.206 25     6.31     7.16  12   
##  1.5    7.22 0.206 25     6.79     7.64   23  
##  3      7.47 0.206 25     7.05     7.89   23  
##  0      7.96 0.206 29     7.54     8.38    3  
## 
## Per = 2:
##  Mate emmean    SE df lower.CL upper.CL .group
##  4.5    6.12 0.206 25     5.70     6.54  1    
##  0      6.23 0.206 29     5.81     6.65  1    
##  3      6.33 0.206 25     5.90     6.75  1    
##  6      6.34 0.206 25     5.91     6.76  1    
##  1.5    6.42 0.206 25     6.00     6.85  1    
## 
## Per = 3:
##  Mate emmean    SE df lower.CL upper.CL .group
##  1.5    6.44 0.206 25     6.02     6.87  1    
##  0      6.55 0.206 29     6.13     6.97  1    
##  3      6.68 0.206 25     6.26     7.10  1    
##  4.5    6.77 0.206 25     6.35     7.19  1    
##  6      6.79 0.206 25     6.36     7.21  1    
## 
## Per = 4:
##  Mate emmean    SE df lower.CL upper.CL .group
##  1.5    5.94 0.206 25     5.52     6.36  1    
##  3      6.02 0.206 25     5.59     6.44  1    
##  4.5    6.04 0.206 25     5.61     6.46  1    
##  0      6.24 0.206 29     5.82     6.66  1    
##  6      6.28 0.206 25     5.86     6.70  1    
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Per = 1:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5     0.7417 0.291 25  2.549  0.1115 
##  0 - 3       0.4914 0.291 25  1.689  0.4583 
##  0 - 4.5     1.2222 0.291 25  4.201  0.0025 
##  0 - 6       1.8686 0.291 25  6.423  <.0001 
##  1.5 - 3    -0.2503 0.291 25 -0.860  0.9086 
##  1.5 - 4.5   0.4805 0.291 25  1.652  0.4803 
##  1.5 - 6     1.1269 0.291 25  3.873  0.0056 
##  3 - 4.5     0.7308 0.291 25  2.512  0.1199 
##  3 - 6       1.3772 0.291 25  4.734  0.0007 
##  4.5 - 6     0.6463 0.291 25  2.222  0.2046 
## 
## Per = 2:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5    -0.1916 0.291 25 -0.658  0.9634 
##  0 - 3      -0.0973 0.291 25 -0.334  0.9971 
##  0 - 4.5     0.1110 0.291 25  0.382  0.9952 
##  0 - 6      -0.1050 0.291 25 -0.361  0.9961 
##  1.5 - 3     0.0943 0.291 25  0.324  0.9974 
##  1.5 - 4.5   0.3026 0.291 25  1.040  0.8344 
##  1.5 - 6     0.0866 0.291 25  0.298  0.9982 
##  3 - 4.5     0.2082 0.291 25  0.716  0.9509 
##  3 - 6      -0.0077 0.291 25 -0.026  1.0000 
##  4.5 - 6    -0.2160 0.291 25 -0.742  0.9443 
## 
## Per = 3:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5     0.1038 0.291 25  0.357  0.9963 
##  0 - 3      -0.1317 0.291 25 -0.453  0.9908 
##  0 - 4.5    -0.2220 0.291 25 -0.763  0.9388 
##  0 - 6      -0.2374 0.291 25 -0.816  0.9233 
##  1.5 - 3    -0.2355 0.291 25 -0.810  0.9252 
##  1.5 - 4.5  -0.3258 0.291 25 -1.120  0.7947 
##  1.5 - 6    -0.3412 0.291 25 -1.173  0.7664 
##  3 - 4.5    -0.0903 0.291 25 -0.310  0.9978 
##  3 - 6      -0.1057 0.291 25 -0.363  0.9960 
##  4.5 - 6    -0.0154 0.291 25 -0.053  1.0000 
## 
## Per = 4:
##  contrast  estimate    SE df t.ratio p.value
##  0 - 1.5     0.2976 0.291 25  1.023  0.8424 
##  0 - 3       0.2198 0.291 25  0.755  0.9408 
##  0 - 4.5     0.2011 0.291 25  0.691  0.9566 
##  0 - 6      -0.0442 0.291 25 -0.152  0.9999 
##  1.5 - 3    -0.0778 0.291 25 -0.268  0.9988 
##  1.5 - 4.5  -0.0965 0.291 25 -0.332  0.9972 
##  1.5 - 6    -0.3418 0.291 25 -1.175  0.7653 
##  3 - 4.5    -0.0187 0.291 25 -0.064  1.0000 
##  3 - 6      -0.2639 0.291 25 -0.907  0.8914 
##  4.5 - 6    -0.2452 0.291 25 -0.843  0.9145 
## 
## P value adjustment: tukey method for comparing a family of 5 estimates
grid <- multcomp::cld(emm)
grid <- grid %>%
    group_by(Per) %>%
    mutate(cld = numbers_to_letters(.group)) %>%
    ungroup()

ggplot(data = grid,
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    geom_text(mapping = aes(label = sprintf("%0.2f %s", emmean, cld)),
              angle = 90,
              vjust = 0,
              nudge_x = -0.05) +
    labs(x = "Níveis de mate",
         y = "???") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

#

Análise da variável camovos

# Modelo de efeitos aleatórios.
mm0 <- lme(camovos ~ Mate + Per + Mate:Per,
           random = ~1 | parcela,
           method = "ML",
           na.action = na.omit,
           data = tb)

# Modelo com correlação AR(1) entre observações.
mm1 <- update(mm0,
              correlation = corAR1(value = 0.1, form = ~per | parcela))

# Estimativa de \rho.
intervals(mm1)$corStruct
##          lower       est.     upper
## Phi -0.2425798 0.05163788 0.3371567
## attr(,"label")
## [1] "Correlation structure:"
# LRT entre modelos aninhados.
anova(mm1, mm0)
##     Model df      AIC       BIC  logLik   Test   L.Ratio p-value
## mm1     1 23 -154.340 -90.22766 100.170                         
## mm0     2 22 -156.228 -94.90321 100.114 1 vs 2 0.1119421  0.7379
# Estimativas dos parâmetros do modelo escolhido.
summary(mm0)
## Linear mixed-effects model fit by maximum likelihood
##  Data: tb 
##        AIC       BIC  logLik
##   -156.228 -94.90321 100.114
## 
## Random effects:
##  Formula: ~1 | parcela
##         (Intercept)   Residual
## StdDev:  0.08285527 0.08668787
## 
## Fixed effects: camovos ~ Mate + Per + Mate:Per 
##                   Value  Std.Error DF  t-value p-value
## (Intercept)   1.9228500 0.05362795 75 35.85537  0.0000
## Mate1.5      -0.1640167 0.07584137 25 -2.16263  0.0403
## Mate3        -0.0544167 0.07584137 25 -0.71751  0.4797
## Mate4.5      -0.2110500 0.07584137 25 -2.78278  0.0101
## Mate6        -0.3778833 0.07584137 25 -4.98255  0.0000
## Per2         -0.4231667 0.05482623 75 -7.71833  0.0000
## Per3         -0.3430667 0.05482623 75 -6.25735  0.0000
## Per4         -0.4148500 0.05482623 75 -7.56663  0.0000
## Mate1.5:Per2  0.2211333 0.07753599 75  2.85201  0.0056
## Mate3:Per2    0.1549333 0.07753599 75  1.99821  0.0493
## Mate4.5:Per2  0.2691000 0.07753599 75  3.47065  0.0009
## Mate6:Per2    0.4851833 0.07753599 75  6.25752  0.0000
## Mate1.5:Per3  0.1639500 0.07753599 75  2.11450  0.0378
## Mate3:Per3    0.1614667 0.07753599 75  2.08247  0.0407
## Mate4.5:Per3  0.3481500 0.07753599 75  4.49017  0.0000
## Mate6:Per3    0.5168500 0.07753599 75  6.66594  0.0000
## Mate1.5:Per4  0.1130833 0.07753599 75  1.45846  0.1489
## Mate3:Per4    0.0579500 0.07753599 75  0.74739  0.4572
## Mate4.5:Per4  0.2339333 0.07753599 75  3.01709  0.0035
## Mate6:Per4    0.4651167 0.07753599 75  5.99872  0.0000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.56143798 -0.60470629  0.03188733  0.58494602  3.38945169 
## 
## Number of Observations: 120
## Number of Groups: 30
# Quadro de testes para os termos de efeito fixo.
anova(mm0)
##             numDF denDF  F-value p-value
## (Intercept)     1    75 7556.111  <.0001
## Mate            4    25    0.457  0.7661
## Per             3    75   37.904  <.0001
## Mate:Per       12    75    6.187  <.0001
# Médias ajustadas para todas as combinações entre níveis.
emm <- emmeans(mm0, specs = ~Per + Mate)
emm
##  Per Mate emmean     SE df lower.CL upper.CL
##  1   0      1.92 0.0536 29     1.81     2.03
##  2   0      1.50 0.0536 29     1.39     1.61
##  3   0      1.58 0.0536 29     1.47     1.69
##  4   0      1.51 0.0536 29     1.40     1.62
##  1   1.5    1.76 0.0536 25     1.65     1.87
##  2   1.5    1.56 0.0536 25     1.45     1.67
##  3   1.5    1.58 0.0536 25     1.47     1.69
##  4   1.5    1.46 0.0536 25     1.35     1.57
##  1   3      1.87 0.0536 25     1.76     1.98
##  2   3      1.60 0.0536 25     1.49     1.71
##  3   3      1.69 0.0536 25     1.58     1.80
##  4   3      1.51 0.0536 25     1.40     1.62
##  1   4.5    1.71 0.0536 25     1.60     1.82
##  2   4.5    1.56 0.0536 25     1.45     1.67
##  3   4.5    1.72 0.0536 25     1.61     1.83
##  4   4.5    1.53 0.0536 25     1.42     1.64
##  1   6      1.54 0.0536 25     1.43     1.66
##  2   6      1.61 0.0536 25     1.50     1.72
##  3   6      1.72 0.0536 25     1.61     1.83
##  4   6      1.60 0.0536 25     1.48     1.71
## 
## d.f. method: containment 
## Confidence level used: 0.95
# Gráfico com intervalos de confiança.
gg1 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Per,
                     y = emmean)) +
    facet_wrap(facets = ~Mate,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Período de postura",
         y = "???") +
    ggtitle(label = NULL,
            subtitle = "Níveis de mate")

# Gráfico com intervalos de confiança.
gg2 <-
ggplot(data = as.data.frame(emm),
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    labs(x = "Níveis de mate",
         y = "???") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

gridExtra::grid.arrange(gg1, gg2, ncol = 1)

# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Per | Mate)
multcomp::cld(emm)
## Mate = 0:
##  Per emmean     SE df lower.CL upper.CL .group
##  2     1.50 0.0536 29     1.39     1.61  1    
##  4     1.51 0.0536 29     1.40     1.62  1    
##  3     1.58 0.0536 29     1.47     1.69  1    
##  1     1.92 0.0536 29     1.81     2.03   2   
## 
## Mate = 1.5:
##  Per emmean     SE df lower.CL upper.CL .group
##  4     1.46 0.0536 25     1.35     1.57  1    
##  2     1.56 0.0536 25     1.45     1.67  1    
##  3     1.58 0.0536 25     1.47     1.69  1    
##  1     1.76 0.0536 25     1.65     1.87   2   
## 
## Mate = 3:
##  Per emmean     SE df lower.CL upper.CL .group
##  4     1.51 0.0536 25     1.40     1.62  1    
##  2     1.60 0.0536 25     1.49     1.71  12   
##  3     1.69 0.0536 25     1.58     1.80   2   
##  1     1.87 0.0536 25     1.76     1.98    3  
## 
## Mate = 4.5:
##  Per emmean     SE df lower.CL upper.CL .group
##  4     1.53 0.0536 25     1.42     1.64  1    
##  2     1.56 0.0536 25     1.45     1.67  1    
##  1     1.71 0.0536 25     1.60     1.82   2   
##  3     1.72 0.0536 25     1.61     1.83   2   
## 
## Mate = 6:
##  Per emmean     SE df lower.CL upper.CL .group
##  1     1.54 0.0536 25     1.43     1.66  1    
##  4     1.60 0.0536 25     1.48     1.71  12   
##  2     1.61 0.0536 25     1.50     1.72  12   
##  3     1.72 0.0536 25     1.61     1.83   2   
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Mate = 0:
##  contrast estimate     SE df t.ratio p.value
##  1 - 2     0.42317 0.0548 75  7.718  <.0001 
##  1 - 3     0.34307 0.0548 75  6.257  <.0001 
##  1 - 4     0.41485 0.0548 75  7.567  <.0001 
##  2 - 3    -0.08010 0.0548 75 -1.461  0.4660 
##  2 - 4    -0.00832 0.0548 75 -0.152  0.9987 
##  3 - 4     0.07178 0.0548 75  1.309  0.5600 
## 
## Mate = 1.5:
##  contrast estimate     SE df t.ratio p.value
##  1 - 2     0.20203 0.0548 75  3.685  0.0024 
##  1 - 3     0.17912 0.0548 75  3.267  0.0087 
##  1 - 4     0.30177 0.0548 75  5.504  <.0001 
##  2 - 3    -0.02292 0.0548 75 -0.418  0.9752 
##  2 - 4     0.09973 0.0548 75  1.819  0.2726 
##  3 - 4     0.12265 0.0548 75  2.237  0.1227 
## 
## Mate = 3:
##  contrast estimate     SE df t.ratio p.value
##  1 - 2     0.26823 0.0548 75  4.892  <.0001 
##  1 - 3     0.18160 0.0548 75  3.312  0.0076 
##  1 - 4     0.35690 0.0548 75  6.510  <.0001 
##  2 - 3    -0.08663 0.0548 75 -1.580  0.3961 
##  2 - 4     0.08867 0.0548 75  1.617  0.3753 
##  3 - 4     0.17530 0.0548 75  3.197  0.0107 
## 
## Mate = 4.5:
##  contrast estimate     SE df t.ratio p.value
##  1 - 2     0.15407 0.0548 75  2.810  0.0314 
##  1 - 3    -0.00508 0.0548 75 -0.093  0.9997 
##  1 - 4     0.18092 0.0548 75  3.300  0.0079 
##  2 - 3    -0.15915 0.0548 75 -2.903  0.0245 
##  2 - 4     0.02685 0.0548 75  0.490  0.9612 
##  3 - 4     0.18600 0.0548 75  3.393  0.0060 
## 
## Mate = 6:
##  contrast estimate     SE df t.ratio p.value
##  1 - 2    -0.06202 0.0548 75 -1.131  0.6716 
##  1 - 3    -0.17378 0.0548 75 -3.170  0.0116 
##  1 - 4    -0.05027 0.0548 75 -0.917  0.7959 
##  2 - 3    -0.11177 0.0548 75 -2.039  0.1833 
##  2 - 4     0.01175 0.0548 75  0.214  0.9965 
##  3 - 4     0.12352 0.0548 75  2.253  0.1186 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
# Médias ajustadas para níveis de um fator fixando níveis dos demais.
emm <- emmeans(mm0, specs = ~Mate | Per)
multcomp::cld(emm)
## Per = 1:
##  Mate emmean     SE df lower.CL upper.CL .group
##  6      1.54 0.0536 25     1.43     1.66  1    
##  4.5    1.71 0.0536 25     1.60     1.82  12   
##  1.5    1.76 0.0536 25     1.65     1.87  12   
##  3      1.87 0.0536 25     1.76     1.98   2   
##  0      1.92 0.0536 29     1.81     2.03   2   
## 
## Per = 2:
##  Mate emmean     SE df lower.CL upper.CL .group
##  0      1.50 0.0536 29     1.39     1.61  1    
##  1.5    1.56 0.0536 25     1.45     1.67  1    
##  4.5    1.56 0.0536 25     1.45     1.67  1    
##  3      1.60 0.0536 25     1.49     1.71  1    
##  6      1.61 0.0536 25     1.50     1.72  1    
## 
## Per = 3:
##  Mate emmean     SE df lower.CL upper.CL .group
##  1.5    1.58 0.0536 25     1.47     1.69  1    
##  0      1.58 0.0536 29     1.47     1.69  1    
##  3      1.69 0.0536 25     1.58     1.80  1    
##  4.5    1.72 0.0536 25     1.61     1.83  1    
##  6      1.72 0.0536 25     1.61     1.83  1    
## 
## Per = 4:
##  Mate emmean     SE df lower.CL upper.CL .group
##  1.5    1.46 0.0536 25     1.35     1.57  1    
##  0      1.51 0.0536 29     1.40     1.62  1    
##  3      1.51 0.0536 25     1.40     1.62  1    
##  4.5    1.53 0.0536 25     1.42     1.64  1    
##  6      1.60 0.0536 25     1.48     1.71  1    
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05
# Comparações múltiplas.
contrast(emm, method = "pairwise")
## Per = 1:
##  contrast   estimate     SE df t.ratio p.value
##  0 - 1.5    1.64e-01 0.0758 25  2.163  0.2264 
##  0 - 3      5.44e-02 0.0758 25  0.718  0.9505 
##  0 - 4.5    2.11e-01 0.0758 25  2.783  0.0694 
##  0 - 6      3.78e-01 0.0758 25  4.983  0.0003 
##  1.5 - 3   -1.10e-01 0.0758 25 -1.445  0.6057 
##  1.5 - 4.5  4.70e-02 0.0758 25  0.620  0.9705 
##  1.5 - 6    2.14e-01 0.0758 25  2.820  0.0642 
##  3 - 4.5    1.57e-01 0.0758 25  2.065  0.2660 
##  3 - 6      3.23e-01 0.0758 25  4.265  0.0021 
##  4.5 - 6    1.67e-01 0.0758 25  2.200  0.2125 
## 
## Per = 2:
##  contrast   estimate     SE df t.ratio p.value
##  0 - 1.5   -5.71e-02 0.0758 25 -0.753  0.9415 
##  0 - 3     -1.01e-01 0.0758 25 -1.325  0.6785 
##  0 - 4.5   -5.80e-02 0.0758 25 -0.765  0.9381 
##  0 - 6     -1.07e-01 0.0758 25 -1.415  0.6243 
##  1.5 - 3   -4.34e-02 0.0758 25 -0.572  0.9779 
##  1.5 - 4.5 -9.33e-04 0.0758 25 -0.012  1.0000 
##  1.5 - 6   -5.02e-02 0.0758 25 -0.662  0.9628 
##  3 - 4.5    4.25e-02 0.0758 25  0.560  0.9796 
##  3 - 6     -6.78e-03 0.0758 25 -0.089  1.0000 
##  4.5 - 6   -4.93e-02 0.0758 25 -0.649  0.9652 
## 
## Per = 3:
##  contrast   estimate     SE df t.ratio p.value
##  0 - 1.5    6.67e-05 0.0758 25  0.001  1.0000 
##  0 - 3     -1.07e-01 0.0758 25 -1.411  0.6263 
##  0 - 4.5   -1.37e-01 0.0758 25 -1.808  0.3915 
##  0 - 6     -1.39e-01 0.0758 25 -1.832  0.3783 
##  1.5 - 3   -1.07e-01 0.0758 25 -1.412  0.6257 
##  1.5 - 4.5 -1.37e-01 0.0758 25 -1.809  0.3911 
##  1.5 - 6   -1.39e-01 0.0758 25 -1.833  0.3779 
##  3 - 4.5   -3.01e-02 0.0758 25 -0.396  0.9945 
##  3 - 6     -3.19e-02 0.0758 25 -0.421  0.9930 
##  4.5 - 6   -1.87e-03 0.0758 25 -0.025  1.0000 
## 
## Per = 4:
##  contrast   estimate     SE df t.ratio p.value
##  0 - 1.5    5.09e-02 0.0758 25  0.672  0.9607 
##  0 - 3     -3.53e-03 0.0758 25 -0.047  1.0000 
##  0 - 4.5   -2.29e-02 0.0758 25 -0.302  0.9981 
##  0 - 6     -8.72e-02 0.0758 25 -1.150  0.7786 
##  1.5 - 3   -5.45e-02 0.0758 25 -0.718  0.9503 
##  1.5 - 4.5 -7.38e-02 0.0758 25 -0.973  0.8646 
##  1.5 - 6   -1.38e-01 0.0758 25 -1.822  0.3840 
##  3 - 4.5   -1.93e-02 0.0758 25 -0.255  0.9990 
##  3 - 6     -8.37e-02 0.0758 25 -1.104  0.8030 
##  4.5 - 6   -6.44e-02 0.0758 25 -0.848  0.9126 
## 
## P value adjustment: tukey method for comparing a family of 5 estimates
grid <- multcomp::cld(emm)
grid <- grid %>%
    group_by(Per) %>%
    mutate(cld = numbers_to_letters(.group)) %>%
    ungroup()

ggplot(data = grid,
       mapping = aes(x = Mate,
                     y = emmean)) +
    facet_wrap(facets = ~Per,
               nrow = 1) +
    geom_point() +
    geom_errorbar(mapping = aes(ymin = lower.CL,
                                ymax = upper.CL),
                  width = 0.1) +
    geom_text(mapping = aes(label = sprintf("%0.2f %s", emmean, cld)),
              angle = 90,
              vjust = 0,
              nudge_x = -0.05) +
    labs(x = "Níveis de mate",
         y = "???") +
    ggtitle(label = NULL,
            subtitle = "Perído de postura")

#

Análise da variável massaovos

# Modelo para um fator.
m0 <- lm(massaovos ~ Mate,
         data = subset(tb, per == 1))

# par(mfrow = c(2, 2))
# plot(m0)
# layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: massaovos
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Mate       4  19.748  4.9371  0.6605 0.6252
## Residuals 25 186.882  7.4753
# Médias ajustadas.
emm <- emmeans(m0, specs = ~Mate)
multcomp::cld(emm)
##  Mate emmean   SE df lower.CL upper.CL .group
##  6      51.0 1.12 25     48.7     53.3  1    
##  4.5    51.3 1.12 25     49.0     53.6  1    
##  3      51.9 1.12 25     49.6     54.2  1    
##  0      52.4 1.12 25     50.1     54.7  1    
##  1.5    53.3 1.12 25     51.0     55.6  1    
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05
# Comparações múltiplas.
contrast(emm, method = "pairwise")
##  contrast  estimate   SE df t.ratio p.value
##  0 - 1.5     -0.882 1.58 25 -0.559  0.9798 
##  0 - 3        0.492 1.58 25  0.311  0.9978 
##  0 - 4.5      1.138 1.58 25  0.721  0.9496 
##  0 - 6        1.368 1.58 25  0.867  0.9064 
##  1.5 - 3      1.374 1.58 25  0.870  0.9051 
##  1.5 - 4.5    2.020 1.58 25  1.280  0.7055 
##  1.5 - 6      2.250 1.58 25  1.425  0.6177 
##  3 - 4.5      0.647 1.58 25  0.410  0.9937 
##  3 - 6        0.876 1.58 25  0.555  0.9803 
##  4.5 - 6      0.230 1.58 25  0.146  0.9999 
## 
## P value adjustment: tukey method for comparing a family of 5 estimates