12.2.1 Visualização gráfica

A seguinte função proporciona uma visualização gráfica de todos os pares de correlação possíveis (scatter-plot)

pairs(datacor)

12.2.2 Estimativa dos coeficientes de correlação

corr <- corr_coef(datacor)
print(corr)
# ---------------------------------------------------------------------------
# Pearson's correlation coefficient
# ---------------------------------------------------------------------------
#       CD    CL    CW    PH    EH    EP    EL    ED
# CD 1.000 0.300 0.484 0.315 0.281 0.175 0.912 0.390
# CL 0.300 1.000 0.738 0.325 0.397 0.391 0.255 0.697
# CW 0.484 0.738 1.000 0.505 0.519 0.425 0.458 0.737
# PH 0.315 0.325 0.505 1.000 0.932 0.638 0.380 0.661
# EH 0.281 0.397 0.519 0.932 1.000 0.870 0.363 0.630
# EP 0.175 0.391 0.425 0.638 0.870 1.000 0.263 0.458
# EL 0.912 0.255 0.458 0.380 0.363 0.263 1.000 0.385
# ED 0.390 0.697 0.737 0.661 0.630 0.458 0.385 1.000
# ---------------------------------------------------------------------------
# p-values for the correlation coefficients
# ---------------------------------------------------------------------------
#          CD       CL       CW       PH       EH       EP       EL       ED
# CD 0.00e+00 1.39e-04 1.54e-10 6.06e-05 3.90e-04 2.88e-02 1.97e-61 4.94e-07
# CL 1.39e-04 0.00e+00 3.92e-28 3.45e-05 2.84e-07 4.55e-07 1.29e-03 4.74e-24
# CW 1.54e-10 3.92e-28 0.00e+00 1.83e-11 3.76e-12 3.25e-08 1.81e-09 5.31e-28
# PH 6.06e-05 3.45e-05 1.83e-11 0.00e+00 1.11e-69 3.10e-19 9.80e-07 5.66e-21
# EH 3.90e-04 2.84e-07 3.76e-12 1.11e-69 0.00e+00 4.80e-49 3.28e-06 1.19e-18
# EP 2.88e-02 4.55e-07 3.25e-08 3.10e-19 4.80e-49 0.00e+00 8.92e-04 1.83e-09
# EL 1.97e-61 1.29e-03 1.81e-09 9.80e-07 3.28e-06 8.92e-04 0.00e+00 6.88e-07
# ED 4.94e-07 4.74e-24 5.31e-28 5.66e-21 1.19e-18 1.83e-09 6.88e-07 0.00e+00

12.2.3 Combinando visualização gráfica e numérica (I)

pairs.panels(datacor)

Na figura acima, os pontos observados são plotados na diagonal inferior. Na diagonal, é apresentada a estimativa da densidade Kernel (linha preta) e um histgrama de cada variável. A diagonal superior contém os coeficientes de correlação.

12.2.4 Combinando visualização gráfica e numérica (II)

A função corr_plot()do pacote metan retorna um gráfico semelhante ao anterior, no entanto possui diversas opções, tais como mudança no tamanho da letra dependendo da magnitude da correlação e indicação de cores para correlações significativas.

corr_plot(datacor)

Figure 12.1: Scatter plot de uma matriz de correlação de Pearson

A função corr_plot() pode ser utilizada com o operador %>%. Em adição, é possível escolher as variáveis a serem plotadas. Para isto, basta digitar o nome das variáveis.

maize %>%
  corr_plot(CD, EL, PERK, NKR, CW,
            shape.point = 19,
            size.point = 2,
            alpha.point = 0.5,
            alpha.diag = 0,
            pan.spacing = 0,
            col.sign = "gray",
            alpha.sign = 0.3,
            axis.labels = TRUE,
            progress = FALSE)

Figure 12.2: Scatter plot de uma matriz de correlação de Pearson

