
FERRAMENTAS GRÁFICAS NO PROCESSO DE SELEÇÃO DE

VARIÁVEIS

Lucas Ragiotto

Dissertação apresentada à Universidade Es-

tadual Paulista “Júlio de Mesquita Filho”

para a obtenção do t́ıtulo de Mestre em Bio-

metria.

BOTUCATU

São Paulo - Brasil

Março – 2019


FERRAMENTAS GRÁFICAS NO PROCESSO DE SELEÇÃO DE

VARIÁVEIS

Lucas Ragiotto

Orientadora: Prof. Dr. Luzia Aparecida Trinca

Dissertação apresentada à Universidade Es-

tadual Paulista “Júlio de Mesquita Filho”

para a obtenção do t́ıtulo de Mestre em Bio-

metria.

BOTUCATU

São Paulo - Brasil

Março – 2019


Agradecimentos

Aos meus pais, Osmar e Andrea, e irmã, Luhanna, por todo o incentivo

e carinho.

A Aline, que todos os dias, sem exceção, me apoiou, ouviu e aconselhou.

À Professora Doutora Luzia Trinca, minha orientadora, pela compre-

ensão, amizade, disposição, incentivo e oportunidade.

À Professora Doutora Berenice Camargo Damasceno e o Professor

Doutor Luciano Barbanti, meus orientadores da graduação, que me aconselharam

e auxiliaram.

Aos amigos do departamento de Bioestat́ıstica, entre eles os fun-

cionários, professores e estudantes.

Às Professoras Doutoras Júlia Maria Pavan Soler e Miriam Harumi

Tsunemi pelas contribuições.

À Fundação de Aperfeiçoamento de Pessoal de Nı́vel Superior - CAPES

pelo apoio financeiro, que foi essencial para o desenvolvimento de toda a pesquisa

realizada durante o mestrado.

Novamente, agradeço por todo suporte. Todos foram essenciais!


Sumário

Página

LISTA DE FIGURAS vi

LISTA DE TABELAS viii

RESUMO x

SUMMARY xii

1 INTRODUÇÃO 1

2 SELEÇÃO DE VARIÁVEIS EM MODELOS DE REGRESSÃO 3

2.1 Modelos Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2 Modelo de regressão Binomial . . . . . . . . . . . . . . . . . . . . . . . . 8

2.3 Diagnósticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4 Seleção de Variáveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.1 Métodos e critérios usuais . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.2 O método fence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.4.3 Ferramentas gráficas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3 MATERIAL E MÉTODOS 28

3.1 Aplicação 1: estudo sobre peso de recém-nascidos (RN) prematuros . . . 28

3.2 Aplicação 2: probabilidade de prenhez em inseminação artificial (IA) de

vacas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32


v

4 RESULTADOS 39

4.1 Aplicação 1: estudo sobre o peso de recém-nascidos (RN) prematuros . . 39

4.2 Aplicação 2: probabilidade de prenhez em inseminação artificial (IA) de

vacas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

5 CONSIDERAÇÕES FINAIS 69

ANEXOS 72

REFERÊNCIAS BIBLIOGRÁFICAS 76


Lista de Figuras

Página

1 Valores da correlação entre as regressoras quantitativas e a resposta peso,

dados da Aplicação 1 (n = 90) . . . . . . . . . . . . . . . . . . . . . . . . 30

2 Valores da correlação entre as regressoras quantitativas e da resposta por

lote, dados da Aplicação 2 (N = 65) . . . . . . . . . . . . . . . . . . . . 38

3 Gráfico de reśıduos do modelo de efeitos principais ajustado para a

variável - peso de RN: envelope quantil-quantil Normal (à esquerda) e

reśıduos em função dos valores ajustados (à direita). . . . . . . . . . . . 40

4 Gráfico de reśıduos do modelo de efeitos principais ajustado para a

variável - ln(peso) de RN: envelope quantil-quantil Normal (à esquerda)

e reśıduos em função dos valores ajustados (à direita). . . . . . . . . . . 40

5 Gráfico para diagnóstico de influência do ajuste do modelo de efeitos

principais para ln(peso) de RN. O raio dos ćırculos são proporcionais à

distância de Cook. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

6 Gráfico de inclusão de variáveis (VIP) em modelos de efeitos principais

para ln(peso) de RN, com destaque da curva de RV. . . . . . . . . . . . . 43

7 Gráficos do valor do componente de perda contra a dimensão do mo-

delo. Modelos de efeitos principais para ln(peso) de RN. Destaque para

as regressoras IG (à esquerda), GR (à direita) e pré natal (inferior). . . . 44

8 Gráfico de probabilidades de seleção de modelos em função do número de

parâmetros com destaque para a regressora pré-natal. Modelos de efeitos

principais para ln(peso) de RN, utilizando bootstrap ponderado. . . . . . 46


vii

9 Gráficos de probabilidades de seleção de modelos usando o método fence.

Modelos de efeitos principais para o ln(peso) de RN. . . . . . . . . . . . 47

10 Gráficos VIP dos termos de interação para os modelos 2 (à esquerda) e 3

(à direita) para o ln(peso) de RN. . . . . . . . . . . . . . . . . . . . . . . 49

11 Retas ajustadas para o ln(peso) de RN segundo o modelo 2. . . . . . . . 51

12 Retas ajustadas para o ln(peso) de RN segundo o modelo 3. . . . . . . . 52

13 Gráficos de diagnóstico para o ajuste do modelo de efeitos principais para

a probabilidade de prenhez sob IA, dados agregados completos (N= 72). 54

14 Gráficos de diagnóstico para o ajuste do modelo de efeitos principais para

a probabilidade de prenhez sob IA, dados agregados reduzidos (N= 65). . 55

15 Gráfico de inclusão de variáveis (VIP) para o modelo linearizado. . . . . 58

16 Gráfico de probabilidade de seleção de modelos em função do número de

parâmetros, com destaque para as regressoras, visualizadas da esquerda

para a direita e de cima a baixo, strr, linr, motr, alhr e progrr, dados

linearizados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

17 Gráfico de probabilidade de seleção de modelos, dados linearizados. . . . 62

18 Gráficos VIP para a interação dos modelos 3 (à esquerda) e 5 (à direita),

dados na Tabela 15, dados linearizados. . . . . . . . . . . . . . . . . . . . 63

19 Gráficos VIP para a interação dos modelos 4 (à esquerda) e 6 (à direita),

dados na Tabela 15, dados linearizados. . . . . . . . . . . . . . . . . . . . 64

20 Diagnóstico do modelo com efeitos da variável de blocos, str e mot, dados

reduzidos (N= 65). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

21 Diagnóstico do modelo com efeitos da variável de blocos, str, mot e lin,

dados reduzidos (N= 65). . . . . . . . . . . . . . . . . . . . . . . . . . . 67


Lista de Tabelas

Página

1 Tabela ANOVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2 Exemplo da ANODEV com um modelo que contém duas regressoras

cont́ınuas e sua interação . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Medidas descritivas das variáveis quantitativas, dados da Aplicação 1

(n = 90) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4 Medidas descritivas para as variáveis binárias, dados da Aplicação 1 (n =

90) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

5 Medidas descritivas das variáveis quantitativas (caracteŕısticas do

sêmem), dados da Aplicação 2 (N = 65) . . . . . . . . . . . . . . . . . . 36

6 Medidas descritivas da variável resposta por lote, dados da Aplicação 2

(N = 65) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

7 ANOVA do modelo de regressão ajustado para o peso de RN . . . . . . . 39

8 ANOVA do modelo de regressão ajustado com y = ln(peso) de RN . . . . 40

9 Modelos selecionados pelos métodos usuais segundo o critério, para o

ln(peso) de RN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

10 Modelos com probabilidade de seleção superior a 0, 20 visualizados na

Figura 8. Modelos de efeitos principais para ln(peso) de RN . . . . . . . 46

11 Modelos selecionados conforme aplicação do mplot para ln(peso) de RN . 48

12 Estimativas dos parâmetros e intervalos de confiança (95%) para os ajus-

tes dos modelos selecionados (Tabela 11) para ln(peso) de RN . . . . . . 49

13 Regressoras presentes nos modelos selecionados segundo os métodos usu-

ais de seleção de variv́eis, dados de prenhez sob IA agregados (N= 65) . . 56


ix

14 Regressoras presentes no modelo linearizado segundo os métodos usuais

para seleção de variáveis, dados de prenhez sob IA . . . . . . . . . . . . . 57

15 Modelos com destaque probabiĺıstico, conforme Figura 16, dados lineari-

zados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

16 Modelos com destaque probabiĺıstico, encontrados no gráfico de probabi-

lidades de seleção de modelos em função do número de parâmetros (não

apresentado) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

17 Parte preditiva dos modelos utilizando-se o mplot, dados linearizados . . 65

18 Regressoras presentes nos modelos finais, feita a seleção de variáveis com

a metologia usual do modelo linearizado . . . . . . . . . . . . . . . . . . 65

19 Estimativas e IC’s (95%) para os parâmetros do ajuste final, dados redu-

zidos (N= 65) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

20 Estimativas da razão de chances e respectivas IC’s (95%), para o modelo

incluindo a variável de blocos, mot, str e lin, para os dados de prenhez

sob IA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68


FERRAMENTAS GRÁFICAS NO PROCESSO DE SELEÇÃO DE

VARIÁVEIS

Autor: LUCAS RAGIOTTO

Orientadora: Prof. Dr. LUZIA APARECIDA TRINCA

RESUMO

Em problemas de regressão, na busca por um modelo parcimonioso, o pesquisa-

dor pode se deparar com adversidades, por exemplo, a existência de colinearidade

entre as regressoras, dificultando a seleção de variáveis. Dessa forma, com a imple-

mentação de ferramentas inspiradas nas propostas de Murray et al. (2013), Müller

& Welsh (2010) e Jiang et al. (2009) no pacote mplot (Tarr et al., 2018) no software

R, pode-se, gráfica e interativamente, estudar em detalhes a estabilidade e a im-

portância de inclusão de covariáveis para a construção de modelos. Neste trabalho,

medidas de estabilidade e probabilidade de inclusão de variáveis foram obtidas pelo

método bootstrap. Medidas resumo de qualidade do ajuste são baseadas no critério

de informação generalizado, que incorpora, como casos particulares, os critérios de

informação de Akaike e o Bayesiano, e reflete a perda (associada ao ajuste de um

modelo simplificado) mais uma penalização à complexidade do modelo. Ao aplicar


xi

a teoria de seleção de variáveis, utilizando as ferramentas gráficas no ajuste de um

modelo de regressão linear Normal e regressão Binomial, foi posśıvel reconhecer seu

potencial e utilidade no processo de formulação de modelos, no qual a incorporação

de conhecimento do especialista da área pode ser feita de maneira natural, já que o

processo não é automático. Isso é mais um diferencial em relação aos métodos usuais

de seleção de variáveis que também foram aplicados aos mesmos conjuntos de dados

para efeito de discussão.


GRAPHICAL TOOLS IN THE PROCESS OF VARIABLE

SELECTION

Author: LUCAS RAGIOTTO

Adviser: Prof. Dr. LUZIA APARECIDA TRINCA

SUMMARY

In regression analysis, the search of a parsimonious model can be difficult due to

collinearities among variables and other problems. Murray et al. (2013), Müller &

Welsh (2010) e Jiang et al. (2009) proposed tools for model stability and variable

inclusion plots that were refined and implemented in the mplot package of Tarr

et al. (2018), which allows interactive graphs and summaries of information relevant

to model building. Stability measures and the probability of variable inclusion are

obtained through bootstrapping. Goodness of fit measures are based on the generali-

zed information criterion, which includes as particular cases the Akaike and Bayesian

information criteria, given by a measure of loss of the fit and a penalization due to

model complexity. Applying the method to fit a Normal linear regression and a

Binomial regression revealed its great potential and usefulness for model building,

allowing expertise knowledge to be incorporated since the selection model is not au-


xiii

tomated. This is a further contrast to the usual selection methods which were also

applied to the same datasets in order to discuss the differences.


1 INTRODUÇÃO

Conhecimento cient́ıfico em diversas áreas é adquirido através da co-

leta e análise estat́ıstica de dados que, frequentemente, tem suporte num modelo.

Conforme bem colocado por Davison (2003), a formulação do modelo envolve sen-

satez, experiência, tentativa e erro. Não raramente o pesquisador pode se deparar

com vários modelos competitivos para o problema em mãos. Um prinćıpio utilizado

na escolha de um modelo é o da parcimônia que favorece modelos simples em de-

trimento de complexos quando as qualidades de seus ajustes são aproximadamente

equivalentes. No caso de modelos de regressão linear múltipla ou linear generalizado,

a aplicação do prinćıpio da parcimônia nem sempre é imediata ou fácil e existem di-

ferentes estratégias de visitas aos posśıveis modelos, assim como vários critérios de

