Exemplo prático

Agora que já temos uma ideia dos procedimentos por trás das análises, vamos utilizar os pacotes asreml-R[butlerASRemlRReferenceManual2018], lme4 [bates_lme4_2015] e lme4breeding [bates_lme4_2015] em bancos de dados reais.

Pacote necessário

library(lme4)
library(lme4breeding)
library(asreml)

Banco de dados

Utilizaremos um banco de dados disponível publicamente neste link, referente ao artigo de Chaves et al. (2022). O banco de dados contém 34 tratamentos avaliados em blocos completos casualizados, com 5 a 10 repetições por tratamento (sim, é desbalanceado), os quais foram avaliados por nove anos consecutivos. Para realizar as análises, escolheremos somente um ano (“Yr5”):

data = read.table("data_cupu.txt", header = TRUE, sep = ";")
data = data[data$Harvests == "Yr5",]
head(data)

Definindo os fatores

Antes de montar o modelo, temos que transformas os vetores-coluna do conjunto de dados em fatores:

data = transform(data,
                 Hybrids = factor(Hybrids),
                 Replicates = factor(Replicates))

str(data)

'data.frame':   250 obs. of  6 variables:
 $ Harvests  : chr  "Yr5" "Yr5" "Yr5" "Yr5" ...
 $ Plots     : int  1171 1172 1173 1174 1175 1176 1177 1178 1179 1181 ...
 $ Replicates: Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 1 ...
 $ Hybrids   : Factor w/ 34 levels "H117","H118",..: 1 1 1 1 1 1 1 1 1 2 ...
 $ yd        : num  19.5 18.6 43.4 39.9 114.1 ...
 $ wb        : int  0 0 0 0 0 0 0 0 0 0 ...

Modelo

Antes de construir o modelo no pacote, é importante tê-lo em mente em termos lineares. No nosso caso, o modelo é o mesmo definido previamente nos exemplos manuais anteriores, com uma única diferença: a dimensão das matrizes e vetores:

\[ \mathbf y = \mathbf X \mathbf b + \mathbf Z \mathbf u + \mathbf e \]

Utilizando a sintaxe do ASReml-R, o modelo pode ser escrito da seguinte forma:

#asreml-R
mod.dbc = asreml(fixed = yd ~ Replicates,
               random = ~Hybrids,
               data = data)

ASReml Version 4.2 15/04/2025 10:21:21
          LogLik        Sigma2     DF     wall
 1     -733.6691      135.1139    240   10:21:21
 2     -733.5223      133.8850    240   10:21:21
 3     -733.3983      132.3498    240   10:21:21
 4     -733.3520      130.9607    240   10:21:21
 5     -733.3511      130.7648    240   10:21:21

A atribuição à um objeto (no nosso caso, mod.dbc1) é imprescindível para obtenção dos resultados. Agora que já ajustamos o modelo, quais os próximos passos?

1. Teste de hipóteses

Para modelos lineares mistos, a significância dos efeitos aleatórios é testadas via LRT (Likelihood Ratio Test):

\[ LRT = (-2 \times LogL_R) - (-2 \times LogL) \] isto é, a diferença entre o logaritmo do ponto máximo da função de verossimilhança residual entre um modelo completo e um modelo reduzido (\(LogL_R\)), isto é, sem o efeito em teste. Já temos o modelo completo, e, portanto, devemos ajustar o modelo reduzido (\(\mathbf y = \mathbf X \mathbf b + \mathbf e\)):

#asreml-R
mod.dbc.red = asreml(fixed = yd ~ Replicates,
               data = data)

Com os dois modelos, calculamos o LRT utilizando a função “lrt.asreml”:

#asreml-R
lrt.asreml(mod.dbc, mod.dbc.red)

Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

                    df LR-statistic Pr(Chisq)   
mod.dbc/mod.dbc.red  1       8.4414  0.001834 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Os valores do teste são comparados com os valores da distribuição \(\chi^2\), com um grau de liberdade e o p-valor desejado (geralmente \(p \le 0.05\) ou \(p \le 0.01\)). O número de graus de liberdade está relacionado com a diferença entre o número de parâmetros de um modelo para outro. Neste caso, a diferença é de somente um parâmetro: o efeito genotípico. Os valores \(\chi^2\) que determinarão se o efeito é significativo ou não são:

Portanto, conclui-se que o efeito genotípico é significativo à 1% de significância, pois \(LRT > \chi^2\) (\(8.44 > 6.63\)). Neste caso, podemos focar nos demais resultados.

2. Componentes de variância

O modelo fornecerá os componentes de variância dos efeitos aleatórios. No nosso exemplo, estes são o efeito genotípico e o efeito residual. O comando para obter os componentes de variância é:

summary(mod.dbc)$varcomp

Perceba que a todo momento, nos referimos à estes valores como “componentes de variância”. De fato, eles compõem a variância total presente no experimento, dada pela variância fenotípica individual(\(\sigma^2_f\)). Portanto, \(\sigma^2_f\) é a soma dos componentes de variância do modelo: \[ \sigma^2_f = \sigma^2_g + \sigma^2_e \]

Esta variância pode ser ponderada pelo número de repetições, gerando a variância fenotípica a nível de médias: \[ \overline{\sigma^2_f} = \sigma^2_g + \frac{\sigma^2_e}{r} \]

em que \(r\) é o número de repetições.

No nosso conjunto de dados, teremos:

compvar = summary(mod.dbc)$varcomp
sig2f = sum(compvar[,1]); sig2f

[1] 151.7376

sig2f_med = compvar["Hybrids","component"]+compvar["units!R","component"]/
  nlevels(data$Replicates);sig2f_med

[1] 34.04925

3. Parâmetros genéticos

Alguns parâmetros genéticos são muito úteis para realizar algumas inferências

• Acurácia \[ r = \sqrt{1- \frac{PEV}{\sigma^2_g}} \]

• Confiabilide

Em que \(PEV\) é a variância do erro de predição.

No nosso conjunto de dados, utilizaremos os seguintes comandos:

pred = predict(mod.dbc, classify = "Hybrids", vcov = T)

PEVi = (pred$pvals$std.error)^2

confi = 1-(PEVi/compvar["Hybrids","component"]) #Confiabilidade por tratamento
acu = sqrt(1-mean(PEVi)/compvar["Hybrids","component"]) #Acurácia

• Herdabilidade

• A nível de indivíduos

\[ H^2 =\frac{\sigma^2_g}{{\sigma^2_f}} \]

• A nível de médias

\[ H^2 =\frac{\sigma^2_g}{\overline{\sigma^2_f}} \]

numrep = nlevels(data$Replicates)

H2in = vpredict(mod.dbc, H2in ~ V1/(V1+V2))
H2med = compvar["Hybrids","component"]/(compvar["Hybrids","component"]+
                                          (compvar["units!R","component"]/numrep))

Em suma, temos:

data.frame(
  "Parâmetro" = c("Acurácia","Herdabilidade (individual)","Herdabilidade (médias)"),
  "Valor" = c(acu, H2in$Estimate, H2med)
)

4. Seleção

Depois realizar o diagnóstico da população em estudo, podemos focar na seleção. O BLUP é a principal unidade de seleção para o melhoramento genético, por ponderar o desempenho do tratamento pela qualidade e quantidade de informações disponíveis (Piepho et al., 2008). Utilizando o ASReml-R, podemos ter o BLUP centralizado na média (média = 0), ou a média BLUP, isto é, o próprio BLUP somado ao intercepto do modelo. Previamente, já realizamos a função para obter a média BLUP:

pred = predict(mod.dbc, classify = "Hybrids")

O seguinte comando fornecerá um vetor com o BLUP centralizado na média:

summary(mod.dbc, coef=T)$coef.random

Unindo os dois e adicionando a confiabilidade por tratamento, teremos:

BLUP = cbind(pred$pvals[,1:2],summary(mod.dbc, coef=T)$coef.random[,1],confi)
colnames(BLUP)[2:4] = c("Média BLUP", "BLUP", "Confiabilidade")
BLUP

Exemplo prático - Blocos incompletos