maize %>%
  corr_plot(CD, EL, PERK, NKR, CW,
            shape.point = 21,
            col.point = "black",
            fill.point = "orange",
            size.point = 2,
            alpha.point = 0.6,
            maxsize = 4,
            minsize = 2,
            smooth = TRUE,
            col.smooth = "black",
            col.sign = "cyan",
            upper = "scatter",
            lower = "corr",
            diag.type = "density",
            col.diag = "cyan",
            pan.spacing = 0,
            lab.position = "bl")

Figure 12.3: Scatter plot de uma matriz de correlação de Pearson

12.2.5 Combinando visualização gráfica e numérica (III)

A função corrplot.mixed() do pacote corrplot³⁶ também é uma boa opção para visualização gráfica, principalmente quando um grande número de combinações está disponível.

• Criando a matrix de correlação utilizando o conjunto de dados datacor.
• 8 variáveis | 28 combinações

cor1 <- cor(datacor)
corrplot.mixed(cor1,
               upper = "ellipse",
               lower = "number",
               number.digits = 2)

12.2.6 Combinando visualização gráfica e numérica (IV)

• Criando a matrix de correlação utilizando o conjunto de dados dataset.
• 14 variáveis | 91 combinações

cor2 <- cor(numeric_var)
pval <- cor.mtest(cor2)$p

A função corrplot() do pacote corrplot permite uma poderosa personalização. Esta função tem a vantagem de apresentar um elevado número de combinações em um gráfico claro e intuitivo.

corrplot(cor2,
         method = "pie",
         p.mat = pval,
         sig.level = 0.05,
         insig = "blank",
         type = "lower",
         diag = F,
         tl.col = "black",
         tl.srt = 45)

Figure 12.4: Gráfico de pizza de uma matriz de correlação de Pearson

corrplot(cor2,
         method = "ellipse",
         p.mat = pval,
         sig.level = 0.05,
         insig = "blank",
         type = "upper",
         diag = F,
         tl.col = "black",
         tl.srt = 45)

12.4.1 Paramétrico

O intervalo de confiança para o coeficiente de correlação pode ser obtido utilizando a função cor.mtest() do pacote corrplot, conforme o seguinte exemplo.

ci_corr <- cor.mtest(datacor)

12.4.2 Não paramétrico

Um estimador não paramétrico do intervalo de confiança do coeficiente de correlação de Pearson foi proposto por Olivoto et al. (2018). Este estimador é baseado no tamanho da amostra e força de associações e pode ser estimado usando a função corr_ci() do pacote metan. É possível estimar o intervalo de confiança declarando o tamanho da amostra (n) e o coeficiente de correlação (r), ou usando um dataframe. O código a seguir calcula o intervalo de confiança para os possíveis pares de correlação entre as variáveis que contém E no nome. Note o benefício do operador %>% neste caso.

maize %>% 
  select(contains("E")) %>%
  corr_ci(verbose = FALSE) %>%
  plot_ci()

12.7.1 Introdução

A análise de trilha vem se destacando na área do melhoramento genético, pois a seleção para melhoria de um caractere desejável que possui difícil mensuração e baixa herdabilidade, pode ser realizada indiretamente por outro caractere, diretamente associado a este, mas que apresente alta herdabilidade e seja de fácil mensuração. Esta técnica é baseada em ideias originalmente desenvolvidas por Sewall Wright (Wright 1921), no entanto desde sua concepção até a consolidação do método, algumas divergências quanto a fidedignidade do método matemático que explica as relações de causa e efeito foram observadas. Em 1922, Henry E. Niles, em um artigo³⁷ publicado na revista Genetics intitulado Correlation, Causation and Wright’s theory of path coefficients, fez uma crítica ao método proposto por Wright, afirmando que a base filosófica do método dos coeficientes de trilha era falha. Niles, testando o método de Wright, evidenciou em alguns de seus resultados coeficientes superiores a |1|, afirmando […]estes resultados são ridículos[…] e que Wright teria de fornecer provas bem mais convincentes do que ele estava apresentando Niles (1922).

