library(metan)
library(rio)

# gerar tabelas html
print_tbl <- function(table, digits = 3, ...){
  knitr::kable(table, booktabs = TRUE, digits = digits, ...)
}

# dados
df_g <- import("http://bit.ly/df_g", setclass = "tbl")
inspect(df_g, verbose = FALSE) %>% print_tbl()

Modelo misto

A função gamem() pode ser usada para analisar experimentos únicos(experimentos unilaterais) usando um modelo de efeito misto de acordo com o seguinte modelo:

$$ y_{ij} = \mu + \alpha_i + \tau_j + \varepsilon_ {ij} $$

onde \(y_ {ij}\) é o valor observado para o \(i\)-ésimo genótipo no \(j\)-ésimo bloco (\(i\) = 1, 2, … \(g\); \(j\) = 1, 2,…, \(r\)); sendo \(g\) e \(r\) o número de genótipos e blocos, respectivamente; \(\alpha_i\) é o efeito aleatório do genótipo \(i\); \(\tau_j\) é o efeito fixo do bloco \(j\); e \(\varepsilon_ {ij}\) é o erro aleatório associado a \(y_{ij}\). Neste exemplo, usaremos os dados de exemplo df_g.

gen_mod <- 
  gamem(df_g,
        gen = GEN,
        rep = BLOCO,
        resp = everything())
## Evaluating trait ALT_PLANT |====                                 | 10% 00:00:00 
Evaluating trait ALT_ESP |========                               | 20% 00:00:00 
Evaluating trait COMPES |============                            | 30% 00:00:00 
Evaluating trait DIAMES |================                        | 40% 00:00:00 
Evaluating trait COMP_SAB |===================                   | 50% 00:00:00 
Evaluating trait DIAM_SAB |=======================               | 60% 00:00:00 
Evaluating trait MGE |==============================             | 70% 00:00:00 
Evaluating trait NFIL |==================================        | 80% 00:00:00 
Evaluating trait MMG |=======================================    | 90% 00:00:00 
Evaluating trait NGE |===========================================| 100% 00:00:01 
## Method: REML/BLUP
## Random effects: GEN
## Fixed effects: REP
## Denominador DF: Satterthwaite's method
## ---------------------------------------------------------------------------
## P-values for Likelihood Ratio Test of the analyzed traits
## ---------------------------------------------------------------------------
##     model ALT_PLANT  ALT_ESP   COMPES   DIAMES COMP_SAB DIAM_SAB      MGE
##  Complete        NA       NA       NA       NA       NA       NA       NA
##  Genotype  2.27e-12 2.36e-13 0.000224 5.95e-07 9.69e-09 0.000311 4.67e-08
##     NFIL     MMG     NGE
##       NA      NA      NA
##  0.00145 8.3e-08 0.00907
## ---------------------------------------------------------------------------
## All variables with significant (p < 0.05) genotype effect

A maneira mais fácil de obter os resultados do modelo acima é usando a função gmd(), ou seu shortcut gmd().

Diagnósticos do modelo

plot(gen_mod, type = "res") # padrão
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

plot(gen_mod, type = "re") # padrão

Detalhes da análise

details <- gmd(gen_mod, "details")
## Class of the model: gamem
## Variable extracted: details
print_tbl(details)

LRT

lrt <- gmd(gen_mod, "lrt") 
## Class of the model: gamem
## Variable extracted: lrt
print_tbl(lrt)

Componentes de variância

vcomp <- gmd(gen_mod, "vcomp")
## Class of the model: gamem
## Variable extracted: vcomp
print_tbl(vcomp)

plot(gen_mod, type = "vcomp")

Parâmetros genéticos

genpar <- gmd(gen_mod, "genpar")
## Class of the model: gamem
## Variable extracted: genpar
print_tbl(genpar)

BLUPs preditos

blupg <- gmd(gen_mod, "blupg")
## Class of the model: gamem
## Variable extracted: blupg
print_tbl(blupg)

# plotar os BLUPS (default)
plot_blup(gen_mod)

# Trait MGE
plot_blup(gen_mod,
          var = "MGE",
          height.err.bar = 0,
          col.shape = c("black", "gray"),
          x.lab = "Massa de grãos por espiga (g)",
          y.lab = "Híbridos de milho")

Modelos mistos - dentro de ambientes

df_ge <- import("http://bit.ly/df_ge", setclass = "tbl")
mod_gen_whithin <- 
    gamem(df_ge,
          gen = GEN,
          rep = BLOCO,
          resp = everything(),
          by = ENV, verbose = FALSE)

gmd(mod_gen_whithin, "lrt") %>%  print_tbl()

gmd(mod_gen_whithin, "vcomp") %>% print_tbl()