medidas de qualidade do ajuste. Devido ao grande número de modelos que precisam

ser avaliados, as estratégias usuais são automatizadas e seguem um critério de or-

denação dos modelos pré-especificado. As estratégias mais populares são os métodos

tipo stepwise e o all subsets, que avalia todos os subconjuntos posśıveis. Embora

úteis, tais estratégias vêm sofrendo cŕıticas, principalmente pela não incorporação

do conhecimento do pesquisador especialista da área de aplicação. Outra cŕıtica é a

instabilidade dos métodos que podem levar a modelos diferentes quando há pequenas

perturbações nos dados. Uma terceira dificuldade é que, embora todas as estratégias

sejam baseadas na parcimônia, cada uma pode levar a um modelo distinto, prin-

cipalmente na presença de colinearidade entre as variáveis regressoras (Tarr et al.,

2018).

Com o intuito de proporcionar melhor entendimento sobre a im-

portância relativa aos diversos modelos estat́ısticos posśıveis para explicar determi-


2

nada caracteŕıstica ou fator, posteriormente chamada de resposta, Müller & Welsh

(2010) e Murray et al. (2013) propuseram os gráficos de inclusão de variáveis e medi-

das resumo da qualidade do ajuste, variando-se a penalidade imposta à complexidade

do modelo. Tarr et al. (2018) implementaram essas ideias e agregaram outras meto-

dologias no pacote mplot do R (R Core Team, 2017), com excelentes recursos para

construção de gráficos dinâmicos, permitindo o estudo de estabilidade de modelos e a

exploração detalhada da importância de cada variável regressora. Sob a metodologia

proposta por estes autores, a obtenção de medidas de estabilidade e probabilidades

de inclusão de variáveis é realizada via reamostragem pelo método bootstrap.

O objetivo dessa dissertação é revisar e explorar o uso desses recursos

gráficos na formulação de modelos de regressão, assim como levantar as posśıveis

limitações do pacote mplot, lançando mão de duas aplicações. Uma delas é no estudo

da relação entre o peso de recém-nascidos prematuros em função de outras covariáveis

relacionadas à gestação e a outra é no estudo da relação entre a probabilidade de

prenhes de vacas inseminadas artificialmente em função de diversas caracteŕısticas

do sêmen. No primeiro caso, usa-se o modelo linear Normal e no segundo um modelo

de regressão loǵıstica, que envolveu diversos desafios como a inclusão de um fator

de blocos, não submetido ao processo de seleção, e inclusão de termos de interação

entre regressoras. Para fins de discussão, os métodos de seleção de variáveis usuais

também foram aplicados em ambos conjuntos de dados.

O trabalho está organizado em caṕıtulos. No Caṕıtulo 2 apresenta-se

uma breve introdução aos modelos de regressão linear e loǵıstica e aos métodos de

seleção de variáveis, incluindo as ferramentas disponibilizadas no mplot. O Caṕıtulo

3 descreve os conjuntos de dados das duas aplicações e apresenta os métodos utili-

zados para a aplicação das ferramentas de seleção. Os resultados são apresentados e

discutidos no Caṕıtulo 4. As considerações finais se encontram no Caṕıtulo 5.


2 SELEÇÃO DE VARIÁVEIS EM MODELOS

DE REGRESSÃO

2.1 Modelos Lineares

A investigação, análise e interpretação de caracteŕısticas tem fator im-

portante nas decisões de diversas áreas de estudo. Devido à evolução tecnológica e

computacional, é posśıvel observar tais caracteŕısticas e coletar suas informações com

certa praticidade. Uma técnica estat́ıstica que auxilia na investigação e modelagem

dos problemas que envolvem essas caracteŕısticas, chamadas de variáveis, é definida

como análise de regressão. Seu objetivo é investigar a contribuição de certas variáveis

na variação de uma variável de interesse, chamada de resposta. Existe uma classe

ampla de modelos de regressão, sendo que a sub-classe dos modelos lineares oferece

grande potencial na condução deste tipo de estudo.

Para a aplicação desta teoria e construção de um modelo matemático

que visa relacionar as variáveis define-se n como o número de realizações do experi-

mento ou amostra, Y o vetor (n× 1) aleatório cujas observações serão consideradas

respostas e xj um vetor (n × 1) da j-ésima variável (j = 1, . . . , k) a qual deseja-se

analisar sua relação com a resposta, chamada de covariável, regressora ou variável

independente. As variáveis regressoras podem ser qualitativas ou quantitativas. Os

coeficientes da regressão ou parâmetros, β0, β1, . . ., βp−1, são constantes estimáveis,

no qual β0 é o intercepto e β1, . . ., βp−1 são os coeficientes das regressoras, alocadas

em um vetor β (p×1), com p = k+1. Cada βj (j = 1, 2, . . . , k) representa a variação

nos valores esperados da resposta conforme há mudança de uma unidade na regres-

sora xj, desde que os valores das outras estejam fixados. Os vetores das covariáveis


4

são organizados em um matriz X (n× p), com n > p, cuja primeira coluna tem seus

valores iguais a 1 no modelo com intercepto. Visto que vários fatores não observados

na amostra ou experimento contribuem na oscilação dos valores da resposta, deve-se

incluir um erro, neste caso denotado por ε (n × 1). Dado os componentes de um

modelo linear, tem-se sua estrutura, definida da forma

Y = Xβ + ε. (1)

Note que o termo “linear” é justificado pela forma que o vetor de parâmetros aparece

no modelo. Quando o modelo contém apenas uma regressora quantitativa ele é

conhecido também como regressão linear simples e com isso os valores de β a serem

estimados são o intercepto e o coeficiente da reta.

Quando o modelo linear especificado em (1) pressupõe que ε ∼

N(0, σ2I) ele é chamado de Normal, no qual I é a matriz identidade (n × n) in-

dicando que o erro se distribui normalmente com média e correlação 0 e variância

constante (homocedasticidade). Correlação nula e normalidade implicam em inde-

pendência. Consequentemente, tem-se

E(Y|X) = Xβ e V ar(Y|X) = σ2I, (2)

então, Y|X ∼ N(Xβ, σ2I).

Nessa perspectiva, deve-se estimar os parâmetros do modelo e para tal

pode-se utilizar de dois métodos que conduzem ao mesmo resultado. Demonstrações e

apresentação dessas metodologias de ajuste do modelo, propriedades dos estimadores

e exploração da qualidade do ajuste estão presentes em inúmeros livros de autores

clássicos como Montgomery et al. (2012), Draper & Smith (2014), Charnet et al.

(2015) e Seber & Lee (2012), por exemplo. Um dos métodos de estimação é o dos

“mı́nimos quadrados” que minimiza a soma dos erros ao quadrado. O outro método

maximiza a função de verossimilhança da amostra y que é dada pelo produto das

funções de densidades Normais de Yi avaliadas nos valores yi (i = 1, . . . , n). Sob as

pressuposições do modelo em (1), para ambos os métodos, o estimador de β é


5

β̂ = (X′X)−1X′Y,

para X de posto completo.

Por conseguinte, é imediato verificar que

E(β̂) = β

e

V ar(β̂) = σ2(X′X)−1,

ou seja, β̂ é não viciado para β e sua variância depende da matriz do modelo X e

da variância do erro σ2.

No entanto, vale notar que a propriedade de não tendenciosidade de β̂

só é válida se a especificação de E(Y|X) em (2) for correta. No caso do modelo real

ser Y = Xβ+ Aθ+ ε mas se na fase de ajuste, o segundo termo (Aθ) for ignorado,

então E(β̂) = (X′X)−1X′(Xβ+Aθ) = β+ (X′X)−1X′Aθ e β̂ só é não viciado para

β se (X′X)−1X′A = 0, ou seja, A e X são ortogonais, o que raramente acontece em

estudos observacionais.

O vetor dos valores ajustados ŷi’s, i = 1, . . . , n, é dado por

ŷ = Xβ̂ = X(X′X)−1X′y = Hy (3)

em que, y é o vetor de observações da variável resposta Y e H = X(X′X)−1X′, de

dimensão (n× n), conhecida como matriz de projeção ortogonal já que ela projeta,

ortogonalmente, y em ŷ. Logo, a diferença entre os vetores de valores ajustados e

observados de Y, definido como reśıduo, é ε̂ = y − ŷ = (I−H)y.

Para a estimação de σ2, utiliza-se a soma dos quadrados dos reśıduos

dada por

n∑
i=1

ε̂2i = SQres =
n∑
i=1

(yi − ŷi)2 = (y − ŷ)′(y − ŷ) = y′(I−H)y

que, também pode ser escrita como

SQres = y′y − β̂′X′y.


6

Sendo assim, a estimativa não viciada de σ2 é dada por

σ̂2 =
SQres

n− p
.

Logo, V̂ ar(β̂) = σ̂2(X′X)−1.

Uma vez estimados os parâmetros, deve-se buscar pela boa adequação

do modelo e identificar as regressoras de importância. Para tal, os pressupostos do

modelo devem ser satisfeitos e a literatura recomenda diversas técnicas exploratórias

e gráficas. Na Seção 2.3, um resumo das técnicas que auxiliam nesse diagnóstico é

apresentado.

Para testar a significância da regressão, ou seja, determinar se existe

relação linear entre a resposta e as covariáveis, utiliza-se o teste de hipótese global

associado à tabela ANOVA, o qual testa a hipótese nula H0 : β1 = .. = βp−1 = 0

contra a hipótese alternativa H1 : βj 6= 0, para algum j, e busca-se evidências de

que ao menos uma regressora é significante na explicação da resposta. Para sua

construção define-se a tabela ANOVA, apresentada na Tabela 1, cujas quantidades

são obtidas por

SQreg =
n∑
i=1

(ŷi − ȳ)2 = β̂′X′y − (
∑n
i=1 yi)

2

n
,

SQres =
n∑
i=1

(yi − ŷi)2 = y′y − β̂′X′y,

SQT =
n∑
i=1

(yi − ȳ)2 = y′y − (
∑n
i=1 yi)

2

n
,

nas quais SQreg ou Soma dos Quadrados da Regressão é a soma de quadrados dos

valores de ŷ, SQT é a soma de quadrados dos valores de y, ambas corrigidas pela

média, e SQres ou Soma dos Quadrados dos Reśıduos é a diferença entre as duas

primeiras. Sob normalidade e sob H0, cada soma de quadrados está associada a uma

distribuição qui-quadrado central com parâmetros espećıficos. Os graus de liberdade

de SQreg, SQres e SQT são p−1, n−p e n−1, respectivamente. Os quadrados médios

são obtidos pela divisão da soma dos quadrados pelos números de graus de liberdade

(G.L.) correspondentes. Assim, a estat́ıstica QMreg

QMres
, sob H0, segue a distribuição F


7

de Snedecor com p − 1 e n − p graus de liberdade. Se seu valor observado, F0,

for superior ao quantil de F1−α; p−1; n−p da distribuição de referência, então, H0 é

rejeitada, ao ńıvel α de significância.

Tabela 1. Tabela ANOVA

Variação Soma dos Quadrados G.L. Quadrado Médio F0

Regressão SQreg p− 1 QMreg QMreg/QMres

Reśıduo SQres n− p QMres

Total SQT n− 1

Testes de hipóteses também podem ser feitos para cada parâmetro, nos

quais a contribuição de uma regressora dado as outras no modelo é avaliada. Estes

se baseiam na estat́ıstica T = β̂j/ÊP (β̂j), que sob H0 segue a distribuição t-Student

com n − p graus de liberdade, em que ÊP (β̂j) é a estimativa do erro-padrão de β̂j,

dado pela raiz quadrada do valor da diagonal principal de V̂ ar(β̂) na posição j + 1

para j = 0, . . . , p.

Intervalos de confiança também são constrúıdos para os parâmetros.

O intervalo com 100(1−α)% de confiança para o j-ésimo coeficiente βj da regressão

(parâmetro) é expresso da forma

β̂j ± |tα/2, n−p| ÊP (β̂j),

com j = 0, . . . , p e tα/2, n−p o quantil α/2 da distribuição t-Student com n− p graus

de liberdade.

O intervalo de confiança de 100(1 − α)% para o valor esperado da

resposta quando as regressoras assumem um ponto x0, E(Y|x0), no espaço das re-

gressoras é dado por

ŷ0 ± |tα/2, n−p|
√
σ̂2x′0(X

′X)−1x0,

com x′0 = [1, x01, x02, . . . , x0k] sendo um ponto particular. O mesmo pode ser calcu-

lado para uma nova observação, y0, logo um intervalo de predição de 100(1− α)% é


8

dado por

ŷ0 ± |tα/2, n−p|
√
σ̂2(1 + x′0(X

′X)−1x0).

Todo ajuste de modelo de regressão deve ser submetido à análise de

reśıduos e diagnóstico afim de validar os pressupostos do modelo. Na Seção 2.3