No ano seguinte, 1923, Sewall Wright em seu artigo³⁸ entitulado The theory of path coefficients: a reply to Niles’s criticism, publicado na mesma revista Genetics, consolida seu método ao concluir que Niles pareceu se basear em conceitos matemáticos incorretos, resultado de uma falha em reconhecer que um coeficiente de trilha não é uma função simétrica de duas variáveis, mas que ele necessariamente tem direção. Este autor conclui seu trabalho afirmando que a análise de trilha não fornece uma fórmula geral para deduzir relações causais a partir do conhecimento das correlações. Ela é, no entanto, dentro de certas limitações, um método de avaliar as consequências lógicas de uma hipótese de relação causal em um sistema de variáveis correlacionadas. Acrescenta ainda que as críticas oferecidas por Niles em nada invalidam a teoria do método ou sua aplicação (Wright 1923). Atualmente, o método estatístico é consolidado, e utilizado mundialmente em diversas áreas da ciência.

12.7.2 Estimativa

A decomposição das correlações lineares em efeitos diretos e indiretos de um conjunto de p-variáveis explicativas é realizada o sistema de equações normais

\[ X'X\hat \beta = X'Y \]

que tem como resolução

\[ \hat \beta = X'X^{-1} X'Y \]

onde \(\hat \beta\) é o vetor dos coeficiente de regressão parcial (\(\hat \beta_1\), \(\hat \beta_2\), \(\hat \beta_3\),…,\(\hat \beta_p\)) para p + 1; \(X'X^{-1}\) é a inversa da matriz de correlação linear entre as variáveis explicativas e \(X'Y\) é a matriz de correlação de cada variável explicativa, com a variável dependente.

Após a estimativa dos coeficientes de regressão (\(\hat \beta_p\)), os efeitos diretos e indiretos do conjunto de p-variáveis explicativas podem ser estimados. Considere o seguinte exemplo, onde um conjunto de variáveis explicativas (a, b, c) são utilizadas para explicar as relações de causa e efeito na resposta de uma variável dependente (y). Após as estimativas dos coeficientes de regressão parcial (\(\hat \beta_1\), \(\hat \beta_2\) e \(\hat \beta_3\)), os efeitos diretos e indiretos de a sobre y são dados por:

\[ r_{a:y} = \hat \beta_1 + \hat \beta_{2_{ra:b}} + \hat \beta_{3_{ra:c}} \]

onde \(r_{a:y}\) é a correlação linear entre a e y, \(\hat \beta_1\) é o efeito direto de a em y; \(\hat \beta_{2_{ra:b}}\) é o efeito indireto de a em y via b e \(\hat \beta_{3_{ra:c}}\) é o efeito indireto de a em y via c. Regressões semelhantes são utilizadas para estimativa dos efeitos de b e c, conforme segue:

\[\begin{gather*} r_{b:y} = \hat \beta_{1_{rb:a}} + \hat \beta_2 + \hat \beta_{3_{rb:c}}\\ r_{c:y} = \hat \beta_{1_{rc:a}} + \hat \beta_{2_{rc:b}} + \hat \beta_3 \end{gather*}\]

Como exemplo numérico, vamos utilizar as variáveis PERK, EH, CDED como variáveis explicativas e a variável KW como dependente, do conjunto de dados maize. Para seleção destas variáveis, a função select() é utilizada.

x <- maize %>% select(PERK, EH, CDED)
y <- maize %>% select(KW)

xx <- cor(x) # Correlação entre as variáveis explicativas
xy <- cor(x, y) # Correlação das explicativas com a dependente
b <- solve(xx) %*% xy # Estimativa dos coeficientes