Nem sempre teremos as melhores condições para implantar o experimento. Por vezes, queremos testar um grande número de tratamentos que não possuem propágulos suficientes para que haja uma quantidade adequada de repetições. Outra situação comum é a limitação da área disponível para o experimento. De fato, a combinação das duas situações são corriqueiras no melhoramento. Neste caso, é comum montar experientos em blocos incompletos. Veja abaixo um comparativo entre o DBC e duas categorias de blocos incompletos, látice e blocos aumentados:

Banco de dados

A seguir, mostraremos um exemplo em um conjunto de dados reais cujos experimentos foram implantados em látice. Os dados estão disponíveis neste link, e foram utilizados no artigo de Dias et al. (2022). O banco de dados contém 36 tratamentos avaliados em látice com 6 blocos e 2 repetições, em quatro locais diferentes. Para realizar as análises, escolheremos somente um local (“E15”):

data.dbi = read.table("data.dbi.txt", header = TRUE, sep=",")
data.dbi = data.dbi[data.dbi$Location == 'E15',]
head(data.dbi)

Definindo os fatores

Modelo

No banco de dados temos os efeitos de repetição e bloco. Geralmente, o efeito de repetição é considerado fixo. Já o efeito de bloco geralmente é aleatório. Quando a recuperação da informação interblocos for interessante, utiliza-se o efeito de blocos como aleatório. Caso contrário, este é tido como fixo. No exemplo, consideraremos este efeito como aleatório. Portanto, nosso modelo será:

\[ \mathbf y = \mathbf X \mathbf r + \mathbf Z_1 \mathbf u + \mathbf Z_2 \mathbf b +\mathbf e \]

em que \(\mathbf r\) é o vetor de efeitos de repetição acompanhado de sua matriz de incidência \(\mathbf X\), \(\mathbf u\) é o vetor de efeitos aleatórios de genótipos, com a matriz de incidência \(\mathbf Z_i\); \(\mathbf b\) é o vetor de efeito aleatórios de blocos dentro de repetições acompanhado de sua matriz de incidência \(\mathbf Z_2\); e \(\mathbf e\) é o vetor dos erros.

No ASReml-R, o modelo é descrito como:

mod.dbi = asreml(fixed = GY ~ Rep,
               random = ~Hybrid + Block,
               data = data.dbi)
mod.dbi = update(mod.dbi)

Uma vez ajustado o modelo, seguiremos os mesmos passos vistos anteriormente para o experimento em DBC:

1. Teste de hipóteses

mod.dbi.red = asreml(fixed = GY ~ Rep,
               random = ~ Block,
               data = data.dbi)

lrt(mod.dbi, mod.dbi.red)

Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

                    df LR-statistic Pr(Chisq)    
mod.dbi/mod.dbi.red  1        11.79 0.0002978 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Componentes de variância

Neste exemplo, estes os efeito aleatórios são os genotípico, de bloco e residual.

summary(mod.dbi)$varcomp

3. Parâmetros genéticos

Primeiro, calcularemos a acurácia e a confiabilidade por tratamento:

compvar = summary(mod.dbi)$varcomp

pred = predict(mod.dbi, classify = "Hybrid", vcov = T)

PEVi = (pred$pvals$std.error)^2

confi = 1-(PEVi/compvar["Hybrid","component"]) #Confiabilidade por tratamento
acu = sqrt(1-mean(PEVi)/compvar["Hybrid","component"]) #Acurácia

[1] 0.8326237

Em seguida, calcularemos a herdabilidade:

H2 = vpredict(mod.dbi, h2 ~ V2/(V1+V2+V3))

    Estimate        SE
h2 0.5349425 0.1206246

4. Seleção

Para obter a média BLUP:

pred = predict(mod.dbi, classify = "Hybrid")

Para o BLUP centralizado na média:

summary(mod.dbi, coef=T)$coef.random[grep("Hybrid", rownames(summary(mod.dbi, coef=T)$coef.random)),]

Unindo os dois e adicionando a confiabilidade por tratamento, teremos:

BLUP = cbind(pred$pvals[,1:2],
             summary(mod.dbi, coef=T)$coef.random[grep("Hybrid", rownames(summary(mod.dbi, coef=T)$coef.random)),][,1],
             confi)
colnames(BLUP)[2:4] = c("Média BLUP", "BLUP", "Confiabilidade")
BLUP

Referências utilizadas nesta apresentação