apresenta-se as principais ferramentas gráficas para esses diagnósticos.

Todavia, o modelo linear Normal nem sempre é uma ferramenta viável

na construção de modelos de predição, pois a variável resposta pode ser caracteri-

zada por observações do tipo binárias ou de contagem e não se adequa às suposições

feitas no modelo em (1). Outra dificuldade encontrada é que a função que relaciona

a resposta esperada e Xβ não se apresenta de forma linear. Uma alternativa de

modelar os dados via modelos lineares é através da aplicação de transformações na

resposta de forma a aproximar à normalidade. Felizmente, na década de 1970, Nel-

der & Wedderburn (1972) introduziriam os modelos lineares generalizados (MLG),

que englobam diversas distribuições para a variável resposta dentro da famı́lia expo-

nencial. Sua teoria pode ser estudada com detalhes em Paula (2013), McCullagh &

Nelder (1989), Hosmer Jr et al. (2013), Chatfield et al. (2010) e Dobson & Barnett

(2018). Neste trabalho, devido à ilustração utilizada exigir um modelo que lida com

uma resposta do tipo binária, apresenta-se, na Seção 2.2, um resumo dos principais

tópicos para ajuste de um modelo deste tipo.

2.2 Modelo de regressão Binomial

É comum, principalmente nas áreas médicas e biológicas, a coleta de

dados cujas variáveis resposta são do tipo categórico. Existem diversas estratégias de

estudo para dados deste tipo, como, por exemplo, análise de tabelas de contingência

e alguns tipos de MLG’s. A acomodação de variáveis explanatórias quantitativas

é mais natural dentro da classe de modelos de regressão generalizado. No caso de

respostas binárias o modelo de regressão Binomial, um caso particular de MLG,

permite modelar a probabilidade de sucesso em um ensaio básico via um preditor.

Para seu ajuste os dados podem se apresentar de duas formas: a resposta na ob-


9

servação i é binária (Yi = 1 ou Yi = 0) ou a resposta é a contagem ou proporção

amostral do número de sucessos observada em um agregado e existem vários agre-

gados. A agregação é posśıvel quando existem repetições do conjunto de valores das

covariáveis. Em ambos os casos, o interesse está em modelar a probabilidade de

sucesso, e as estimativas dos parâmetros são equivalentes embora algumas medidas

de diagnóstico do ajuste possam não ter interpretações idênticas. Na exposição que

segue não se fará distinção entre os casos. A suposição é a de que os ensaios que

geram as variáveis binárias são independentes, mas não identicamente distribúıdas,

já que a probabilidade de sucesso depende das regressoras, assim como o parâmetro

de cada agregado que depende do número de repetições.

Da mesma maneira como estudado em modelos lineares, o modelo de

regressão Binomial tem o objetivo de modelar/estimar o valor esperado da variável

resposta dadas as observações das regressoras. Nesse caso a esperança da variável

binária se encontra entre os valores [0, 1] e a variância depende da sua esperança,

violando duas das suposições no modelo em (1) (o espaço do preditor linear é o

conjunto dos números reais e homogeneidade de variâncias). Para que o preditor

linear Xβ componha o modelo para E(Y ), sem que os valores da estimação dos

parâmetros sejam limitados, pode-se utilizar a função loǵıstica. A função loǵıstica

transforma o preditor no intervalo requerido, dado que sua forma é sigmoidal com

asśıntotas em 0 e 1, conforme explicado em Giolo (2017). Todavia, existem outras

funções que satisfazem tais necessidades para este tipo de dados, como a probito

e complemento log-log (Collett, 2002; Hosmer Jr et al., 2013; Agresti, 2012; Paula,

2013). Como a função loǵıstica é a mais popular, apenas ela será considerada na

sequência desta exposição.

Para uma variável aleatória Y , dado um conjunto de regressoras x,

tem-se que E(Y |x) = P (Y = 1|x) = π(x). Na regressão loǵıstica, π(x) é modelado

por

π(x) =
exp(β′x)

1 + exp(β′x)
. (4)

Agora β, o vetor de coeficientes da regressão, continua de dimensão p e x é um


10

vetor coluna cujos valores são 1 e demais valores observados para as k regressoras.

Enquanto π(x) fornece a probabilidade de ocorrer determinado evento (Y = 1|x),

1 − π(x) torna-se seu complemento, ou seja, a probabilidade de não ocorrência de

tal evento (Y = 0|x). Consequentemente, a variância de Y |x é π(x)(1 − π(x)),

explicitando a dependência da variância em x.

A razão entre as probabilidades de sucesso e fracasso é definida como

chance. Ao aplicar o logaritmo na chance, obtém-se a função logito, que está direta-

mente associada ao preditor linear, ou seja,

ln

[
π(x)

1− π(x)

]
= β′x. (5)

A equação (5), no contexto dos modelos lineares generalizados, é vista como uma

função de ligação canônica associada ao modelo Binomial.

A estimação dos parâmetros do modelo de regressão Binomial é feita

pelo método de máxima verossimilhança. A função de verossimilhança para o modelo

Binomial é o produto de funções de probabilidades Binomiais dada por

L(β; y) =
N∏
i=1

 ni

yi

 πyii (1− πi)ni−yi

em que πi = π(xi), N é o número de grupos ou agregados e ni (i = 1, . . . , N)

o número de repetições do agregado i. Nota-se que as probabilidades de sucesso πi

dependem dos β’s pela equação (4), então a função de verossimilhança é reescrita em

função do preditor linear β′xi. Com a maximização de L(β; y), ou equivalentemente

lnL(β; y), obtém-se os estimadores dos parâmetros. Nesse sentido, tem-se que

lnL(β; y) =
∑
i

ln

 ni

yi

+ yi ln(πi) + (ni − yi) ln(1− πi)

 =

=
∑
i

ln

 ni

yi

+ yi ln
(

πi
1− πi

)
+ ni ln(1− πi)

 =

=
∑
i

ln

 ni

yi

+ yi ηi − ni ln(1 + eηi)




11

com ηi = ln
[

πi
1−πi

]
(5)
= β′xi =

∑k
j=0 βjxji e x0i = 1 para todo valor de i. A maxi-

mização ocorre quando as derivadas de lnL(β; y) são igualadas a zero, em relação a

cada βj, fornecendo o estimador do parâmetro j. Devido à presença de covariáveis

esse processo de maximização não tem solução fechada, sendo então necessário a uti-

lização de um método iterativo, por exemplo, o algoritmo conhecido como método

de Score de Fisher, para obter β̂ (Collett, 2002). Desenvolvendo essa metodologia,

chega-se à expressão usual ao do método de mı́nimos quadrados reponderados. Para

a função loǵıstica, os valores da variável dependente zi = ηi+(yi−niπi)/[niπi(1−πi)]

são ajustadas pelas k regressoras em cada iteração, usando os pesos da diagonal prin-

cipal da matriz V (N ×N), expressos como

υi = niπi(1− πi), (6)

com i = 1, . . . , N . Nota-se que os elementos da diagonal de V correspondem à

variância de uma variável com distribuição Binomial (ni, πi). Portanto, após um

certo número de iterações, ao ocorrer a covergência, obtém-se os valores de β̂ (Collett,

2002).

Uma medida resumo da qualidade do ajuste de um determinado modelo

M é dada pela função desvio (Giolo, 2017), definida por

DM = 2 {lnLS − lnLM}

na qual LS é a verossimilhança do modelo saturado, ou seja, que possui tantos

parâmetros quantas observações existirem e LM é a verossimilhança do modeloM que

apresenta menor número de parâmetros do que S, ambas avaliadas nos estimadores

de máximas verossimilhanças. Fazendo uma relação com o modelo linear Normal,

naquele caso, DM coincide com a SQres(M) (Soma dos Quadrados dos Reśıduos

segundo o modelo M).

Conforme Collett (2002) esclarece, o estimador da função desvio DM ,

por ser obtida da estat́ıstica ln da razão de verossimilhanças, segue, assintotica-

mente, a distribuição qui-quadrado com N − pM graus de liberdade (pM = número


12

de parâmetros do modelo M) no caso de dados agregados. Assintoticamente quer

dizer que ni → ∞, ou seja, o número de repetições por agregado deve ser muito

grande (i = 1, 2, . . . , N) para satisfazer a condição. Assim, muitos pesquisadores

utilizam DM como uma medida da falta de ajuste do modelo M .

No caso de dados binários Collett (2002) mostra que a função desvio

não é um critério adequado como medida de falta de ajuste já que ela depende das

observações apenas através dos estimadores π̂(xi), e que não é posśıvel associá-la à

distribuição qui-quadrado.

Por meio da função desvio, é posśıvel generalizar a tabela ANOVA

(Tabela 1) para os GLM’s que agora é chamada de tabela ANODEV (Tabela 2).

Sua utilidade está na avaliação do(s) efeito(s) de uma (ou mais) regressora(s) na

presença de outra(s), sequencialmente, ou seja, partindo-se do modelo mais simples

para algum mais complicado. Ilustra-se, na Tabela 2, o caso no qual um efeito é

avaliado dado as outras regressoras presentes no modelo anterior. Sua contribuição é

medida pela diferença entre os valores da função desvio dos dois modelos. No caso de

modelos encaixados, ou seja, o mais simples é um caso particular do mais complexo,

pode-se usar o teste da razão de verossimilhanças (Collett, 2002).

Tabela 2. Exemplo da ANODEV com um modelo que contém duas regressoras

cont́ınuas e sua interação

Modelo Diferença G.L. Função desvio Diferença

Nulo DN

X1 1 DX1 DN −DX1

X2|X1 1 DX1,X2 DX1 −DX1,X2

X1 ×X2|X1, X2 1 DX1,X2,X1×X2 DX1,X2 −DX1,X2,X1×X2

Assim, para testar as hipóteses H0 : βj = 0 contra H1 : βj 6= 0 usa-

se o resultado de que, sob H0, a diferença entre as duas funções desvio (Tabela

2) tem distribuição qui-quadrado com 1 grau de liberdade (χ2
1). Logo, se o valor

observado da diferença for maior do que o quantil χ2
1,1−α, ao ńıvel de significância α,


13

a hipótese nula é rejeitada. O teste pode ser aplicado para um grupo de parâmetros,

simultaneamente, tal que o número de graus de liberdade é dado pela diferença no

número de parâmetros dos dois modelos.

Para a construção de intervalos de confiança (IC’s), seja o estimador

da variância do estimador β dada por V̂ ar(β̂) = (X′V̂X)−1. Para cada β̂j (j =

0, 1, . . . , k) em particular, seu erro-padrão é a raiz quadrada do (j+1)-ésimo elemento

da diagonal principal de V̂ ar(β̂), denotado como ÊP (β̂j). Assim, o IC para β̂j,

baseado na estat́ıstica de Wald (Hosmer Jr et al., 2013), é expresso como

β̂j ± z1−α/2 × ÊP (β̂j), (7)

em que z1−α/2 é o quantil de ordem 1 − α/2 da distribuição Normal padrão. O IC

para ln
[

π(x)
1−π(x)

]
= β′x é obtido através de

β̂′x± z1−α/2 × ÊP (β̂′x),

em que ÊP (β̂′x) = [x′(X′V̂X)−1x]1/2 com x′ = (1, x1, . . . , xp) sendo o vetor de

valores das regressoras fixados. Por conseguinte, os IC’s aproximados para π(x) são

encontrados via transformação inversa.

Uma maneira de interpretar os parâmetros através do ajuste é com

a utilização da razão de chances ajustadas. Esta medida de associação em análise

de tabelas de contingência é chamada de razão de chances (odds ratio). Para cada

regressora xj pode-se calcular a razão de chances que mede o efeito do aumento de

uma unidade sobre a chance de sucesso, dado que os valores das outras regressoras

estão fixadas. Logo, o estimador de razão de chances ajustada (ÔR) é dado por

ÔR = exp(β̂).

Seja x-j o vetor de valores das regressoras exceto a regressora xj e β-j o vetor β

excluindo-se βj. Da equação (5) tem-se que

π(xj|x-j)

1− π(xj|x-j)
= exp(β′-jx-j + βjxj)


14

e
π(xj + 1|x-j)

1− π(xj + 1|x-j)
= exp(β′-jx-j + βj + βjxj)

tal que

OR(xj + 1|x-j) =

π(xj+1|x-j)
1−π(xj+1|x-j)
π(xj |x-j)

1−π(xj |x-j)

= exp(βj)

cujo estimador é

ÔR(xj + 1|x-j) = exp(β̂j).

Dado o IC para βj (equação 7), e usando a transformação inversa é posśıvel obter

intervalos de confiança aproximados para OR. Como OR é uma razão de chances, se

o valor de 1 estiver inclúıdo no intervalo significa que xj não tem efeito significante.

Interpretações similares podem ser feitas para regressoras que são qualitativas. Para