PERK_KW_DIR <- b[1] # Direto de PERK em KW
PERK_KW_EH <- b[2] * xx[1,2] # Indireto de PERK em KW via EH
PERK_KW_CDED <- b[3] * xx[1,3] # Indireto de PERK em KW via CDED

EH_KW_PERK <- b[1] * xx[2,1] # Indireto de EH em KW via PERK
EH_KW_DIR <- b[2] # Direto de EH em KW
EH_KW_CDED <- b[3] * xx[2,3] # Indireto de EH em KW via CDED

CDED_KW_PERK <- b[1] * xx[3,1] # Indireto de CDED em KW via PERK
CDED_KW_EH <- b[2] * xx[3,2] # Indireto de CDED em KW via EH
CDED_KW_DIR <- b[3] # Direto de CDED em KW

# Coeficientes de trilha (direto na diagonal, indireto fora da diagonal)
coeff <- matrix(c(PERK_KW_DIR, PERK_KW_EH, PERK_KW_CDED,
                  EH_KW_PERK, EH_KW_DIR, EH_KW_CDED,
                  CDED_KW_PERK, CDED_KW_EH, CDED_KW_DIR),
                ncol = 3)
rownames(coeff) <- colnames(coeff) <- c("PERK", "EH", "CDED")
coeff
#             PERK          EH        CDED
# PERK -0.10410823 0.002222624  0.05948542
# EH   -0.01473329 0.690110850 -0.04548514
# CDED  0.09200901 0.010613425 -0.16102929

Utilizando os conhecimentos acumulados até agora, estes mesmos coeficientes podem ser estimados de maneira mais “elegantemente”, utilizando o código abaixo.

n <- ncol(xx)
Coeff <- data.frame(xx)
for (i in 1:n) {
  for (j in 1:n) {
    Coeff[j, i] <- b[j] * xx[j, i]
  }
}
rownames(coeff) <- colnames(coeff) <- c("PERK", "EH", "CDED")
Coeff
#             PERK          EH        CDED
# PERK -0.10410823 0.002222624  0.05948542
# EH   -0.01473329 0.690110850 -0.04548514
# CDED  0.09200901 0.010613425 -0.16102929

12.7.3 Multicolinearidade

Embora esta análise revele associações de causa e efeito, sua estimativa é baseada em princípios de regressão múltipla. Assim, as estimativas dos parâmetros podem ser enviesadas devido a natureza complexa dos dados, em que a resposta da variável dependente está ligada a um grande número de variáveis explicativas, que são muitas vezes correlacionadas ou multicolineares entre si (Graham 2003). Assim, sempre que duas supostas variáveis explicativas se apresentam altamente associadas, é difícil estimar as relações de cada variável explicativa individualmente, uma vez que vários parâmetros resolvem o sistema de equações normais. A esta particularidade é atribuída o nome de multicolinearidade (Blalock 1963).

Os principais meios utilizados para identificar o grau de multicolinearidade em uma matriz de variáveis explicativas são os seguintes:

• Número de condição (CN): O número de condição é calculado pela razão entre o maior e menor autovalor (\(\lambda\)) da matriz de correlação \(X'X\), de acordo com a expressão

\[ {\rm{NC = }}\frac{{\lambda_{\rm{Max}}}}{{\lambda_{\rm{Min}}}} \]

O grau de multicolinearidade é considerado fraco se NC \(\leq\) 100, moderado se 100 \(\leq\) NC \(\leq\) 1000 e severo quando NC > 1000.

• Determinante da matriz \(X'X\) (D): O determinante da cada matriz de correlação é estimado pelo produto de seus respectivos autovalores, para \(\lambda_j > 0\), de acordo com a expressão

\[ \mathop D\nolimits_{{\boldsymbol{X'X}}} {\rm{ = }}\prod\limits_{j = 1}^p {\lambda j} \]

Um determinantes muito próximo a zero indica dependência linear entre as características explicativas, indicando problemas graves de multicolinearidade.

• Fator de inflação de variância (VIF): Os (VIFs) são utilizados para medir o quanto a variância dos coeficientes de regressão estimados (\(\hat \beta_k\)) foi inflada em comparação a quando os caracteres explicativos não são linearmente associados. A estimativa do VIF do k-ésimo elemento de \(\hat \beta\) é dada pela soma dos quocientes do quadrado de cada componente do autovetor pelo seu respectivo autovalor associado, de acordo com a expressão

\[ \mathop {{\rm{VIF}}}\nolimits_{\mathop \beta \nolimits_k } = \left( {\frac{{\mathop {{\rm{(AV}}}\nolimits_{KC1} {)^2}}}{{\lambda 1}} + \frac{{{{(\mathop {{\rm{AV}}}\nolimits_{KC2} )}^2}}}{{\lambda 2}} + ... + \frac{{{{(\mathop {{\rm{AV}}}\nolimits_{KCp} )}^2}}}{{\lambda p}}} \right) \]

onde \(\mathop {{\rm{VIF}}}\nolimits_{\mathop \beta \nolimits_k }\) é o fator de inflação de variância o k-ésimo elemento de \(\beta\) para k = 1, 2, …, p; \(\mathop {{\rm{EV}}}\nolimits_{KC{\rm{1}}}\) é o componente do k-ésimo autovetor para k = 1, 2, …, p e C = 1, 2, …, p; e \(\lambda\) é o autovalor associado ao respectivo autovetor para \(\lambda\) = 1, 2, …, p. Os VIFs também podem ser considerados como os elementos da diagonal da inversa da matriz \(X'X\). Considera-se que a presença de VIFs maiores que 10 é um indicativo de multicolinearidade .

cor_data <- cor(datacor)
n <- nrow(datacor)
cold1 <- colindiag(cor_data, n = n)
print(cold1)
# Severe multicollinearity in the matrix! Pay attention on the variables listed bellow
# CN = 1830.678
# Matrix determinant: 2.32e-05 
# Largest correlation: PH x EH = 0.932 
# Smallest correlation: CD x EP = 0.175 
# Number of VIFs > 10: 3 
# Number of correlations with r >= |0.8|: 3 
# Variables with largest weight in the last eigenvalues: 
# EH > PH > EP > CL > ED > CW > CD > EL

cold2 <- colindiag(datacor)
print(cold2)
# Severe multicollinearity in the matrix! Pay attention on the variables listed bellow
# CN = 1830.678
# Matrix determinant: 2.32e-05 
# Largest correlation: PH x EH = 0.932 
# Smallest correlation: CD x EP = 0.175 
# Number of VIFs > 10: 3 
# Number of correlations with r >= |0.8|: 3 
# Variables with largest weight in the last eigenvalues: 
# EH > PH > EP > CL > ED > CW > CD > EL

cold3 <- colindiag(data_ge2, by = ENV)
print(cold3)
# # A tibble: 4 x 2
#   ENV   data      
#   <fct> <list>    
# 1 A1    <colindig>
# 2 A2    <colindig>
# 3 A3    <colindig>
# 4 A4    <colindig>

non_collinear_vars(data_ge2)
#          Parameter                                       values
# 1       Predictors                                           10
# 2              VIF                                         7.16
# 3 Condition Number                                       56.797
# 4      Determinant                                 0.0008810515
# 5         Selected PERK, EP, CDED, NKR, PH, NR, TKW, EL, CD, ED
# 6          Removed                          EH, CL, CW, KW, NKE
non_collinear_vars(data_ge2, EH, CL, CW, KW, NKE, max_vif = 5)
#          Parameter          values
# 1       Predictors               4
# 2              VIF           2.934
# 3 Condition Number          11.248
# 4      Determinant    0.2400583901
# 5         Selected NKE, EH, CL, CW
# 6          Removed              KW

12.7.4 Métodos para ajustar a multicolinearidade

Embora os problemas relacionados a multicolinearidade se apresente como uma dificuldade na estimativa de coeficientes de trilha, algumas medidas podem ser tomadas visando mitigar seus efeitos indesejáveis, quando esta for detectada pelos métodos acima citados. Sabe-se atualmente, que a exclusão das variáveis responsáveis por inflar a variância de um coeficiente de regressão é um dos métodos mais indicados para reduzir a multicolinearidade em matrizes de variáveis explicativas (T. Olivoto, Souza, et al. 2017). A identificação destas variáveis, no entanto, pode ser uma tarefa difícil. Recentemente, T. Olivoto, Nardino, et al. (2017) propuseram a utilização de procedimentos stepwise juntamente com análise de trilha sequencial visando identificar um conjunto de variáveis com alto poder explicativo, mas que não se apresentem altamente correlacionadas. Quando a exclusão das variáveis responsáveis não é um procedimento considerado pelo pesquisador, por exemplo, devido a um número reduzido de variáveis explicativa, ou pela importância em conhecer seus efeitos, uma terceira opção é realizar a análise de trilha com todas as variáveis explicativas, porém com a inclusão de um pequeno valor nos elementos da diagonal \(X'X\), conhecida como regressão em crista³⁹. Este procedimento, no entanto superestima os efeitos diretos, principalmente daquelas variáveis com alto VIF. (T. Olivoto, Souza, et al. 2017).

12.7.5 Análise tradicional

Esta sessão está focada principalmente em três objetivos: (i) diagnóstico da multicolinearidade; (ii) seleção de variáveis preditoras; e (iii) estimação dos coeficientes de trilha. Embora esta seja a sequência correta a ser seguida para a estimativa dos coeficientes de trilha, utilizaremos somente a função path_coeff(), que possibilita todas estas abordagens. Para isto, o conjunto de dados maize será utilizado.

pathtodas <- 
  maize %>%
  path_coeff(KW, verbose = FALSE)

Com a função acima, os coeficientes de trilha foram estimados considerando a variável KW (massa de grãos por espiga) como dependente, e todas as outras variáveis numéricas do conjunto de dados como explicativas. A função summary() pode ser utilizada para resumir os resultados da análise. Note que os coeficientes foram estimados considerando todos os níveis dos fatores ENV, GEN, e REP.

De acordo com o NC, VIF e determinante da matriz, a multicolinearidade na matriz das variáveis explicativa é severa. Por exemplo, foram observados oito VIFs > 10 e o determinante da matriz foi de \(8.619 \times 10^{-12}\). A análise dos autovalores-autovetores (pathtodas$weightvar) indicou que, em ordem de importância, as variáveis que mais contribuem para a multicolinearidadeque são: CL > ED > CDED > EH > CW > PH > NKE > EP > TKW > PERK > NR > EL > NKR > CD. Conforme discutido, temos basicamente duas opções para contornar os problemas da elevada multicolineridade em nossos dados. Excluir as variáveis responsáveis pela multicolinearidade, ou manter todas as variáveis e incluir um fator de correção na diagonal da matrix . Vamos começar pela última opção.

12.7.6 Incluindo um fator de correção (k)

A variância dos coeficientes de regressão é reduzida quando quando um valor k para \(0 < k \leq 1\) é incluído na diagonal da matriz de correlação \(X'X\). Com esta técnica, a estimativa dos coeficientes de regressão é dado por:

\[ \hat \beta = (X'X+Ik)^{-1} X'Y \]

A escolha da magnitude de k, no entanto, deve ser cautelosa. Sugere-se que o valor a ser incluído seja aquele menor possível que estabilize os coeficientes de regressão (\(\beta_p\)). Felizmente, grande parte deste trabalho já foi realizado quando utilizamos a função path_coeff() na sessão anterior. Um conjunto de estimativas de \(\beta_p\) foi estimado com 101 valores de k, (\(k = 0, 0.01, ..., 1\)) e representado graficamente.

pathtodas$plot

Figure 12.5: Valores de beta obtidos com 101 valores de k

print()

O gráfico acima nos proporciona uma interpretação sobre qual é o valor de k mais indicado a ser utilizado. Para fins didátidos, escolheremos, por enquanto, o valor de k igual a 0.04, valor no qual os coeficientes de regressão da maioria das variávies se estabiliza. Para incluir este valor de correção, utilizaremos novamente a função path_coeff(), no entanto, agora, incluiremos o argumento correction = 0.04.

pathtodas_k <- 
  maize %>% 
  path_coeff(KW, correction = 0.04, verbose = FALSE)

Com a inclusão do fator de correção (k = 0.04) na diagonal de \(X'X\), a multicolinearidade foi classificada como moderada (NC = 141). Inevitavelmente, temos duas opções para obtermos menores níveis de multicollinearidade. A primeira é aumentar o valor de k, digamos, para 0.1. Isto iria reduzir ainda mais o nível de multicolinearidade em nossa matriz, no entanto, o viés na estimativa dos coeficientes aumentaria. A segunda (e mais razoável) opção, é a exclusão das variáveis que mais causam problemas de multicolinearidade. Por exemplo, podemos considerar as variáveis com maior peso nos últimos autovalores, ou aquelas com maior VIF. As variáveis CL e EL apresentam alta correlação (veja a seção estimativas dos coeficientes de correlação), assim poderíamos manter apenas uma destas variáveis. A mesma interpretação pode ser considerada para CDED. Esta é uma co-variável (razão do diâmetro do sabugo e diâmetro da espiga, CDED = CD/ED). Vamos considerar então a exclusão destas variáveis.

12.7.7 Excluindo variáveis

O ajuste do novo modelo excluindo estas variáveis é facilmente realizado. Para isto, iremos utilizar dois argumentos da função path_coeff() não vistos até agora: pred e exclude. As variáveis informadas em pred podem ser as variáveis preditoras (default) ou as variáveis a serem excluídas, se exclude = TRUE. Vamos ao exemplo.

path_exclude <- 
  maize %>% 
  path_coeff(resp = KW,
             pred = c(PERK, EH, CDED),
             exclude = TRUE,
             verbose = FALSE)

Abaixo, um resumo das três abordagens realizadas até agora, utilizando o resumo apresentado na tabela abaixo.

Conforme também observado por Hoerl and Kennard (1970) e T. Olivoto, Souza, et al. (2017), a exclusão de variáveis responsáveis pela multicolinearidade foi mais eficiente que a inclusão do k, proporcionando menores níveis de multicolinearidade e maior precisão do modelo (maior \(R^2\) e menor residual). Os níveis de multicolinearidade ao excluir as variáveis ainda preocupam. Vimos que tanto a identificação das variáveis responsáveis pela multicolinearidae quanto o ajuste do modelo declarando preditores específicos é um procedimento relativamente simples utilização a função path_coeff(). Mas, e se algum procedimento estatistico-computacional facilitasse ainda mais essa tarefa? Vamos, a partir de agora, considerar isso.

T. Olivoto, Nardino, et al. (2017) sugeriram a utilização de regressões stepwise para seleção de um conjunto de preditores com minima multicolinearidade em análise de trilha. Esta opção está disponível na função path_coeff(). Baseado em um algoritmo heurístico iterativo executado pelo argumento brutstep = TRUE, um conjunto de preditores com mínima multicolienaridade é selecionado com base nos valores de VIF. Posteriormente, uma série de regressões stepwise são ajustadas. A primeira regressão stepwise é ajustada considerando \(p-1\) variáveis preditoras selecionadas, sendo p o número de variáveis selecionadas no processo iterativo. O segundo modelo ajusta uma regressão considerando \(p-2\) variávies selecionadas, e assim por diante até o último modelo, que considera apenas duas variáveis selecionadas. Vamos ao exemplo.

12.7.8 Seleção de variáveis em análise de trilha

path_step <- 
  maize %>% 
  path_coeff(resp = KW,
             brutstep = TRUE)
# --------------------------------------------------------------------------
# The algorithm has selected a set of 10 predictors with largest VIF = 7.16. 
# Selected predictors: PERK EP CDED NKR PH NR TKW EL CD ED 
# A forward stepwise-based selection procedure will fit 8 models.
# --------------------------------------------------------------------------
# Adjusting the model 1 with 9 predictors (12.5% concluded)
# Adjusting the model 2 with 8 predictors (25% concluded)
# Adjusting the model 3 with 7 predictors (37.5% concluded)
# Adjusting the model 4 with 6 predictors (50% concluded)
# Adjusting the model 5 with 5 predictors (62.5% concluded)
# Adjusting the model 6 with 4 predictors (75% concluded)
# Adjusting the model 7 with 3 predictors (87.5% concluded)
# Adjusting the model 8 with 2 predictors (100% concluded)
# Done!
# --------------------------------------------------------------------------
# Summary of the adjusted models 
# --------------------------------------------------------------------------
#    Model  AIC Numpred    CN Determinant    R2 Residual maxVIF
#  MODEL_1 1099       9 50.99     0.00216 0.945   0.0545   6.96
#  MODEL_2 1098       8 45.48     0.00396 0.945   0.0550   6.80
#  MODEL_3 1097       7 37.46     0.02618 0.944   0.0555   6.42
#  MODEL_4 1103       6 35.70     0.03481 0.942   0.0582   5.71
#  MODEL_5 1116       5 26.39     0.08049 0.936   0.0642   5.71
#  MODEL_6 1129       4 21.07     0.17146 0.930   0.0705   5.71
#  MODEL_7 1170       3  4.97     0.55639 0.907   0.0933   1.78
#  MODEL_8 1232       2  1.57     0.95068 0.860   0.1402   1.05
# --------------------------------------------------------------------------

Três objetos são criados por esta função: Summary, Models e Selectedpred. O objeto Summary contém um resumo do procedimento, listando o número do modelo, o valor do AIC, o diagnóstico da multicolinearidade e os valores de R2 e residual. O objeto Models, contém todos os modelos ajustados, e o objeto Selectedpred, contém o nome das variávies preditoras selecionadas no processo iterativo. Podemos notar que o algorítmo selecionou um conjunto com 10 preditores (PERK, EP, CDED, NKR, PH, NR, TKW, EL, CD, ED) que apresenta multicolinearidade em niveis aceitáveis. Assim, qualquer um destes modelos poderia ser utilizado sem maiores problemas em relação à isto. O procedimento stepwise realizado com diferentes números de variáveis selecionados também permite a seleção de um modelo mais parcimonio, tarefa que ficará a critério do pesquisador.

12.7.9 Análise de trilha para cada nível de um fator

Em alguns casos, é de interesse estimar os coeficientes de trilha para cada nível de um fator, por exemplo, para cada ambiente em um ensaio multi-ambiente. Utilizando o argumento by, isto é facilmente realizado. . Para fins didáticos vamos estimar os coeficientes para cada ambiente do conjunto de dados maize.

path_levels <- 
  maize %>%
  group_by(ENV) %>%
  path_coeff(resp = KW,
             pred = c(TKW, NKE, PERK))
# Weak multicollinearity. 
# Condition Number = 4.697
# You will probably have path coefficients close to being unbiased. 
# Weak multicollinearity. 
# Condition Number = 2.296
# You will probably have path coefficients close to being unbiased. 
# Weak multicollinearity. 
# Condition Number = 4.244
# You will probably have path coefficients close to being unbiased. 
# Weak multicollinearity. 
# Condition Number = 1.884
# You will probably have path coefficients close to being unbiased.