mais detalhes sobre estas e outras interpretações no ajuste de regressão Binomial

e análises de tabelas de contingência veja Giolo (2017) e Hosmer Jr et al. (2013).

Na Seção 2.3 apresentam-se os principais diagnósticos para o modelo de regressão

loǵıstica e o linear.

2.3 Diagnósticos

A busca pela boa adequação do modelo gera diversas teorias e conse-

quentemente ferramentas, muitas vezes gráficas, para avaliar e compreender qual o

peso de cada observação ou um conjunto delas na inferência dos parâmetros. Nesse

sentido, a identificação de pontos que se destacam, dado um componente avaliativo, é

de grande importância para determinar seu impacto nos aspectos gerais do modelo.

Na teoria dos MLG, cada função de probabilidade considerada em estudo, possui

maneiras diferentes de mensurar o impacto das observações feitas e também avaliar

a qualidade do ajuste, consequentemente, existem diferentes formas de construir me-

didas de avaliação dessas observações. Visto que as aplicações feitas neste trabalho

foram constrúıdas com base nos modelos de regressão linear e Binomial, apresenta-se

os principais diagnósticos para os mesmos.


15

Ao considerar um modelo linear Normal deve-se, inicialmente, verificar

se os pressupostos feitos para ε da equação (1) (basicamente normalidade, homoce-

dasticidade e não correlação) são satisfeitos, com isso considera-se sua estimativa

ε̂, denominada por reśıduo bruto. Nesse sentido, definem-se dois outros tipos de

reśıduos com base nos reśıduos brutos (Montgomery et al., 2012), o padronizado

dado por

ri =
ε̂i

S
√

1− hii
,

com S =
√
σ̂2 e hii sendo o i-ésimo valor da diagonal da matriz H (equação 3) e o

reśıduo studentizado, definido como

ε̂∗(−i) =
ε̂i

S(−i)
√

1− hii
,

com S(−i) sendo S resultante do modelo ajustado sem a observação i. Uma vez que

as pré-suposições do modelo estão de acordo, este reśıduo segue uma distribuição t

de Student com n− p− 1 graus de liberdade.

Definidos os reśıduos, pode-se verificar o pressuposto da normalidade

com um gráfico dos reśıduos padronizados contra os quantis teóricos de uma Normal

padrão ou os reśıduos studentizados contra os quantis teóricos da distribuição t. Sua

construção é feita ao plotar os reśıduos em ordem crescente contra os quantis teóricos

de uma Normal ou de uma distribuição t de Student. Para melhor comparação

entre os elementos que compõem o gráfico, constrói-se, a partir de simulações, uma

banda de confiança denominada envelope (Atkinson, 1981). Nesse gráfico, deve-se

esperar que os pontos estejam dentro dessa faixa criada, pois a presença de uma

certa quantidade de pontos fora dessa banda pode indicar inadequação do modelo

ou presença de observações estranhas/at́ıpicas ou ainda de posśıvel influência no

ajuste.

Para verificação do pressuposto de variância constante (homocedasti-

cidade), deve-se plotar os reśıduos contra os valores ajustados (ŷi’s). Nele, deve-se

buscar por pontos espalhados pelo gráfico sem observar tendências e com faixa de va-

riação aproximadamente constante. Uma forma de amenizar tendências crescentes de


16

variabilidade é utilizando alguma transformação da variável resposta. Uma famı́lia

de transformações bastante conhecida é a famı́lia Box-Cox (Box & Cox, 1964), mais

detalhes sobre ela e outras transformações podem ser estudadas em Montgomery

et al. (2012) e Faraway (2016). Falta de ajuste ou necessidade de inclusão de outros

termos no modelo podem ser investigados plotando-se os reśıduos em função de cada

regressora. O padrão esperado é o mesmo descrito acima.

Três são os componentes de pontos do banco de dados que podem

comprometer o modelo. São eles, os pontos at́ıpicos/extremos (outlier), ponto de

alavanca (Hoaglin & Welsch, 1978) e pontos de influência. Um ponto outlier é

aquele que apresenta alto reśıduo, provocado usualmente por valores extremos na

resposta. Para a identificação desses pontos pode-se utilizar o gráfico envelope e

identificar os valores fora da faixa simulada. Outra maneira é plotar os valores dos

reśıduos padronizados ou studentizados contra a ordem das observações, buscando

por pontos fora do intervalo −2 e +2, visto que uma vez padronizados, espera-se que

95% deles estejam dentro da faixa referida, segundo as propriedades da densidade

Normal padrão.

Um ponto de alavanca é aquele que se destaca em relação ao centro

do espaço das regressoras. Para identificá-los plotam-se os pontos de hii, que são os

valores da diagonal principal da matriz H, contra a ordem das observações, buscando

por valores maiores do que 2p/n, conforme recomendado por Paula (2013) e Faraway

(2016). Em geral, deve-se avaliar o motivo pelo qual os pontos se destacam e têm

seu comportamento diferente dos demais, sendo assim, busca-se por alternativas que

possam tratá-los ou até mesmo reavaliar o modelo constrúıdo.

Além de observações at́ıpicas, certos conjuntos de dados não raramente

podem apresentar pontos indesejáveis de alta influência no ajuste. Um ponto é dito

influente quando o ajuste ou estimativa de um ou mais parâmetros da regressão sofre

alteração destacada quando o ponto é removido da análise. Assim, a presença de

pontos influentes podem levar a estimativas viciadas e pode, inclusive, comprometer

a qualidade do ajuste e consequentemente afetar o processo de seleção de variáveis.


17

É importante esclarecer que um ponto considerado outlier não necessariamente tem

influência no modelo. Nesse sentido, utiliza-se de ferramentas para a identificação de

posśıveis pontos influentes. A maneira mais simples é utilizando as medidas de hii

(pontos de alavanca), já que as propriedades da matriz H estabelece que no modelo

“ideal”os valores de hii deveriam ser todos iguais a p/n (Montgomery et al., 2012).

Entre as medidas mais conhecidas para identificar pontos influentes

são os DFBETAS e DFFITS. O DFBETAS é uma medida que indica quanto cada

coeficiente da regressão muda se a i-ésima observação for removida dos dados e o

modelo reajustado, enquanto os DFFITS mede a influência dos valores ajustados.

Detalhes sobre essas estat́ısticas podem ser consultadas em Montgomery et al. (2012)

e Faraway (2016). No sentido de resumir a informação contida nos DFBETAS, pois

a mesma reflete qual a influência do ponto em cada coeficiente, utiliza-se a distância

de Cook (Cook, 1977), definida da forma

r2i
p
× hii

1− hii
.

com i = 1, . . . , n. Em geral, valores que se distanciam dos demais devem ser ana-

lisados, mas como referência busca-se pelos pontos próximos de 1, logo, um gráfico

dessas medidas contra cada observação pode auxiliar nesse procedimento.

Para um modelo de regressão loǵıstica, as medidas de diagnóstico pro-

postas para o modelo de regressão linear Normal sofrem adaptações, e diversos tipos

de reśıduos são considerados. Tomando dados agregados e agora yi sendo o número

de sucessos no agregado i, o reśıduo componente do desvio (Paula, 2013) é dado por:

di = sinal(yi − ŷi)
{

2yi ln

(
yi
ŷi

)
+ 2(ni − yi) ln

(
ni − yi
ni − ŷi

)}1/2

,

em que ŷi é a estimativa do número de sucessos e ni o número de realizações repetidas

no grupo i. O reśıduo de Pearson é dado por

δi =
yi − ni π̂i√
ni π̂i (1 − π̂i)

,

com i = 1, .., N e π̂i = π̂i(x) é a estimativa da proporção esperada em cada lote i.


18

Para padronizá-los, considera-se a matriz

H = V1/2X(X′VX)−1X′V1/2 (8)

sendo V (N×N) a matriz de pesos utilizados no ajuste do modelo (conforme equação

6). Sendo assim, ao dividir di e δi por
√

(1− hii), com hii sendo os elementos da

diagonal principal da matriz H, obtém-se os reśıduos componente do desvio padro-

nizados (rdi) e os reśıduos de Pearson padronizados (rδi), respectivamente. Outro

reśıduo que auxilia na identificação de pontos at́ıpicos é chamado de reśıduo da

verossimilhança (Collett, 2002) e definido da forma

rLi
= sinal(yi − ŷi)

√
hii r2δi + (1− hii) r2di .

Note que rLi
é uma combinação de r2δi e r2di . Se hii for pequeno então rLi

é similar a

r2di .

Para a identificação de outliers no modelo de regressão Binomial, dado

que os reśıduos componente do desvio se distribuem normalmente (Giolo, 2017; Col-

lett, 2002), um gráfico envelope pode ser constrúıdo. Logo, busca-se por pontos

fora da faixa simulada. Outra maneira é plotar rdi e/ou rLi
contra a ordem das

observações que deve conter, em sua maioria, pontos no intervalo de −2 e +2, pelos

mesmos motivos citados no modelo linear Normal.

Na busca por pontos influentes, uma das formas de medir, aproxima-

damente, a influência da i-ésima observação no ajuste geral do modelo, é avaliar

o valor de r2Li
, eles medem a diminuição no valor do desvio do modelo quando a

observação i é removida dos dados e, neste caso, é importante avaliar valores que

se destacam dos demais. Outra medida que auxilia na busca desses pontos, é pelos

elementos da diagonal principal da matriz H (equação (8)), os hii’s. São chamados

de pontos de alavanca os valores que se destacarem no gráfico de hii contra os ı́ndices

das observações. Neste caso, procura-se por pontos que sejam maiores do que 2p/N

(Collett, 2002).

Da mesma forma vista em um modelo linear Normal, a distância de

Cook também é uma ferramenta útil na identificação de pontos influentes para o


19

modelo de regressão Binomial (Collett, 2002), sendo definida da forma

ψi =
r2δihii

p(1− hii)
, i = 1, . . . , N.

Novamente, busca-se por pontos que se destacam em um gráfico de ψi contra a ordem

das observações. Para efeitos de comparação os mesmos devem ter valores pequenos,

enquanto valores próximos de 1 devem ser avaliados. De forma geral, é importante

tentar descobrir os motivos que levam pontos dos dados possúırem valores fora do

padrão, embora na prática isso nem sempre é viável, para que se possa tratá-los.

Mudanças na distribuição da variável resposta, inclusão de termos das regressoras,

entre outras possibilidades, podem ser utilizadas com o intuito de diminuir os pro-

blemas diagnosticados no ajuste. Outras maneiras de medir a influência, bem como

análise da adequação e procedimentos para posśıveis transformações nas variáveis,

podem ser estudados com detalhes em Collett (2002) e Paula (2013).

2.4 Seleção de Variáveis

2.4.1 Métodos e critérios usuais

Os prinćıpios da modelagem estat́ıstica preconizam o ajuste de um

modelo parcimonioso aos dados observados (Montgomery et al., 2012; Draper &

Smith, 2014), ou seja, um modelo com número reduzido de parâmetros de forma que

apenas as variáveis preditoras de importância sejam inclúıdas, porém mantendo-se a

utilidade prática do modelo. Com isso, a interpretação e o uso posterior do modelo

parcimonioso torna-se relevante pela sua simplicidade, exigindo menos tempo e gasto

financeiro com posśıveis predições, aspectos muitas vezes levados em consideração

por pesquisadores e proprietários dos dados. A busca por um modelo que se adéque

às exigências mencionadas pode ser realizado via teoria de seleção de variáveis. A

bibliografia especializada apresenta várias técnicas e critérios clássicos que podem ser

utilizados para selecionar variáveis na modelagem estat́ıstica, seja no modelo linear

Normal, generalizado ou outros envolvendo diversas regressoras.


20

A dificuldade encontrada no processo de seleção deve-se à existência,

na grande maioria, de colinearidade entre as regressoras causando instabilidade no

processo de seleção. Esse problema é mais frequente nos estudos observacionais já

que nos experimentos a ortogonalidade entre as regressoras pode ser obtida via pla-

nejamento, o que simplifica a seleção dos efeitos significantes devido à independência

entre as estimativas dos diversos efeitos. Os métodos revisados neste trabalho são

aplicáveis apenas aos casos nos quais o número de observações (n) é maior do que o

número de parâmetros posśıveis no modelo (p). No entanto, o desenvolvimento tec-

nológico da atualidade permite mensuração de um número muito grande de variáveis

numa mesma unidade ou indiv́ıduo, gerando o caso de n << p, como é o caso de

estudos epidemiológicos genéticos (Wasserman & Roeder, 2009). Nesse contexto,

o problema de seleção de variáveis é muito complexo, exigindo alternativas de es-

timação que incorporam penalidades às estimativas dos parâmetros, forçando-as a

se aproximarem do valor nulo (Tibshirani, 1996).

A busca por um modelo parcimonioso leva à seguinte problemática,

poucas regressoras no ajuste pode conduzir a estimativas e predições tendenciosas e

um grande número delas conduz à instabilidade e consequente aumento da imprecisão

na estimação e na predição de novas observações, além de poder resultar em modelos

superajustados. Um modelo superajustado representa os dados quase que fielmente,

explicando até parte da variabilidade aleatória intŕınseca, o que não é desejável desde

que dados amostrais ou experimentais, pela sua natureza, apresentam imperfeições.

Assim, modelos superajustados são perigosos para novas previsões.

A colinearidade também atrapalha o processo de seleção de variáveis

levando à instabilidade na escolha do modelo final. Além do mais, existe uma vari-

edade de métodos e critérios para a seleção. É sabido, por exemplo, que variações

dos métodos do tipo stepwise nem sempre convergirão para o mesmo modelo, e,

mesmo fixando-se o método, pequenas perturbações nos dados também podem con-

duzir a modelos distintos. Outro aspecto é a existência de diferentes critérios de

ordenação da qualidade do ajuste que podem levar a certa insegurança e dificuldade


21

para escolher o melhor modelo num grande conjunto de possibilidades.

Para introdução dos critérios de seleção é importante estabelecer a

notação que será utilizada. Assume-se n observações independentes da variável Y

e uma matriz X n × p (p < n) de posto completo p cujas colunas estão indexadas

pelos elementos do conjunto α∗ = {1, 2, .., p}. Seja α um subconjunto qualquer de

α∗ tal que contrói-se a matriz Xα de dimensão n × pα e de posto pα com pα ≤ p.

A matriz Xα é a parte preditora do modelo α e βα são os parâmetros de regressão.

O problema de seleção é o de descobrir α que resulta no modelo parcimonioso que

explica Y em função de Xα. Os critérios de seleção recorrentes são o critério de

informação de Akaike (AIC), proposto por Akaike (1973), e o critério de informação

Bayesiano (BIC), proposto por Schwarz et al. (1978). Visto que AIC e BIC se ba-

seiam na mesma equação para impor alguma penalidade à complexidade do modelo,

Konishi & Kitagawa (1996) ampliou essa teoria introduzindo o critério de informação

generalizado, dado pela equação

GIC(α, λ) = Q̂(α) + λ pα (9)

sendo Q̂(α) o componente que mede a falta de ajuste do modelo α, por exemplo, soma

de quadrados dos reśıduos ou −2 ln(βα, φ; y), em que φ é o parâmetro de dispersão

(φ = 1 para binomial), cuja complexidade é medida por pα (posto do modelo α)

e λ é o multiplicador da penalidade. Quando Q̂(α) = −2 ln(βα, φ; y) e λ = 2 ou

λ = ln(n), obtém-se o AIC e BIC, respectivamente. Na comparação de dois ou mais

modelos, sob os critérios de informação, aquele com valor menor é prefeŕıvel já que

Q̂(α), de certa forma, mede o distanciamento do modelo α ao “verdadeiro” modelo.

Para definições de outros critérios veja Faraway (2016).

Para proceder à uma seleção utilizando as abordagens usuais é ne-

cessário decidir pelo critério e pelo procedimento de visitação dos diversos modelos.

O procedimento de visitação dos modelos é referido como método de seleção. Uma

classe de métodos para seleção de regressoras bastante usual é chamada de stepwise,

que inclui três abordagens: forward, backward e stepwise (Montgomery et al., 2012).

A forward adiciona regressoras ao modelo uma a uma, partindo-se de um modelo


22

nulo, a backward exclui regressoras uma a uma ao considerar, inicialmente, modelo

completo, ou seja, aquele com número máximo de regressoras dispońıveis, enquanto

a stepwise combina os dois procedimentos, começando pelo modelo completo. No

procedimento forward, uma vez que determinada variável é inclúıda, ela segue no

modelo até o final do processo. Já no backward, uma vez que a variável é exclúıda

do modelo, ela não volta a ser testada.

A inclusão ou exclusão das variáveis é feita seguindo algum critério, que

informa se a regressora é importante no modelo que está sendo constrúıdo. O critério

pode ser um teste de significância do efeito da regressora, por exemplo t, F , Wald,

ou R2 (Montgomery et al., 2012; Charnet et al., 2015), ou critérios de informação

como AIC e BIC (Faraway, 2016), entre outros. No caso do critério ser um teste,

recomenda-se que ao ńıvel de significância estipulado não seja muito pequeno para

dar maior chance de investigação de efeito de regressoras na presença de outras (veja

Box & Draper (1987)). Caso o ńıvel de significância for muito exigente (pequeno) a

entrada de variv́eis (forward) pode ser rara, assim como a sáıda (backward) muito

frequente. Combinando, por exemplo, o método forward e o critério AIC, é feito a

comparação dos valores AIC para cada modelo posśıvel a cada passo. Dessa maneira,

ao iniciar, a variável selecionada será aquela cujo ajuste produzir o menor AIC entre

todos, inclusive do modelo nulo. O procedimento se repete até que o menor valor de

AIC seja maior do que o valor do modelo constrúıdo até aquele passo.

Além dos métodos tipo stepwise para seleção de modelos pode-se

também proceder uma busca exaustiva na qual todos os subconjuntos de regres-

soras posśıveis de um modelo completo são visitados e comparados conforme um

critério pré estipulado, por exemplo, o valor de AIC. Considerando que β0 esteja

contido em todos os modelos, para um total de p− 1 regressoras, existem 2p−1 mo-

delos posśıveis para estimar e examinar. Portanto, o método é viável nos estudos

com poucas variáveis independentes mas torna-se custoso conforme p aumenta, em-

bora, com a evolução dos métodos e equipamentos computacionais tal custo tende a

descrescer.


23

2.4.2 O método fence

Jiang et al. (2008) introduziram o método fence, ou cerca, para auxiliar

na seleção de modelos dentro da classe de modelos mistos. Um modelo é dito ser de

efeitos mistos quando, além da parte fixa que compõe o modelo de regressão usual,

ele inclui também outros efeitos, além dos erros, que são supostos variarem aleato-

riamente. Tal classe de modelos encontra aplicações nos problemas com dados de

medidas repetidas ou que exibem algum tipo de agrupamento entre as observações,

devido à forma de amostragem. Referências clássicas sobre esses modelos são, por

exemplo, Pinheiro & Bates (2000) e Littell et al. (2006). Motivados pela dificul-

dade da definição ou interpretação de critérios baseados na equação (9) no contexto

de modelos mistos, juntamente com a evolução computacional, Jiang et al. (2008)

introduziram o método fence ou cerca, posteriormente simplificado em Jiang et al.

(2009), com variações em Jiang (2014).

O método fence aplica a ideia de cercar subgrupos de modelos com

propriedades similiares. O processo é baseado na condição

Q̂(α)− Q̂(α∗) ≤ c, (10)

no qual Q̂(α) é uma medida de falta de ajuste do modelo α em análise, Q̂(α∗) é

uma medida de falta de ajuste do modelo completo α∗, que contém o efeito de todas

as regressoras e c é uma constante que determina um ponto de separação ou de

corte. Vale ressaltar que os modelos candidatos são aqueles encaixados em α∗, ou

seja, seus parâmetros formam um subconjunto dos parâmetros de α∗. Jiang et al.

(2009) apresenta o método em detalhes e Jiang (2014) faz uma revisão indicando

extensões e variações do método, inclusive para problemas com p >> n, nos quais o

método fence foi aplicado previamente para reduzir a dimensionalidade previamente

às técnicas de regularização.


24

2.4.3 Ferramentas gráficas

Recentemente Tarr et al. (2018) desenvolveram uma ferramenta com-

putacional e gráfica para assessorar o pesquisador no processo de seleção de variáveis,

influenciado fortemente nos trabalhos de Murray et al. (2013), Müller & Welsh (2010)

e Jiang et al. (2009), disponibilizada no pacote mplot (Tarr et al., 2018) para o soft-

ware R (R Core Team, 2017). Esse pacote permite a construção de gráficos que se

baseiam na teoria do critério de informação generalizado, no método fence e estudos

de estabilidade de modelos. Para estabilidade é utilizado o método de reamostra-

gem de bootstrap ponderado exponencialmente que consiste na reponderação das

observações usando pesos simulados de uma função de densidade exponencial com

média 1 (Murray et al., 2013). Existem diversas variações do método bootstrap, inclu-

sive com pesos, conforme revistos em Gotwalt et al. (2018). Com base nos resultados

de simulações de Jin et al. (2001) conclui-se que o uso de pesos dados por realizações

de uma variável aleatória com média e variância iguais a um oferece estimadores com

boas propriedades. Mais detalhes do método são dados em Murray et al. (2013). Em

Müller & Welsh (2010) define-se que um procedimento/método de seleção de modelo

é considerado instável quando, para uma pequena variação na penalidade (λ) na

equação (9), pode-se selecionar modelos de diferentes dimensões.

O mplot oferece ferramentas para o aux́ılio no processo de seleção de

modelos no caso linear e linear generalizado. Dois comandos principais calculam as

medidas para a investigação da contribuição das variáveis e estabilidade. Uma vez

criado os objetos, pelas informações obtidas nos dois comandos, é posśıvel a cons-

trução de cinco tipos de gráficos: o gráfico para exploração do espaço dos modelos,

o gráfico de estabilidade, o gráfico de inclusão de variável e os dois gráficos baseados

no método fence.

O gráfico para exploração do espaço dos modelos é um diagrama do

componente de perda, −2 ln(βα, φ; y), em função da dimensão do modelo. Para

cada dimensão, todos os modelos posśıveis, incluindo o intercepto, são ajustados.

No gráfico é permitido o uso de uma legenda para destacar os modelos que in-


25

cluem/excluem determinada covariável, proporcionando uma visualização informa-

tiva da contribuição da covariável ou ajuste, ou seja, na diminuição de−2 ln(βα, φ; y).

Uma outra versão desse gráfico permite o estudo de estabilidade de modelos reali-

zado por reamostragens bootstrap. Ele concede a proporção de vezes que o modelo,

para dada dimensão, apresentou o valor mı́nimo do componente de perda entre as

repetições realizadas. Cada ponto no gráfico, que plota −2 ln(βα, φ; y) versus a

dimensão do modelo, é representado por um ćırculo com diâmetro proporcional à

frequência de vezes que o modelo foi selecionado como sendo o melhor nas reamos-

tras da dimensão considerada.

O gráfico de inclusão de variáveis (VIP) também é constrúıdo com base

nas informações do processo de reamostragem. Para isso, toma-se B reamostras con-

siderando um valor de λ, que varia em pequenos intervalos, e calcula-se a proporção

de vezes que esta variável estava presente no modelo final, sendo esse o que tem

menor valor no critério de informação generalizado (equação (9)) para a penalidade

λ definida inicialmente. O gráfico VIP é um diagrama de linhas traçadas pelas pro-

babilidades em função de λ. Espera-se que para valores pequenos, a maioria das

regressoras estejam presentes nos ajustes e conforme há o aumento da penalidade,

poucas pertençam a ele. Nesse processo é posśıvel incluir uma variável redundante

(RV), simulada a partir de uma distribuição Normal padrão. Seu objetivo é fornecer

uma curva de base em VIP que serve de base de comparação para as outras variáveis.

Curvas próximas ou abaixo da RV podem ter sido inclúıdas por acaso e indicam re-

gressoras de baixa relevância no modelo. Ao considerar os dois primeiros gráficos

mencionados, os modelos que contêm RV não são válidos.

Outras alternativas para o estudo de estabilidade de modelos são ofe-

recidos pelo mplot, através da aplicação do método fence simplificado (Jiang et al.,

2009), que foi proposto inicialmente para os modelos mistos, mas adaptado para

GLM’s por Tarr et al. (2018). Ilustrando o caso do modelo linear Normal, primei-

ramente escolha o conjunto de valores de c e ajuste o modelo completo α∗, obtendo

Q̂(α∗). Para cada valor de c segue-se os passos:


26

1. Reamostre (B vezes) parametricamente ε ∼ N(0, σ̂2
α∗), em que σ̂2

α∗ é o quadrado

médio dos reśıduos uma vez ajustado o modelo α∗ .

Para cada reamostra, identifique o modelo de menor dimensão α̂(c) que satisfaz

a condição na equação (10). Caso exista mais de um modelo do mesmo tamanho

dentro da cerca, utilize o que tem menor Q̂(α̂(c)).

Calcule a probabilidade emṕırica de selecionar o modelo α, dada por p∗(α, c) =

P ∗(α̂(c) = α) = #(α̂(c)=α)
B

;

2. Obtenha a probabilidade de seleção global p∗(c) para dado c, que é dada por

p∗(c, α) máximo. Nesse passo, selecione α que apresenta a maior probabilidade

para dado c;

Após esse processo, faça o gráfico de p∗(c) em relação a c e procure

por picos e/ou regiões de estabilidade, ou seja, que mantêm o mesmo modelo numa

faixa de valores de c.

Uma variação do procedimento é, no passo 1, ao invés de usar apenas o

“melhor”modelo, considerar todos do mesmo tamanho e ponderar sua probabilidade

por 1/k, tal que k é o número de modelos de tamanho α dentro da cerca. Ao avaliar

os dois posśıveis gráficos constrúıdos, a preferência é dada ao modelo que pertencer ao

pico com menor c, ou seja, prefere-se o modelo cuja probabilidade de inclusão dentro

da cerca é alta e cuja distância ao modelo completo é pequena. Nessa variação, serão

encontradas probabilidades menores, entretanto pode-se destacar novos modelos que

sofreram pouco com a aplicação da ponderação. Logo, ao construir esses gráficos,

aplica-se o conceito de estabilidade ao buscar por intervalos na abscissa que contêm

um único ajuste selecionado. Sendo assim, é de interesse verificar a existência de

picos e intervalos de c que contêm um único ajuste destacado.

Na geração de todos os gráficos é necessário informar o número de

reamostras (B) que devem ser feitas. Murray et al. (2013) optaram por B = 1000.

Tarr et al. (2018) mostraram uma previsão do tempo gasto em cada função principal

ao fixar B = 50 e variar o número de regressoras e recomendaram o uso de B = 150


27

e B = 200. Neste trabalho, optou-se pela variação de B entre 150 e 1000.

Para aumentar a praticidade das ferramentas propostas, uma das

funções do mplot conduz a uma página no navegador de internet que possibilita

a interação do usuário com todos os gráficos gerados pelas funções. Portanto, ao uti-

lizar essas ferramentas, o pesquisador será capaz de identificar modelos e variáveis

que possam contribuir significativamente na explicação da resposta. Para exemplos

desse recurso, veja Tarr et al. (2018).


3 MATERIAL E MÉTODOS

Para a aplicação das ferramentas que assessoram o processo de seleção

de variáveis dispońıveis no mplot, utilizaram-se dois conjuntos de dados, o primeiro

ilustra a seleção de variáveis no modelo linear Normal e o segundo no modelo de

regressão loǵıstica. Em ambas as aplicações, a necessidade de inclusão de interações

entre regressoras obedece o prinćıpio da hierarquia entre efeitos, muito utilizado em

análise de dados e planejamento de experimentos e bem explicado em Wu & Hamada

(2009). Assim, primeiramente fez-se a busca pelos efeitos principais e posteriormente

para as interações entre os mesmos. Decidiu-se pelo uso desse prinćıpio uma vez que

o mplot não distingue entre efeito principal e interação, ou seja, a função trata, por

exemplo, o termo de interação entre duas regressoras como se fosse outra regressora,

o que não faz sentido na modelagem. Outro motivo para o uso deste prinćıpio é que,

para dados observacionais com muitas regressoras, a inclusão de por exemplo todos

efeitos principais e interações pode levar ao problema de colinearidade entre colunas

da matriz X e até ao problema n < p. Nas seções 3.1 e 3.2 são apresentadas algumas

informações sobre os dados para as duas aplicações. Todas as análises e cálculos

foram realizados no programa R (R Core Team, 2017).

3.1 Aplicação 1: estudo sobre peso de recém-nascidos (RN)

prematuros

O primeiro conjunto foi obtido de colaborações entre pesquisadores do

Departamento de Bioestat́ıstica (IBB, Unesp) e Faculdade de Medicina de Botucatu

(Unesp) e está relacionado ao estudo de Prigenzi et al. (2008), que investiga fato-


29

res de risco associados à mortalidade de recém-nascidos. Uma parcela dos dados

está exposta no Anexo 3. Neste trabalho, estudou-se a relação do peso de recém-

nascidos prematuros e outras caracteŕısticas da mãe, do peŕıodo gestacional e do

recém-nascido. A amostra contém 90 observações e treze variáveis com todos os

dados completos, sendo cinco quantitativas (idade da mãe, número de consultas,

número de gestações, número de partos, idade gestacional e peso do recém-nascido)

e sete qualitativas binárias (hábito de fumar da mãe - sim ou não; realização do

pré-natal - sim ou não; mãe com hipertensão arterial sistêmica - sim ou não; uti-

lização de corticoide pela mãe - sim ou não; gestação única - sim ou não; alguma

outra doença durante a gestação - sim ou não; sexo do recém-nascido - masculino ou

feminino). Nas Tabelas 3 e 4 são apresentadas algumas medidas descritivas de cada

caracteŕıstica.

Tabela 3. Medidas descritivas das variáveis quantitativas, dados da Aplicação 1 (n =

90)

Variável Código Mı́n. Mediana Média Máx. D.P.

Peso dos recém-nascidos (g) peso 635 1290 1318 2450 322

Idade da mãe (anos) idade mae 14 26 26,47 44 6,6

Número de consultas n consultas 0 4 4,4 14 2,5

Número de gestações n gest 1 2 2,36 7 1,4

Número de partos n partos 0 1 0,97 5 1,2

Idade gestacional (semanas) IG 24 30 30 36 2,4


30

Tabela 4. Medidas descritivas para as variáveis binárias, dados da Aplicação 1 (n =

90)

Variável Código Categoria Número %

Hábito de fumar fumo Sim 21 23,33

Realização do pré-natal pre natal Sim 81 90,00

Hipertensão HAS Sim 13 14,44

Uso de corticoide corticoide Sim 60 66,66

Gestação única GU Sim 75 83,33

Doença na gestação doenca na ges Sim 64 71,11

Sexo do recém-nascido GR Feminino 48 53,33

Figura 1 - Valores da correlação entre as regressoras quantitativas e a resposta peso,

dados da Aplicação 1 (n = 90)


31

Nota-se pela Figura 1 que, analisando as correlações lineares (Pear-

son) entre pares de variáveis, a regressora mais correlacionada com o peso de RN é

IG (idade gestacional), enquanto que as demais regressoras apresentam baix́ıssimo

coeficiente de correlação com a resposta. Há alta colinearidade entre as covariáveis

n gest e n partos, o que é esperado, pois para ocorrência de um parto é necessário

que ocorra a gestação (número de partos é sempre menor ou igual ao número de

gestações).

O objetivo nesse estudo é investigar as relações entre o peso e as demais

variáveis, explorando alguns métodos de seleção de variáveis e recursos oferecidos

no pacote mplot. Para isso, iniciou-se com um ajuste de um modelo de regressão

múltipla com erros Normais aplicando a função lm e seguido de algumas análises de

diagnóstico preliminares. Para verificar o pressuposto da normalidade, utilizou-se a

função qqPlot do pacote car (Fox & Weisberg, 2011), que resultou na Figura 3 e

4. Dessas análises, concluiu-se a necessidade de transformar a variável resposta, o

que foi feito seguindo a metodologia de Box-Cox (Box & Cox, 1964) e aplicando-se

a função box.Cox do pacote car. Para explorar os métodos de seleção de variáveis

aplicou-se os procedimentos usuais de stepwise utilizando os critérios de AIC e BIC

pela função step do R.

Utilizando os recursos do mplot, foram identificados alguns modelos

potenciais para o problema, sendo que em alguns deles abriu-se a possibilidade de

investigação de termos de interação entre variáveis regressoras. Seguindo o prinćıpio

de hierarquia entre os efeitos, no mplot, os termos de interações foram avaliados

aplicando-se o processo de seleção em um modelo que contém em sua resposta os

reśıduos de um ajuste feito com os efeitos principais selecionados e, na parte preditiva,

as posśıveis interações entre esses efeitos. Nota-se que esse procedimento é feito para

todos os modelos potenciais.

Em suma, os passos tomados até a obtenção dos modelos potenciais

foram:

1) ajuste um modelo linear Normal com todas as regressoras;


32

2) análise de diagnóstico indicou heterogeneidade de variâncias, o que foi resolvido

com a transformação Box-Cox;

3) aplicação dos procedimentos de seleção de variáveis para a variável resposta

peso transformada e efeitos principais das regressoras;

4) aplicação dos procedimentos de seleção de variáveis nos efeitos de interação das

variáveis selecionadas no item 3.

Nos passos 3 e 4 tantos os procedimentos automáticos quanto o mplot foram aplica-

dos. No Anexo 1 são apresentados os códigos utilizados para construção e geração

das ferramentas gráficas.

3.2 Aplicação 2: probabilidade de prenhez em inseminação

artificial (IA) de vacas

O segundo conjunto de dados, se refere a um estudo das relações entre

a probabilidade de prenhez de vacas inseminadas artificialmente (IA) e as diversas

caracteŕısticas do sêmen utilizado, obtido a partir das colaborações do Departamento

de Bioestat́ıstica (IBB, Unesp) e pesquisadores da área de veterinária. No estudo

um total de 3657 vacas foram IA com sêmen de touros também distintos, e classi-

ficadas em prenhez (sucesso) ou não (fracasso). Além das caracteŕısticas do sêmen

que englobam treze covariáveis, com observações completas, o estudo apresenta ou-

tras variáveis ambientais e de manejo que devem ser levadas em consideração como

as fazendas, os inseminadores, as datas da realização da IA e outros fatores. Essas

informações foram condensadas em uma variável categórica com quatro grupos, in-

clúıda na modelagem, como uma variável de blocagem. Em geral, mais de uma vaca

foi inseminada com sêmen da mesma amostra, havendo repetições de condições. Cada

condição é definida como uma combinação espećıfica dos valores das caracteŕısticas

do sêmen e do fator bloco, que será referida como lote. Assim, a modelagem pode

proceder a partir de dados agregados ou binários. Na forma agregada, o conjunto


33

de dados ficou composto por 72 lotes cujas respostas são em termos de números de

vacas no lote e número de sucessos. Os lotes são desbalanceados quanto ao número

de vacas. Na forma binária, a resposta de cada vaca é 0 (fracasso) ou 1 (sucesso).

O tipo de variável resposta leva ao modelo de regressão Binomial capaz

de modelar a probabilidade de prenhez em função das diversas covariáveis. Tanto a

versão para dados agregados, que supõe variáveis binomiais independentes aos lotes,

cada uma com parâmetro n espećıfico, quanto a versão para dados binários, que supõe

variáveis Bernoulli independentes, levam às mesmas estimativas de efeitos, uma vez

que as covariáveis são as mesmas. As probabilidades de sucesso são vistas como

funções dos efeitos das covariáveis, inclusive a de blocagem. Entretanto, para fins de

diagnóstico do ajuste, a forma agregada se apresenta mais vantajosa (Collett, 2002).

Neste trabalho, as duas formas foram utilizadas, cada qual escolhida para resolver os

problemas de implementação computacional enfrentados nas diversas fases, conforme

será esclarecido no decorrer da descrição dos métodos.

Para a construção do modelo a estratégia foi, inicialmente, a partir

de uma análise preliminar, explorar a presença de observações influentes e/ou ex-

tremas e, num primeiro momento, remover aquelas observações com indicações de

altamente espúrias ou influentes, antes da aplicação da metodologia de seleção de

variáveis. Após a seleção do modelo (ou dos modelos que se mostraram mais pro-

missores), procedeu-se às técnicas de diagnóstico do ajuste e posśıveis interpretações

dos resultados.

Assim, inicialmente foi realizada uma análise preliminar ajustando um

modelo com todos os efeitos principais das covariáveis, na versão de dados agre-

gados (72 lotes), com função de ligação logito, caracterizando então o modelo de

regressão loǵıstica. Para o diagnóstico do ajuste utilizaram-se os recursos do pacote

hnp (Moral et al., 2017). Essa inspeção mostrou 7 lotes extremamente influentes

no ajuste e decidiu-se trabalhar com apenas N = 65 lotes, totalizando nv = 2399

vacas inseminadas para a aplicação da teoria de seleção de variáveis, visto que o

ajuste com esses lotes não destaca pontos influentes e parece satisfazer às pressu-


34

posições do modelo loǵıstico, conforme brevemente descritos no Cápitulo 2, Seção

2.3. Poderia ser também de interesse verificar o quanto a presença desses pontos,

considerados at́ıpicos e/ou influentes, afetam ou dificultam os processos de seleção,

porém esse procedimento não foi realizado neste trabalho. Análises descritivas das

caracteŕısticas do sêmen dos 65 lotes estão apresentadas nas Tabelas 5 e 6.

No caso da regressão Binomial, devido à utilização do método bootstrap

para o cálculo das várias medidas pelo pacote mplot é necessário que a variável

resposta seja limitada em [0; 1], podendo ser a proporção de sucessos e declarando-se

ni (número total de vacas no lote i) como o peso da observação i (i = 1, 2, . . . , 65) ou

o uso dos dados na forma expandida. Esta é uma limitação do pacote que, felizmente,

pode ser contornada devido à equivalência entre os ajustes.

Para tal, declara-se na função glm do R a opção weight = w, na qual

w é o vetor coluna cujos elementos são os ni’s. É posśıvel também trabalhar com

os dados expandidos, ou seja, com uma observação para cada vaca inseminada com

resposta binária (0 ou 1) pois, como cada lote recebeu apenas um tipo de sêmem

e, é sabido, o total de vacas inseminadas e prenhas, é posśıvel replicar as linhas

necessárias para as covariáveis. Outra obrigatoriedade do pacote utilizado é que a

função de ligação utilizada, quando considerado um modelo Binomial, deve ser a

loǵıstica, visto que existem outras possibilidades para a mesma (Collett, 2002).

Nesta aplicação, existe um fator de delineamento do estudo, o fator

de blocagem que agrega os efeitos de fazenda, inseminadores, época e outros, con-

trolando a heterogeneidade externa que não deve sofrer o processo de seleção. O

mplot também não é capaz de diferenciar esse tipo de variável das demais sujeitas

ao processo de seleção. Assim, foi necessário buscar uma alternativa para aplicar o

procedimento mplot. Ao invés de utilizar a variável prenhez (binária) ou a proporção

de sucessos, utilizaram-se os reśıduos de um modelo loǵıstico ajustado apenas para o

fator “bloco”. Para isso, foi utilizado a estratégia de linearização do modelo loǵıstico

e aplicação do método de mı́nimos quadrados ponderados. Esse método produz es-

timativas equivalentes às do modelo loǵıstico usual (Hosmer et al., 1989). Os passos


35

são:

1. ajuste uma regressão loǵıstica, com todos os efeitos principais, para dados

binários (yi = 0; 1 com i = 1, 2, . . . , nv), em que nv é o número de vacas

inseminadas;

2. extraia as probabilidades estimadas, π̂i;

3. crie uma nova variável dependente, dada por

zi = ln

(
π̂i

1− π̂i

)
+

yi − π̂i
π̂i (1− π̂i)

,

com a segunda parcela definida como os reśıduos tipo working ;

4. considere os pesos como

vi = π̂i (1− π̂i);

5. ajuste um modelo linear Normal com ponderação no qual o vetor resposta é

z = (z1, z2, . . . , znv), as regressoras são as indicadoras de blocos e os pesos dado

pelo vetor v = (v1, v2, . . . , vnv). Coloque os reśıduos desse ajuste no vetor r;

6. ajuste um modelo linear Normal similar ao passo 5 porém com cada regressora,

digamos x, como resposta e coloque os reśıduos desse ajuste no vetor xr. Repita

para todas as regressoras que deseja explorar, nomeando o vetor de reśıduos

com o nome da regressora acrescido da letra r. A ideia nesse passo é remover

de cada regressora aquele componente que é explicado pelo fator bloco. Note

que, a ordem da regressora utilizada nesse passo não afetará o resultado final,

pois tem-se um ajuste separado para cada regressora;

7. ajuste um modelo linear Normal ponderado de r em função dos reśıduos obtidos

no passo 6;

8. aplique os comandos do mplot para gerar os gráficos.


36

No caso da aplicação, as regressoras que sofrerão os ajustes no passo 6 são aquelas

apresentadas na Tabela 5.

Após encontrar alguns modelos candidatos potenciais para explicar

a probabilidade de sucesso na inseminação, explorou-se a inclusão de termos de

interação de primeira ordem entre as regressoras pertinentes para cada um deles.

Novamente, a estratégia é aplicando-se o mesmo procedimento utilizado na Aplicação

1, baseado no prinćıpio da hierarquia entre efeitos.

Tabela 5. Medidas descritivas das variáveis quantitativas (caracteŕısticas do

sêmem), dados da Aplicação 2 (N = 65)

Variável Código Mı́n. Mediana Média Máx. D.P.

Membrana ı́ntegra (%) mint 4,66 50,96 47,61 74.98 16,23

Atividade mitocondrial (%) mitp 19,81 47,39 46,66 72,61 12,05

Motilidade do sêmem (%) mot 4,30 39,20 35,95 67,90 14,59

Defeitos totais (%) deft 0,06 0,15 0,16 0,40 0,08

Defeitos maiores (%) defm 0,02 0,08 0,10 0,28 0,07

Velocidade linear (µm/s) vsl 54,93 75,11 77,30 105,72 10,57

Velocidade curviĺınea (ou total) vcl 101,50 147,20 152,20 228,50 26

Velocidade média da traje-

tória do deslocamento (µm/s)
vap 64,13 89,04 91,28 128,03 12,80

Amplitude do deslocamento

lateral da cabeça (µm)
alh 4,33 6,52 6,84 9,74 1,31

Motilidade progressiva progr 2,60 29,10 26,74 48,20 11,07

Linearidade =(vsl/vcl)×100 lin 37,83 51,42 51,35 64,41 5,56

Retilinearidade =(vsl/vap)×100 str 73,23 84,05 83,51 90,56 4,36


37

Tabela 6. Medidas descritivas da variável resposta por lote, dados da Aplicação 2

(N = 65)

Variável (por lote) Mı́n. Mediana Média Máx. D.P.

Número de vacas (ni) 10 20 36,91 138 32,08

Número de vacas prenhas (yi) 2 11 17,26 65 15,32

Proporção de vacas prenhas (p̃i) 0,11 0,47 0,49 0,92 0,14

Em suma, o roteiro de análise dos dados da Aplicação 2 foi:

1) ajuste de um modelo de regressão binomial com função de ligação logito

(N=72);

2) análise diagnóstico indicou 7 lotes altamente influentes optando-se pela ex-

clusão destes;

3) construção do modelo linearizado incluindo o efeito da regressora de blocos;

4) aplicação dos procedimentos de seleção de variáveis no modelo de efeitos prin-

cipais (N=65);

5) aplicação dos procedimentos de seleção de variáveis nos efeitos de interação das

variáveis selecionadas no item 3.

Nos passos 4 e 5 tantos os procedimentos automáticos quanto o mplot foram apli-

cados. No Anexo 2 são apresentados os códigos utilizados para contrução e geração

das ferramentas gráficas.

Nota-se pela Figura 2 que, analisando as correlações lineares (Pearson)

entre pares de variáveis, as únicas regressoras que apontam algum destaque com a

proporção de prenhes são lin e str, com estimativas negativas, embora os coefi-

cientes estimados sejam baixos em valor absoluto. Essas duas covariáveis mostram

alt́ıssima correlação entre si. Vários outros coeficientes são altos, tanto positivos

quanto negativos, indicando problemas de colinearidade entre as regressoras, o que

torna a seleção de variáveis mais complicada.


38

Figura 2 - Valores da correlação entre as regressoras quantitativas e da resposta por

lote, dados da Aplicação 2 (N = 65)


4 RESULTADOS

4.1 Aplicação 1: estudo sobre o peso de recém-nascidos

(RN) prematuros

Nessa aplicação pretende-se explorar as técnicas de seleção de variáveis

para modelos explicativos do peso de recém-nascidos (RN) em função de variáveis

relacionadas à gestação. Considerou-se, inicialmente, um modelo Gaussiano com

apenas os efeitos principais das posśıveis regressoras. Embora a análise de variância

(Tabela 7) mostre a significância da regressão, ou seja, traga evidência de que ao

menos uma regressora possui relação linear com a resposta (F0,95;12;77 =1,9), nota-

se graficamente, na Figura 3, à direita, que o pressuposto de variância constante

não é satisfeito, pois verifica-se um aumento da dispersão dos reśıduos conforme

o peso estimado aumenta. O gráfico envelope (à esquerda) também confirma a

existência de reśıduos positivos ultrapassando a banda. Consequentemente, utilizou-

se a transformação Box-Cox, que indicou a transformação logaŕıtmica para o peso,

ou seja, y = ln(peso). Propõe-se então um novo modelo, cujas informações sobre o

ajuste são apresentadas na Tabela 8 e na Figura 4.

Tabela 7. ANOVA do modelo de regressão ajustado para o peso de RN

Fonte de variação S.Q. G.L. Q.M. F0

Regressão 4562050 12 380171 6,3

Reśıduo 4669544 77 60643

Total 9231594 89


40

−2 −1 0 1 2

−
2

−
1

0
1

2

Quantis Normal

R
es

íd
uo

s 
pa

dr
on

iz
ad

os

600 800 1000 1200 1400 1600 1800

−
2

−
1

0
1

2

Peso ajustado

R
es

íd
uo

s 
pa

dr
on

iz
ad

os

Figura 3 - Gráfico de reśıduos do modelo de efeitos principais ajustado para a variável

- peso de RN: envelope quantil-quantil Normal (à esquerda) e reśıduos em

função dos valores ajustados (à direita).

Tabela 8. ANOVA do modelo de regressão ajustado com y = ln(peso) de RN

Fonte de variação S.Q. G.L. Q.M. F0

Regressão 0,0509 12 0,0042 6,8

Reśıduo 0,0479 77 0,0006

Total 0,0988 89

−2 −1 0 1 2

−
2

−
1

0
1

2

Quantis Normal

R
es

íd
uo

s 
pa

dr
on

iz
ad

os

6.6 6.8 7.0 7.2 7.4

−
0.

4
−

0.
2

0.
0

0.
2

Y Ajustado

R
es

íd
uo

s 
pa

dr
on

iz
ad

os

Figura 4 - Gráfico de reśıduos do modelo de efeitos principais ajustado para a variável

- ln(peso) de RN: envelope quantil-quantil Normal (à esquerda) e reśıduos

em função dos valores ajustados (à direita).


41

O gráfico na Figura 5 relaciona os reśıduos com os valores indicadores

de alavanca, os hii’s, elementos da diagonal da matriz H. A área dos ćırculos é

proporcional aos valores das distâncias de Cook, permitindo a exploração de posśıveis

pontos influentes no ajuste. Dois pontos, 22 e 67, apresentam hii > 2p/n ≈0,29. Em

particular, a observação 22, é de um bebê cuja mãe não realizou o pré-natal e, dentre

essas mães, foi a que apresentou o maior número de gestações, 5, e apenas 3 partos. Já

a observação 67 apresenta o maior número de consultas entre todas, 14. Entretanto,

em termos da distância de Cook essas observações não se destacam. Verifica-se que

três outros pontos, 29, 12 e 87, apresentam destaque em termos de distância de

Cook, porém seus valores são 0,120; 0,123 e 0,135, respectivamente, considerados

ainda baixos e provavelmente não muito influentes no ajuste. Os outros pontos não

se destacaram o suficiente para que se faça um posśıvel tratamento.

0.05 0.10 0.15 0.20 0.25 0.30 0.35

−
2

−
1

0
1

2

Valores h

R
es

íd
uo

s 
S

tu
de

nt
iz

ad
os

12
29

87

Figura 5 - Gráfico para diagnóstico de influência do ajuste do modelo de efeitos

principais para ln(peso) de RN. O raio dos ćırculos são proporcionais à

distância de Cook.


42

Aplicando os métodos stepwise, backward e forward para seleção de

variáveis, obteve-se os mesmos modelos quando o critério foi mantido. Os modelos

selecionados são indicados na Tabela 9, mostrando que, ao considerar um valor de λ

maior, pois ln 90 > 2, o ajuste encontrado contém uma regressora a menos.

Tabela 9. Modelos selecionados pelos métodos usuais segundo o critério, para o

ln(peso) de RN

Critério Modelos

AIC (λ = 2) E(Y ) = β0 + β1 IG+ β2 GR1 + β3 pre natal1

BIC (λ = ln 90) E(Y ) = β0 + β1 IG+ β2 GR1

Considerando a metodologia do pacote mplot, verificam-se diferentes

possibilidades de modelos ao analisar estabilidade e observar quais covariáveis se des-

tacam. A Figura 6 apresenta o gráfico de inclusão de variáveis que mostra a regressora

idade gestacional (IG) sendo inclúıda com probabilidade 1, independentemente da

penalização (λ). A variável sexo do RN (GR) também se destaca com altas proba-

bilidades de inclusão para uma grande faixa de valores de λ. Com probabilidades

que decrescem com rapidez, porém bem maiores do que as da variável de referência

(RV), aparece a regressora binária pre natal. Duas outras variáveis que apresentam

probabilidades altas para uma penalidade baixa e que são ligeiramente maiores que

as de RV, são as binárias gestação única (GU) e doença hipertensiva (HAS), porém as

mesmas não seriam selecionadas em favor da simplicidade do modelo.

É posśıvel destacar no gráfico da Figura 6 qualquer regressora que o

pesquisador desejar, neste caso optou-se pela RV que auxilia na escolha de variáveis,

como mencionado. Outro destaque, é a indicação feita para as penalidades quando

considerados os critéros AIC (λ = 2) e BIC (λ = ln 90).


43

Figura 6 - Gráfico de inclusão de variáveis (VIP) em modelos de efeitos principais

para ln(peso) de RN, com destaque da curva de RV.

A Figura 7 apresenta gráficos de estabilidade dos modelos nos quais

para cada tamanho de modelo (número de regressoras mais intercepto) uma medida

de perda ou falta de ajuste é apresentada. Enquanto o ponto azul, para 1 parâmetro,

apresenta o modelo apenas com intercepto, o ponto vermelho, com 14 parâmetros,

mostra um ajuste composto pelos efeitos principais das 12 regressoras mais o in-

tercepto e a variável redundante (RV). Destacam-se, nos gráficos, os modelos que

contêm, ou não, determinada variável de interesse, neste caso analisou-se com três

regressoras que possuem maiores probabilidades na Figura 6. Todavia, é posśıvel

destacar qualquer covariável de interesse. Nota-se no primeiro gráfico que a inclusão

da regressora IG tem grande impacto nos valores do componente de perda, mesmo no

modelo em que ela se apresenta sozinha. Verifica-se também que, dado IG, ao incluir

a covariável GR, ocorre a diminuição substancial dos valores no segundo gráfico. Já

para a regressora pre natal não se observa clara separação entre os modelos, logo

deve-se buscar pelo aux́ılio de outras ferramentas para investigar sua importância.


44

Figura 7 - Gráficos do valor do componente de perda contra a dimensão do modelo.

Modelos de efeitos principais para ln(peso) de RN. Destaque para as re-

gressoras IG (à esquerda), GR (à direita) e pré natal (inferior).


45

A Figura 8 complementa a Figura 7 ao relacionar o valor do compo-

nente de perda ao número de parâmetros do modelo, com a informação adicional

sobre a probabilidade de seleção dos modelos com determinada dimensão represen-

tado pelo tamanho do ćırculo. Na figura também pode-se destacar uma variável

que, no caso mostrado, foi a variável pré natal. Assim, o modelo nulo e o com-

pleto apresentam probabilidade 1 pois, por conterem o número mı́nimo e o máximo

de parâmetros, respectivamente, eles sempre serão selecionados na reamostragem de

seu tamanho. Note que modelos sem o intercepto não são permitidos. Modelos de

destaque são aqueles com probabilidades altas e valores baixos da função de perda e

número de parâmetros. As probabilidades dos três modelos em destaque, no canto

inferior esquerdo do gráfico, estão apresentadas na Tabela 10. Nota-se que mesmo

com pouco destaque, o modelo que contém pré-natal foi selecionado em 28% das

vezes, dado a dimensão quatro, tornando sua inclusão uma nova opção para a mo-

delagem. No tipo de gráfico desta figura, pode-se destacar qualquer regressora de

interesse uma a uma, entretando, decidiu-se mostrar apenas a pré natal, pois ela

não apresentou de forma evidente sua importância nas outras figuras analisadas.

Por se tratar de uma ferramenta interativa, quando um ćırculo é sele-

cionado, como ilustrado no gráfico da Figura 8, uma caixa com algumas informações

sobre o modelo apontado é disponibilizada. Dentre as informações tem-se o modelo,

o valor de −2∗Log-likelihood (LL), a variável de destaque (var.ident) e a probabili-

dade de seleção feita pelo método de reamostragem (prob). A letra k indica o número

de parâmetros do modelo, mas não apresenta com exatidão (devido ao efeito jitter

usado no gráfico para possibilitar a identificação de modelos do mesmo tamanho),

logo deve ser aproximado para um número inteiro mais próximo.


46

Figura 8 - Gráfico de probabilidades de seleção de modelos em função do número de

parâmetros com destaque para a regressora pré-natal. Modelos de efeitos

principais para ln(peso) de RN, utilizando bootstrap ponderado.

Tabela 10. Modelos com probabilidade de seleção superior a 0, 20 visualizados na

Figura 8. Modelos de efeitos principais para ln(peso) de RN

Variável presente Probabilidade

IG 1,0

IG e GR 0,77

IG, GR e pre natal 0,28


47

Figura 9 - Gráficos de probabilidades de seleção de modelos usando o método fence.

Modelos de efeitos principais para o ln(peso) de RN.

A Figura 9 mostra o gráfico para avaliar a estabilidade de modelos

usando o método fence. Nesse gráfico, deve-se procurar regiões em c que destaquem

apenas um modelo e pontos com altas probabilidades (picos). O gráfico superior

é a versão quando apenas o melhor modelo é considerado no cálculo de p∗(α, c) e

o inferior quando todos são considerados e as probabilidades são ponderadas. Am-

bos apresentam padrão similar, pois as maiores probabilidades indicam o ajuste que


48

contem apenas a regressora IG. Outro modelo de destaque é o que contem as regres-

soras IG e GR, pois se destacam num intervalo de c com probabilidades consideráveis.

Nos mesmos gráficos, os outros modelos que são visualizados não se destacaram em

nenhum intervalo e suas probabilidades são pequenas, logo são descartados.

A Tabela 11 apresenta as regressoras dos modelos selecionados para

investigação mais detalhada. Para os modelos 2 e 3 existe a possibilidade de inves-

tigação de inclusão de termos de interação. Para esta investigação, ajusta-se um

modelo cuja variável resposta é formada pelos reśıduos do ajuste do modelo com os

efeitos principais, das regressoras de interesse, e a parte preditiva contém somente as

interações, conforme indicado pelos próprios autores do pacote, Tarr et al. (2018).

Tabela 11. Modelos selecionados conforme aplicação do mplot para ln(peso) de RN

Modelo Variável presente

1 IG

2 IG e GR

3 IG, GR e pré-natal

Os gráficos de inclusão de variáveis para as interações nos modelos 2

e 3 estão apresentados na Figura 10, indicando que a inclusão de interações não é

importante em nenhum deles. No gráfico, nota-se que para os dois modelos, as in-

terações apresentam curvas de probabilidades baixas que caem rapidamente com a

penalidade estando próximas ou abaixo da curva de RV. A exploração das demais

ferramentas, para estudar a relevância das interações, não revelou resultados diferen-

tes. Portanto, busca-se pelas adequações dos modelos que contêm apenas os efeitos

principais das regressoras indicadas na Tabela 11. Estimativas pontuais e intervalos

com 95% de confiança são apresentados na Tabela 12.

Para fins de exploração, aplicou-se também a seleção de termos de

interação de primeira ordem nos modelos selecionados pelos métodos automáticos

(stepwise) usando os critérios AIC e BIC (Tabela 9). Em nenhum dos casos, termos

de interação foram selecionados, indicando como modelos finais os mesmos modelos


49

2 e 3 da Tabela 11.

Figura 10 - Gráficos VIP dos termos de interação para os modelos 2 (à esquerda) e

3 (à direita) para o ln(peso) de RN.

Tabela 12. Estimativas dos parâmetros e intervalos de confiança (95%) para os ajus-

tes dos modelos selecionados (Tabela 11) para ln(peso) de RN

Estimativa (I.C.)

Regressora Modelo 1 Modelo 2 Modelo 3

Intercepto 5,228 (4,739; 5,717) 5,228 (4,759; 5,695) 5,231 (4,765; 5,696)

IG 0,064 (0,048; 0,080) 0,066 (0,050; 0,082) 0,064 (0,048; 0,079)

GR (feminino) - -0,112 (-0,186; -0,038) -0,114 (-0,188; -0,041)

pre natal (sim) - - 0,089 (-0,037; 0,215)

Análises de diagnóstico dos ajustes dos três modelos não evidenciaram

problemas com os pressupostos e não foram encontrados pontos de influência con-

sideráveis. Logo, obteve-se as estimativas dos parâmetros dos modelos finais e os

respectivos intervalos de confiança, feitos pela função confint no R, apresentados na

Tabela 12. Dessa forma, mantendo-se fixas as demais variáveis do modelo, conforme

a idade gestacional cresce há um aumento esperado do peso. O efeito esperado de

sexo é negativo, indicando que em média bebês do sexo feminino apresentam menor

peso ao nascer. Avaliando o efeito de pré-natal, tem-se que, ao manter as outras


50

regressoras fixas, espera-se um aumento do peso. Deve-se notar, entretanto, que

o IC para o efeito dessa variável, engloba o valor 0, indicando evidência fraca de

efeito dessa variável. Deve-se destacar também que, das 90 observações, apenas 9

não fizeram o pré-natal, uma posśıvel explicação para a falta de poder na detecção

do efeito desse fator. Do ponto de vista prático, essa evidência fraca é importante já

que alerta aos benef́ıcios do pré-natal, principalmente em gestações de risco.

Avaliando o modelo 1, nota-se que para cada semana a mais de

gestação, o peso estimado do RN ficará 6, 6% maior do que o peso da semana ante-

rior. Para os outros modelos essa porcentagem é similar, já que suas estimativas são

muito próximas.

Estimando a diferença entre sexos de RN em que y = ln(peso) no modelo 2, tem-se

que

ŷ(GR = 0)− ŷ(GR = 1) =

= (5, 228 + 0, 066× IG− 0, 112× 0)− (5, 228 + 0, 066× IG− 0, 112× 1) = 0, 112.

Como ŷ = ̂ln(peso), fazendo-se a tranformação tem-se p̂eso = 1, 119.

Logo, a estimativa da diferença esperada entre os pesos dos RN do sexo masculino e

feminino é constante, não importando os valores de IG, indicando que o RN do sexo

masculino tem peso esperado de 11, 9% maior do que o sexo feminino.

Para o modelo 3 a estimativa da diferença da resposta entre o sexo

masculino e feminino, considerando fixas as covariáveis IG e pré natal, tem-se que

ŷ(GR = 0)− ŷ(GR = 1) =

= (5, 231 + 0, 064× IG− 0, 114× 0 + 0, 089× pre natal)−

(5, 231 + 0, 064× IG− 0, 114× 1 + 0, 089× pre natal) = 0, 114.

Aplicando-se a transformação inversa tem-se p̂eso = 1, 121. Logo,

estima-se uma diferença constante em que o valor esperado do peso do RN do sexo

masculino é 12, 1% maior do que o sexo feminino. Considerando a realização ou não

do pré-natal dado as outras covariáveis fixas, a diferença entre ŷ é dado por

ŷ(pre natal = 1)− ŷ(pre natal = 0) =


51

= (5, 231 + 0, 064× IG− 0, 114×GR + 0, 089× 1)−

(5, 231 + 0, 064× IG− 0, 114×GR + 0, 089× 0) = 0, 089.

Novamente, aplicando a inversa, tem-se que p̂eso = 1, 93. Logo, a

estimativa da diferença dos pesos é constante, sendo que o RN cuja mãe realizou o

pré-natal, apresenta uma estimativa de peso esperado de 9, 3% maior do que os que

não realizaram.

A Figura 11, apresenta os valores de ln(peso) contra os valores de IG

e as retas ajustadas segundo o modelo 2.

24 26 28 30 32 34 36

6
.4

6
.6

6
.8

7
.0

7
.2

7
.4

7
.6

7
.8

IG (semanas)

lo
g

(p
e

s
o

)

Masculino
Feminino

Figura 11 - Retas ajustadas para o ln(peso) de RN segundo o modelo 2.

Na Figura 12, tem-se os valores do ln(peso) contra a idade gestacio-

nal, com as retas ajustadas segundo o modelo 3. Nota-se que um recém-nascido

do sexo masculino ao não realizar o pré-natal, tem seu ln(peso) esperado, e conse-

quentemente, seu peso esperado, próximo de um recém-nascido do sexo feminino que

realizou o pré-natal.


52

24 26 28 30 32 34 36

6
.4

6
.6

6
.8

7
.0

7
.2

7
.4

7
.6

7
.8

IG (semanas)

lo
g

(p
e

s
o

)

Masc. c/ pré−natal
Fem. c/ pré−natal
Masc. s/ pré−natal
Fem. s/ pré−natal

Figura 12 - Retas ajustadas para o ln(peso) de RN segundo o modelo 3.


53

4.2 Aplicação 2: probabilidade de prenhez em inseminação

artificial (IA) de vacas

Para estudar a relação entre a probabilidade de prenhez de vacas e

as caracteŕısticas do sêmen do touro, considerou-se, inicialmente, um modelo de re-

gressão Binomial com função de ligação loǵıstica, e apenas os efeitos principais no

preditor linear. Aqui destaca-se que o estudo apresenta um número razoável de co-

variáveis (12), além do fator de blocos, não sendo viável a exploração de um modelo

muito complexo, por exemplo, incluindo termos de interação entre as covariáveis.

Assim, antes de proceder com os métodos de seleção, realizou-se uma análise de di-

agnóstico preliminar na tentativa de explorar a presença de pontos problemáticos

tanto no ajuste quanto no processo de seleção já que estes podem influenciar na

seleção. Segundo Faraway (2016), na existência de observações influentes, estas de-

vem ser, pelo menos temporariamente, removidas do conjunto de dados antes de

aplicar qualquer estratégia de seleção de variáveis. Outro ponto a se destacar é que,

para regressão Binomial, ap