Métodos Estatísticos Avançados em Bioinformática
Profa. Dra. Elisângela Lizzi & Profa. Dra. Glaucia Bressan
INTRODUÇÃO
A disciplina Métodos Estatísticos Avançados em Bioinformática, concebida e ministrada pelas Profa. Dra. Elisângela Lizzi e Profa. Dra. Glaucia Bressan, nasce da convergência entre a estatística e a matemática computacional, áreas fundamentais para o enfrentamento dos desafios impostos pelos dados biológicos contemporâneos. Sua proposta emerge da necessidade de ir além de abordagens meramente operacionais, frequentemente centradas na aplicação automática de ferramentas, para promover um entendimento aprofundado dos fundamentos teóricos, das suposições dos modelos e das implicações inferenciais dos métodos utilizados em bioinformática.
Nesse contexto, trata-se de uma disciplina em nível de doutorado, estruturada para formar pesquisadores capazes de compreender, criticar e desenvolver metodologias, e não apenas utilizá-las como “caixas-pretas”. O material didático é inédito e foi integralmente concebido para esta disciplina, refletindo uma curadoria rigorosa de conteúdos e aplicações atuais. Espera-se que a experiência seja intelectualmente exigente e, ao mesmo tempo, altamente proveitosa, contribuindo para a formação analítica e científica dos participantes.
Com carinho científico, Elisangela e Glaucia!
1.1. Apresentação geral da área e objetivos
A bioinformática contemporânea demanda abordagens estatísticas não triviais e capazes de lidar com dados de alta dimensionalidade, heterogeneidade e complexidade estrutural. Esta disciplina tem como objetivo fornecer fundamentos teóricos e práticos de métodos estatísticos avançados aplicados à análise de dados biológicos, com ênfase em inferência, modelagem e validação de resultados. Busca-se capacitar os estudantes a compreender, implementar e interpretar técnicas modernas amplamente utilizadas em genômica, transcriptômica, metagenômica e áreas correlatas, todos ancoradas na linguagem computacional R.
1.2. Estratégias de operacionalização dos conteúdos
Os conteúdos serão desenvolvidos por meio de uma abordagem teórico-prática, integrando exposições conceituais, estudos dirigidos, implementação computacional e análise de estudos de caso em quatro módulos distintos, sendo: Técnicas de Permutação e Reamostragem;AMOVA – Análise Molecular de Variância;Modelos Lineares Generalizados (GLM) e Redes Bayesianas – Fundamentos e Algoritmos. As atividades incluirão resolução de problemas, simulações, aplicação em bases de dados reais e crítica de literatura científica, promovendo a consolidação do aprendizado de forma aplicada. Pois ao final teremos, as aplicações práticas e discussão de literatura científica.
1.3. Integração entre estatística, computação e aplicações biológicas
A disciplina enfatiza a natureza interdisciplinar da bioinformática, articulando conceitos estatísticos e matemáticos com algoritmos computacionais e problemas biológicos reais. Serão exploradas estratégias de modelagem que conectam teoria estatística à interpretação biológica, permitindo ao estudante compreender tanto os aspectos metodológicos quanto as implicações práticas dos resultados obtidos.
1.4. Visão geral das ferramentas (uso do R)
Será adotada a linguagem R como principal ferramenta computacional, considerando sua ampla utilização na análise de dados biológicos e estatísticos. Serão abordados pacotes específicos para cada módulo, promovendo autonomia na implementação e reprodutibilidade das análises.
MÓDULO 2: Técnicas de Permutação e Reamostragem
A reamostragem baseia-se em calcular estimativas a partir de repetidas amostragens dentro de uma única amostra observada. Esta abordagem é fundamental quando a distribuição amostral teórica de um estimador é desconhecida ou difícil de definir (ex: mediana, coeficientes não lineares) Nesse ponto temos limitações da Estatística Clássica Dependência de suposições como normalidade e variâncias iguais (ANOVA, Regressão). Dificuldade em amostras pequenas, onde o Teorema do Limite Central (TLC) pode não ser válido. Incerteza na definição da distribuição de estimadores complexos.
2.1. Fundamentos de testes não-paramétricos
Os testes não paramétricos constituem uma classe de métodos estatísticos que não assumem uma forma funcional específica para a distribuição dos dados (por exemplo, distribuição normal conhecido como normalidade dos dados). Em vez disso, baseiam-se em propriedades invariantes dos dados, como ordenação(ranqueamento) ou reamostragem, sendo particularmente úteis em contextos com amostras pequenas, distribuições assimétricas ou presença de outliers, muito comum em dados de saúde e biológicos. Diferente dos testes clássicos, os testes de permutação não dependem de distribuições teóricas de probabilidade (distribuição de poisson, distribuição normal, distribuição t-student e muitas outras). Eles constroem a distribuição de referência “embaralhando” os dados observados sob a hipótese nula (\(H_0\)).
# Dados: Área classificada em duas situações (ex: 1 vs 2 imagens) [cite: 83]
situacao_1 <- c(70, 51, 60, 90, 43, 15, 25, 103)
mean(situacao_1)## [1] 57.125## [1] 63.125## [1] 6# Simulação Monte Carlo [cite: 82]
combinados <- c(situacao_1, situacao_2)
n_perm <- 5000 #numero de replicações
diffs_perm <- replicate(n_perm, {
shuffled <- sample(combinados) #embaralhando
mean(shuffled[9:16]) - mean(shuffled[1:8])
})
p_valor <- sum(abs(diffs_perm) >= abs(obs_diff)) / n_perm
cat("Diferença Observada:", obs_diff, "\np-valor:", p_valor)## Diferença Observada: 6
## p-valor: 0.72542.2. Testes de aleatorização e testes de permutação
Diferentemente dos testes paramétricos clássicos, os testes de permutação não dependem de distribuições de probabilidade teóricas conhecidas . Sob a hipótese nula (\(H_0\)), assume-se que os rótulos dos grupos são intercambiáveis. A distribuição nula da estatística de teste é então construída empiricamente por meio da geração de todas (ou um grande número de) permutações possíveis dos dados observados. Esse procedimento permite calcular valores-p exatos ou aproximados, mantendo forte aderência aos dados e reduzindo dependência de suposições teóricas e probabilísticas.
Os testes de aleatorização são conceitualmente relacionados, sendo frequentemente utilizados em experimentos com alocação aleatória, onde a inferência é baseada no mecanismo de randomização do delineamento experimental.
2.3. Métodos de reamostragem e estimação de incerteza
Os métodos de reamostragem, como o bootstrap, consistem na geração repetida de amostras a partir dos dados observados, tipicamente com reposição. Essas amostras são utilizadas para aproximar a distribuição amostral de estimadores, permitindo a estimação de medidas de incerteza, como erro padrão, viés e intervalos de confiança. Tais abordagens são amplamente aplicáveis, inclusive em situações em que soluções analíticas são inviáveis ou inexistentes.
2.4. Jackknife: princípios e aplicações
O método jackknife baseia-se na reamostragem sistemática por exclusão de observações individuais (leave-one-out). Para um conjunto de dados com \(n\) observações, são geradas \(n\) subamostras, cada uma omitindo um elemento. A variabilidade entre as estimativas obtidas nessas subamostras é utilizada para estimar o viés e o erro padrão do estimador original. Matemáticamente definido da seguinte forma:
i-ésimo pseudovalor: \(x_{(i)}^{*} = n\hat{\theta} - (n-1)\hat{\theta}_{(i)}\)
Estimador Jackknife: \(\hat{\theta}_{jk} = \frac{1}{n}\sum_{i=1}^{n}x_{(i)}^{*}\)
Variância de Jackknife: \(Var_{jk}(\hat{\theta}) \cong \frac{n-1}{n}\sum_{i=1}^{n}(\hat{\theta}_{(i)} - \hat{\theta}_{(\cdot)})^{2}\)
Fonte: Rennó, CM (2011)[Jackknife, Bootstrap e outros métodos de reamostragem- INPE].
O jackknife é particularmente útil para estimadores suaves e possui menor custo computacional em comparação ao bootstrap, embora com limitações em cenários mais complexos.
dados <- c(2.2, 3.5, 3.4, 6.7, 6.2, 8.2, 9.2, 7.9, 9.0, 10.1)
n <- length(dados)
# Função para média geométrica
mean_geom <- function(x) exp(mean(log(x)))
# Replicação Jackknife
jack_reps <- sapply(1:n, function(i) mean_geom(dados[-i]))
variancia_jk <- ((n-1)/n) * sum((jack_reps - mean(jack_reps))^2)
cat("Média Geométrica:", mean_geom(dados), "\nVariância Jackknife:", variancia_jk)## Média Geométrica: 5.984431
## Variância Jackknife: 1.0118872.5. Bootstrap: variantes, propriedades e usos em bioinformática
O bootstrap é um método de reamostragem baseado na geração de múltiplas amostras com reposição a partir do conjunto de dados original. Seu objetivo principal é aproximar a distribuição amostral de um estimador, especialmente em situações onde a forma analítica dessa distribuição é desconhecida ou intratável.
Entre as principais variantes, destacam-se:
Bootstrap não paramétrico: reamostragem direta dos dados observados, sem de distribuição. Bootstrap paramétrico: amostras são geradas a partir de uma distribuição ajustada aos dados. Bootstrap estratificado: preserva a estrutura de subgrupos ou classes durante a reamostragem. Bootstrap BCa (bias-corrected and accelerated): corrige viés e assimetria na construção de intervalos de confiança.
Do ponto de vista teórico, o bootstrap apresenta consistência assintótica sob condições gerais e permite estimar erro padrão, viés e intervalos de confiança com elevada flexibilidade. No entanto, sua performance pode ser limitada em amostras muito pequenas ou em estatísticas altamente sensíveis a outliers.
Figura 2: Esquema de reamostragem Bootstrap.Fonte: Rennó, CM (2011)[Jackknife, Bootstrap e outros métodos de reamostragem- INPE].
Diferente do Jackknife, o Bootstrap realiza o sorteio aleatório com reposição, permitindo um número muito maior de replicações.Temos as seguintes formulações matemáticas:
Estimativa Bootstrap: \(\hat{\theta}_{b} = \frac{1}{B}\sum_{j=1}^{B}\theta_{(j)}\)
Variância Bootstrap: \(Var_{b}(\hat{\theta}) = \frac{1}{B-1}\sum_{j=1}^{B}(\hat{\theta}_{(j)} - \bar{\theta}_{b})^{2}\)
Na bioinformática, o bootstrap é amplamente utilizado, por exemplo, na avaliação de robustez de árvores filogenéticas (via suporte bootstrap), na inferência de redes gênicas e na estimativa de estabilidade de assinaturas moleculares em dados ômicos (genômica, transcriptômica, proteômica), frequentemente caracterizados por alta dimensionalidade e baixa razão amostra/variável.
func_boot <- function(data, indices) mean_geom(data[indices])
boot_res <- boot(data = dados, statistic = func_boot, R = 2000)
print(boot_res)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = dados, statistic = func_boot, R = 2000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 5.984431 0.07788026 0.9122543
2.6. Validação cruzada: estratégias, particionamento e avaliação de modelos
A validação cruzada (cross-validation) é um conjunto de técnicas utilizadas para avaliar o desempenho preditivo de modelos estatísticos e de aprendizado de máquina, especialmente no contexto de generalização para dados não observados.
A ideia central consiste em particionar os dados em subconjuntos de treino e teste de forma sistemática. As principais estratégias incluem:
Holdout: divisão simples entre treino e teste (ex.: 70/30). k-fold cross-validation: os dados são divididos em k subconjuntos; o modelo é treinado em k-1 folds e testado no fold restante, repetindo o processo k vezes. Leave-one-out (LOOCV): caso extremo de k-fold, onde k = n. Validação cruzada estratificada: preserva a proporção de classes em cada partição, essencial em dados desbalanceados. Validação cruzada aninhada (nested CV): utilizada para seleção de hiperparâmetros, evitando viés otimista na avaliação do modelo.
A escolha da estratégia impacta diretamente o viés e a variância da estimativa de desempenho. Em geral, há um trade-off entre custo computacional e precisão da estimativa.
As métricas de avaliação variam conforme o problema (classificação, regressão, sobrevivência), incluindo acurácia, AUC, F1-score, erro quadrático médio, entre outras.
# Exemplo de regressão do material (X e Y) [cite: 92, 94]
df <- data.frame(
x = c(1.2, 1.9, 2.8, 4.3, 5.5, 6.0, 7.2, 7.8, 9.1, 10.3, 11.1, 13.7, 14.7, 2.5, 3.6, 5.6, 7.8, 10.1, 12.0, 12.4),
y = c(16.4, 13.3, 18.4, 21.4, 27.7, 19.6, 23.0, 27.2, 25.8, 32.1, 34.0, 44.2, 41.0, 12.8, 22.1, 23.3, 24.7, 31.9, 38.0, 39.2)
)
ctrl <- trainControl(method = "cv", number = 5)
modelo <- train(y ~ x, data = df, method = "lm", trControl = ctrl)
print(modelo)## Linear Regression
##
## 20 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 17, 17, 15, 15, 16
## Resampling results:
##
## RMSE Rsquared MAE
## 3.000461 0.9469342 2.691638
##
## Tuning parameter 'intercept' was held constant at a value of TRUEVamos analisar caso- a - caso:
A.Holdout (Divisão Simples)É a forma mais básica de validação. A amostra é dividida aleatoriamente em dois subconjuntos, sendo treinamento (\(D_{train}\)) e teste/validação (\(D_{test}\)). Nogeral usa-se a proporção comum: \(70/30\) ou \(80/20\), baseando- se na distribuição dos dados (distribuição de probabilidade de pareto). Temos que a fórmula do erro: \(Err_{holdout} = \frac{1}{n_{test}} \sum_{i \in D_{test}} L(y_i, \hat{f}(x_i))\)
Onde \(L\) é a função de perda (ex: Erro Quadrático Médio) e \(\hat{f}\) é o modelo treinado em \(D_{train}\).
B. k-fold Cross-ValidationOs dados são divididos em \(k\) partições (folds) de tamanhos aproximadamente iguais.Vamos entender o processo, o modelo é treinado em \(k-1\) folds e testado no fold restante. E, isso se repete \(k\) vezes. A formulação matemática é : \(CV_{(k)} = \frac{1}{k} \sum_{i=1}^{k} MSE_i\)
C. Leave-one-out (LOOCV)É o caso extremo do k-fold, onde \(k = n\) (número total de observações).O processo é o seguinte: em cada iteração, apenas uma observação é deixada de fora para teste.A fórmulação matemática é: \(CV_{(n)} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_{(i)})^2\)Onde \(\hat{y}_{(i)}\) é a previsão para a observação \(i\) usando o modelo treinado sem ela.
D. Validação Cruzada Estratificada Essencial para dados desbalanceados (ex: estudos de doenças raras ou mortalidade específica).O principio norteador é garantir que cada fold tenha a mesma proporção de classes (ou distribuição da variável resposta) que o dataset original. Se a prevalência de uma doença é \(10\%\), cada fold terá exatamente \(10\%\) de casos positivos.
E. Validação Cruzada Aninhada (Nested CV)Utilizada para seleção de hiperparâmetros (ajuste de modelo) e estimativa de erro de generalização simultaneamente.Trabalha-se com o loop externo, onde se estima o erro do modelo e o loop interno, onde se realiza a busca pelos melhores hiperparâmetros.O objetivo é evitar o “viés otimista”, onde o modelo parece melhor do que é porque os mesmos dados foram usados para ajustar e testar.
# Carregando pacotes necessários
library(caret)
library(dplyr)
# Criando o dataset do material (Slides 10)
df <- data.frame(
x = c(1.2, 1.9, 2.8, 4.3, 5.5, 6.0, 7.2, 7.8, 9.1, 10.3,
11.1, 13.7, 14.7, 2.5, 3.6, 5.6, 7.8, 10.1, 12.0, 12.4),
y = c(16.4, 13.3, 18.4, 21.4, 27.7, 19.6, 23.0, 27.2, 25.8, 32.1,
34.0, 44.2, 41.0, 12.8, 22.1, 23.3, 24.7, 31.9, 38.0, 39.2)
)
# 1. Holdout (70/30)
set.seed(123)
index <- createDataPartition(df$y, p = 0.7, list = FALSE)
train_set <- df[index, ]
test_set <- df[-index, ]
model_holdout <- lm(y ~ x, data = train_set)
pred_holdout <- predict(model_holdout, test_set)
rmse_holdout <- RMSE(pred_holdout, test_set$y)
# 2. k-fold (k=5)
ctrl_kfold <- trainControl(method = "cv", number = 5)
model_kfold <- train(y ~ x, data = df, method = "lm", trControl = ctrl_kfold)
# 3. Leave-One-Out (LOOCV)
ctrl_loocv <- trainControl(method = "LOOCV")
model_loocv <- train(y ~ x, data = df, method = "lm", trControl = ctrl_loocv)
# 4. Estratificada (Exemplo com classificação simulada)
# Para regressão, o caret faz a estratificação baseada em quantis da variável y
ctrl_strat <- trainControl(method = "cv", number = 5)
model_strat <- train(y ~ x, data = df, method = "lm", trControl = ctrl_strat)
# Comparação de Resultados
resultados <- data.frame(
Metodo = c("Holdout", "k-fold (5)", "LOOCV"),
RMSE = c(rmse_holdout, model_kfold$results$RMSE, model_loocv$results$RMSE)
)
print(resultados)## Metodo RMSE
## 1 Holdout 2.966490
## 2 k-fold (5) 2.936746
## 3 LOOCV 3.000624Dica: O LOOCV costuma ser computacionalmente caro em datasets grandes, mas para o tamanho de amostra do exemplo (n=20), ele é muito estável e ideal. Para os dados de alta dimensionalidade (2.7 da ementa), a Validação Cruzada Aninhada será sua melhor aliada para garantir que a seleção de variáveis não vicie os resultados.
2.7. Aplicações computacionais em dados de alta dimensionalidade
Dados de alta dimensionalidade, caracterizados por um número de variáveis (\(p\)) muito superior ao número de observações (\(n\)), são comuns em áreas como bioinformática, epidemiologia molecular e ciência de dados. Esse cenário introduz desafios específicos, como sobreajuste (overfitting), instabilidade de modelos e dificuldades na interpretação.
Técnicas de reamostragem e validação são fundamentais nesse contexto, pois permitem:
Avaliar estabilidade de modelos e seleções de variáveis (ex.: seleção de genes relevantes); Reduzir overfitting, por meio de validação cruzada rigorosa; Estimar incerteza em cenários sem soluções analíticas robustas; Comparar modelos preditivos de forma mais confiável.
Além disso, abordagens computacionais modernas frequentemente combinam reamostragem com métodos de regularização (como LASSO e Ridge), aprendizado de máquina (Random forests, SVM) e técnicas de redução de dimensionalidade (PCA, t-SNE, UMAP).
Em bioinformática, essas estratégias são essenciais para análise de dados de sequenciamento de nova geração (NGS), estudos de expressão gênica, análise de variantes genéticas e integração de dados multi-ômicos, onde a robustez estatística e a reprodutibilidade são requisitos críticos.
Vamos aos exemplos:
- Simulando um cenário p >> n (dados tipo expressão gênica)
set.seed(123)
n <- 50 # número de amostras
p <- 1000 # número de variáveis (genes)
# matriz de preditores
X <- matrix(rnorm(n * p), nrow = n, ncol = p)
# apenas alguns genes realmente relevantes
beta <- c(rep(2, 10), rep(0, p - 10))
# variável resposta (regressão)
y <- X %*% beta + rnorm(n)
# versão para classificação (caso binário)
y_class <- as.factor(ifelse(y > median(y), "Doente", "Controle"))2)Redução de dimensionalidade (PCA), frequentementes usado em dados ômicos para explorar estrutura global.
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 5.50391 5.35638 5.29642 5.23227 5.19847 5.10144 5.0695
## Proportion of Variance 0.03029 0.02869 0.02805 0.02738 0.02702 0.02602 0.0257
## Cumulative Proportion 0.03029 0.05898 0.08704 0.11441 0.14144 0.16746 0.1932
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 5.05816 5.02599 5.00730 4.92785 4.91010 4.88355 4.82976
## Proportion of Variance 0.02558 0.02526 0.02507 0.02428 0.02411 0.02385 0.02333
## Cumulative Proportion 0.21875 0.24401 0.26908 0.29336 0.31747 0.34132 0.36465
## PC15 PC16 PC17 PC18 PC19 PC20 PC21
## Standard deviation 4.80558 4.77036 4.74849 4.69543 4.65956 4.61020 4.5939
## Proportion of Variance 0.02309 0.02276 0.02255 0.02205 0.02171 0.02125 0.0211
## Cumulative Proportion 0.38774 0.41050 0.43305 0.45509 0.47680 0.49806 0.5192
## PC22 PC23 PC24 PC25 PC26 PC27 PC28
## Standard deviation 4.56690 4.54509 4.5386 4.46552 4.45415 4.43642 4.39044
## Proportion of Variance 0.02086 0.02066 0.0206 0.01994 0.01984 0.01968 0.01928
## Cumulative Proportion 0.54002 0.56068 0.5813 0.60122 0.62106 0.64074 0.66001
## PC29 PC30 PC31 PC32 PC33 PC34 PC35
## Standard deviation 4.34322 4.31905 4.27409 4.26476 4.23175 4.1948 4.15706
## Proportion of Variance 0.01886 0.01865 0.01827 0.01819 0.01791 0.0176 0.01728
## Cumulative Proportion 0.67888 0.69753 0.71580 0.73399 0.75190 0.7695 0.78677
## PC36 PC37 PC38 PC39 PC40 PC41 PC42
## Standard deviation 4.10989 4.08520 4.05406 4.02783 4.01132 3.98323 3.92975
## Proportion of Variance 0.01689 0.01669 0.01644 0.01622 0.01609 0.01587 0.01544
## Cumulative Proportion 0.80367 0.82035 0.83679 0.85301 0.86910 0.88497 0.90041
## PC43 PC44 PC45 PC46 PC47 PC48 PC49
## Standard deviation 3.90651 3.85856 3.82445 3.78895 3.71912 3.68330 3.61347
## Proportion of Variance 0.01526 0.01489 0.01463 0.01436 0.01383 0.01357 0.01306
## Cumulative Proportion 0.91567 0.93056 0.94519 0.95954 0.97338 0.98694 1.00000
## PC50
## Standard deviation 5.041e-15
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00# Visualização dos dois primeiros componentes
plot(pca$x[,1], pca$x[,2], col = y_class,
xlab = "PC1", ylab = "PC2", pch = 19)
legend("topright", legend = levels(y_class),
col = 1:2, pch = 19)
- Estabilidade de seleção de variáveis (Bootstrap) :Aqui você avalia quais variáveis aparecem consistentemente.
library(glmnet)
n_boot <- 100
selected_vars <- matrix(0, nrow = n_boot, ncol = p)
for (i in 1:n_boot) {
idx <- sample(1:n, replace = TRUE)
X_boot <- X[idx, ]
y_boot <- y[idx]
model <- cv.glmnet(X_boot, y_boot, alpha = 1)
coef_vec <- coef(model, s = "lambda.min")[-1]
selected_vars[i, ] <- as.numeric(coef_vec != 0)
}
# Frequência de seleção
selection_freq <- colMeans(selected_vars)
selection_freq## [1] 0.02 0.00 0.02 0.77 0.15 0.00 0.45 0.30 0.19 0.13 0.00 0.00 0.38 0.01
## [15] 0.00 0.01 0.07 0.01 0.03 0.01 0.00 0.00 0.00 0.00 0.00 0.01 0.03 0.00
## [29] 0.00 0.01 0.00 0.00 0.11 0.48 0.00 0.00 0.14 0.00 0.01 0.10 0.03 0.00
## [43] 0.02 0.01 0.06 0.00 0.02 0.00 0.00 0.03 0.09 0.00 0.04 0.01 0.00 0.00
## [57] 0.01 0.04 0.02 0.00 0.00 0.61 0.00 0.01 0.00 0.00 0.01 0.00 0.01 0.03
## [71] 0.03 0.00 0.01 0.00 0.01 0.01 0.01 0.07 0.01 0.00 0.00 0.00 0.21 0.04
## [85] 0.00 0.00 0.02 0.00 0.02 0.02 0.02 0.00 0.05 0.09 0.00 0.13 0.00 0.01
## [99] 0.00 0.01 0.00 0.00 0.00 0.00 0.01 0.01 0.00 0.01 0.12 0.00 0.00 0.00
## [113] 0.00 0.01 0.00 0.00 0.04 0.02 0.01 0.04 0.00 0.01 0.02 0.01 0.00 0.00
## [127] 0.00 0.03 0.00 0.01 0.00 0.09 0.01 0.00 0.00 0.04 0.00 0.00 0.00 0.00
## [141] 0.10 0.01 0.00 0.00 0.00 0.03 0.05 0.01 0.00 0.00 0.00 0.04 0.01 0.01
## [155] 0.00 0.05 0.10 0.02 0.03 0.02 0.00 0.01 0.00 0.00 0.03 0.00 0.02 0.00
## [169] 0.05 0.07 0.02 0.02 0.07 0.00 0.00 0.02 0.00 0.12 0.00 0.00 0.00 0.00
## [183] 0.00 0.02 0.00 0.01 0.00 0.00 0.02 0.02 0.00 0.01 0.25 0.02 0.00 0.00
## [197] 0.00 0.58 0.00 0.01 0.06 0.00 0.00 0.00 0.00 0.04 0.07 0.02 0.00 0.00
## [211] 0.00 0.04 0.00 0.00 0.00 0.02 0.01 0.00 0.00 0.03 0.01 0.00 0.00 0.00
## [225] 0.00 0.02 0.00 0.02 0.00 0.03 0.00 0.00 0.00 0.04 0.00 0.04 0.00 0.00
## [239] 0.01 0.02 0.00 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.07 0.00 0.00
## [253] 0.00 0.01 0.02 0.00 0.00 0.00 0.02 0.02 0.07 0.01 0.04 0.00 0.01 0.05
## [267] 0.00 0.08 0.03 0.05 0.00 0.00 0.00 0.00 0.00 0.03 0.00 0.01 0.01 0.00
## [281] 0.02 0.01 0.00 0.01 0.01 0.00 0.01 0.00 0.00 0.00 0.04 0.02 0.02 0.00
## [295] 0.05 0.01 0.00 0.00 0.32 0.00 0.00 0.01 0.13 0.01 0.00 0.02 0.00 0.00
## [309] 0.00 0.00 0.00 0.01 0.03 0.00 0.01 0.02 0.00 0.24 0.00 0.00 0.00 0.00
## [323] 0.00 0.00 0.01 0.00 0.00 0.02 0.14 0.00 0.08 0.09 0.00 0.00 0.02 0.16
## [337] 0.00 0.07 0.00 0.08 0.00 0.00 0.01 0.00 0.01 0.00 0.00 0.00 0.01 0.00
## [351] 0.15 0.01 0.00 0.00 0.00 0.10 0.00 0.00 0.01 0.00 0.00 0.07 0.15 0.00
## [365] 0.00 0.01 0.00 0.01 0.02 0.00 0.07 0.00 0.01 0.00 0.00 0.00 0.04 0.01
## [379] 0.00 0.04 0.11 0.00 0.00 0.00 0.01 0.00 0.03 0.01 0.17 0.00 0.20 0.00
## [393] 0.04 0.03 0.02 0.00 0.00 0.02 0.02 0.00 0.00 0.01 0.01 0.00 0.19 0.11
## [407] 0.00 0.02 0.02 0.00 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
## [421] 0.01 0.04 0.02 0.00 0.02 0.00 0.04 0.03 0.00 0.00 0.00 0.00 0.10 0.00
## [435] 0.01 0.00 0.00 0.00 0.06 0.00 0.29 0.12 0.00 0.01 0.00 0.06 0.00 0.02
## [449] 0.02 0.08 0.02 0.00 0.00 0.08 0.00 0.00 0.00 0.00 0.02 0.00 0.00 0.00
## [463] 0.02 0.00 0.27 0.02 0.00 0.42 0.02 0.03 0.04 0.00 0.08 0.10 0.02 0.00
## [477] 0.00 0.00 0.01 0.00 0.05 0.02 0.01 0.00 0.10 0.00 0.01 0.00 0.01 0.03
## [491] 0.04 0.00 0.01 0.01 0.00 0.00 0.00 0.03 0.00 0.00 0.01 0.00 0.00 0.00
## [505] 0.00 0.00 0.00 0.10 0.00 0.01 0.00 0.00 0.01 0.07 0.03 0.03 0.06 0.13
## [519] 0.00 0.00 0.01 0.00 0.04 0.00 0.10 0.01 0.01 0.00 0.00 0.02 0.00 0.00
## [533] 0.00 0.01 0.04 0.00 0.00 0.00 0.05 0.00 0.01 0.01 0.00 0.04 0.01 0.01
## [547] 0.05 0.01 0.02 0.00 0.00 0.00 0.03 0.00 0.00 0.01 0.00 0.00 0.00 0.00
## [561] 0.02 0.00 0.00 0.02 0.00 0.10 0.01 0.08 0.00 0.02 0.02 0.01 0.00 0.00
## [575] 0.01 0.05 0.02 0.01 0.09 0.02 0.01 0.00 0.00 0.01 0.03 0.00 0.08 0.02
## [589] 0.17 0.00 0.01 0.00 0.13 0.03 0.01 0.00 0.00 0.00 0.07 0.00 0.05 0.00
## [603] 0.00 0.02 0.06 0.05 0.01 0.00 0.02 0.06 0.00 0.03 0.01 0.00 0.07 0.00
## [617] 0.01 0.00 0.00 0.01 0.00 0.11 0.00 0.00 0.03 0.00 0.03 0.01 0.07 0.00
## [631] 0.01 0.01 0.00 0.01 0.03 0.03 0.00 0.00 0.00 0.02 0.00 0.01 0.06 0.01
## [645] 0.00 0.00 0.00 0.00 0.01 0.00 0.01 0.00 0.00 0.01 0.00 0.00 0.01 0.04
## [659] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.06 0.00 0.00 0.00 0.01 0.00
## [673] 0.00 0.01 0.44 0.00 0.00 0.00 0.00 0.00 0.05 0.12 0.02 0.00 0.00 0.00
## [687] 0.02 0.30 0.04 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.11
## [701] 0.04 0.00 0.01 0.02 0.00 0.04 0.03 0.01 0.03 0.00 0.00 0.01 0.02 0.00
## [715] 0.00 0.03 0.00 0.00 0.01 0.05 0.00 0.00 0.00 0.00 0.00 0.08 0.00 0.02
## [729] 0.00 0.00 0.01 0.00 0.02 0.02 0.00 0.00 0.05 0.01 0.02 0.11 0.14 0.06
## [743] 0.00 0.04 0.00 0.00 0.01 0.00 0.04 0.09 0.05 0.01 0.00 0.04 0.04 0.00
## [757] 0.01 0.00 0.02 0.00 0.05 0.01 0.01 0.16 0.05 0.27 0.00 0.02 0.29 0.19
## [771] 0.01 0.00 0.00 0.16 0.00 0.01 0.04 0.00 0.04 0.00 0.02 0.02 0.46 0.00
## [785] 0.02 0.11 0.00 0.00 0.01 0.00 0.00 0.20 0.01 0.01 0.00 0.00 0.00 0.00
## [799] 0.00 0.00 0.26 0.02 0.05 0.00 0.00 0.01 0.01 0.00 0.04 0.00 0.00 0.00
## [813] 0.00 0.00 0.00 0.00 0.03 0.00 0.00 0.00 0.00 0.01 0.06 0.00 0.02 0.00
## [827] 0.02 0.05 0.00 0.01 0.00 0.00 0.06 0.01 0.00 0.00 0.00 0.02 0.32 0.00
## [841] 0.10 0.16 0.01 0.01 0.00 0.01 0.00 0.00 0.02 0.06 0.02 0.16 0.00 0.00
## [855] 0.00 0.00 0.00 0.04 0.01 0.00 0.02 0.00 0.01 0.00 0.00 0.00 0.01 0.03
## [869] 0.16 0.00 0.02 0.00 0.04 0.00 0.39 0.07 0.04 0.00 0.00 0.00 0.00 0.00
## [883] 0.01 0.00 0.04 0.00 0.01 0.00 0.00 0.07 0.02 0.00 0.00 0.00 0.00 0.00
## [897] 0.00 0.00 0.00 0.03 0.00 0.15 0.00 0.00 0.00 0.00 0.01 0.03 0.00 0.00
## [911] 0.01 0.00 0.04 0.00 0.01 0.09 0.00 0.02 0.01 0.05 0.00 0.00 0.06 0.03
## [925] 0.03 0.00 0.00 0.01 0.02 0.06 0.10 0.03 0.03 0.00 0.00 0.01 0.00 0.00
## [939] 0.00 0.00 0.01 0.01 0.00 0.05 0.01 0.00 0.03 0.00 0.01 0.01 0.18 0.00
## [953] 0.01 0.06 0.01 0.01 0.00 0.00 0.00 0.00 0.01 0.03 0.00 0.02 0.00 0.00
## [967] 0.01 0.03 0.03 0.02 0.00 0.02 0.00 0.00 0.01 0.06 0.00 0.01 0.00 0.00
## [981] 0.00 0.00 0.19 0.09 0.21 0.02 0.01 0.00 0.02 0.02 0.00 0.01 0.09 0.01
## [995] 0.00 0.00 0.02 0.15 0.01 0.01## [1] 4 62## [1] 2- SVM (máquinas de vetor de suporte) :Muito usado em classificação de dados genômicos.
library(e1071)
svm_model <- svm(X, y_class, kernel = "linear")
pred_svm <- predict(svm_model, X)
table(pred_svm, y_class)## y_class
## pred_svm Controle Doente
## Controle 25 0
## Doente 0 25Resumo das aplicações (Módulo 02: Técnicas de Permutação e Reamostragem)
As técnicas de permutação, reamostragem e validação cruzada possuem aplicações amplas e estratégicas, especialmente em contextos onde os métodos estatísticos clássicos apresentam limitações. De forma integrada, elas permitem realizar inferência estatística mais robusta, avaliar modelos preditivos e lidar com dados complexos — como os de alta dimensionalidade.
No âmbito inferencial, os testes de permutação e aleatorização são aplicados para comparar grupos sem depender de suposições paramétricas, sendo particularmente úteis em dados biológicos e de saúde, frequentemente não normais ou com amostras pequenas. Já os métodos de reamostragem como bootstrap e jackknife são empregados para estimar incerteza (erro padrão, viés e intervalos de confiança) de estimadores complexos, inclusive quando não há soluções analíticas disponíveis.
Em bioinformática e áreas correlatas, o bootstrap tem papel central na avaliação de robustez de inferências, como no suporte de árvores filogenéticas, estabilidade de redes gênicas e reprodutibilidade de assinaturas moleculares. O jackknife, por sua vez, é utilizado em cenários com menor custo computacional, principalmente para avaliar a influência de observações individuais.
No contexto de modelagem preditiva, a validação cruzada é fundamental para estimar o desempenho de modelos em dados não observados, auxiliando na prevenção de sobreajuste (overfitting). Estratégias como k-fold, LOOCV e validação aninhada são amplamente aplicadas em aprendizado de máquina, especialmente na seleção de modelos e ajuste de hiperparâmetros.
Em cenários de alta dimensionalidade (p >> n) — típicos de dados ômicos — essas técnicas tornam-se indispensáveis. Elas permitem avaliar a estabilidade de seleção de variáveis (por exemplo, genes relevantes), melhorar a generalização de modelos e garantir maior confiabilidade nos resultados. Frequentemente, são combinadas com métodos modernos como regularização (LASSO, Ridge), algoritmos de aprendizado de máquina (Random Forest, SVM) e técnicas de redução de dimensionalidade (PCA).
Em síntese, essas abordagens constituem um conjunto essencial de ferramentas para análise estatística contemporânea, proporcionando maior robustez, flexibilidade e confiabilidade na inferência e na modelagem de dados complexos.
Referências bibliográficas do módulo 02
EFRON, Bradley; TIBSHIRANI, Robert J. An introduction to the bootstrap. New York: Chapman & Hall, 1993.
DAVISON, Anthony C.; HINKLEY, David V. Bootstrap methods and their application. Cambridge: Cambridge University Press, 1997.
GOOD, Phillip I. Permutation tests: a practical guide to resampling methods for testing hypotheses. 3. ed. New York: Springer, 2005.
WASSERMAN, Larry. All of nonparametric statistics. New York: Springer, 2006.
HASTIE, Trevor; TIBSHIRANI, Robert; FRIEDMAN, Jerome. The elements of statistical learning: data mining, inference, and prediction. 2. ed. New York: Springer, 2009.
JAMES, Gareth; WITTEN, Daniela; HASTIE, Trevor; TIBSHIRANI, Robert. An introduction to statistical learning: with applications in R. 2. ed. New York: Springer, 2021.
BISHOP, Christopher M. Pattern recognition and machine learning. New York: Springer, 2006.
GENTLEMAN, Robert et al. Bioinformatics and computational biology solutions using R and Bioconductor. New York: Springer, 2005.
DATTA, Somnath; NETTLETON, Daniel S. Statistical analysis of next generation sequencing data. New York: Springer, 2014.
SPEED, Terry (ed.). Statistical analysis of gene expression data. New York: Chapman & Hall/CRC, 2003.
VENABLES, William N.; RIPLEY, Brian D. Modern applied statistics with S. 4. ed. New York: Springer, 2002.
KUHN, Max; JOHNSON, Kjell. Applied predictive modeling. New York: Springer, 2013.
DAVINACK, Andrew. Molecular ecology & evolution: an introduction. 2024. Disponível sob licença Creative Commons Attribution-NonCommercial 4.0 International.
MÓDULO 03: Analise Molecular De Variância (AMOVA)
3.1. Conceitos básicos de biologia molecular e variância molecular
Biologia molecular
A biologia molecular explora a estrutura e a função das moléculas que compõem as células vivas, particularmente o DNA, o RNA e as proteínas, que são os blocos de construção da vida. Este campo explica como essas moléculas interagem para realizar os processos essenciais para a sobrevivência e reprodução dos organismos. Na ecologia molecular, os pesquisadores aplicam esse conhecimento molecular fundamental para abordar questões sobre a dinâmica ecológica e evolutiva das populações.
Por exemplo, ao compreender a replicação do DNA e a expressão gênica, os ecólogos moleculares podem investigar variações genéticas dentro e entre populações, ajudando a desvendar os mecanismos moleculares da adaptação a diferentes condições ambientais. Essa abordagem é crucial para estudar como as espécies evoluem em resposta às pressões ecológicas e para identificar a base genética de características que contribuem para a sobrevivência e o sucesso reprodutivo em diversos ambientes.
Além disso, técnicas de biologia molecular, como o sequenciamento de DNA e a PCR (reação em cadeia da polimerase), são ferramentas indispensáveis na ecologia molecular. Essas técnicas permitem que os cientistas examinem a composição genética de organismos de diferentes ecossistemas, fornecendo informações sobre a estrutura populacional, a diversidade genética e a história evolutiva. Por exemplo, a capacidade de amplificar regiões específicas do DNA por meio da PCR possibilita o estudo detalhado de marcadores genéticos em indivíduos e populações, facilitando investigações sobre fluxo gênico, deriva genética e pressões seletivas.
Ou seja, um sólido conhecimento de biologia molecular é necessário para interpretar com precisão os resultados de estudos de ecologia molecular. Compreender as fontes de variação genética, como mutações, duplicação gênica e transferência horizontal de genes, permite que os ecólogos cheguem a conclusões fundamentadas sobre a estrutura genética das populações e suas respostas adaptativas às mudanças ambientais.
Variancia molecular
A Análise de Variância Molecular (AMOVA) é um método estatístico amplamente utilizado em ecologia molecular para quantificar a variação genética dentro e entre populações. Desenvolvida por Laurent Excoffier(Excoffier et al., 1992) no início da década de 1990, a AMOVA permite que pesquisadores particionem a diversidade genética em múltiplos níveis hierárquicos, como dentro de indivíduos, entre indivíduos dentro de populações e entre as próprias populações. É importante ter uma compreensão básica dessa análise, pois a utilizaremos em laboratório. A AMOVA baseia-se na premissa de que a variação genética pode ser decomposta em componentes correspondentes a diferentes níveis de organização biológica. Comparando sequências genéticas ou frequências alélicas entre esses níveis, a AMOVA calcula a proporção da variação genética total atribuível a cada nível. Essa partição é crucial para a compreensão da distribuição da diversidade genética e dos processos que impulsionam a estrutura populacional.
3.2. Estruturas hierárquicas de dados genéticos
Os dados de entrada da AMOVA, normalmente utiliza dados moleculares, como sequências de DNA, genótipos de microssatélites ou polimorfismos de nucleotídeo único (SNPs). Os dados de entrada são organizados de acordo com os níveis da estrutura hierárquica que está sendo estudada.
3.3. Matrizes de distância: construção e escolha do critério
Os cálculos de distância da AMOVA, baseia-se em calcular as distâncias genéticas ou dissimilaridades entre todos os pares de indivíduos ou haplótipos, frequentemente utilizando medidas como a estatística F, que se baseia no conceito de estatística F de Wright.
A AMOVA não trabalha diretamente com os dados brutos, mas com uma matriz de distâncias genéticas entre indivíduos ou haplótipos.
Sejam \(i\) e \(j\) dois indivíduos, a distância genética é definida como:
\[ d_{ij} \]
Dependendo do tipo de dado, pode-se utilizar:
Distância Euclidiana
\[ d_{ij} = \sum_{k=1}^{L} (x_{ik} - x_{jk})^2 \]
Distância de Hamming (para sequências)
\[ d_{ij} = \sum_{k=1}^{L} \delta(x_{ik}, x_{jk}) \]
onde:
\[ \delta(x_{ik}, x_{jk}) = \begin{cases} 0, & \text{se } x_{ik} = x_{jk} \\ 1, & \text{caso contrário} \end{cases} \]
Matriz de distância
\[ \mathbf{D} = \begin{bmatrix} 0 & d_{12} & \cdots & d_{1n} \\ d_{21} & 0 & \cdots & d_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ d_{n1} & d_{n2} & \cdots & 0 \end{bmatrix} \]
3.4. Partição da variância e estatísticas derivadas (Φ-statistics)
A variância genética total é particionada em componentes. Para uma hierarquia simples de dois níveis (dentro e entre populações, analogo a anova atradicional), a AMOVA calcula o FST, análogo ao FST de Wright, para medir a diferenciação genética entre populações. Para estruturas mais complexas, a AMOVA pode calcular o FSC (entre grupos dentro das populações) e o FCT (entre grupos).
A AMOVA decompõe a variância total em componentes hierárquicos.
Soma de quadrados total
\[ SS_T = \frac{1}{2n} \sum_{i=1}^{n} \sum_{j=1}^{n} d_{ij}^2 \]
Componentes de variância
- Variância entre populações:
\[ \sigma^2_{\text{entre}} \]
- Variância dentro das populações:
\[ \sigma^2_{\text{dentro}} \]
Variância total
\[ \sigma^2_T = \sigma^2_{\text{entre}} + \sigma^2_{\text{dentro}} \]
Deste modo, temos o que podemos chamar de Estatística ΦST (análogo ao FST)
\[ \Phi_{ST} = \frac{\sigma^2_{\text{entre}}}{\sigma^2_T} \]
Com poderamos a interpretação:
- \(\Phi_{ST} \approx 0\): baixa
diferenciação genética
- \(\Phi_{ST} \approx 1\): alta diferenciação genética
Assim, obtém-se uma estrutura com três níveis hierárquicos, sendo:
- Entre grupos:
\[ \Phi_{CT} = \frac{\sigma^2_{\text{grupos}}}{\sigma^2_T} \]
- Entre populações dentro dos grupos:
\[ \Phi_{SC} = \frac{\sigma^2_{\text{populações}}}{\sigma^2_{\text{populações}} + \sigma^2_{\text{dentro}}} \]
- Total:
\[ \Phi_{ST} = \frac{\sigma^2_{\text{grupos}} + \sigma^2_{\text{populações}}}{\sigma^2_T} \]
3.5. Formulação do modelo estatístico da AMOVA
Em termos de formulação do modelo estatístico, A AMOVA pode ser interpretada como uma extensão da ANOVA tradicional, porém é baseada em distâncias.
Temos que o modelo geral é dado por: \[ d_{ij}^2 = \mu + \alpha_g + \beta_{p(g)} + \epsilon_{ij} \]
onde:
- \(\mu\): média global
- \(\alpha_g\): efeito do grupo
- \(\beta_{p(g)}\): efeito da
população dentro do grupo
- \(\epsilon_{ij}\): erro residual
Com isso obtemos os Quadrados médios, sendo:
\[ MS = \frac{SS}{df} \]
- Quadrado médio entre populações:
\[ MS_{\text{entre}} \]
- Quadrado médio dentro:
\[ MS_{\text{dentro}} \]
e, por fim temos as estimativas dos componentes de variância: \[ \sigma^2_{\text{entre}} = \frac{MS_{\text{entre}} - MS_{\text{dentro}}}{n} \]
\[ \sigma^2_{\text{dentro}} = MS_{\text{dentro}} \]
3.6. Teste de hipóteses e significância
Para avaliar a significância dos componentes da variância, a AMOVA normalmente utiliza testes de permutação, realocando aleatoriamente indivíduos ou haplótipos para diferentes grupos e recalculando os componentes da variância para criar uma distribuição contra a qual os valores observados podem ser testados.
Em termos de (ho), ou seja Hipótese nula, têm-se:
\[ H_0: \sigma^2_{\text{entre}} = 0 \]
Ou seja, não há estrutura genética entre populações.
Quando se baseia em testes por permutação, temos que:
\[ p = \frac{\#(\Phi^* \geq \Phi_{\text{obs}}) + 1}{N_{\text{perm}} + 1} \]
onde:
- \(\Phi^*\): valores obtidos por
permutação
- \(\Phi_{\text{obs}}\): valor
observado
- \(N_{\text{perm}}\): número de permutações
3.7. Métodos de validação e abordagens por permutação
Reatribuição aleatória de indivíduos entre populações para gerar a distribuição nula.
- Bootstrap
\[ \hat{\theta}^{*} = \text{estatística recalculada em amostras reamostradas} \]
- Intervalo de confiança
\[ IC_{(1-\alpha)} = [\theta_{\alpha/2}, \theta_{1-\alpha/2}] \]
- Conexão com o FST clássico
\[ F_{ST} = \frac{H_T - H_S}{H_T} \]
onde:
- \(H_T\): heterozigosidade
total
- \(H_S\): heterozigosidade média dentro das populações
Observação didática:A AMOVA generaliza a ANOVA tradicional substituindo valores individuais por distâncias genéticas, permitindo a análise direta de dados moleculares sem assumir normalidade.
3.8. Aplicações em genética de populações e filogeografia
A filogeografia é uma disciplina científica que combina princípios da filogenética e da genética de populações para compreender a distribuição geográfica de linhagens genéticas em escalas espaciais e temporais . Ela explora como eventos históricos, como glaciações, formação de rios e soerguimento de montanhas, combinados com fatores ecológicos, moldaram a estrutura genética e a distribuição das espécies. Ao rastrear os padrões geográficos e históricos da variação genética, os pesquisadores podem inferir eventos demográficos passados, como expansões, contrações e migrações populacionais. Essas inferências dependem da análise de marcadores genéticos, como o DNA mitocondrial ou microssatélites, para revelar a história evolutiva das populações. Por exemplo, ao estudar uma espécie distribuída por diferentes regiões geográficas, a filogeografia ajuda a identificar linhagens genéticas distintas que podem corresponder a barreiras históricas ou períodos de isolamento . Essa informação aprimora nossa compreensão dos processos evolutivos, conectando dados genéticos com mudanças climáticas e geológicas passadas, fornecendo, em última análise, insights sobre padrões de biodiversidade e orientando esforços de conservação. Nesse exmplo, vamos trabalhar com os dados Aeut, são informações de microssatélites, ou seja são regiões do DNA que funcionam como ‘contadores de repetições’. Cada indivíduo tem um número diferente dessas repetições, e essa variação é usada como marcador genético para estudar diversidade e estrutura populacional.
################################################################
# EXEMPLO 1 — AMOVA SIMPLES (VERSÃO FUNCIONAL COM pacote PEGAS)
#Visão geral do que o código faz
#1-Carregando dados genéticos reais
#2-Definindo a estrutura hierárquica
#3-Preparando o objeto para análise
#4-Executando a AMOVA
#5-Testando significância com permutação
############################################################
# 1. Instalar e carregar pacotes necessários
#if(!require(poppr)) install.packages("poppr")
library(poppr)
# 2. Carregar dados de exemplo do próprio R (Aphanomyces euteiches - fungo)
#Aeut = Aphanomyces euteiches (patógeno de plantas) são dados de microssatélites
# com estrutura populacional real
data(Aeut)
# 3. Definir a hierarquia (Estrutura: Pop -> Subpop)
# O conjunto de dados original tem colunas de estratificação
strata(Aeut) <- other(Aeut)$population_hierarchy[-1]
# 4. Converter para genclone para análise
## Converte o objeto para um formato que: reconhece clonalidade;é compatível com poppr.amova()
agc <- as.genclone(Aeut)
# 5. Executar a AMOVA
# Hierarquia: ~Pop/Subpop (Populações dentro de Grupos)
amova_result <- poppr.amova(agc, ~Pop/Subpop)
# 6. Ver o resultado
print(amova_result)## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Pop 1 1051.2345 1051.234516
## Between samples Within Pop 16 273.4575 17.091091
## Within samples 169 576.5059 3.411277
## Total 186 1901.1979 10.221494
##
## $componentsofcovariance
## Sigma %
## Variations Between Pop 11.063446 70.006786
## Variations Between samples Within Pop 1.328667 8.407483
## Variations Within samples 3.411277 21.585732
## Total variations 15.803391 100.000000
##
## $statphi
## Phi
## Phi-samples-total 0.7841427
## Phi-samples-Pop 0.2803128
## Phi-Pop-total 0.7000679# 7. Testar a significância (permutações)
## Embaralha indivíduos entre grupos com isso recalcula AMOVA 999 vezes, depois ria distribuição nula.
signif_result <- randtest(amova_result, nrepet = 999)
plot(signif_result)
## class: krandtest lightkrandtest
## Monte-Carlo tests
## Call: randtest.amova(xtest = amova_result, nrepet = 999)
##
## Number of tests: 3
##
## Adjustment method for multiple comparisons: none
## Permutation number: 999
## Test Obs Std.Obs Alter Pvalue
## 1 Variations within samples 3.411277 -29.92366 less 0.001
## 2 Variations between samples 1.328667 21.00331 greater 0.001
## 3 Variations between Pop 11.063446 10.28113 greater 0.001- Resultados Principais dos Componentes do Resultado em R* Variância entre grupos/populações: Indica a diferenciação genética. ou que os Índices de fixação que medem a estrutura populacional. Significância (P-value): Calculada com randtest, verifica se a estrutura é real.
A poppr.amova funciona melhor com objetos genind ou genclone e pode usar tanto ade4 quanto pegas internamente para realizar os cálculos de variância.
Entendendo passo a passo as Interpretação dos resultados da AMOVA
Sobre os dados analisados: Esse conjunto de dados representa 187 indivíduos de um patógeno, analisados em 56 regiões do DNA. Os dados são do tipo presença/ausência, ou seja, indicam se um fragmento genético está presente ou não,ss dados Aeut são marcadores AFLP (presença/ausência). Esses indivíduos estão organizados em duas populações, coletadas em regiões diferentes. A partir disso, usamos a AMOVA para avaliar se a variação genética está mais dentro das populações ou entre elas.
## /// GENIND OBJECT /////////
##
## // 187 individuals; 56 loci; 56 alleles; size: 72.7 Kb
##
## // Basic content
## @tab: 187 x 56 matrix of allele counts
## @loc.n.all: number of alleles per locus (range: 56-56)
## @ploidy: ploidy of each individual (range: 2-2)
## @type: PA
## @call: old2new_genind(object = x, donor = new(class(x)))
##
## // Optional content
## @pop: population of each individual (group size range: 90-97)
## @strata: a data frame with 2 columns ( Pop, Subpop )
## @other: a list containing: population_hierarchyEm termos de tradução simples: 187 indivíduos → 187 fungos diferentes 56 loci → 56 “pontos do DNA” analisados 56 alelos → basicamente um marcador por locus
Ou seja, cada indivíduo foi analisado em 56 regiões do DNA. Cada região pode ter variação, e isso nos permite comparar geneticamente os indivíduos.
Nesse ponto: @tab: 187 x 56 matrix of allele counts Signofica que temos uma matriz, onde linhas = indivíduos, e colunas = loci Logo, cada valor indica presença/ausência do marcador.
Depois:@type: PA Que significa PA = Presence/Absence (presença/ausência).
Então, Isso NÃO é microssatélite clássico (com tamanhos tipo 150/152), na verdade estamos trabalhando dados AFLP (dominantes)
E, o que isso significa biologicamente? Para cada locus:
1 → fragmento presente 0 → fragmento ausente
Em vez de contar repetições (como nos microssatélites), aqui estamos verificando se um pedaço do DNA aparece ou não.
Por fim, @ploidy: 2 Significa: organismo diploide, ou seja, tem duas cópias do DNA.
Em termos @pop: population of each individual (group size range: 90-97) Tem-se 2 populações, aproximadamente 90–97 indivíduos em cada. Os indivíduos foram coletados em dois locais diferentes — e queremos saber se eles são geneticamente diferentes.
Existe a hierarquia: @other: population_hierarchy Trazendo informações extras, como:subpopulações que serão os níveis hierárquicos usados na AMOVA.
Dados de Aphanomyces euteiches, um patógeno que causa podridão radicular. A coleta dos dados acnteceu em: Oregon (EUA) e Washington (EUA). A ideia geral é a interpretação ecológica, então estamos testando se:“Os fungos desses dois locais são geneticamente diferentes?”
Por critérios didáticos, é importante definir nesse exemplo que Microssatélite vs AFLP, ambos são marcadores genéticos, ou seja, formas de “ler diferenças no DNA”,mas eles funcionam de maneiras bem diferentes, conforme tabela abaixo.
| Característica | Microssatélites (SSR) | AFLP |
|---|---|---|
| Tipo de dado | Número de repetições | Presença / ausência (0/1) |
| Tipo de marcador | Codominante | Dominante |
| Informação genética | Alta (distingue heterozigotos) | Menor (não distingue heterozigotos) |
| Forma dos dados | Ex: 150/154 | Ex: 1 ou 0 |
| Interpretação | Quantitativa | Binária |
| Resolução genética | Alta | Moderada |
| Complexidade | Maior | Mais simples |
| Custo/tempo | Mais caro | Mais rápido e barato |
| Uso comum | Estrutura fina, parentesco | Diversidade geral, estrutura |
3.8.1. Partição da variância genética
A análise de variância molecular (AMOVA) revelou como a variação genética está distribuída entre os diferentes níveis hierárquicos analisados.
Observa-se que:
- 70,01% da variação genética ocorre entre populações (Between
Pop)
- 8,41% ocorre entre amostras dentro das populações (Between
samples Within Pop)
- 21,59% ocorre dentro das amostras (Within samples)
Esses resultados indicam que a maior parte da diversidade genética está associada às diferenças entre populações.
3.8.2. Interpretação biológica da estrutura genética
A predominância da variação entre populações sugere que:
- Existe forte estruturação genética
- As populações são geneticamente diferenciadas entre
si
- O fluxo gênico entre populações é limitado
Por outro lado, a variação dentro das amostras (21,59%) indica que ainda há diversidade genética interna, mas em menor proporção quando comparada à diferenciação entre populações.
3.8.3. Estatísticas Φ (Phi)
As estatísticas Φ, análogas aos índices F de Wright, reforçam essa interpretação:
- ΦST = 0,784 → indica alta diferenciação
genética total
- ΦPop-total = 0,700 → evidencia forte diferenciação
entre populações
- Φsamples-Pop = 0,280 → indica diferenciação moderada entre amostras dentro das populações
Valores elevados de ΦST (> 0,25) são geralmente interpretados como evidência de estrutura genética acentuada.
3.8.4. Teste de significância (permutações)
Os testes de permutação apresentaram:
- p = 0,001 para todos os níveis de variação
Isso indica que os componentes de variância observados são estatisticamente significativos e altamente improváveis de ocorrer ao acaso.
3.8.5. Interpretação dos gráficos de permutação
Os histogramas representam a distribuição dos valores esperados sob a hipótese nula (ausência de estrutura genética), enquanto o ponto preto indica o valor observado.
- Para variação dentro das amostras, o valor
observado é menor que o esperado ao acaso
- Para variação entre amostras e entre populações, os valores observados são significativamente maiores que os esperados
Isso reforça que a estrutura genética detectada é real e não resultado de variação aleatória, veja resumo, conforme figura 01.
Figura 1: Síntese do módulo e estudo de casoFonte: Obtido pelo notebook LM-2026].
3.8.6. Conclusão geral
A AMOVA evidencia que a maior parte da variação genética está concentrada entre populações, indicando forte diferenciação genética e baixa conectividade entre elas.
Esse padrão é consistente com cenários de:
- isolamento geográfico
- barreiras ao fluxo gênico
- história evolutiva diferenciada entre populações
Portanto, conclui-se que as populações analisadas apresentam uma estrutura genética bem definida e estatisticamente suportada.
Referências bibliográficas do módulo 03
Davinack, A. (2024). Molecular Ecology & Evolution: An Introduction. Open educational resource licensed under Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0).
Excoffier, L., Smouse, P. E., & Quattro, J. M. (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics, 131(2), 479–491.
Wright, S. (1951). The genetical structure of populations. Annals of Eugenics, 15(1), 323–354.
Nei, M. (1987). Molecular Evolutionary Genetics. New York: Columbia University Press.
Hedrick, P. W. (2011). Genetics of populations (4th ed.). Sudbury, MA: Jones & Bartlett Learning.
Avise, J. C. (2000). Phylogeography: The history and formation of species. Cambridge, MA: Harvard University Press.
Peakall, R., & Smouse, P. E. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research. Bioinformatics, 28(19), 2537–2539.
Kamvar, Z. N., Tabima, J. F., & Grünwald, N. J. (2014). Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ, 2, e281.
Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics, 24(11), 1403–1405.
Excoffier, L., & Lischer, H. E. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources, 10(3), 564–567.
Selkoe, K. A., & Toonen, R. J. (2006). Microsatellites for ecologists: a practical guide to using and evaluating microsatellite markers. Ecology Letters, 9(5), 615–629.
Vos, P., Hogers, R., Bleeker, M., Reijans, M., van de Lee, T., Hornes, M., et al. (1995). AFLP: a new technique for DNA fingerprinting. Nucleic Acids Research, 23(21), 4407–4414.
Lowe, A., Harris, S., & Ashton, P. (2004). Ecological genetics: design, analysis, and application. Oxford: Blackwell Publishing.
MÓDULO 04. Modelos Lineares Generalizados (GLM)
Notas das autoras/professoras: Referência básica: Cordeiro, Demétrio (2008)
Os Modelos Lineares Generalizados constituem uma das principais ferramentas da Estatística moderna para análise de dados, especialmente em situações nas quais as suposições do modelo linear clássico não são adequadas. Em muitos problemas reais, a variável resposta não apresenta distribuição normal, podendo assumir a forma de contagens, proporções ou dados assimétricos, além de frequentemente apresentar variância dependente da média. Nesses casos, o uso de modelos lineares tradicionais torna-se limitado ou inadequado, exigindo abordagens mais flexíveis. Por exemplo, são apropriados quando se deseja modelar o número de ocorrências de um evento (dados de contagem), a probabilidade de sucesso em um experimento (dados binários ou proporcionais) ou ainda variáveis contínuas com distribuição não normal.
A aplicação de um Modelo Linear Generalizado requer a definição de três componentes fundamentais. Primeiramente, é necessário escolher uma distribuição adequada para a variável resposta, geralmente pertencente à família exponencial, como as distribuições normal, binomial ou de Poisson. Em seguida, deve-se especificar o componente sistemático do modelo, representado por um preditor linear que relaciona as variáveis explicativas aos parâmetros do modelo. Por fim, é preciso definir uma função de ligação, responsável por conectar a média da variável resposta ao preditor linear.
Dessa forma, os GLMs oferecem uma estrutura unificada e flexível, capaz de acomodar diferentes tipos de dados e relações, tornando-se uma ferramenta indispensável na análise estatística contemporânea.
4.1. Revisão do conceito de modelos lineares e suas limitações
Diversas distribuições de probabilidade amplamente utilizadas na Estatística podem ser reunidas em uma classe denominada família exponencial de distribuições. Essa classe inclui distribuições como normal, binomial, Poisson, gama e binomial negativa, entre outras.
A importância da família exponencial decorre de suas propriedades matemáticas, que permitem uma formulação unificada de diversos modelos estatísticos. Em particular, essa classe desempenha papel central na teoria dos Modelos Lineares Generalizados (MLG), possibilitando a modelagem de diferentes tipos de dados dentro de um mesmo arcabouço teórico.
Além disso, a família exponencial apresenta propriedades importantes, como a existência de estatísticas suficientes de dimensão reduzida e relações simples entre momentos e parâmetros, o que facilita tanto a inferência quanto a estimação.
Família Exponencial Uniparamétrica
Uma distribuição pertence à família exponencial uniparamétrica se sua função densidade (ou de probabilidade) pode ser escrita na forma:
\[ f(x; \theta) = h(x) \exp \left[ \eta(\theta) \, t(x) - b(\theta) \right] \]
onde: - \(\theta\) é o parâmetro; - \(\eta(\theta)\) é o parâmetro canônico; - \(t(x)\) é uma estatística suficiente; - \(b(\theta)\) é uma função de normalização; - \(h(x)\) é uma função que não depende de \(\theta\).
Uma forma especial importante é a forma canônica, dada por:
\[ f(x; \theta) = h(x) \exp \left[ \theta x - b(\theta) \right] \]
Nessa parametrização: - \(\theta\) é chamado de parâmetro canônico; - a função \(b(\theta)\) está diretamente relacionada aos momentos da distribuição.
As propriedades fundamentais incluem:
\[ E(X) = b'(\theta) \]
\[ \text{Var}(X) = b''(\theta) \]
Componente Aleatório
Na formulação dos Modelos Lineares Generalizados, o componente aleatório é definido a partir da família exponencial, incorporando um parâmetro de dispersão \(\phi > 0\):
\[ f(y; \theta, \phi) = \exp \left\{ \frac{1}{\phi} [y\theta - b(\theta)] + c(y, \phi) \right\} \]
Nesse contexto:
- \(\mu = E(Y) = b'(\theta)\)
- \(\text{Var}(Y) = \phi \, b''(\theta)\)
Define-se ainda a função de variância:
\[ V(\mu) = b''(\theta) \]
O parâmetro \(\phi\) controla a dispersão dos dados: - \(\phi = 1\) em distribuições como Poisson e binomial; - \(\phi \neq 1\) em distribuições como normal e gama.
Função Geradora de Momentos
A função geradora de momentos (FGM) de uma variável aleatória \(Y\) na família exponencial é dada por:
\[ M(t) = E(e^{tY}) \]
Para a forma geral da família exponencial, a FGM pode ser expressa como:
\[ M(t) = \exp \left\{ \frac{1}{\phi} \left[ b(\theta + \phi t) - b(\theta) \right] \right\} \]
A função geradora de cumulantes é definida como:
\[ K(t) = \log M(t) \]
e permite obter os cumulantes da distribuição por diferenciação sucessiva.
Em particular: - o primeiro cumulante é a média; - o segundo cumulante é a variância.
Essa propriedade é essencial para o desenvolvimento teórico dos MLG, especialmente em resultados assintóticos.
Estatística Suficiente
Uma estatística \(T(Y)\) é dita suficiente para o parâmetro \(\theta\) se contém toda a informação sobre \(\theta\) presente na amostra.
Para distribuições da família exponencial, existe uma estatística suficiente de dimensão reduzida, geralmente da forma:
\[ T = \sum_{i=1}^{n} t(Y_i) \]
Essa propriedade decorre do teorema da fatoração de Neyman-Fisher, que estabelece que a função de verossimilhança pode ser decomposta em uma parte que depende apenas de \(T\) e outra independente de \(\theta\).
Essa característica é extremamente importante, pois reduz a complexidade da inferência estatística.
Família Exponencial Multiparamétrica
A família exponencial pode ser generalizada para múltiplos parâmetros, sendo definida por:
\[ f(x; \boldsymbol{\theta}) = h(x) \exp \left[ \sum_{i=1}^{k} \eta_i(\boldsymbol{\theta}) t_i(x) - b(\boldsymbol{\theta}) \right] \]
onde: - \(\boldsymbol{\theta}\) é um vetor de parâmetros; - \(t_i(x)\) são estatísticas suficientes; - \(\eta_i(\boldsymbol{\theta})\) são parâmetros canônicos.
Na forma canônica:
\[ f(x; \boldsymbol{\theta}) = h(x) \exp \left[ \sum_{i=1}^{k} \theta_i t_i(x) - b(\boldsymbol{\theta}) \right] \]
As propriedades dos momentos são generalizadas como:
\[ E[T_i(X)] = \frac{\partial b(\boldsymbol{\theta})}{\partial \theta_i} \]
\[ \text{Cov}(T_i, T_j) = \frac{\partial^2 b(\boldsymbol{\theta})}{\partial \theta_i \partial \theta_j} \]
Essa generalização permite modelar distribuições mais complexas, como a multinomial e a normal com parâmetros desconhecidos, ampliando o escopo da modelagem estatística.
4.2. Família exponencial e função de ligação
A modelagem estatística desempenha um papel central na análise de dados em diversas áreas do conhecimento, como agricultura, economia, medicina e engenharia. Um dos principais desafios consiste na seleção de modelos que sejam suficientemente simples, mas capazes de representar adequadamente os dados observados.
Nesse contexto, os Modelos Lineares Generalizados (MLG), introduzidos por Nelder e Wedderburn (1972), constituem uma estrutura unificadora que abrange diferentes técnicas estatísticas tradicionalmente tratadas de forma separada. Esses modelos estendem os modelos clássicos de regressão ao permitir maior flexibilidade na distribuição da variável resposta e na relação entre variáveis.
Os MLG são definidos a partir de três componentes principais:
Componente aleatório: a variável resposta segue uma distribuição pertencente à família exponencial, incluindo distribuições como normal, binomial, Poisson e gama.
Componente sistemático: corresponde a uma combinação linear das variáveis explicativas, usualmente representada por:
\[ \eta = X\beta \]
Função de ligação: estabelece a relação entre a média da variável resposta e o preditor linear, podendo assumir diferentes formas, como identidade, logarítmica ou logística.
No modelo clássico de regressão linear, assume-se que a variável resposta segue uma distribuição normal com média igual ao preditor linear e variância constante. Nesse caso, a função de ligação é a identidade.
Entretanto, essa estrutura nem sempre é adequada, especialmente quando os dados são discretos, assimétricos ou apresentam variância não constante. Os MLG surgem, portanto, como uma generalização natural, permitindo o uso de diferentes distribuições e funções de ligação, ampliando significativamente o escopo da modelagem estatística.
Diversos métodos estatísticos podem ser interpretados como casos particulares dos MLG, incluindo regressão linear, regressão logística, modelos de Poisson e análise de variância. Essa unificação torna os MLG uma ferramenta fundamental na Estatística moderna.
Definição dos Modelos Lineares Generalizados
Os Modelos Lineares Generalizados (MLG) constituem uma extensão dos modelos lineares clássicos, permitindo a modelagem de variáveis resposta que seguem distribuições pertencentes à família exponencial. Essa estrutura amplia significativamente a aplicabilidade dos modelos estatísticos, especialmente em situações onde as suposições de normalidade e homocedasticidade não são adequadas.
Formalmente, um MLG é definido por três componentes fundamentais:
Componente Aleatório
Seja \(Y_1, Y_2, \ldots, Y_n\) um conjunto de variáveis aleatórias independentes, em que cada \(Y_i\) segue uma distribuição pertencente à família exponencial, com função densidade (ou de probabilidade) dada por:
\[ f(y_i; \theta_i, \phi) = \exp \left\{ \frac{1}{\phi} \left[ y_i \theta_i - b(\theta_i) \right] + c(y_i, \phi) \right\} \]
onde: - \(\theta_i\) é o parâmetro canônico; - \(\phi > 0\) é o parâmetro de dispersão; - \(b(\cdot)\) e \(c(\cdot)\) são funções conhecidas.
As propriedades dessa família implicam que:
\[ E(Y_i) = \mu_i = b'(\theta_i) \]
\[ \text{Var}(Y_i) = \phi \, b''(\theta_i) = \phi \, V(\mu_i) \]
em que \(V(\mu_i)\) é a função de variância.
Componente Sistemático
O componente sistemático é definido por um preditor linear:
\[ \eta_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} \]
ou, de forma matricial:
\[ \eta = X\beta \]
onde: - \(X\) é a matriz de delineamento (dimensão \(n \times p\)); - \(\beta\) é o vetor de parâmetros desconhecidos; - \(\eta_i\) é o preditor linear associado à observação \(i\).
Função de Ligação
A relação entre a média da variável resposta \(\mu_i\) e o preditor linear \(\eta_i\) é estabelecida por meio de uma função de ligação \(g(\cdot)\), tal que:
\[ g(\mu_i) = \eta_i \]
ou, equivalentemente:
\[ \mu_i = g^{-1}(\eta_i) \]
A escolha da função de ligação depende da natureza dos dados e da distribuição assumida. Exemplos comuns incluem: - Função identidade: \(g(\mu) = \mu\) - Função logarítmica: \(g(\mu) = \log(\mu)\) - Função logística (logit): \(g(\mu) = \log\left( \frac{\mu}{1 - \mu} \right)\)
Estrutura Geral do Modelo
Portanto, um Modelo Linear Generalizado é completamente especificado pela escolha de: 1. Uma distribuição da família exponencial para \(Y_i\); 2. Um preditor linear \(\eta_i = X\beta\); 3. Uma função de ligação \(g(\cdot)\) que relaciona \(\mu_i\) e \(\eta_i\).
Essa estrutura permite modelar diferentes tipos de dados, como: - Dados contínuos (distribuição normal, gama); - Dados de contagem (distribuição de Poisson); - Dados binários ou proporcionais (distribuição binomial).
Exemplo em R: Ajuste de um Modelo Linear Generalizado com
glm()
A função glm() do R permite ajustar Modelos Lineares
Generalizados de forma simples, especificando a distribuição da variável
resposta e a função de ligação.
Neste exemplo, consideramos uma variável resposta binária (0 = não, 1 = sim), modelada em função de uma variável explicativa \(x\).
# Gerando dados simulados
set.seed(123)
n <- 100
x <- runif(n, 0, 10)
# Probabilidade verdadeira (função logística)
p <- 1 / (1 + exp(-( -2 + 0.5 * x )))
# Variável resposta binária
y <- rbinom(n, size = 1, prob = p)
# Ajuste do modelo logístico
modelo_logistico <- glm(y ~ x, family = binomial(link = "logit"))
# Resumo do modelo
summary(modelo_logistico)##
## Call:
## glm(formula = y ~ x, family = binomial(link = "logit"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0006 0.5346 -3.742 0.000182 ***
## x 0.5867 0.1194 4.915 8.86e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 130.684 on 99 degrees of freedom
## Residual deviance: 91.716 on 98 degrees of freedom
## AIC: 95.716
##
## Number of Fisher Scoring iterations: 5Observação
Os MLG generalizam o modelo linear clássico ao permitir: - distribuições não normais para a variável resposta; - relações não lineares entre a média e o preditor linear; - variâncias dependentes da média.
Essas características tornam os MLG uma ferramenta extremamente flexível e amplamente utilizada na análise de dados em diversas áreas do conhecimento.
4.3. Modelo de regressão binomial e a extensão multinomial
O modelo binomial é outro caso particular fundamental dos Modelos Lineares Generalizados, sendo especialmente adequado para a análise de dados em que a variável resposta representa o número de sucessos em um número fixo de tentativas independentes. Esse tipo de modelo é amplamente utilizado em situações envolvendo proporções, probabilidades ou respostas dicotômicas (sucesso/fracasso, presença/ausência, sim/não).
Seja \(Y\) uma variável aleatória com distribuição binomial com parâmetros \(m\) (número de tentativas) e \(\pi\) (probabilidade de sucesso em cada tentativa). A função de probabilidade é dada por
\[ P(Y = y) = \binom{m}{y} \pi^y (1 - \pi)^{m - y}, \quad y = 0,1,\ldots,m. \]
O valor esperado e a variância de \(Y\) são, respectivamente,
\[ E(Y) = m\pi \]
e
\[ \mathrm{Var}(Y) = m\pi(1 - \pi). \]
Quando se trabalha com proporções, é comum considerar a variável \(Y/m\), cujo valor esperado é \(\pi\) e cuja variância é \(\pi(1 - \pi)/m\). Assim, a variabilidade depende da média, característica típica dos modelos da família exponencial.
No contexto dos Modelos Lineares Generalizados, a variável resposta é modelada por meio de sua média \(\mu = E(Y) = m\pi\), sendo a probabilidade \(\pi\) relacionada ao preditor linear por meio de uma função de ligação apropriada.
A função de ligação canônica para o modelo binomial é a função logito, definida por
\[ g(\pi) = \log\left(\frac{\pi}{1 - \pi}\right). \]
Dessa forma, o modelo pode ser escrito como
\[ \log\left(\frac{\pi_i}{1 - \pi_i}\right) = \eta_i = x_i^T \beta, \]
onde \(x_i^T\) representa o vetor de covariáveis da i-ésima observação e \(\beta\) é o vetor de parâmetros desconhecidos.
A inversa da função de ligação fornece
\[ \pi_i = \frac{\exp(\eta_i)}{1 + \exp(\eta_i)}, \]
o que garante que as probabilidades estimadas estejam sempre no intervalo \((0,1)\).
O modelo logístico, obtido a partir da função de ligação logito, é amplamente utilizado na prática devido à sua interpretação conveniente. Os coeficientes \(\beta\) podem ser interpretados em termos de log-odds, isto é, logaritmo da razão de chances. Mais especificamente, uma variação unitária em uma covariável está associada a uma mudança aditiva no logaritmo da razão de chances.
Além da função logito, outras funções de ligação podem ser utilizadas, como a função probito,
\[ g(\pi) = \Phi^{-1}(\pi), \]
onde \(\Phi\) é a função de distribuição acumulada da normal padrão, e a função complemento log-log,
\[ g(\pi) = \log[-\log(1 - \pi)]. \]
Essas funções alternativas são úteis em diferentes contextos, dependendo da natureza dos dados e da interpretação desejada.
Em aplicações práticas, os dados binomiais podem surgir de duas formas principais:
- dados binomiais agrupados, em que se observa o número de sucessos em um total conhecido de tentativas;
- dados binários individuais, em que cada observação assume valor 0 ou 1.
No primeiro caso, a resposta é frequentemente especificada na forma \((y_i, m_i - y_i)\), enquanto no segundo caso utiliza-se diretamente a variável binária.
O modelo binomial é amplamente aplicado em diversas áreas, incluindo estudos epidemiológicos, experimentos biológicos, análises de sobrevivência simplificadas e problemas de classificação. Em experimentos dose–resposta, por exemplo, é comum modelar a probabilidade de sucesso (como morte ou resposta a um tratamento) em função da dose aplicada.
Assim como no modelo de Poisson, a adequação do modelo binomial pode ser avaliada por meio da análise do desvio e de estatísticas baseadas nos resíduos, como o qui-quadrado de Pearson. Essas ferramentas permitem verificar se o modelo ajustado descreve adequadamente os dados observados.
Em síntese, o modelo binomial desempenha um papel central na modelagem de dados categóricos e de proporções. Sua flexibilidade, aliada à possibilidade de diferentes funções de ligação, torna-o uma ferramenta essencial dentro do arcabouço dos Modelos Lineares Generalizados.
Exemplo: Cypermethrin
O exemplo do Cypermethrin ilustra a aplicação de Modelos Lineares Generalizados com distribuição binomial, sendo um caso clássico de modelagem de dados de proporções (ou respostas binárias agrupadas).
Nesse contexto, a variável resposta corresponde ao número de insetos mortos em um total fixo de indivíduos expostos a diferentes doses de um inseticida. Assim, cada observação pode ser interpretada como o resultado de um experimento binomial, no qual se modela a probabilidade de morte dos insetos em função de variáveis explicativas.
O principal objetivo da análise é investigar como a probabilidade de mortalidade varia com a dose aplicada e verificar se essa relação difere entre os sexos dos insetos. Para isso, são ajustados diferentes modelos, variando:
-a forma de inclusão da dose (como fator ou como variável quantitativa); -a presença ou não do efeito de sexo; -a possibilidade de interação entre dose e sexo.
Um aspecto central desse exemplo é a utilização da função de ligação logito, que transforma a probabilidade de sucesso (morte) em uma escala linear, permitindo o uso de técnicas de regressão. Essa transformação é fundamental, pois a relação entre dose e probabilidade geralmente apresenta formato sigmoide, típico de experimentos do tipo dose–resposta.
Além disso, o exemplo evidencia a importância da escolha da estrutura do modelo. Ao tratar a dose como variável quantitativa (por meio do logaritmo da dose), obtém-se um modelo mais parsimonioso e interpretável, capaz de capturar a tendência dos dados com menor número de parâmetros.
A comparação entre modelos, realizada por meio da análise de desvio, permite avaliar a significância dos efeitos incluídos e identificar a especificação mais adequada. Em particular, a análise pode indicar se há diferenças entre os sexos ou se as curvas dose–resposta são paralelas.
Por fim, os gráficos com curvas ajustadas desempenham papel fundamental na interpretação dos resultados, permitindo visualizar o comportamento da probabilidade de mortalidade em função da dose e comparar diretamente os grupos analisados.
# ==========================================
# Exemplo - Cypermethrin
# ==========================================
# Dados
y <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
sex <- factor(rep(c("M", "F"), c(6, 6)))
ldose <- rep(0:5, 2)
dose <- 2^ldose
dose <- factor(dose)
Cyper.dat <- data.frame(sex, dose, ldose, y)
# Visualização inicial
plot(
ldose, y / 20,
pch = c(rep("*", 6), rep("+", 6)),
col = c(rep("green", 6), rep("red", 6)),
xlab = "log(dose)",
ylab = "Proportion killed"
)
# Resposta binomial agrupada
resp <- cbind(y, 20 - y)
# Modelos com dose tratada como fator
mod1 <- glm(resp ~ 1, family = binomial, data = Cyper.dat)
mod2 <- glm(resp ~ dose, family = binomial, data = Cyper.dat)
mod3 <- glm(resp ~ sex, family = binomial, data = Cyper.dat)
mod4 <- glm(resp ~ dose + sex, family = binomial, data = Cyper.dat)
anova(mod1, mod2, mod4, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: resp ~ 1
## Model 2: resp ~ dose
## Model 3: resp ~ dose + sex
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 11 124.876
## 2 6 15.152 5 109.72 < 2.2e-16 ***
## 3 5 5.013 1 10.14 0.001451 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table
##
## Model 1: resp ~ 1
## Model 2: resp ~ sex
## Model 3: resp ~ dose + sex
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 11 124.876
## 2 10 118.799 1 6.077 0.0137 *
## 3 5 5.013 5 113.786 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelos com log(dose) como covariável quantitativa
mod5 <- glm(resp ~ ldose, family = binomial, data = Cyper.dat)
mod6 <- glm(resp ~ sex + ldose - 1, family = binomial, data = Cyper.dat)
mod7 <- glm(resp ~ ldose / sex, family = binomial, data = Cyper.dat)
mod8 <- glm(resp ~ ldose * sex, family = binomial, data = Cyper.dat)
anova(mod1, mod5, mod6, mod8, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: resp ~ 1
## Model 2: resp ~ ldose
## Model 3: resp ~ sex + ldose - 1
## Model 4: resp ~ ldose * sex
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 11 124.876
## 2 10 16.984 1 107.892 < 2.2e-16 ***
## 3 9 6.757 1 10.227 0.001384 **
## 4 8 4.994 1 1.763 0.184209
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table
##
## Model 1: resp ~ 1
## Model 2: resp ~ ldose
## Model 3: resp ~ ldose/sex
## Model 4: resp ~ ldose * sex
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 11 124.876
## 2 10 16.984 1 107.892 < 2.2e-16 ***
## 3 9 5.044 1 11.940 0.0005495 ***
## 4 8 4.994 1 0.051 0.8221411
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## glm(formula = resp ~ sex + ldose - 1, family = binomial, data = Cyper.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## sexF -3.4732 0.4685 -7.413 1.23e-13 ***
## sexM -2.3724 0.3855 -6.154 7.56e-10 ***
## ldose 1.0642 0.1311 8.119 4.70e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 126.2269 on 12 degrees of freedom
## Residual deviance: 6.7571 on 9 degrees of freedom
## AIC: 42.867
##
## Number of Fisher Scoring iterations: 4# Gráfico com curvas ajustadas
plot(
c(1, 32), c(0, 1),
type = "n",
xlab = "log(dose)",
ylab = "Proportions",
log = "x"
)
points(
2^ldose, y / 20,
pch = c(rep("*", 6), rep("+", 6)),
col = c(rep("green", 6), rep("red", 6))
)
ld <- seq(0, 5, 0.1)
lines(
2^ld,
predict(
mod6,
newdata = data.frame(
ldose = ld,
sex = factor(rep("M", length(ld)), levels = levels(Cyper.dat$sex))
),
type = "response"
),
col = "green"
)
lines(
2^ld,
predict(
mod6,
newdata = data.frame(
ldose = ld,
sex = factor(rep("F", length(ld)), levels = levels(Cyper.dat$sex))
),
type = "response"
),
col = "red"
)
Exemplo: Tribolium (múltiplos tratamentos)
O exemplo do Tribolium castaneum representa uma aplicação mais avançada de Modelos Lineares Generalizados com distribuição binomial, sendo particularmente relevante para a análise de experimentos dose–resposta com múltiplos tratamentos.
Nesse contexto, a variável resposta corresponde ao número de insetos mortos em um total fixo de indivíduos, caracterizando dados binomiais agrupados. O estudo considera diferentes níveis de dose e três tipos de inseticidas (DDT, BHC e a combinação DDT+BHC), permitindo investigar simultaneamente o efeito da dose e do tipo de tratamento sobre a probabilidade de mortalidade.
O principal objetivo da análise é avaliar:
-se a mortalidade depende da dose aplicada; -se há diferenças entre os inseticidas; -e se a relação entre dose e resposta é semelhante (ou não) entre os tratamentos.
Para isso, são ajustados diversos modelos com diferentes estruturas, incluindo:
-modelos com dose tratada como fator (mais flexíveis, porém menos parsimoniosos); -modelos com o logaritmo da dose como covariável quantitativa (mais interpretáveis); -modelos com e sem interação entre dose e inseticida.
Um aspecto central desse exemplo é a análise da estrutura das curvas dose–resposta. Quando o modelo inclui apenas efeitos aditivos (sem interação), assume-se que as curvas são paralelas na escala do logito, ou seja, que o efeito da dose é semelhante para todos os inseticidas. Por outro lado, a inclusão de termos de interação permite verificar se as inclinações diferem entre os tratamentos, indicando respostas distintas à dose.
A comparação entre modelos, realizada por meio da análise de desvio e testes qui-quadrado, possibilita identificar a estrutura mais adequada para descrever os dados. Esse processo evidencia um dos principais aspectos dos MLG: a construção do modelo é incremental, guiada tanto por critérios estatísticos quanto pela interpretação científica do problema.
Os gráficos nas escalas do logito e da proporção desempenham papel fundamental na análise, permitindo visualizar o ajuste dos modelos e interpretar o comportamento das curvas para cada inseticida. Em particular, esses gráficos facilitam a compreensão de diferenças entre tratamentos e da adequação das suposições do modelo.
Por fim, esse exemplo destaca a flexibilidade dos MLG na modelagem de dados complexos e reforça a importância da escolha adequada da estrutura do modelo, especialmente quando múltiplos fatores e possíveis interações estão envolvidos.
# ==========================================
# Exemplo 7.3 - Tribolium
# Programa R organizado para R 4.5.3
# ==========================================
# Dados
dose <- c(2.00, 2.64, 3.48, 4.59, 6.06, 8.00)
dose <- rep(dose, 3)
d <- factor(dose)
ldose <- log(dose)
y <- c(3, 5, 19, 19, 24, 35,
2, 14, 20, 27, 41, 40,
28, 37, 46, 48, 48, 50)
m <- c(50, 49, 47, 50, 49, 50,
50, 49, 50, 50, 50, 50,
50, 50, 50, 50, 50, 50)
insecticid <- factor(c(rep("DDT", 6),
rep("BHC", 6),
rep("DDT+BHC", 6)))
Tribolium.dat <- data.frame(dose, y, m, insecticid, d, ldose)
# Resposta binomial agrupada
resp <- cbind(y, m - y)
# Proporções observadas
pe <- y / m
# ==========================================
# Gráficos iniciais
# ==========================================
plot(
dose, pe,
pch = c(rep("*", 6), rep("+", 6), rep("-", 6)),
col = c(rep("green", 6), rep("red", 6), rep("blue", 6)),
xlab = "Dose",
ylab = "Proporções de insetos mortos"
)
plot(
ldose, pe,
pch = c(rep("*", 6), rep("+", 6), rep("-", 6)),
col = c(rep("green", 6), rep("red", 6), rep("blue", 6)),
xlab = "Log(Dose)",
ylab = "Proporções de insetos mortos"
)
# ==========================================
# Ajuste dos modelos
# ==========================================
# Modelo nulo
mod1 <- glm(resp ~ 1, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod1), df.residual(mod1))## [1] 0## [1] 347.1374## [1] 0# Modelo com dose como fator
mod2 <- glm(resp ~ d, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod2), df.residual(mod2))## [1] 0## [1] 218.9461## [1] 0# Modelo com inseticida
mod3 <- glm(resp ~ insecticid, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod3), df.residual(mod3))## [1] 0## [1] 215.0503## [1] 0# Modelo com inseticida e dose como fator
mod4 <- glm(resp ~ insecticid + d, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod4), df.residual(mod4))## [1] 0.2316562##
## Call:
## glm(formula = resp ~ insecticid + d, family = binomial, data = Tribolium.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4187 0.2980 -8.116 4.83e-16 ***
## insecticidDDT -0.6937 0.1959 -3.541 0.000399 ***
## insecticidDDT+BHC 2.5616 0.2580 9.928 < 2e-16 ***
## d2.64 1.1490 0.3247 3.539 0.000402 ***
## d3.48 2.3395 0.3330 7.025 2.13e-12 ***
## d4.59 2.6384 0.3351 7.873 3.46e-15 ***
## d6.06 3.4017 0.3489 9.750 < 2e-16 ***
## d8 3.9400 0.3657 10.774 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 413.648 on 17 degrees of freedom
## Residual deviance: 12.859 on 10 degrees of freedom
## AIC: 92.33
##
## Number of Fisher Scoring iterations: 4## [1] 11.7971## [1] 0.2988656## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: resp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 17 413.65
## insecticid 2 178.93 15 234.71 < 2.2e-16 ***
## d 5 221.85 10 12.86 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo com log(dose) e efeito do inseticida sem intercepto global
mod5 <- glm(resp ~ insecticid + ldose - 1, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod5), df.residual(mod5))## [1] 0.09462334##
## Call:
## glm(formula = resp ~ insecticid + ldose - 1, family = binomial,
## data = Tribolium.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## insecticidBHC -3.8425 0.3327 -11.550 < 2e-16 ***
## insecticidDDT -4.5553 0.3611 -12.613 < 2e-16 ***
## insecticidDDT+BHC -1.4248 0.2851 -4.998 5.79e-07 ***
## ldose 2.6958 0.2157 12.498 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 429.268 on 18 degrees of freedom
## Residual deviance: 21.282 on 14 degrees of freedom
## AIC: 92.753
##
## Number of Fisher Scoring iterations: 4## [1] 20.32354## [1] 0.1202665## Analysis of Deviance Table
##
## Model 1: resp ~ 1
## Model 2: resp ~ d
## Model 3: resp ~ insecticid + d
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 17 413.65
## 2 12 242.64 5 171.00 < 2.2e-16 ***
## 3 10 12.86 2 229.78 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table
##
## Model 1: resp ~ 1
## Model 2: resp ~ insecticid
## Model 3: resp ~ insecticid + d
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 17 413.65
## 2 15 234.71 2 178.93 < 2.2e-16 ***
## 3 10 12.86 5 221.85 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo com estrutura aninhada
mod6 <- glm(resp ~ ldose / insecticid, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod6), df.residual(mod6))## [1] 0.03751499##
## Call:
## glm(formula = resp ~ ldose/insecticid, family = binomial, data = Tribolium.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5619 0.3034 -11.738 < 2e-16 ***
## ldose 2.5227 0.2190 11.517 < 2e-16 ***
## ldose:insecticidDDT -0.4608 0.1257 -3.667 0.000245 ***
## ldose:insecticidDDT+BHC 2.2581 0.2449 9.222 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 413.648 on 17 degrees of freedom
## Residual deviance: 24.712 on 14 degrees of freedom
## AIC: 96.183
##
## Number of Fisher Scoring iterations: 5# Modelo apenas com log(dose)
mod7 <- glm(resp ~ ldose, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod7), df.residual(mod7))## [1] 0## [1] 219.84## [1] 0# Modelo com inclinações específicas
mod8 <- glm(resp ~ insecticid * ldose - insecticid, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod8), df.residual(mod8))## [1] 0.03751499## [1] 28.0365## [1] 0.01407002# Modelo aditivo com inseticida e log(dose)
mod9 <- glm(resp ~ insecticid + ldose, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod9), df.residual(mod9))## [1] 0.09462334## [1] 20.32354## [1] 0.1202665##
## Call:
## glm(formula = resp ~ insecticid + ldose, family = binomial, data = Tribolium.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.8425 0.3327 -11.550 < 2e-16 ***
## insecticidDDT -0.7128 0.1981 -3.598 0.00032 ***
## insecticidDDT+BHC 2.4177 0.2381 10.154 < 2e-16 ***
## ldose 2.6958 0.2157 12.498 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 413.648 on 17 degrees of freedom
## Residual deviance: 21.282 on 14 degrees of freedom
## AIC: 92.753
##
## Number of Fisher Scoring iterations: 4## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: resp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 17 413.65
## insecticid 2 178.93 15 234.71 < 2.2e-16 ***
## ldose 1 213.43 14 21.28 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo com interação completa
mod10 <- glm(resp ~ insecticid * ldose, family = binomial, data = Tribolium.dat)
1 - pchisq(deviance(mod10), df.residual(mod10))## [1] 0.1190809##
## Call:
## glm(formula = resp ~ insecticid * ldose, family = binomial, data = Tribolium.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.0428 0.4972 -8.131 4.27e-16 ***
## insecticidDDT 0.1278 0.7118 0.180 0.8575
## insecticidDDT+BHC 1.9221 0.7722 2.489 0.0128 *
## ldose 2.8381 0.3392 8.367 < 2e-16 ***
## insecticidDDT:ldose -0.5602 0.4680 -1.197 0.2313
## insecticidDDT+BHC:ldose 0.5503 0.6662 0.826 0.4088
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 413.65 on 17 degrees of freedom
## Residual deviance: 17.89 on 12 degrees of freedom
## AIC: 93.361
##
## Number of Fisher Scoring iterations: 4## [1] 17.61094## [1] 0.1280245# ==========================================
# Gráfico na escala do logito
# ==========================================
plot(
c(1.8, 8), c(-3, 3.5),
type = "n",
xlab = "dose",
ylab = "Logit(proporções)",
log = "x"
)
points(
dose, log(pe / (1 - pe)),
pch = c(rep("*", 6), rep("+", 6), rep("-", 6)),
col = c(rep("green", 6), rep("red", 6), rep("blue", 6))
)
ld <- seq(log(2), log(8), 0.005)
lines(
exp(ld),
predict(
mod5,
newdata = data.frame(
ldose = ld,
insecticid = factor(rep("DDT", length(ld)),
levels = levels(Tribolium.dat$insecticid))
),
type = "link"
),
col = "green"
)
lines(
exp(ld),
predict(
mod5,
newdata = data.frame(
ldose = ld,
insecticid = factor(rep("BHC", length(ld)),
levels = levels(Tribolium.dat$insecticid))
),
type = "link"
),
col = "red"
)
lines(
exp(ld),
predict(
mod5,
newdata = data.frame(
ldose = ld,
insecticid = factor(rep("DDT+BHC", length(ld)),
levels = levels(Tribolium.dat$insecticid))
),
type = "link"
),
col = "blue"
)
# ==========================================
# Gráfico na escala original das proporções
# ==========================================
plot(
c(1.8, 8), c(0, 1),
type = "n",
xlab = "dose",
ylab = "proporções",
log = "x"
)
points(
dose, pe,
pch = c(rep("*", 6), rep("+", 6), rep("-", 6)),
col = c(rep("green", 6), rep("red", 6), rep("blue", 6))
)
ld <- seq(log(2), log(8), 0.005)
lines(
exp(ld),
predict(
mod5,
newdata = data.frame(
ldose = ld,
insecticid = factor(rep("DDT", length(ld)),
levels = levels(Tribolium.dat$insecticid))
),
type = "response"
),
col = "green"
)
lines(
exp(ld),
predict(
mod5,
newdata = data.frame(
ldose = ld,
insecticid = factor(rep("BHC", length(ld)),
levels = levels(Tribolium.dat$insecticid))
),
type = "response"
),
col = "red"
)
lines(
exp(ld),
predict(
mod5,
newdata = data.frame(
ldose = ld,
insecticid = factor(rep("DDT+BHC", length(ld)),
levels = levels(Tribolium.dat$insecticid))
),
type = "response"
),
col = "blue"
)
4.4. Modelo de regressão de Poisson para dados de contagem
O modelo de Poisson é um caso particular importante dos Modelos Lineares Generalizados, sendo especialmente adequado para a análise de dados na forma de contagens. Esse tipo de modelo é utilizado quando a variável resposta assume valores inteiros não negativos, como números de ocorrências, eventos, indivíduos ou unidades observadas em um determinado intervalo de tempo, área ou volume.
Seja \(Y\) uma variável aleatória com distribuição de Poisson de parâmetro \(\mu > 0\). Sua função de probabilidade é dada por
\[ P(Y = y) = \frac{e^{-\mu}\mu^y}{y!}, \quad y = 0,1,2,\ldots \]
No modelo de Poisson, o valor esperado e a variância coincidem, isto é,
\[ E(Y) = \mu \]
e
\[ \mathrm{Var}(Y) = \mu. \]
Essa propriedade caracteriza a distribuição e mostra que a variabilidade cresce com a média. Além disso, a distribuição de Poisson pertence à família exponencial e, por isso, pode ser incorporada naturalmente ao contexto dos Modelos Lineares Generalizados.
A parametrização canônica do modelo é obtida escrevendo-se
\[ \theta = \log(\mu), \]
de modo que a função de ligação canônica é a função logarítmica. Assim, em um modelo de regressão de Poisson, tem-se
\[ \eta_i = \log(\mu_i), \]
em que \(\eta_i\) é o preditor linear associado à i-ésima observação. Como consequência,
\[ \mu_i = \exp(\eta_i), \]
o que garante que as médias ajustadas sejam sempre positivas, como deve ocorrer em dados de contagem.
Uma forma geral do modelo de Poisson pode ser escrita como
\[ \log(\mu_i) = x_i^T \beta, \]
em que \(x_i^T\) representa o vetor de covariáveis da i-ésima observação e \(\beta\) é o vetor de parâmetros a serem estimados. Essa formulação permite interpretar os coeficientes em termos multiplicativos sobre a média, já que incrementos lineares no preditor correspondem a efeitos exponenciais sobre \(\mu_i\).
O modelo de Poisson é particularmente útil quando as contagens dependem de uma medida de exposição conhecida, como tempo de observação, área amostrada ou volume analisado. Nessas situações, utiliza-se frequentemente uma variável offset, que entra no modelo como um termo conhecido e fixo. Por exemplo, se \(\nu_i\) representa uma unidade de exposição conhecida, pode-se escrever
\[ \mu_i = \lambda \nu_i, \]
onde \(\lambda\) é uma taxa de ocorrência. Aplicando o logaritmo, obtém-se
\[ \log(\mu_i) = \log(\lambda) + \log(\nu_i), \]
em que \(\log(\nu_i)\) entra no modelo como offset.
Esse tipo de formulação é bastante comum em aplicações biológicas. Em estudos de diluição em série, por exemplo, o número esperado de organismos em uma subamostra pode ser modelado como proporcional ao volume analisado. Quando se trabalha com presença ou ausência do organismo em subamostras independentes, também surgem modelos relacionados, com outras funções de ligação, mostrando a flexibilidade do arcabouço dos MLG. :contentReferenceoaicite:3
Do ponto de vista assintótico, quando \(\mu\) é grande, a distribuição de Poisson pode ser aproximada por uma distribuição normal com média \(\mu\) e variância \(\mu\). Essa propriedade é útil em algumas aproximações analíticas e reforça a conexão entre o modelo de Poisson e outras técnicas estatísticas clássicas. :contentReferenceoaicite:4
Em síntese, o modelo de Poisson desempenha papel central na análise estatística de contagens. Sua importância decorre não apenas de sua interpretação natural em muitos problemas práticos, mas também de sua estrutura matemática simples, que o torna um dos modelos especiais mais relevantes dentro da teoria dos Modelos Lineares Generalizados.
Exemplo: contagem de bactérias
O exemplo de contagem de bactérias ilustra a aplicação de Modelos Lineares Generalizados com distribuição de Poisson, sendo um caso típico de modelagem de dados de contagem ao longo do tempo.
Nesse contexto, a variável resposta corresponde ao número de bactérias observadas em diferentes instantes de tempo. Como se trata de contagens (valores inteiros não negativos), a suposição de normalidade não é adequada, sendo mais apropriado assumir uma distribuição de Poisson, na qual a média está relacionada à variância.
O objetivo principal do exemplo é investigar como o número de bactérias varia ao longo do tempo e verificar qual forma funcional melhor descreve essa relação. Para isso, são comparados dois modelos:
um modelo linear simples, considerando o tempo original como variável explicativa; um modelo utilizando o logaritmo do tempo, com o intuito de capturar possíveis relações não lineares.
Essa comparação evidencia um ponto fundamental na modelagem com MLG: a escolha da forma funcional das covariáveis pode ser tão importante quanto a escolha da distribuição da variável resposta. Transformações como o logaritmo frequentemente permitem linearizar relações originalmente não lineares, resultando em melhor ajuste do modelo.
Além disso, o uso da função de ligação logarítmica (característica da distribuição de Poisson) garante que as predições da média sejam sempre positivas, o que é coerente com a natureza dos dados de contagem.
Por fim, o exemplo destaca a importância da análise de ajuste, realizada por meio da comparação de desvios (via teste qui-quadrado), permitindo selecionar o modelo mais adequado para descrever o fenômeno observado.
# ==========================================
# Exemplo - Contagem de bactérias
# ==========================================
# Dados
tim <- c(0, 1, 2, 6, 12)
count <- c(31, 26, 19, 15, 20)
# Transformações
ltime <- log(tim + 0.1)
lcount <- log(count)
# Data frame
bacteria.dat <- data.frame(tim, count, ltime, lcount)
# Visualização inicial
par(mfrow = c(1, 2))
plot(tim, count, xlab = "Time in months", ylab = "Counts")
plot(ltime, lcount, xlab = "Log(time in months)", ylab = "Log(counts)")
par(mfrow = c(1, 1))
# ==========================================
# Modelo 1: Poisson com tempo original
# ==========================================
mod1 <- glm(count ~ tim, family = poisson)
anova(mod1, test = "Chisq")## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 4 7.0672
## tim 1 2.5423 3 4.5249 0.1108# ==========================================
# Modelo 2: Poisson com log do tempo
# ==========================================
mod2 <- glm(count ~ ltime, family = poisson)
anova(mod2, test = "Chisq")## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 4 7.0672
## ltime 1 5.2334 3 1.8338 0.02216 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ==========================================
# Gráfico com ajuste do modelo
# ==========================================
plot(c(0, 12), c(15, 31), type = "n",
xlab = "Time in months", ylab = "Counts")
points(tim, count, pch = "*")
# Sequência para predição
x <- seq(0, 12, 0.1)
lp <- predict(mod2, newdata = data.frame(ltime = log(x + 0.1)))
# Função inversa do log (Poisson → link log)
yhat <- exp(lp)
# Curva ajustada
lines(x, yhat, lty = 1, col = "blue")
Exemplo: Armadilhas adesivas
O exemplo das armadilhas adesivas ilustra a aplicação de Modelos Lineares Generalizados com distribuição de Poisson na análise de dados de contagem organizados em tabelas de contingência. Trata-se de um caso clássico de modelagem log-linear, no qual o interesse está em investigar a associação entre duas variáveis categóricas.
Nesse contexto, a variável resposta corresponde ao número de insetos capturados, classificados de acordo com duas características: a cor da armadilha e o sexo dos insetos. Assim, os dados são organizados em uma tabela 2×2, na qual cada célula contém uma contagem observada.
O principal objetivo da análise é verificar se existe associação entre as variáveis consideradas — isto é, se a captura de insetos depende simultaneamente da cor da armadilha e do sexo, ou se essas variáveis são independentes.
Para isso, são ajustados dois modelos principais: -um modelo saturado, que inclui o termo de interação entre as variáveis e, portanto, reproduz exatamente os dados observados; -um modelo de independência, que considera apenas os efeitos principais, assumindo que não há interação entre as variáveis.
A comparação entre esses modelos é realizada por meio da análise de desvio, permitindo avaliar se a inclusão do termo de interação melhora significativamente o ajuste. Caso o modelo de independência apresente ajuste inadequado, isso indica evidência de associação entre as variáveis.
Um aspecto importante desse exemplo é a interpretação do coeficiente de interação no modelo saturado. Esse coeficiente está diretamente relacionado ao logaritmo da razão de chances (odds ratio), uma medida amplamente utilizada para quantificar a associação entre variáveis categóricas em tabelas de contingência.
Além disso, o uso da distribuição de Poisson é justificado pela natureza dos dados (contagens), e a função de ligação logarítmica garante a coerência das estimativas, assegurando valores positivos para as médias.
Por fim, esse exemplo evidencia a conexão entre os Modelos Lineares Generalizados e métodos clássicos de análise de tabelas de contingência, mostrando como testes de independência podem ser interpretados dentro do framework dos MLG. Essa abordagem unificada amplia a capacidade analítica e fornece uma base mais geral para a modelagem estatística de dados categóricos.
# ==========================================
# Exemplo - Armadilhas adesivas
# ==========================================
# Dados: tabela 2x2
y <- c(246, 17, 458, 32)
armcor <- factor(c(1, 1, 2, 2))
sexo <- factor(c(1, 2, 1, 2))
count.dat <- data.frame(armcor, sexo, y)
# Visualizar os dados
count.dat## armcor sexo y
## 1 1 1 246
## 2 1 2 17
## 3 2 1 458
## 4 2 2 32# ==========================================
# Razão de chances observada
# ==========================================
odds_ratio_obs <- 246 * 32 / (17 * 458)
odds_ratio_obs## [1] 1.011045# ==========================================
# Modelo log-linear saturado
# ==========================================
mod1 <- glm(y ~ armcor * sexo, family = poisson, data = count.dat)
# Estatística de Pearson residual
sum(residuals(mod1, type = "pearson")^2)## [1] 1.797633e-27## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 3 750.89
## armcor 1 69.51 2 681.38 <2e-16 ***
## sexo 1 681.38 1 0.00 <2e-16 ***
## armcor:sexo 1 0.00 0 0.00 0.9718
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## glm(formula = y ~ armcor * sexo, family = poisson, data = count.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.50533 0.06376 86.348 < 2e-16 ***
## armcor2 0.62154 0.07905 7.863 3.75e-15 ***
## sexo2 -2.67212 0.25078 -10.655 < 2e-16 ***
## armcor2:sexo2 0.01098 0.31036 0.035 0.972
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 7.5089e+02 on 3 degrees of freedom
## Residual deviance: -2.0872e-14 on 0 degrees of freedom
## AIC: 33.299
##
## Number of Fisher Scoring iterations: 3# O modelo saturado reproduz os dados
# O coeficiente de interação corresponde ao log da razão de chances ajustada
exp(mod1$coef[4])## armcor2:sexo2
## 1.011045# ==========================================
# Modelo de independência (efeitos principais)
# ==========================================
mod2 <- glm(y ~ armcor + sexo, family = poisson, data = count.dat)
# Estatística de Pearson residual
sum(residuals(mod2, type = "pearson")^2)## [1] 0.00125277## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 3 750.89
## armcor 1 69.51 2 681.38 < 2.2e-16 ***
## sexo 1 681.38 1 0.00 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] 0.9717516##
## Call:
## glm(formula = y ~ armcor + sexo, family = poisson, data = count.dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.50487 0.06241 88.21 < 2e-16 ***
## armcor2 0.62225 0.07644 8.14 3.94e-16 ***
## sexo2 -2.66496 0.14775 -18.04 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 750.888715 on 3 degrees of freedom
## Residual deviance: 0.001254 on 1 degrees of freedom
## AIC: 31.3
##
## Number of Fisher Scoring iterations: 3# ==========================================
# Comparação entre os modelos
# ==========================================
anova(mod2, mod1, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ armcor + sexo
## Model 2: y ~ armcor * sexo
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 1 0.001254
## 2 0 0.000000 1 0.001254 0.97184.5. Estimação dos parâmetros e interpretação dos coeficientes
Função de Verossimilhança
Nos Modelos Lineares Generalizados (MLG), a estimação dos parâmetros é realizada por meio do método da máxima verossimilhança. Considera-se um conjunto de variáveis aleatórias independentes \(Y_1, Y_2, \ldots, Y_n\), cada uma com distribuição pertencente à família exponencial.
A função de verossimilhança é obtida a partir da distribuição conjunta das observações e pode ser escrita como o produto das funções de densidade individuais. Em termos logarítmicos, a log-verossimilhança assume a forma:
\[ \ell(\theta, \phi) = \sum_{i=1}^{n} \left[ \frac{y_i \theta_i - b(\theta_i)}{\phi} + c(y_i, \phi) \right] \]
onde: - \(\theta_i\) é o parâmetro canônico; - \(\phi\) é o parâmetro de dispersão; - \(b(\cdot)\) e \(c(\cdot)\) são funções conhecidas.
A média da variável resposta é dada por:
\[ \mu_i = E(Y_i) = b'(\theta_i) \]
e a variância é dada por:
\[ \mathrm{Var}(Y_i) = \phi \, b''(\theta_i) \]
O objetivo é estimar os parâmetros \(\beta\) que aparecem no preditor linear, por meio da maximização da função de verossimilhança.
Função Escore e Informação
Para obter os estimadores de máxima verossimilhança, derivamos a log-verossimilhança em relação aos parâmetros \(\beta\), obtendo a chamada função escore.
A função escore pode ser escrita como:
\[ U(\beta) = \frac{\partial \ell}{\partial \beta} \]
Nos MLG, essa expressão pode ser reescrita em termos da média \(\mu_i\) e do preditor linear \(\eta_i\), resultando em:
\[ U(\beta) = X^T W (y - \mu) \]
onde: - \(X\) é a matriz de delineamento; - \(y\) é o vetor de observações; - \(\mu\) é o vetor de médias; - \(W\) é uma matriz diagonal de pesos.
A matriz de pesos depende da variância da resposta e da função de ligação, sendo dada por:
\[ W_i = \frac{1}{\mathrm{Var}(Y_i)} \left( \frac{\partial \mu_i}{\partial \eta_i} \right)^2 \]
A matriz de informação de Fisher é definida como:
\[ I(\beta) = -E\left[ \frac{\partial^2 \ell}{\partial \beta \partial \beta^T} \right] \]
e, nos GLM, assume a forma:
\[ I(\beta) = X^T W X \]
Essa matriz desempenha papel fundamental na estimação e na obtenção de propriedades assintóticas dos estimadores.
Método de Estimação: Algoritmo de Fisher Scoring
Como as equações de verossimilhança não possuem solução analítica em geral, utiliza-se um método iterativo para estimar os parâmetros. Um dos métodos mais utilizados é o Fisher Scoring, que é uma versão do método de Newton-Raphson.
A atualização dos parâmetros é dada por:
\[ \beta^{(t+1)} = \beta^{(t)} + I(\beta^{(t)})^{-1} U(\beta^{(t)}) \]
Substituindo as expressões da função escore e da matriz de informação, obtém-se:
\[ \beta^{(t+1)} = \beta^{(t)} + (X^T W X)^{-1} X^T W (y - \mu) \]
Esse procedimento pode ser interpretado como um problema de mínimos quadrados ponderados, no qual se define uma variável ajustada:
\[ z = \eta + \frac{y - \mu}{\frac{\partial \mu}{\partial \eta}} \]
Assim, em cada iteração, ajusta-se um modelo de regressão linear ponderada da forma:
\[ z = X\beta \]
com pesos dados por:
\[ W_i = \frac{1}{\mathrm{Var}(Y_i)} \left( \frac{\partial \mu_i}{\partial \eta_i} \right)^2 \]
Esse algoritmo é conhecido como Iteratively Reweighted Least Squares (IRLS).
O processo iterativo continua até que haja convergência, isto é, até que as estimativas dos parâmetros não sofram variações significativas entre iterações sucessivas.
Método de Iteração: Mínimos Quadrados Ponderados Iterativos (IRLS)
Para resolver as equações de verossimilhança nos MLG, utiliza-se um procedimento iterativo conhecido como Iteratively Reweighted Least Squares (IRLS), ou Mínimos Quadrados Ponderados Iterativos.
A ideia central desse método é aproximar o problema de máxima verossimilhança por uma sequência de problemas de mínimos quadrados ponderados.
Em cada iteração, define-se uma variável auxiliar \(z_i\), chamada de variável dependente ajustada, dada por:
\[ z_i = \eta_i + \frac{y_i - \mu_i}{\frac{\partial \mu_i}{\partial \eta_i}} \]
Além disso, são definidos pesos \(w_i\), que dependem da variância da resposta e da função de ligação:
\[ w_i = \left( \frac{\partial \mu_i}{\partial \eta_i} \right)^2 \Big/ \mathrm{Var}(Y_i) \]
Com essas definições, o problema de estimação passa a ser resolver, em cada passo, um modelo de regressão linear ponderada da forma:
\[ z = X\beta \]
em que cada observação é ponderada por \(w_i\).
O algoritmo segue os seguintes passos:
- Inicializar valores para \(\beta\);
- Calcular \(\eta_i\) e \(\mu_i\);
- Construir a variável ajustada \(z_i\);
- Calcular os pesos \(w_i\);
- Ajustar um modelo de mínimos quadrados ponderados;
- Atualizar \(\beta\);
- Repetir até convergência.
Esse procedimento converge, sob condições gerais, para o estimador de máxima verossimilhança.
Comentário
O método IRLS é uma das ideias mais importantes dos Modelos Lineares Generalizados, pois conecta dois mundos:
- a máxima verossimilhança (mais geral e teórica);
- os mínimos quadrados (mais intuitivos e computacionalmente simples).
Na prática, funções como glm() no R utilizam
internamente esse algoritmo, permitindo ao usuário ajustar modelos
complexos de forma transparente.
Assim, embora o processo matemático envolva iterações, o uso computacional torna a aplicação dos MLG bastante acessível.
4.6. Resíduos, medidas de ajuste e diagnóstico do modelo
Introdução
Após o ajuste de um Modelo Linear Generalizado (GLM), é fundamental avaliar a qualidade do ajuste obtido. Essa avaliação permite verificar se o modelo especificado é adequado para descrever os dados observados e se as suposições adotadas são razoáveis.
Diferentemente do modelo linear clássico, em que os resíduos apresentam propriedades bem conhecidas (como normalidade e variância constante), nos MLG a análise de resíduos é mais complexa, devido à dependência entre média e variância e à possibilidade de distribuições não normais para a variável resposta.
Nesse contexto, os resíduos desempenham um papel central no diagnóstico do modelo, sendo utilizados para:
- identificar observações discrepantes;
- avaliar a adequação da função de ligação;
- verificar a presença de padrões não explicados pelo modelo;
- detectar possíveis problemas de especificação.
A análise de resíduos em MLG busca, portanto, construir medidas que permitam comparar valores observados e ajustados de forma apropriada, levando em consideração a estrutura probabilística do modelo.
Tipos de Resíduos
Nos Modelos Lineares Generalizados, diferentes tipos de resíduos podem ser definidos, cada um com propriedades específicas e utilidades distintas na análise diagnóstica.
Resíduos brutos
Os resíduos mais simples são definidos como a diferença entre o valor observado e o valor ajustado:
\[ e_i = y_i - \mu_i \]
onde: - \(y_i\) é o valor observado; - \(\mu_i\) é o valor ajustado pelo modelo.
Entretanto, esses resíduos não são adequados para análise direta, pois não levam em conta a variabilidade da resposta, que depende da média.
Resíduos padronizados pela variância
Para contornar essa limitação, consideram-se resíduos que incorporam a variância da resposta. Um exemplo importante é o resíduo de Pearson:
\[ r_i^P = \frac{y_i - \mu_i}{\sqrt{\mathrm{Var}(Y_i)}} \]
Esse resíduo ajusta a diferença entre observado e esperado pela variabilidade do modelo, tornando os valores mais comparáveis entre si.
A soma dos quadrados desses resíduos é conhecida como estatística de Pearson:
\[ X^2 = \sum_{i=1}^{n} (r_i^P)^2 \]
Essa quantidade é frequentemente utilizada para avaliar a qualidade do ajuste do modelo.
Resíduos de desvio
Outro tipo importante são os resíduos de desvio, baseados na contribuição individual de cada observação para o desvio do modelo.
O desvio é uma medida de discrepância entre o modelo ajustado e o modelo saturado (que ajusta perfeitamente os dados). Os resíduos de desvio são definidos de forma a refletir essa contribuição individual, sendo dados por:
\[ r_i^D = \text{sign}(y_i - \mu_i) \sqrt{d_i} \]
onde \(d_i\) representa a contribuição da i-ésima observação para o desvio total.
Esses resíduos apresentam propriedades mais próximas da normalidade, sendo amplamente utilizados em análises gráficas e diagnósticas.
Comentário
A existência de diferentes tipos de resíduos em MLG reflete a complexidade desses modelos em comparação com a regressão linear clássica. Não existe um único resíduo “ideal”, sendo necessário escolher o tipo mais apropriado de acordo com o objetivo da análise.
Em geral:
- resíduos de Pearson são úteis para avaliar a variabilidade global;
- resíduos de desvio são preferidos para análise gráfica e identificação de pontos discrepantes.
Na prática, softwares como o R permitem extrair facilmente esses
resíduos por meio da função residuals() com diferentes
opções, facilitando a análise diagnóstica dos modelos ajustados.
Análise de resíduos e diagnóstico para o modelo clássico de regressão: Resíduos Padronizados e Influência das Observações
Nos Modelos Lineares Generalizados, a análise de resíduos simples pode não ser suficiente para identificar observações influentes ou avaliar adequadamente o ajuste do modelo. Isso ocorre porque diferentes observações podem ter diferentes níveis de influência na estimação dos parâmetros.
Para lidar com esse problema, introduzem-se conceitos adicionais, como padronização dos resíduos e medidas de alavancagem (leverage).
Matriz de projeção e alavancagem (leverage)
No contexto dos MLG, a matriz de projeção (ou matriz “hat”) generalizada é dada por:
\[ H = W^{1/2} X (X^T W X)^{-1} X^T W^{1/2} \]
onde: - \(X\) é a matriz de delineamento; - \(W\) é a matriz de pesos do modelo.
Os elementos diagonais dessa matriz, denotados por \(h_{ii}\), são chamados de alavancagens (leverages).
A alavancagem mede o grau de influência de uma observação na estimação dos valores ajustados. Observações com valores elevados de \(h_{ii}\) têm maior potencial de influenciar o modelo.
Resíduos padronizados
Para levar em conta a variabilidade e a influência de cada observação, define-se o resíduo de Pearson padronizado como:
\[ r_i^{*} = \frac{y_i - \mu_i}{\sqrt{\mathrm{Var}(Y_i)(1 - h_{ii})}} \]
Esse resíduo ajusta o resíduo de Pearson pela alavancagem, tornando-o mais adequado para identificar observações discrepantes.
Valores elevados de \(|r_i^{*}|\) indicam possíveis outliers.
Resíduos de desvio padronizados
De forma análoga, os resíduos de desvio também podem ser padronizados:
\[ r_i^{(D)*} = \frac{r_i^{(D)}}{\sqrt{1 - h_{ii}}} \]
Esses resíduos apresentam boas propriedades assintóticas e são frequentemente utilizados em análises gráficas.
Resíduos studentizados
Os resíduos padronizados ainda dependem de estimativas da variância obtidas com base em todos os dados. Para melhorar a identificação de observações influentes, utilizam-se os resíduos studentizados, que consideram a variabilidade estimada sem a observação em análise.
Embora a formulação exata seja mais complexa nos GLM, a ideia central é avaliar o impacto da exclusão de cada observação na estimação do modelo.
Influência das observações
Além dos resíduos, é importante avaliar o impacto de cada observação sobre os parâmetros estimados. Observações influentes são aquelas cuja remoção provoca mudanças significativas nas estimativas.
Uma medida amplamente utilizada para avaliar essa influência é a distância de Cook, definida, em termos gerais, como uma função da variação nos parâmetros estimados ao se remover uma observação.
Valores elevados da distância de Cook indicam observações potencialmente influentes.
Interpretação conjunta
A análise de resíduos e influência deve ser feita de forma integrada:
- resíduos grandes indicam discrepâncias entre valores observados e ajustados;
- alavancagem alta indica potencial de influência estrutural;
- combinação de alta alavancagem e grande resíduo indica observações altamente influentes.
Comentário
Nem todas as observações têm o mesmo peso na construção do modelo.
Enquanto alguns pontos apenas apresentam desvios naturais, outros podem exercer forte influência sobre o ajuste, distorcendo as estimativas. Por isso, a análise de resíduos deve sempre considerar:
- o tamanho do resíduo;
- a alavancagem da observação;
- o impacto global na estimação.
Na prática, essas medidas são utilizadas conjuntamente em gráficos de diagnóstico, permitindo identificar padrões, outliers e possíveis problemas de especificação do modelo.
Funções como residuals(), hatvalues() e
cooks.distance() no R facilitam a obtenção dessas medidas,
tornando a análise de diagnóstico uma etapa acessível e essencial no uso
de Modelos Lineares Generalizados.
Análise de resíduos e diagnóstico para modelos lineares generalizados
A análise gráfica desempenha um papel central na avaliação da adequação dos Modelos Lineares Generalizados (MLG). Embora medidas numéricas, como resíduos e estatísticas de ajuste, forneçam informações importantes, os gráficos permitem identificar padrões, tendências e anomalias que dificilmente seriam percebidos apenas por meio de valores numéricos.
Nesta seção, são apresentados os principais gráficos utilizados na análise de diagnóstico em MLG, com ênfase na interpretação dos resíduos e na identificação de observações influentes.
Gráfico de resíduos versus valores ajustados
Um dos gráficos mais utilizados é o gráfico dos resíduos em função dos valores ajustados:
\[ r_i \text{ vs } \hat{\mu}_i \]
Esse gráfico permite avaliar:
- a adequação da função de ligação;
- a presença de padrões sistemáticos não capturados pelo modelo;
- a existência de heterogeneidade não explicada.
Em um modelo bem ajustado, espera-se que os resíduos estejam distribuídos de forma aleatória em torno de zero, sem padrões estruturais.
A presença de tendências (como curvas ou funis) pode indicar problemas de especificação do modelo.
Gráfico de resíduos versus covariáveis
Outra análise importante consiste em representar os resíduos em função das variáveis explicativas:
\[ r_i \text{ vs } x_{ij} \]
Esse gráfico permite identificar:
- relações não lineares entre a resposta e as covariáveis;
- necessidade de transformações nas variáveis explicativas;
- possíveis variáveis omitidas.
Se padrões sistemáticos forem observados, isso sugere que o modelo não está capturando adequadamente a relação entre as variáveis.
Gráfico de alavancagem (leverge)
A análise da alavancagem pode ser feita por meio de gráficos dos valores \(h_{ii}\), que representam a influência estrutural de cada observação.
Observações com valores elevados de alavancagem:
- estão associadas a combinações incomuns de covariáveis;
- têm maior potencial de influenciar o ajuste do modelo.
No entanto, alta alavancagem por si só não indica problema — é necessário considerar também os resíduos.
Gráfico de resíduos padronizados versus alavancagem (leverage)
Uma ferramenta poderosa de diagnóstico é o gráfico que combina resíduos padronizados e alavancagem:
\[ r_i^* \text{ vs } h_{ii} \]
Esse gráfico permite identificar:
- observações com alto resíduo e alta alavancagem;
- pontos potencialmente influentes.
Frequentemente, adicionam-se curvas de nível associadas à distância de Cook, facilitando a identificação de observações que exercem grande influência no modelo.
Gráfico de probabilidade normal
Embora os MLG não assumam normalidade da variável resposta, os resíduos padronizados podem ser analisados por meio de gráficos de probabilidade normal.
Nesse tipo de gráfico, avalia-se se os resíduos seguem aproximadamente uma distribuição simétrica em torno de zero.
Desvios significativos da linha de referência podem indicar:
- presença de outliers;
- inadequação do modelo;
- distribuição incorreta para a variável resposta.
Gráfico half-normal
Uma alternativa ao gráfico de probabilidade normal é o gráfico half-normal, que é particularmente útil na análise de resíduos em MLG.
Nesse gráfico, os valores absolutos dos resíduos ordenados são comparados com quantis de uma distribuição half-normal.
A principal vantagem desse gráfico é facilitar a identificação de observações discrepantes, especialmente quando acompanhado de envelopes simulados.
Observações que se afastam significativamente da faixa esperada indicam possíveis problemas no ajuste do modelo.
Interpretação dos gráficos
A análise de diagnóstico deve ser realizada de forma conjunta, considerando diferentes gráficos e medidas. Em particular:
- resíduos aleatórios indicam bom ajuste;
- padrões sistemáticos sugerem inadequação do modelo;
- pontos extremos podem indicar outliers ou observações influentes;
- alta alavancagem combinada com grandes resíduos indica forte influência.
# ============================================================
# RESÍDUOS E DIAGNÓSTICOS EM GLM
# Modelo: Regressão de Poisson
# ============================================================
rm(list = ls())
cat("\014") # limpa console no RStudioset.seed(123)
# ============================================================
# 1. GERANDO DADOS SIMULADOS
# ============================================================
n <- 80
x1 <- runif(n, 0, 10)
x2 <- rnorm(n, mean = 5, sd = 1.5)
# Verdadeiro preditor linear
eta <- 1.0 + 0.18 * x1 - 0.12 * x2
# Média verdadeira
mu <- exp(eta)
# Resposta Poisson
y <- rpois(n, lambda = mu)
# Inserindo alguns pontos problemáticos para fins didáticos
y[c(10, 35)] <- y[c(10, 35)] + c(15, 20) # possíveis outliers
x1[70] <- 15 # ponto de alta alavancagem
y[70] <- y[70] + 10
dados <- data.frame(id = 1:n, y = y, x1 = x1, x2 = x2)
# Visualização rápida dos dados
head(dados)## id y x1 x2
## 1 1 3 2.875775 3.957940
## 2 2 7 7.883051 4.688124
## 3 3 2 4.089769 3.101905
## 4 4 6 8.830174 8.253434
## 5 5 5 9.404673 6.811943
## 6 6 3 0.455565 3.315337## id y x1 x2
## Min. : 1.00 Min. : 0.000 Min. : 0.006248 Min. :1.536
## 1st Qu.:20.75 1st Qu.: 2.000 1st Qu.: 2.610014 1st Qu.:4.089
## Median :40.50 Median : 3.000 Median : 4.765563 Median :4.925
## Mean :40.50 Mean : 4.675 Mean : 5.146869 Mean :5.001
## 3rd Qu.:60.25 3rd Qu.: 6.000 3rd Qu.: 7.659209 3rd Qu.:5.790
## Max. :80.00 Max. :21.000 Max. :15.000000 Max. :8.281# ============================================================
# 2. AJUSTE DO MODELO
# ============================================================
modelo <- glm(y ~ x1 + x2, family = poisson(link = "log"), data = dados)
summary(modelo)##
## Call:
## glm(formula = y ~ x1 + x2, family = poisson(link = "log"), data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.35946 0.21037 6.462 1.03e-10 ***
## x1 0.13282 0.01643 8.083 6.34e-16 ***
## x2 -0.11930 0.03765 -3.169 0.00153 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 258.46 on 79 degrees of freedom
## Residual deviance: 186.03 on 77 degrees of freedom
## AIC: 436.34
##
## Number of Fisher Scoring iterations: 5## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 79 258.46
## x1 1 62.368 78 196.09 2.85e-15 ***
## x2 1 10.062 77 186.03 0.001513 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ============================================================
# 3. VALORES AJUSTADOS E RESÍDUOS
# ============================================================
ajustados <- fitted(modelo)
eta_hat <- predict(modelo, type = "link")
res_bruto <- residuals(modelo, type = "response")
res_pearson <- residuals(modelo, type = "pearson")
res_desvio <- residuals(modelo, type = "deviance")
res_working <- residuals(modelo, type = "working")
# Resíduos padronizados / studentizados
res_std_pearson <- rstandard(modelo, type = "pearson")
res_std_dev <- rstandard(modelo, type = "deviance")
res_stud <- rstudent(modelo)
# Medidas de influência
lev <- hatvalues(modelo)
cook <- cooks.distance(modelo)
dff <- dffits(modelo)
# Tabela-resumo
diagnostico <- data.frame(
id = dados$id,
y = dados$y,
ajustados = ajustados,
res_bruto = res_bruto,
res_pearson = res_pearson,
res_desvio = res_desvio,
res_std_pearson = res_std_pearson,
res_std_dev = res_std_dev,
res_stud = res_stud,
leverage = lev,
cook = cook,
dffits = dff
)
head(diagnostico, 10)## id y ajustados res_bruto res_pearson res_desvio res_std_pearson
## 1 1 3 3.558094 -0.5580940 -0.2958684 -0.3041573 -0.2995418
## 2 2 7 6.341859 0.6581409 0.2613428 0.2570075 0.2641076
## 3 3 2 4.630170 -2.6301705 -1.2223214 -1.3793313 -1.2454296
## 4 4 6 4.700248 1.2997524 0.5995153 0.5746535 0.6307769
## 5 5 5 6.024857 -1.0248568 -0.4175320 -0.4302971 -0.4312595
## 6 6 3 2.785535 0.2144645 0.1284994 0.1269012 0.1312256
## 7 7 3 4.648146 -1.6481456 -0.7644617 -0.8180189 -0.7702484
## 8 8 14 7.627477 6.3725229 2.3073897 2.0638216 2.3509671
## 9 9 7 3.879778 3.1202222 1.5840982 1.4217620 1.6017317
## 10 10 19 3.992176 15.0078242 7.5112619 5.4100138 7.5660730
## res_std_dev res_stud leverage cook dffits
## 1 -0.3079336 -0.3077317 0.02437599 0.0007472628 -0.03111889
## 2 0.2597264 0.2598184 0.02082677 0.0004945415 0.02421544
## 3 -1.4054078 -1.3998501 0.03676456 0.0197339502 -0.17643515
## 4 0.6046187 0.6071965 0.09666479 0.0141922009 0.12654202
## 5 -0.4444442 -0.4436297 0.06264894 0.0041435019 -0.07348016
## 6 0.1295936 0.1296611 0.04111888 0.0002461459 0.01715370
## 7 -0.8242110 -0.8234293 0.01496909 0.0030052865 -0.06506118
## 8 2.1027990 2.1124296 0.03672829 0.0702463425 0.26562186
## 9 1.4375885 1.4413829 0.02189690 0.0191450370 0.13825339
## 10 5.4494917 5.4858596 0.01443618 0.2795035687 0.45985664# ============================================================
# 4. IDENTIFICAÇÃO DE PONTOS SUSPEITOS
# ============================================================
# Regras práticas
limite_res <- 2
limite_lev <- 2 * length(coef(modelo)) / n
limite_cook <- 4 / n
limite_dff <- 2 * sqrt(length(coef(modelo)) / n)
cat("\nLimites práticos:\n")##
## Limites práticos:## |resíduo studentizado| > 2## leverage > 0.075## Cook > 0.05## |DFFITS| > 0.3873suspeitos <- subset(
diagnostico,
abs(res_stud) > limite_res |
leverage > limite_lev |
cook > limite_cook |
abs(dffits) > limite_dff
)
suspeitos[order(-abs(suspeitos$res_stud)), ]## id y ajustados res_bruto res_pearson res_desvio res_std_pearson
## 35 35 21 2.506183 18.4938167 11.6820788 7.2314830 11.8859322
## 10 10 19 3.992176 15.0078242 7.5112619 5.4100138 7.5660730
## 44 44 0 3.119051 -3.1190513 -1.7660836 -2.4976194 -1.7852313
## 65 65 2 7.502516 -5.5025159 -2.0088978 -2.3909552 -2.0538263
## 59 59 2 7.344230 -5.3442295 -1.9720230 -2.3420906 -2.0052034
## 46 46 0 2.430103 -2.4301031 -1.5588788 -2.2045875 -1.5790418
## 31 31 3 8.413994 -5.4139936 -1.8664513 -2.1541318 -1.9153703
## 8 8 14 7.627477 6.3725229 2.3073897 2.0638216 2.3509671
## 28 28 1 4.676705 -3.6767050 -1.7001561 -2.0659677 -1.7116223
## 34 34 13 7.003372 5.9966276 2.2659664 2.0221872 2.3022420
## 24 24 16 9.638858 6.3611419 2.0489072 1.8694573 2.1424690
## 32 32 6 10.746777 -4.7467773 -1.4479706 -1.5809482 -1.6365869
## 4 4 6 4.700248 1.2997524 0.5995153 0.5746535 0.6307769
## 70 70 14 13.340037 0.6599631 0.1806929 0.1792329 0.2216403
## 68 68 8 8.503167 -0.5031669 -0.1725527 -0.1742978 -0.1814501
## res_std_dev res_stud leverage cook dffits
## 35 7.3576731 7.5563725 0.03400747 1.657852359 1.04793782
## 10 5.4494917 5.4858596 0.01443618 0.279503569 0.45985664
## 44 -2.5246983 -2.5111953 0.02133615 0.023160625 -0.24245879
## 65 -2.4444284 -2.4288267 0.04327248 0.063595991 -0.33774947
## 59 -2.3814975 -2.3700957 0.03282042 0.045481269 -0.28477918
## 46 -2.2331024 -2.2188907 0.02537527 0.021639104 -0.23345973
## 31 -2.2105908 -2.1966535 0.05042818 0.064942602 -0.32997258
## 8 2.1027990 2.1124296 0.03672829 0.070246342 0.26562186
## 28 -2.0799009 -2.0754138 0.01335314 0.013216498 -0.15648733
## 34 2.0545602 2.0627544 0.03126499 0.057020896 0.23864152
## 24 1.9548247 1.9715536 0.08543312 0.142928433 0.38586606
## 32 -1.7868865 -1.7553333 0.21721709 0.247747646 -0.60687479
## 4 0.6046187 0.6071965 0.09666479 0.014192201 0.12654202
## 70 0.2198495 0.2204517 0.33536268 0.008262401 0.09983050
## 68 -0.1832852 -0.1831105 0.09566564 0.001160968 -0.03810628# ============================================================
# 5. MEDIDAS GLOBAIS DE AJUSTE
# ============================================================
cat("\n--- Medidas globais ---\n")##
## --- Medidas globais ---## Desvio residual: 186.0291## GL residual: 77## p-valor (desvio): 5.211109e-11## Qui-quadrado de Pearson: 281.833## p-valor (Pearson): 0# ============================================================
# 6. GRÁFICOS DE DIAGNÓSTICO
# ============================================================
par(mfrow = c(3, 3), mar = c(4, 4, 2, 1))
# (1) Observado vs Ajustado
plot(ajustados, dados$y,
xlab = "Valores ajustados", ylab = "Valores observados",
main = "Observado vs Ajustado", pch = 19)
abline(0, 1, col = 2, lty = 2)
# (2) Resíduo de desvio vs ajustado
plot(ajustados, res_desvio,
xlab = "Valores ajustados", ylab = "Resíduo de desvio",
main = "Resíduo de desvio vs Ajustado", pch = 19)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = c(2, 1, 2))
# (3) Resíduo de Pearson vs ajustado
plot(ajustados, res_pearson,
xlab = "Valores ajustados", ylab = "Resíduo de Pearson",
main = "Pearson vs Ajustado", pch = 19)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = c(2, 1, 2))
# (4) Resíduo studentizado vs índice
plot(res_stud,
xlab = "Índice", ylab = "Resíduo studentizado",
main = "Resíduo studentizado", pch = 19)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = c(2, 1, 2))
text(which(abs(res_stud) > 2), res_stud[abs(res_stud) > 2],
labels = which(abs(res_stud) > 2), pos = 3, cex = 0.8)
# (5) Resíduos vs x1
plot(dados$x1, res_desvio,
xlab = "x1", ylab = "Resíduo de desvio",
main = "Resíduo de desvio vs x1", pch = 19)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = c(2, 1, 2))
# (6) Resíduos vs x2
plot(dados$x2, res_desvio,
xlab = "x2", ylab = "Resíduo de desvio",
main = "Resíduo de desvio vs x2", pch = 19)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = c(2, 1, 2))
# (7) Alavancagem
plot(lev,
xlab = "Índice", ylab = "Leverage",
main = "Alavancagem", pch = 19)
abline(h = limite_lev, col = 2, lty = 2)
text(which(lev > limite_lev), lev[lev > limite_lev],
labels = which(lev > limite_lev), pos = 3, cex = 0.8)
# (8) Distância de Cook
plot(cook,
xlab = "Índice", ylab = "Cook's distance",
main = "Distância de Cook", pch = 19)
abline(h = limite_cook, col = 2, lty = 2)
text(which(cook > limite_cook), cook[cook > limite_cook],
labels = which(cook > limite_cook), pos = 3, cex = 0.8)
# (9) DFFITS
plot(dff,
xlab = "Índice", ylab = "DFFITS",
main = "DFFITS", pch = 19)
abline(h = c(-limite_dff, limite_dff), col = 2, lty = 2)
text(which(abs(dff) > limite_dff), dff[abs(dff) > limite_dff],
labels = which(abs(dff) > limite_dff), pos = 3, cex = 0.8)
# ============================================================
# 7. QQ-PLOT DOS RESÍDUOS DE DESVIO
# ============================================================
par(mfrow = c(1, 1))
qqnorm(res_desvio, main = "QQ-plot dos resíduos de desvio", pch = 19)
qqline(res_desvio, col = 2, lty = 2)
# ============================================================
# 8. RESÍDUOS PADRONIZADOS VS ALAVANCAGEM
# ============================================================
par(mfrow = c(1, 1))
plot(lev, res_std_dev,
xlab = "Alavancagem",
ylab = "Resíduo de desvio padronizado",
main = "Resíduo padronizado vs Alavancagem",
pch = 19)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = c(2, 1, 2))
abline(v = limite_lev, lty = 2, col = 4)
# Curvas aproximadas de Cook
p <- length(coef(modelo))
h_seq <- seq(0.001, max(lev) * 1.1, length.out = 200)
for (D in c(0.5, 1)) {
curva <- sqrt(D * p * (1 - h_seq) / h_seq)
lines(h_seq, curva, lty = 3, col = "darkgreen")
lines(h_seq, -curva, lty = 3, col = "darkgreen")
}
text(lev, res_std_dev, labels = ifelse(cook > limite_cook, dados$id, ""),
pos = 3, cex = 0.8)
# ============================================================
# 9. HALF-NORMAL PLOT COM ENVELOPE SIMULADO
# ============================================================
# Função auxiliar para half-normal plot com envelope
half_normal_envelope <- function(modelo, nsim = 99, seed = 321) {
set.seed(seed)
n <- length(modelo$y)
mu_hat <- fitted(modelo)
# resíduos observados
r_obs <- sort(abs(residuals(modelo, type = "deviance")))
# quantis half-normal teóricos
p_i <- (1:n - 0.375) / (n + 0.25)
hn_q <- qnorm((p_i + 1) / 2)
# simulações
sim_res <- matrix(NA, nrow = n, ncol = nsim)
for (j in 1:nsim) {
y_sim <- rpois(n, lambda = mu_hat)
mod_sim <- glm(y_sim ~ x1 + x2,
family = poisson(link = "log"),
data = dados)
sim_res[, j] <- sort(abs(residuals(mod_sim, type = "deviance")))
}
env_inf <- apply(sim_res, 1, quantile, probs = 0.025)
env_med <- apply(sim_res, 1, quantile, probs = 0.50)
env_sup <- apply(sim_res, 1, quantile, probs = 0.975)
plot(hn_q, r_obs,
xlab = "Quantis half-normal teóricos",
ylab = "|Resíduo de desvio| ordenado",
main = "Half-normal plot com envelope",
pch = 19)
lines(hn_q, env_inf, lty = 2, col = 2)
lines(hn_q, env_med, lty = 1, col = 4)
lines(hn_q, env_sup, lty = 2, col = 2)
}
half_normal_envelope(modelo)
# ============================================================
# 10. GRÁFICOS PADRÃO DO R PARA glm()
# ============================================================
# Alguns métodos plot() para glm podem variar um pouco
# entre versões, mas vale mostrar também:
par(mfrow = c(2, 2))
plot(modelo)
# ============================================================
# 11. COMPARAÇÃO: MODELO COMPLETO VS MODELO SEM PONTOS SUSPEITOS
# ============================================================
ids_remover <- suspeitos$id
dados_limpos <- subset(dados, !(id %in% ids_remover))
modelo_limpo <- glm(y ~ x1 + x2, family = poisson(link = "log"), data = dados_limpos)
cat("\n--- Comparação dos coeficientes ---\n")##
## --- Comparação dos coeficientes ---coef_comp <- cbind(
completo = coef(modelo),
sem_suspeitos = coef(modelo_limpo)
)
print(round(coef_comp, 4))## completo sem_suspeitos
## (Intercept) 1.3595 0.9273
## x1 0.1328 0.1904
## x2 -0.1193 -0.1167##
## --- Resumo modelo original ---##
## Call:
## glm(formula = y ~ x1 + x2, family = poisson(link = "log"), data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.35946 0.21037 6.462 1.03e-10 ***
## x1 0.13282 0.01643 8.083 6.34e-16 ***
## x2 -0.11930 0.03765 -3.169 0.00153 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 258.46 on 79 degrees of freedom
## Residual deviance: 186.03 on 77 degrees of freedom
## AIC: 436.34
##
## Number of Fisher Scoring iterations: 5##
## --- Resumo modelo sem suspeitos ---##
## Call:
## glm(formula = y ~ x1 + x2, family = poisson(link = "log"), data = dados_limpos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.92734 0.28698 3.231 0.00123 **
## x1 0.19044 0.02542 7.493 6.73e-14 ***
## x2 -0.11669 0.05312 -2.197 0.02803 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 111.832 on 64 degrees of freedom
## Residual deviance: 50.959 on 62 degrees of freedom
## AIC: 252.11
##
## Number of Fisher Scoring iterations: 5# ============================================================
# 12. COMENTÁRIOS FINAIS PARA A AULA
# ============================================================
cat("\nINTERPRETAÇÃO DIDÁTICA:\n")##
## INTERPRETAÇÃO DIDÁTICA:## 1. Resíduos grandes sugerem observações mal ajustadas.## 2. Alavancagem alta indica pontos com grande potencial de influência.## 3. Cook e DFFITS ajudam a identificar observações influentes.## 4. O QQ-plot e o half-normal ajudam a avaliar a adequação global do ajuste.## 5. Comparar o modelo com e sem pontos suspeitos mostra o impacto dessas observações.Comentário – Exemplo de Resíduos e Diagnósticos em GLM
O presente exemplo ilustra, de forma integrada, as principais ferramentas de análise de resíduos e diagnóstico em Modelos Lineares Generalizados, considerando um modelo de regressão de Poisson ajustado a dados simulados.
Inicialmente, os dados foram gerados de acordo com uma estrutura conhecida, permitindo controlar o comportamento da variável resposta e, ao mesmo tempo, introduzir artificialmente observações discrepantes e pontos de alta alavancagem. Essa estratégia tem caráter didático, pois possibilita observar, de maneira clara, como diferentes medidas de diagnóstico reagem à presença de problemas no modelo.
Após o ajuste do modelo por meio da função glm(), foram
obtidos diferentes tipos de resíduos, incluindo resíduos brutos, de
Pearson e de desvio. Esses resíduos fornecem informações complementares
sobre o ajuste do modelo, sendo os resíduos de desvio e os padronizados
os mais adequados para análise gráfica.
A análise de influência foi realizada por meio de medidas como alavancagem, distância de Cook e DFFITS. Observações com valores elevados nessas métricas indicam potencial impacto na estimação dos parâmetros, especialmente quando combinadas com resíduos elevados. A identificação dessas observações é fundamental para avaliar a robustez do modelo.
Os gráficos de diagnóstico desempenham papel central na interpretação dos resultados. O gráfico de resíduos versus valores ajustados permite verificar a presença de padrões sistemáticos, enquanto os gráficos de resíduos em função das covariáveis ajudam a identificar possíveis relações não capturadas pelo modelo. O QQ-plot e o gráfico half-normal fornecem uma avaliação global do ajuste, destacando possíveis desvios em relação ao comportamento esperado dos resíduos.
O gráfico de resíduos padronizados versus leverage, complementado por curvas associadas à distância de Cook, permite identificar observações simultaneamente discrepantes e influentes, sendo uma das ferramentas mais completas para diagnóstico.
Por fim, a comparação entre o modelo ajustado com todos os dados e o modelo obtido após a remoção de observações suspeitas evidencia o impacto dessas observações na estimação dos parâmetros. Essa etapa reforça a importância da análise de diagnóstico como parte essencial do processo de modelagem, permitindo construir modelos mais confiáveis e interpretáveis.
Em síntese, este exemplo demonstra que a análise de resíduos e diagnósticos não deve ser tratada como uma etapa opcional, mas sim como parte integrante da modelagem estatística, especialmente em contextos envolvendo Modelos Lineares Generalizados.
4.7. Aplicação utilizando Reamostragem
Aplicação: Validação Cruzada e Bootstrap em Regressão
Fonte: James et al. (2023)
Introdução
Nesta aplicação, serão apresentados dois importantes métodos de reamostragem amplamente utilizados em Estatística e Aprendizado de Máquina: a validação cruzada e o bootstrap. Esses métodos são úteis para avaliar o desempenho preditivo de modelos e para estimar a variabilidade de estimadores quando fórmulas analíticas nem sempre são suficientes ou adequadas.
Usaremos o conjunto de dados Auto, disponível no pacote
ISLR2, para investigar a relação entre o consumo de
combustível do automóvel (mpg) e sua potência
(horsepower). Inicialmente, será aplicada a abordagem de
conjunto de validação, seguida da validação
cruzada leave-one-out (LOOCV) e da validação cruzada
k-fold. Em seguida, será utilizado o bootstrap
para estimar o erro-padrão dos coeficientes de modelos de regressão
linear.
Essa aplicação permite compreender, na prática, como métodos de reamostragem auxiliam na comparação entre modelos com diferentes graus de complexidade, além de fornecerem estimativas mais robustas da incerteza associada aos parâmetros estimados. A estrutura da aplicação segue o exemplo apresentado no material-base, especialmente nas seções sobre Validation Set Approach, LOOCV, k-Fold Cross-Validation e Bootstrap.
## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : int 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : int 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : int 3504 3693 3436 3433 3449 4341 4354 4312 4425 3850 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : int 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : int 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
## - attr(*, "na.action")= 'omit' Named int [1:5] 33 127 331 337 355
## ..- attr(*, "names")= chr [1:5] "33" "127" "331" "337" ...## mpg cylinders displacement horsepower weight
## Min. : 9.00 Min. :3.000 Min. : 68.0 Min. : 46.0 Min. :1613
## 1st Qu.:17.00 1st Qu.:4.000 1st Qu.:105.0 1st Qu.: 75.0 1st Qu.:2225
## Median :22.75 Median :4.000 Median :151.0 Median : 93.5 Median :2804
## Mean :23.45 Mean :5.472 Mean :194.4 Mean :104.5 Mean :2978
## 3rd Qu.:29.00 3rd Qu.:8.000 3rd Qu.:275.8 3rd Qu.:126.0 3rd Qu.:3615
## Max. :46.60 Max. :8.000 Max. :455.0 Max. :230.0 Max. :5140
##
## acceleration year origin name
## Min. : 8.00 Min. :70.00 Min. :1.000 amc matador : 5
## 1st Qu.:13.78 1st Qu.:73.00 1st Qu.:1.000 ford pinto : 5
## Median :15.50 Median :76.00 Median :1.000 toyota corolla : 5
## Mean :15.54 Mean :75.98 Mean :1.577 amc gremlin : 4
## 3rd Qu.:17.02 3rd Qu.:79.00 3rd Qu.:2.000 amc hornet : 4
## Max. :24.80 Max. :82.00 Max. :3.000 chevrolet chevette: 4
## (Other) :3651. Abordagem do conjunto de validação
Dividir os dados em dois subconjuntos: um para treinamento e outro para validação. Ajustamos o modelo no conjunto de treinamento e avaliamos o erro quadrático médio (MSE) no conjunto de validação.
set.seed(1)
train <- sample(1:nrow(Auto), nrow(Auto)/2)
lm.fit1 <- lm(mpg ~ horsepower, data = Auto, subset = train)
pred1 <- predict(lm.fit1, newdata = Auto[-train, ])
mse1 <- mean((Auto$mpg[-train] - pred1)^2)
mse1## [1] 23.26601Modelos polinomiais
Agora ajustamos modelos mais flexíveis, com termos polinomiais de grau 2 e 3.
lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data = Auto, subset = train)
pred2 <- predict(lm.fit2, newdata = Auto[-train, ])
mse2 <- mean((Auto$mpg[-train] - pred2)^2)
mse2## [1] 18.71646lm.fit3 <- lm(mpg ~ poly(horsepower, 3), data = Auto, subset = train)
pred3 <- predict(lm.fit3, newdata = Auto[-train, ])
mse3 <- mean((Auto$mpg[-train] - pred3)^2)
mse3## [1] 18.79401Podemos repetir o procedimento com outra partição dos dados para observar a variabilidade dos resultados.
set.seed(2)
train2 <- sample(1:nrow(Auto), nrow(Auto)/2)
lm.fit1b <- lm(mpg ~ horsepower, data = Auto, subset = train2)
pred1b <- predict(lm.fit1b, newdata = Auto[-train2, ])
mse1b <- mean((Auto$mpg[-train2] - pred1b)^2)
lm.fit2b <- lm(mpg ~ poly(horsepower, 2), data = Auto, subset = train2)
pred2b <- predict(lm.fit2b, newdata = Auto[-train2, ])
mse2b <- mean((Auto$mpg[-train2] - pred2b)^2)
lm.fit3b <- lm(mpg ~ poly(horsepower, 3), data = Auto, subset = train2)
pred3b <- predict(lm.fit3b, newdata = Auto[-train2, ])
mse3b <- mean((Auto$mpg[-train2] - pred3b)^2)
c(MSE_linear = mse1b,
MSE_quadratico = mse2b,
MSE_cubico = mse3b)## MSE_linear MSE_quadratico MSE_cubico
## 25.72651 20.43036 20.38533Observa-se que os valores de MSE podem mudar quando alteramos a divisão entre treinamento e validação. Ainda assim, é possível perceber que modelos polinomiais tendem a apresentar desempenho superior ao modelo puramente linear, sugerindo que a relação entre mpg e horsepower não é estritamente linear.
Desta forma, os erros variam conforme a partição dos dados. Modelos polinomiais tendem a melhorar o ajuste.
2. LOOCV
A validação cruzada leave-one-out utiliza uma observação por vez como validação e ajusta o modelo com todas as demais. Esse procedimento é repetido para todas as observações.
glm.fit <- glm(mpg ~ horsepower, data = Auto)
loocv.err1 <- cv.glm(Auto, glm.fit)$delta[1]
loocv.err1## [1] 24.23151#calculamos os erros LOOCV para modelos polinomiais de grau 1 até 10.
cv.error <- rep(NA, 10)
for (i in 1:10) {
glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)
cv.error[i] <- cv.glm(Auto, glm.fit)$delta[1]
}
cv.error## [1] 24.23151 19.24821 19.33498 19.42443 19.03321 18.97864 18.83305 18.96115
## [9] 19.06863 19.49093plot(1:10, cv.error, type = "b", pch = 19,
xlab = "Grau do polinômio",
ylab = "Erro LOOCV",
main = "LOOCV para modelos polinomiais")
Comentário
O LOOCV permite uma estimativa mais estável do erro de teste do que a simples divisão treino-validação, pois utiliza praticamente toda a amostra para ajuste em cada iteração. Em geral, observa-se uma redução acentuada do erro ao passar do modelo linear para o quadrático, enquanto aumentos adicionais no grau do polinômio produzem pouca ou nenhuma melhora substancial. Isso sugere que um modelo quadrático pode ser suficiente para capturar a estrutura principal da relação entre as variáveis.
3. Validação cruzada k-fold
Agora utilizamos validação cruzada com k = 10, uma escolha bastante comum na prática.
set.seed(17)
cv.error.10 <- rep(NA, 10)
for (i in 1:10) {
glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)
cv.error.10[i] <- cv.glm(Auto, glm.fit, K = 10)$delta[1]
}
cv.error.10## [1] 24.27207 19.26909 19.34805 19.29496 19.03198 18.89781 19.12061 19.14666
## [9] 18.87013 20.95520Visualização gráfica
#plot(1:10, cv.error.10, type="b", pch=19)
plot(1:10, cv.error.10, type = "b", pch = 19,
xlab = "Grau do polinômio",
ylab = "Erro 10-fold CV",
main = "10-fold Cross-Validation para modelos polinomiais")
Comentário
A validação cruzada 10-fold costuma apresentar custo computacional menor do que o LOOCV, mantendo boa qualidade na estimativa do erro preditivo. Os resultados geralmente são muito próximos dos obtidos com LOOCV, reforçando a conclusão de que o modelo quadrático oferece um bom equilíbrio entre simplicidade e desempenho. Assim, o uso de termos polinomiais de ordem muito alta não parece trazer vantagem prática relevante neste exemplo.
4. Bootstrap para estimar a variabilidade dos coeficientes
Agora utilizaremos o bootstrap para estimar o erro-padrão dos coeficientes do modelo linear simples.
Primeiro, definimos uma função que retorna os coeficientes do modelo ajustado em uma amostra bootstrap.
boot.fn <- function(data, index) {
coef(lm(mpg ~ horsepower, data = data, subset = index))
}
#Coeficientes estimados com a amostra original:
boot.fn(Auto, 1:nrow(Auto))## (Intercept) horsepower
## 39.9358610 -0.1578447Gerando duas amostras bootstrap manualmente:
## (Intercept) horsepower
## 40.3404517 -0.1634868## (Intercept) horsepower
## 40.1186906 -0.1577063Agora aplicamos o bootstrap com 1000 repetições.
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Auto, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 39.9358610 0.0553942585 0.843931305
## t2* -0.1578447 -0.0006285291 0.007367396Comparando com os erros-padrão usuais do modelo linear:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.9358610 0.717498656 55.65984 1.220362e-187
## horsepower -0.1578447 0.006445501 -24.48914 7.031989e-81Comentário
Bootstrap fornece estimativas robustas de erro-padrão. O Bootstrap fornece estimativas empíricas do erro-padrão dos coeficientes a partir de reamostragens com reposição. Ao comparar os resultados com os erros-padrão obtidos pela fórmula usual do modelo linear, podemos verificar que os valores são semelhantes, mas não idênticos. Isso ocorre porque o bootstrap não depende tão fortemente das hipóteses teóricas do modelo, podendo refletir melhor a variabilidade real quando certas suposições clássicas são questionáveis. Conforme discutido no material-base, isso é particularmente útil quando há indícios de não linearidade ou quando as condições do modelo linear são apenas aproximadamente atendidas.
5. Bootstrap no modelo quadrático
Agora repetimos a ideia para um modelo quadrático.
boot.fn2 <- function(data, index) {
coef(lm(mpg ~ horsepower + I(horsepower^2), data = data, subset = index))
}
set.seed(1)
boot.res2 <- boot(Auto, boot.fn2, R = 1000)
boot.res2##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Auto, statistic = boot.fn2, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 56.900099702 3.511640e-02 2.0300222526
## t2* -0.466189630 -7.080834e-04 0.0324241984
## t3* 0.001230536 2.840324e-06 0.0001172164Comparando com os erros-padrão usuais do modelo quadrático:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.900099702 1.8004268063 31.60367 1.740911e-109
## horsepower -0.466189630 0.0311246171 -14.97816 2.289429e-40
## I(horsepower^2) 0.001230536 0.0001220759 10.08009 2.196340e-21Comentário
No modelo quadrático, a comparação entre os erros-padrão via bootstrap e os erros-padrão usuais tende a apresentar concordância ainda melhor. Isso sugere que o ajuste quadrático representa de forma mais adequada a relação entre mpg e horsepower, reduzindo parte da inadequação estrutural observada no modelo linear simples. Em termos práticos, isso reforça a interpretação de que a associação entre consumo e potência é curvilínea, e não apenas linear.
Esta aplicação mostrou como métodos de reamostragem podem ser empregados para avaliar modelos estatísticos de maneira mais confiável. A abordagem do conjunto de validação evidenciou que diferentes partições dos dados podem produzir diferentes estimativas de erro, o que motiva o uso de métodos mais estáveis, como LOOCV e k-fold cross-validation.
Os resultados indicaram que o modelo quadrático é mais adequado do que o modelo linear simples para descrever a relação entre mpg e horsepower, enquanto modelos de ordem mais alta não trouxeram ganhos substanciais. Além disso, o bootstrap mostrou-se uma ferramenta útil para estimar a variabilidade dos coeficientes de regressão, especialmente quando se deseja uma alternativa menos dependente das hipóteses clássicas.
Portanto, esta aplicação ilustra bem a relevância dos métodos de reamostragem tanto para seleção de modelos quanto para inferência estatística, constituindo uma etapa fundamental na análise de dados moderna.
Referências bibliográficas do módulo 04
Cordeiro, G. M. (2008). Modelos Lineares Generalizados e Extensões. 400p.Disponível em https://www2.ufjf.br/clecio_ferreira/wp-content/uploads/sites/234/2013/05/Livro-Gauss-e-Clarice.pdf.
James, G., Witten, D., Hastie, T., & Tibshirani, R. (2023). An introduction to statistical learning: With applications in R (2nd ed.). Springer.
Nelder, J. A.; Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society, A, 135, 370-384.
Módulo 05. Redes Bayesianas – Fundamentos e Algoritmos
Uma rede Bayesiana pode ser considerada como uma forma de modelo gráfico probabilístico usado para representar conhecimento sobre um domínio de dados.
É composta por uma estrutura de rede, que consiste em um grafo acíclico direcionado, e um conjunto de tabelas de probabilidade. Os nós da estrutura da rede representam as variáveis e os arcos entre os nós representam relações de dependência entre as variáveis correspondentes.
Um arco começando em um nó X e terminando em um nó Y estabelece X como pai de Y e Y como filho de X. Uma rede Bayesiana pode ser utilizada como um classificador, calculando a probabilidade condicional de um nó, chamado nó classe, dados os valores das probabilidades dos outros nós.
5.1. Conceitos fundamentais
As redes Bayesianas, também conhecidas como modelos gráficos probabilísticos, permitem representar de forma concisa as distribuições de probabilidade conjuntas de um conjunto de variáveis aleatórias \(X = \{X_1, X_2, ..., X_p\}\) por meio de um grafo acíclico direcionado (DAG) \(G = (V, A)\). Cada nó \(v_i \in V\) corresponde a uma variável aleatória \(X_i\).
Estrutura do Grafo e Fatoração de Probabilidade
A correspondência entre a estrutura do grafo e a distribuição de probabilidade é baseada na ausência de uma aresta específica que represente independência condicional. Essa correspondência é formalmente conhecida como independência (Pearl, 1988) e é definida da seguinte forma.
Definição 5.1 (Mapas). Um grafo \(G\) é um mapa de independência (I-map) da distribuição de probabilidade \(P\) se existe uma correspondência um-para-um entre as independências de \(P\) e os nós de \(G\), tal que, para todos os subconjuntos disjuntos \(A, B, C\) de variáveis em \(X\), temos:
\[ A \perp B \mid C \Leftarrow A \perp_G B \mid C \tag{5.1} \]
De forma semelhante, \(G\) é um mapa de dependência (D-map) de \(P\) se:
\[ A \perp B \mid C \Rightarrow A \perp_G B \mid C \tag{5.2} \]
O grafo \(G\) é dito um mapa perfeito de \(P\) se ele for simultaneamente um D-map e um I-map, isto é:
\[ A \perp B \mid C \Longleftrightarrow A \perp_G B \mid C \tag{5.3} \]
Nesse caso, diz-se que \(P\) é isomórfica ou fiel a \(G\).
A correspondência entre a estrutura do DAG \(G\) e as relações de independência condicional que ele representa é esclarecida pelo critério de d-separação, como descrito a seguir.
Definição 5.2 (d-separação). Sejam \(A, B\) e \(C\) subconjuntos disjuntos de nós em um DAG \(G\). Diz-se que \(A\) está d-separado de \(B\) por \(C\), denotado \(A \perp_G B \mid C\), se, ao longo de qualquer sequência de arcos entre um nó em \(A\) e um nó em \(B\), existe um nó \(v\) que satisfaz uma das seguintes condições:
- \(v\) possui arestas convergentes (ou seja, existem duas arestas apontando para \(v\) a partir de nós adjacentes ao caminho) e nem \(v\) nem seus descendentes pertencem a \(C\);
- \(v\) está em \(C\) e não possui arestas convergentes.
A propriedade de Markov das redes Bayesianas, que decorre diretamente da d-separação, permite que a distribuição de probabilidade conjunta das variáveis aleatórias em \(X\) seja fatorada como o produto das distribuições condicionais de probabilidade (as distribuições locais associadas a cada variável \(X_i\)). Essa é uma aplicação direta da regra da cadeia (Korb e Nicholson, 2010). No caso de variáveis discretas, a fatoração da função de distribuição conjunta é dada por:
\[ P_X(X) = \prod_{i=1}^{p} P(X_i \mid \Pi_{X_i}) \tag{5.4} \]
onde \(\Pi_{X_i}\) representa o conjunto de pais de \(X_i\).
No caso de variáveis contínuas, a fatoração da densidade conjunta é:
\[ f_X(X) = \prod_{i=1}^{p} f(X_i \mid \Pi_{X_i}) \tag{5.5} \]
Resultados semelhantes também são válidos para distribuições de probabilidade mistas, ou seja, aquelas que envolvem tanto variáveis discretas quanto contínuas.
Conexões Fundamentais
Existem três configurações fundamentais entre três nós em uma rede Bayesiana:
- Conexão serial: \(A \rightarrow C \rightarrow B\)
- Conexão divergente: \(C \rightarrow A\) e \(C \rightarrow B\)
- Conexão convergente (v-structure): \(A \rightarrow C \leftarrow B\)
Cada uma dessas estruturas implica diferentes relações de independência condicional:
- Em conexões serial e divergente, \(A\) e \(B\) tornam-se independentes dado \(C\)
- Em conexões convergentes, \(A\) e \(B\) são independentes, mas tornam-se dependentes quando condicionado em \(C\).
Considere as conexões fundamentais (Jensen, 2001). Existem três possíveis configurações de três nós e duas arestas. Na conexão convergente (estrutura em V), o nó \(C\) recebe arestas de \(A\) e \(B\), o que viola ambas as condições da Definição 5.2. Portanto, pode-se concluir que \(C\) não d-separa \(A\) e \(B\). Isso implica que \(A\) e \(B\) são dependentes dado \(C\), e, como \(\Pi_A = \emptyset\), \(\Pi_B = \emptyset\) e \(\Pi_C = \{A, B\}\), temos:
\[ P(A,B,C) = P(C \mid A,B)P(A)P(B) \tag{5.6} \]
A partir da propriedade de Markov introduzida na Equação (5.4), observa-se que \(A\) e \(B\) são incondicionalmente independentes. Por outro lado, \(A\) e \(B\) tornam-se dependentes ao se condicionar em \(C\).
Nas conexões serial e divergente, as condições da Definição 5.2 são satisfeitas. Para a conexão serial, temos \(\Pi_A = \emptyset\), \(\Pi_B = \{C\}\) e \(\Pi_C = \{A\}\), portanto:
\[ P(A,B,C) = P(B \mid C)P(C \mid A)P(A) \tag{5.7} \]
Para a conexão divergente, temos \(\Pi_A = \{C\}\), \(\Pi_B = \{C\}\) e \(\Pi_C = \emptyset\), logo:
\[ P(A,B,C) = P(A \mid C)P(B \mid C)P(C) \tag{5.8} \]
Estruturas Equivalentes
Diferentes DAGs podem representar a mesma distribuição de probabilidade se implicarem as mesmas independências condicionais. Esses grafos são chamados de estruturas equivalentes de Markov.
Observa-se que conexões seriais e divergentes resultam em fatorações equivalentes. Cada uma pode ser obtida a partir da outra por aplicações repetidas do teorema de Bayes. Tais estruturas probabilisticamente equivalentes são chamadas de estruturas equivalentes de Markov.
Como a equivalência é simétrica e transitiva, cada conjunto de estruturas equivalentes forma uma classe de equivalência. Para identificar qualquer DAG em uma classe de equivalência, é necessário identificar todos os arcos que pertencem a pelo menos uma estrutura em V (Chickering, 1995).
Essas classes são frequentemente representadas por grafos parcialmente direcionados completos (CPDAGs), que contêm apenas arestas direcionadas que não participam de estruturas em V e arestas não direcionadas para as demais. Tais arestas são chamadas de completas.
Markov Blankets
Markov Blanket de um nó \(A\) é o conjunto mínimo de nós que o torna independente do restante da rede. Esse conjunto inclui:
- os pais de \(A\)
- os filhos de \(A\)
- os outros pais dos filhos de \(A\)
Definição O Markov Blanket de um nó \(A \in V\) é o menor subconjunto \(S \subset V\) tal que:
\[ A \perp V \setminus \{A, S\} \mid S \tag{5.9} \]
Os Markov Blankets facilitam a comparação entre redes Bayesianas e grafos não direcionados. Um grafo Bayesiano pode ser transformado em um grafo não direcionado equivalente por meio do processo de moralização, que consiste em:
- conectar nós não adjacentes que participam de uma estrutura em V;
- remover a direção das arestas.
# ============================================================
# REDES BAYESIANAS
#Fonte: Bayesian Networks in R: with Applications in Systems Biology
# Radhakrishnan Nagarajan; Marco Scutari; Sophie Lèbre
# Versão comentada para R 4.5.3
# Glaucia Bressan e Elisangela Lizzi
# PPGAB - UTFPR
# ============================================================
# ------------------------------------------------------------
# 0) INSTALAÇÃO E CARREGAMENTO DE PACOTES
# ------------------------------------------------------------
# Pacote principal da seção:
# install.packages("bnlearn")
library(bnlearn)
# Pacote opcional para gráficos
# Em alguns sistemas, Rgraphviz precisa ser instalado via Bioconductor:
# install.packages("BiocManager")
# BiocManager::install("Rgraphviz")
# Se quiser usar os gráficos depois:
# library(Rgraphviz)
# ------------------------------------------------------------
# 1) CARREGANDO O CONJUNTO DE DADOS marks
# ------------------------------------------------------------
# dataset marks, disponível em bnlearn. O marks é um dataset: sintético, contínuo, gerado por uma rede Bayesiana gaussiana
#com estrutura causal conhecida
data(marks)
# Estrutura dos dados
str(marks)## 'data.frame': 88 obs. of 5 variables:
## $ MECH: num 77 63 75 55 63 53 51 59 62 64 ...
## $ VECT: num 82 78 73 72 63 61 67 70 60 72 ...
## $ ALG : num 67 80 71 63 65 72 65 68 58 60 ...
## $ ANL : num 67 70 66 70 70 64 65 62 62 62 ...
## $ STAT: num 81 81 81 68 63 73 68 56 70 45 ...## MECH VECT ALG ANL STAT
## 1 77 82 67 67 81
## 2 63 78 80 70 81
## 3 75 73 71 66 81
## 4 55 72 63 70 68
## 5 63 63 65 70 63
## 6 53 61 72 64 73# ------------------------------------------------------------
# 2) CRIANDO E MANIPULANDO ESTRUTURAS DE REDE
# ------------------------------------------------------------
# 2.1 Criando um grafo vazio com os nomes das variáveis
ug <- empty.graph(names(marks))
ug##
## Random/Generated Bayesian network
##
## model:
## [MECH][VECT][ALG][ANL][STAT]
## nodes: 5
## arcs: 0
## undirected arcs: 0
## directed arcs: 0
## average markov blanket size: 0.00
## average neighbourhood size: 0.00
## average branching factor: 0.00
##
## generation algorithm: Empty# 2.2 Adicionando os arcos do grafo não direcionado original
# em bnlearn, arestas não direcionadas são representadas por
# dois arcos em sentidos opostos.
arcs(ug) <- matrix(c(
"MECH", "VECT",
"VECT", "MECH",
"MECH", "ALG",
"ALG", "MECH",
"VECT", "ALG",
"ALG", "VECT",
"ALG", "ANL",
"ANL", "ALG",
"ALG", "STAT",
"STAT", "ALG",
"ANL", "STAT",
"STAT", "ANL"
), byrow = TRUE, ncol = 2,
dimnames = list(NULL, c("from", "to")))
ug##
## Random/Generated Bayesian network
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 6
## undirected arcs: 6
## directed arcs: 0
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 0.00
##
## generation algorithm: Empty# 2.3 Criando o DAG equivalente
dag <- empty.graph(names(marks))
arcs(dag) <- matrix(c(
"VECT", "MECH",
"ALG", "MECH",
"ALG", "VECT",
"ANL", "ALG",
"STAT", "ALG",
"STAT", "ANL"
), byrow = TRUE, ncol = 2,
dimnames = list(NULL, c("from", "to")))
dag##
## Random/Generated Bayesian network
##
## model:
## [STAT][ANL|STAT][ALG|ANL:STAT][VECT|ALG][MECH|VECT:ALG]
## nodes: 5
## arcs: 6
## undirected arcs: 0
## directed arcs: 6
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 1.20
##
## generation algorithm: Empty## [1] "MECH" "VECT" "ALG" "ANL" "STAT"## from to
## [1,] "VECT" "MECH"
## [2,] "ALG" "MECH"
## [3,] "ALG" "VECT"
## [4,] "ANL" "ALG"
## [5,] "STAT" "ALG"
## [6,] "STAT" "ANL"## [1] 6# ------------------------------------------------------------
# 3) CRIANDO A MESMA REDE A PARTIR DE MATRIZ DE ADJACÊNCIA
# ------------------------------------------------------------
# Ordem das variáveis:
# MECH, VECT, ALG, ANL, STAT
ord <- c("MECH", "VECT", "ALG", "ANL", "STAT")
mat <- matrix(c(
0, 0, 0, 0, 0,
1, 0, 0, 0, 0,
1, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 1, 1, 0
), nrow = 5, byrow = TRUE,
dimnames = list(ord, ord))
mat## MECH VECT ALG ANL STAT
## MECH 0 0 0 0 0
## VECT 1 0 0 0 0
## ALG 1 1 0 0 0
## ANL 0 0 1 0 0
## STAT 0 0 1 1 0dag2 <- empty.graph(ord)
amat(dag2) <- mat
# Verificando se dag e dag2 são iguais
all.equal(dag, dag2)## [1] TRUE# ------------------------------------------------------------
# 4) CRIANDO UMA REDE POR MODIFICAÇÕES SUCESSIVAS
# ------------------------------------------------------------
dag3 <- empty.graph(ord)
dag3 <- set.arc(dag3, "VECT", "MECH")
dag3 <- set.arc(dag3, "ALG", "MECH")
dag3 <- set.arc(dag3, "ANL", "ALG")
dag3 <- set.arc(dag3, "STAT", "ALG")
dag3 <- set.arc(dag3, "STAT", "ANL")
dag3##
## Random/Generated Bayesian network
##
## model:
## [VECT][STAT][ANL|STAT][ALG|ANL:STAT][MECH|VECT:ALG]
## nodes: 5
## arcs: 5
## undirected arcs: 0
## directed arcs: 5
## average markov blanket size: 2.40
## average neighbourhood size: 2.00
## average branching factor: 1.00
##
## generation algorithm: Empty## [1] "Different number of directed/undirected arcs"## [1] TRUE# ------------------------------------------------------------
# 5) INVESTIGANDO PROPRIEDADES DA REDE
# ------------------------------------------------------------
# 5.1 Ordem topológica dos nós
node.ordering(dag)## [1] "STAT" "ANL" "ALG" "VECT" "MECH"#A ordem topológica apresenta uma sequência válida para percorrer os nós respeitando a direção dos arcos.
#Isso significa que nenhum nó aparece antes de seus pais na estrutura.
# 5.2 Vizinhança de ALG
nbr(dag, "ALG")## [1] "MECH" "VECT" "ANL" "STAT"## [1] "ALG" "STAT"#O Markov Blanket do nó ANL é formado por ALG e STAT.
#Isso significa que, conhecendo essas duas variáveis, ANL torna-se condicionalmente independente de todas as demais variáveis da rede.
# 5.4 Verificando a simetria do Markov Blanket
"ANL" %in% mb(dag, "ALG")## [1] TRUE## [1] TRUE# 5.5 Filhos, pais e outros pais dos filhos de VECT
child <- children(dag, "VECT")
par <- parents(dag, "VECT")
o.par <- unlist(lapply(child, parents, x = dag))
child## [1] "MECH"## [1] "ALG"## [1] "VECT" "ALG"## [1] "MECH" "ALG"## [1] "MECH" "ALG"# ------------------------------------------------------------
# 6) SCORE E EQUIVALÊNCIA DE REDES
# ------------------------------------------------------------
# 6.1 Score log-likelihood gaussiano
bnlearn::score(dag, data = marks, type = "loglik-g")## [1] -1695.589#para garantir que o score utilizado é o do pacote bnlearn
#Esse valor representa o log-likelihood gaussiano da rede.
#Ele mede quão bem a estrutura proposta ajusta-se aos dados observados.
#Em geral, valores maiores (menos negativos) indicam melhor ajuste.
# 6.2 Criando uma rede equivalente ao inverter o arco STAT -> ANL
dag.eq <- reverse.arc(dag, "STAT", "ANL")
dag.eq##
## Random/Generated Bayesian network
##
## model:
## [ANL][STAT|ANL][ALG|ANL:STAT][VECT|ALG][MECH|VECT:ALG]
## nodes: 5
## arcs: 6
## undirected arcs: 0
## directed arcs: 6
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 1.20
##
## generation algorithm: Empty## [1] -1695.589# ------------------------------------------------------------
# 7) V-STRUCTURES E GRAFO MORAL
# ------------------------------------------------------------
# V-structures presentes em dag
vstructs(dag)## X Z Y## X Z Y# Observação:
# Em algumas versões antigas aparecia algo como vstructs(..., moral = TRUE),
# mas isso não funciona nas versões atuais do bnlearn.
# Para trabalhar com moralização, use moral() separadamente.
moral(dag)##
## Random/Generated Bayesian network
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 6
## undirected arcs: 6
## directed arcs: 0
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 0.00
##
## generation algorithm: Empty##
## Random/Generated Bayesian network
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 6
## undirected arcs: 6
## directed arcs: 0
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 0.00
##
## generation algorithm: Empty# ------------------------------------------------------------
# 8) PLOTANDO AS REDES
# Requer Rgraphviz
# ------------------------------------------------------------
if (requireNamespace("Rgraphviz", quietly = TRUE)) {
library(Rgraphviz)
hl2 <- list(arcs = vstructs(dag2, arcs = TRUE))
hl3 <- list(arcs = vstructs(dag3, arcs = TRUE))
#rodar uma linha por vez, sao 4 figuras
graphviz.plot(dag2, highlight = hl2, layout = "fdp", main = "dag2")
graphviz.plot(dag3, highlight = hl3, layout = "fdp", main = "dag3")
graphviz.plot(cpdag(dag2), highlight = hl2, layout = "fdp", main = "cpdag(dag2)")
graphviz.plot(cpdag(dag3), highlight = hl3, layout = "fdp", main = "cpdag(dag3)")
} else {
message("Rgraphviz não está instalado. Bloco de gráficos foi pulado.")
}
5.2. Modelagem de Redes Bayesianas Estáticas
O processo de construção de uma rede Bayesiana a partir de dados é conhecido como aprendizado e envolve duas etapas principais:
- aprendizado da estrutura
- aprendizado dos parâmetros
A primeira etapa, chamada aprendizado da estrutura, consiste em identificar a estrutura do grafo que melhor representa as dependências entre os dados. Idealmente, deve-se minimizar a divergência entre a distribuição obtida e a distribuição verdadeira no espaço de probabilidade.
Diversas abordagens têm sido propostas, podendo ser classificadas em três categorias principais:
- métodos baseados em restrições;
- métodos baseados em pontuação;
- métodos híbridos.
A segunda etapa, chamada aprendizado de parâmetros, consiste em estimar os parâmetros da distribuição global a partir das distribuições locais associadas à estrutura previamente obtida.
Aprendizado de Estrutura
O objetivo é identificar a estrutura do grafo que melhor representa as dependências entre as variáveis nos dados. Idealmente, busca-se uma distribuição que se aproxime da verdadeira distribuição dos dados.
Existem três abordagens principais:
- métodos baseados em restrições
- métodos baseados em pontuação
- métodos híbridos
5.3. Algoritmos de Aprendizagem Baseados em Restrições Constraint-Based)
Esses métodos utilizam testes de independência condicional para determinar a estrutura do grafo. O algoritmo Inductive Causation (IC) é um exemplo clássico dessa abordagem (Verma and Pearl, 1991).
Inicialmente, constrói-se um grafo completo. Em seguida, as arestas são removidas com base em testes de independência condicional. Posteriormente, identificam-se estruturas em V e orientam-se as arestas restantes.
O primeiro passo identifica quais pares de variáveis estão conectados por um arco, independentemente de sua direção. Essas variáveis não podem ser independentes dado qualquer outro subconjunto de variáveis, pois não podem ser d-separadas. Esse passo também pode ser visto como um procedimento de seleção regressiva (backward selection), iniciando a partir do modelo saturado com um grafo completo e realizando sua poda com base em testes estatísticos de independência condicional.
O segundo passo trata da identificação das v-structures entre todos os pares de nós não adjacentes A e B que possuem um vizinho comum C. Por definição, v-structures são a única conexão fundamental na qual os dois nós não adjacentes não são independentes condicionalmente ao terceiro. Portanto, se existir um subconjunto de nós que contenha C e que d-separe A e B, então os três nós fazem parte de uma v-structure centrada em C. Essa condição pode ser verificada realizando um teste de independência condicional para A e B contra todos os possíveis subconjuntos de seus vizinhos comuns que incluam C. Ao final do segundo passo, tanto o esqueleto quanto as v-structures da rede são conhecidos, de modo que a classe de equivalência à qual a rede bayesiana pertence fica unicamente identificada.
O terceiro e último passo do algoritmo IC identifica os arcos compelidos (compelled arcs) e os orienta recursivamente para obter o DAG parcialmente completo (Completed Partially Directed Acyclic Graph – CPDAG) que descreve a classe de equivalência identificada pelos passos anteriores.
Variações desse algoritmo frequentemente utilizam o conceito de Markov Blanket para reduzir a complexidade computacional (número de testes necessários).
Algoritmos Baseados em Restrição para Aprendizagem de Estrutura
PC: primeira aplicação prática do algoritmo IC (Spirtes et al., 2001), um procedimento de seleção regressiva (backward selection) a partir do grafo saturado.
Grow-Shrink (GS): baseado no algoritmo Grow-Shrink de cobertor de Markov (Markov blanket) (Margaritis, 2003), uma abordagem simples de detecção do cobertor de Markov por seleção progressiva (forward selection).
Incremental Association (IAMB): baseado no algoritmo Incremental Association de cobertor de Markov (Tsamardinos et al., 2003), um esquema de seleção em duas fases baseado em uma seleção progressiva seguida de uma seleção regressiva.
Fast Incremental Association (Fast-IAMB): uma variante do IAMB que utiliza seleção progressiva especulativa passo a passo (speculative stepwise forward selection) para reduzir o número de testes de independência condicional (Yaramakala e Margaritis, 2005).
Interleaved Incremental Association (Inter-IAMB): outra variante do IAMB que utiliza seleção progressiva passo a passo (forward stepwise selection) (Tsamardinos et al., 2003) para evitar falsos positivos na fase de detecção do cobertor de Markov.
# ------------------------------------------------------------
# 9) APRENDIZADO DE ESTRUTURA COM bnlearn
#
# ------------------------------------------------------------
# 9.1 Grow-Shrink
bn.gs <- gs(marks)
bn.gs##
## Bayesian network learned via Constraint-based methods
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 6
## undirected arcs: 6
## directed arcs: 0
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 0.00
##
## learning algorithm: Grow-Shrink
## conditional independence test: Pearson's Correlation
## alpha threshold: 0.05
## tests used in the learning procedure: 80##
## Bayesian network learned via Constraint-based methods
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 6
## undirected arcs: 6
## directed arcs: 0
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 0.00
##
## learning algorithm: IAMB
## conditional independence test: Pearson's Correlation
## alpha threshold: 0.05
## tests used in the learning procedure: 97## [1] TRUE##
## Bayesian network learned via Constraint-based methods
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 6
## undirected arcs: 6
## directed arcs: 0
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 0.00
##
## learning algorithm: Inter-IAMB
## conditional independence test: Pearson's Correlation
## alpha threshold: 0.05
## tests used in the learning procedure: 100## [1] TRUE##
## Bayesian network learned via Constraint-based methods
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 6
## undirected arcs: 6
## directed arcs: 0
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 0.00
##
## learning algorithm: IAMB
## conditional independence test: Pearson's Correlation (MC)
## alpha threshold: 0.05
## permutations: 5000
## tests used in the learning procedure: 101## [1] TRUE##
## Bayesian network learned via Constraint-based methods
##
## model:
## [partially directed graph]
## nodes: 5
## arcs: 6
## undirected arcs: 1
## directed arcs: 5
## average markov blanket size: 3.20
## average neighbourhood size: 2.40
## average branching factor: 1.00
##
## learning algorithm: PC (Stable)
## conditional independence test: Pearson's Correlation
## alpha threshold: 0.05
## tests used in the learning procedure: 60##
## Bayesian network learned via Score-based methods
##
## model:
## [MECH][VECT|MECH][ALG|MECH:VECT][ANL|ALG][STAT|ALG:ANL]
## nodes: 5
## arcs: 6
## undirected arcs: 0
## directed arcs: 6
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 1.20
##
## learning algorithm: Hill-Climbing
## score: BIC (Gauss.)
## penalization coefficient: 2.238668
## tests used in the learning procedure: 34
## optimized: TRUE# ------------------------------------------------------------
# 10) COMPARANDO SCORES DAS REDES APRENDIDAS
# ------------------------------------------------------------
#bnlearn::score(bn.gs, data = marks, type = "bic-g")
#bnlearn::score(bn.hc, data = marks, type = "bic-g")
#bnlearn::score(bn.pc, data = marks, type = "bic-g")
#Esses valores permitem comparar quantitativamente diferentes estruturas aprendidas.
#Como o BIC penaliza modelos excessivamente complexos, ele busca equilíbrio entre ajuste aos dados e simplicidade estrutural.
# ------------------------------------------------------------
# 11) PLOTANDO AS REDES APRENDIDAS
# ------------------------------------------------------------
#if (requireNamespace("Rgraphviz", quietly = TRUE)) {
library(Rgraphviz)
#plotar separadamente
graphviz.plot(bn.gs, main = "Grow-Shrink")
#} else {
# message("Rgraphviz não está instalado. Gráficos das redes aprendidas foram pulados.")
#}
# ------------------------------------------------------------
# 12) BLOCO OPCIONAL: pcalg
# ------------------------------------------------------------
# Este bloco é opcional porque costuma dar mais problema de dependência.
if (requireNamespace("pcalg", quietly = TRUE) &&
requireNamespace("graph", quietly = TRUE)) {
library(pcalg)
library(graph)
suff <- list(C = cor(marks), n = nrow(marks))
pc.fit <- pcalg::pc(
suffStat = suff,
indepTest = gaussCItest,
alpha = 0.05,
labels = colnames(marks)
)
pc.fit
} else {
message("pcalg/graph não estão instalados. Bloco opcional do pcalg foi pulado.")
}## Object of class 'pcAlgo', from Call:
## pcalg::pc(suffStat = suff, indepTest = gaussCItest, alpha = 0.05,
## labels = colnames(marks))
## Number of undirected edges: 1
## Number of directed edges: 5
## Total number of edges: 65.4. Algoritmos de Aprendizagem Baseados em Pontuação (Score-Based)
Os algoritmos baseados em pontuação atribuem uma medida de qualidade a cada estrutura candidata e procuram maximizar essa pontuação.
Exemplos incluem:
- Greedy search, como hill-climbing ou Busca Tabu (Bouckaert, 1995): exploram o espaço de busca partindo de uma estrutura de rede (usualmente o grafo vazio) e adicionando, deletando ou invertendo um arco por vez até que o score não possa mais ser melhorado.
Passos do algoritmo Hill Climbing
Escolher uma estrutura inicial de rede \(G\), geralmente um grafo vazio, embora outra estrutura inicial também possa ser utilizada.
Calcular a pontuação (score) da estrutura inicial \(G\) com base em uma função de avaliação apropriada.
Armazenar essa pontuação como a melhor pontuação atual (maxscore).
Gerar todas as modificações locais possíveis na rede atual que não produzam ciclos, incluindo:
- adição de um arco;
- remoção de um arco;
- inversão da direção de um arco.
- adição de um arco;
Calcular a pontuação de cada rede modificada gerada no passo anterior.
Substituir a estrutura atual pela estrutura modificada caso alguma apresente pontuação superior à da rede corrente.
Repetir os passos 4 a 6 até que nenhuma modificação produza melhora na pontuação, retornando então a melhor estrutura encontrada.
algoritmos genéticos (Larranaga et al., 1997): imitam a evolução natural através da seleção iterativa dos modelos “mais aptos” e da hibridização de suas características. O espaço de busca é explorado através dos operadores estocásticos de cruzamento (que combina a estrutura de duas redes) e mutação (que introduz alterações aleatórias).
simulated annealing (Bouckaert, 1995): realiza uma busca local estocástica, aceitando mudanças que aumentam a pontuação da rede e, ao mesmo tempo, permitindo mudanças que a diminuem com uma probabilidade inversamente proporcional à diminuição da pontuação.
As funções de pontuação mais comuns incluem:
- BDe (Bayesian Dirichlet equivalent);
- BIC (Bayesian Information Criterion).
5.5. Algoritmos Híbridos e heurísticas combinadas
Os algoritmos híbridos combinam abordagens baseadas em restrições e em pontuação.
Um exemplo é o algoritmo Max-Min Hill-Climbing (MMHC), de Tsamardinos et al. (2006), que utiliza heurísticas para restringir o espaço de busca e, posteriormente, aplica otimização baseada em pontuação.
Escolha de Distribuições, Testes de independência condicional e Scores de rede
A modelagem depende do tipo de dados:
- Dados discretos: utilizam tabelas de probabilidade condicional
- Dados contínuos: assumem frequentemente distribuição normal
Testes comuns de independência incluem:
- Informação mútua
- Teste qui-quadrado de Pearson
- Testes baseados em correlação parcial
Em princípio, há muitas escolhas possíveis tanto para as funções de distribuição global quanto local, dependendo da natureza dos dados e dos objetivos da análise. No entanto, a literatura tem se concentrado principalmente em dois casos:
Variáveis multinomiais: adequadas para conjuntos de dados discretos/categóricos e, com frequência, também para dados discretizados. Tanto a distribuição global quanto as distribuições locais são multinomiais, e estas últimas são representadas como tabelas de probabilidade condicional (CPTs). Por esse motivo, esses modelos são frequentemente chamados de redes Bayesianas discretas, em contraste com as redes contínuas ou mistas.
Variáveis normais multivariadas: quando essa representação é utilizada para conjuntos de dados contínuos, as distribuições locais associadas aos nós são distribuições normais multivariadas, em que algumas distribuições locais são lineares em variáveis ligadas por restrições lineares. As distribuições locais são, na verdade, modelos lineares com os pais desempenhando o papel de variáveis explicativas. Essas redes são conhecidas como redes Bayesianas Gaussianas (Geiger e Heckerman, 1994; Neapolitan, 2003).
Outras suposições distribucionais exigem algoritmos de aprendizado ad hoc ou apresentam sérias limitações devido à dificuldade de especificar as funções de distribuição no espaço de parâmetros. Por exemplo, modelos para dados não exatos, como os apresentados em Bøttcher (2004), são inconvenientes, já que a escolha das distribuições é o nó.
A escolha de um conjunto particular de distribuições globais e locais pode, às vezes, ser menos crucial do que a escolha dos testes de independência condicional a serem usados para aprender a estrutura da rede Bayesiana.
Os testes de independência condicional e os escores de rede para dados discretos são funções das frequências \(n_{ijk}\) na tabela de contingência das observações para as variáveis aleatórias \(X\) e \(Y\), e de todas as configurações das variáveis de condicionamento \(Z\). Dois testes de independência condicional comumente utilizados são os seguintes:
- Informação mútua (Cover e Thomas, 2006), uma medida de distância da teoria da informação definida como
\[ \mathrm{MI}(X, Y \mid Z) = \sum_{i=1}^{R} \sum_{j=1}^{C} \sum_{k=1}^{L} \frac{n_{ijk}}{n} \log \frac{n_{ijk} n_{++k}}{n_{i+k} n_{+jk}}. \tag{5.10} \]
Ela é proporcional ao teste da razão de verossimilhanças \(G^2\) (eles diferem por um fator \(2n\), em que \(n\) é o tamanho da amostra), e está relacionada ao desvio dos modelos testados.
- O clássico teste \(\chi^2\) de Pearson para tabelas de contingência,
\[ X^2(X, Y \mid Z) = \sum_{i=1}^{R} \sum_{j=1}^{C} \sum_{k=1}^{L} \frac{(n_{ijk} - m_{ijk})^2}{m_{ijk}}, \qquad \text{onde } m_{ijk} = \frac{n_{i+k} n_{+jk}}{n_{++k}}. \tag{5.11} \]
Em ambos os casos, a hipótese nula de independência pode ser testada tanto pela distribuição assintótica \(\chi^2_{(R-1)(C-1)L}\) quanto pela abordagem de permutação de Monte Carlo descrita em Edwards (2000). Outras escolhas possíveis são o teste exato de Fisher e os estimadores por encolhimento (shrinkage estimators) para a informação mútua definidos por Hausser e Strimmer (2009) e estudados em Scutari e Brogini (2012).
Os scores de rede comumente encontrados na literatura incluem os seguintes:
O score bayesiano Dirichlet equivalente (BDe), a densidade a posteriori associada a uma distribuição a priori uniforme tanto sobre o espaço das estruturas de rede quanto sobre os parâmetros de cada distribuição local (Heckerman et al., 1995).
O Critério de Informação Bayesiano (BIC), o score de razão de verossimilhança definido como
\[ \mathrm{BIC} = \sum_{i=1}^{n} \log P_X\!\left(X_i \mid \Pi_{X_i}\right) - \frac{d}{2}\log n, \tag{5.12} \]
em que \(d\) é o número de parâmetros da distribuição global. É universalmente equivalente ao princípio do comprimento mínimo de descrição (minimum description length — MDL) proposto por Rissanen (2007), ainda que tenha uma interpretação completamente diferente. O BIC converge assintoticamente para o escore BDe a posteriori.
Diz-se que essas funções de score são equivalentes por score, uma vez que atribuem o mesmo escore a redes equivalentes de Markov. Elas também são decomponíveis em componentes correspondentes a cada nó, o que representa uma vantagem computacional significativa ao aprender a estrutura da rede: as únicas partes que precisam ser computadas são aquelas que diferem entre as redes comparadas.
No caso contínuo, os testes de independência condicional e os escores de rede são funções dos coeficientes de correlação parcial \(\rho_{XY \mid Z}\) de \(X\) e \(Y\), dado \(Z\). Dois testes de independência condicional comumente utilizados são os seguintes:
- O teste t exato para o coeficiente de correlação parcial de Pearson, definido como
\[ t(X, Y \mid Z) = \rho_{XY \mid Z} \sqrt{ \frac{n - 2 - |Z|}{1 - \rho_{XY \mid Z}^{2}} }, \tag{5.13} \]
e distribuído como uma t de Student com \(n - |Z| - 2\) graus de liberdade.
- O teste Z de Fisher, uma transformação do coeficiente de correlação linear com distribuição normal assintótica, definido como
\[ Z(X, Y \mid Z) = \sqrt{n - |Z| - 3} \cdot \frac{1}{2} \log \frac{1 + \rho_{XY \mid Z}}{1 - \rho_{XY \mid Z}}. \tag{5.14} \]
Ambos os testes também podem ser realizados usando abordagens de permutação de Monte Carlo, como as descritas em Legendre (2000). Outras escolhas possíveis são o teste de informação mútua definido em Kullback (1968), que é proporcional ao teste correspondente da razão de verossimilhanças logarítmicas, ou os estimadores por encolhimento desenvolvidos por Schäfer e Strimmer (2005).
Os scores de rede comumente usados são, novamente, o BIC, desta vez definido como
\[ \mathrm{BIC} = \sum_{i=1}^{n} \log f_X\!\left(X_i \mid \Pi_{X_i}\right) - \frac{d}{2}\log n \tag{5.15} \]
e o escore bayesiano Gaussiano equivalente (BGe), a densidade a posteriori da rede associada a uma distribuição a priori uniforme tanto sobre o espaço das estruturas de rede quanto sobre os parâmetros das distribuições locais (Geiger e Heckerman, 1994).
Aprendizado de Parâmetros
Após definir a estrutura, os parâmetros das distribuições condicionais são estimados. Isso pode ser feito por:
- máxima verossimilhança
- métodos bayesianos
# ------------------------------------------------------------
# 13) APRENDIZADO DE PARÂMETROS
# ------------------------------------------------------------
# 13.1 Ajustando os parâmetros da rede bn
bn.hc <- hc(marks) #usar uma rede que já sai direcionada
#deixamos gs(), iamb() e pc.stable() mais para a parte de aprendizado de estrutura.
fitted <- bn.fit(bn.hc, data = marks)
fitted$ALG##
## Parameters of node ALG (Gaussian distribution)
##
## Conditional density: ALG | MECH + VECT
## Coefficients:
## (Intercept) MECH VECT
## 25.3619809 0.1833755 0.3577122
## Standard deviation of the residuals: 8.080725# 13.2 Substituindo manualmente os parâmetros do nó ALG
#parents(bn.hc, "ALG") para descobrir os pais
fitted$ALG <- list(
coef = c("(Intercept)" = 25,
"MECH" = 0.5,
"VECT" = 0.25),
sd = 6.5
)
fitted$ALG##
## Parameters of node ALG (Gaussian distribution)
##
## Conditional density: ALG | MECH + VECT
## Coefficients:
## (Intercept) MECH VECT
## 25.00 0.50 0.25
## Standard deviation of the residuals: 6.5# 13.3 Criando listas manuais de parâmetros para outros nós
#######################################################
dag <- empty.graph(c("MECH", "VECT", "ALG", "ANL", "STAT"))
arcs(dag) <- matrix(c(
"VECT", "MECH",
"ALG", "MECH",
"ALG", "VECT",
"ANL", "ALG",
"STAT", "ALG",
"STAT", "ANL"
), byrow = TRUE, ncol = 2,
dimnames = list(NULL, c("from", "to")))
MECH.par <- list(
coef = c("(Intercept)" = -10,
"VECT" = 0.5,
"ALG" = 0.6),
sd = 13
)
VECT.par <- list(
coef = c("(Intercept)" = 10,
"ALG" = 0.5),
sd = 13
)
ALG.par <- list(
coef = c("(Intercept)" = 25,
"ANL" = 0.5,
"STAT" = 0.25),
sd = 6.5
)
ANL.par <- list(
coef = c("(Intercept)" = 10,
"STAT" = 0.5),
sd = 13
)
STAT.par <- list(
coef = c("(Intercept)" = 50),
sd = 13
)
dist <- list(
MECH = MECH.par,
VECT = VECT.par,
ALG = ALG.par,
ANL = ANL.par,
STAT = STAT.par
)
fitted2 <- custom.fit(dag, dist = dist)
fitted2##
## Bayesian network parameters
##
## Parameters of node MECH (Gaussian distribution)
##
## Conditional density: MECH | VECT + ALG
## Coefficients:
## (Intercept) VECT ALG
## -10.0 0.5 0.6
## Standard deviation of the residuals: 13
##
## Parameters of node VECT (Gaussian distribution)
##
## Conditional density: VECT | ALG
## Coefficients:
## (Intercept) ALG
## 10.0 0.5
## Standard deviation of the residuals: 13
##
## Parameters of node ALG (Gaussian distribution)
##
## Conditional density: ALG | ANL + STAT
## Coefficients:
## (Intercept) ANL STAT
## 25.00 0.50 0.25
## Standard deviation of the residuals: 6.5
##
## Parameters of node ANL (Gaussian distribution)
##
## Conditional density: ANL | STAT
## Coefficients:
## (Intercept) STAT
## 10.0 0.5
## Standard deviation of the residuals: 13
##
## Parameters of node STAT (Gaussian distribution)
##
## Conditional density: STAT
## Coefficients:
## (Intercept)
## 50
## Standard deviation of the residuals: 13Discretização
Quando há variáveis contínuas, pode-se aplicar discretização para facilitar o aprendizado. Métodos incluem:
- discretização baseada em conhecimento prévio
- métodos baseados em frequência
- métodos baseados em entropia
Essa etapa envolve um trade-off entre precisão e eficiência computacional.
# ------------------------------------------------------------
# 14) DISCRETIZAÇÃO
#
# ------------------------------------------------------------
# Discretizando marks em duas classes
dmarks <- discretize(marks, method = "interval", breaks = 2)
str(dmarks)## 'data.frame': 88 obs. of 5 variables:
## $ MECH: Factor w/ 2 levels "[0,38.5]","(38.5,77]": 2 2 2 2 2 2 2 2 2 2 ...
## $ VECT: Factor w/ 2 levels "[9,45.5]","(45.5,82]": 2 2 2 2 2 2 2 2 2 2 ...
## $ ALG : Factor w/ 2 levels "[15,47.5]","(47.5,80]": 2 2 2 2 2 2 2 2 2 2 ...
## $ ANL : Factor w/ 2 levels "[9,39.5]","(39.5,70]": 2 2 2 2 2 2 2 2 2 2 ...
## $ STAT: Factor w/ 2 levels "[9,45]","(45,81]": 2 2 2 2 2 2 2 2 2 1 ...## MECH VECT ALG ANL STAT
## 1 (38.5,77] (45.5,82] (47.5,80] (39.5,70] (45,81]
## 2 (38.5,77] (45.5,82] (47.5,80] (39.5,70] (45,81]
## 3 (38.5,77] (45.5,82] (47.5,80] (39.5,70] (45,81]
## 4 (38.5,77] (45.5,82] (47.5,80] (39.5,70] (45,81]
## 5 (38.5,77] (45.5,82] (47.5,80] (39.5,70] (45,81]
## 6 (38.5,77] (45.5,82] (47.5,80] (39.5,70] (45,81]##
## Bayesian network learned via Constraint-based methods
##
## model:
## [undirected graph]
## nodes: 5
## arcs: 2
## undirected arcs: 2
## directed arcs: 0
## average markov blanket size: 0.80
## average neighbourhood size: 0.80
## average branching factor: 0.00
##
## learning algorithm: Grow-Shrink
## conditional independence test: Mutual Information (disc.)
## alpha threshold: 0.05
## tests used in the learning procedure: 36##
## Bayesian network learned via Score-based methods
##
## model:
## [MECH][VECT|MECH][ALG|VECT][ANL|ALG][STAT|ALG]
## nodes: 5
## arcs: 4
## undirected arcs: 0
## directed arcs: 4
## average markov blanket size: 1.60
## average neighbourhood size: 1.60
## average branching factor: 0.80
##
## learning algorithm: Hill-Climbing
## score: BIC (disc.)
## penalization coefficient: 2.238668
## tests used in the learning procedure: 26
## optimized: TRUE## [1] "Different number of directed/undirected arcs"# Ajustando parâmetros para a rede discreta
fitted.disc <- bn.fit(bn.dhc, data = dmarks)
fitted.disc$ALG##
## Parameters of node ALG (multinomial distribution)
##
## Conditional probability table:
##
## VECT
## ALG [9,45.5] (45.5,82]
## [15,47.5] 0.5806452 0.2280702
## (47.5,80] 0.4193548 0.77192985.6. Causalidade de Pearl
Nesta seção, introduzimos o conceito de causalidade em redes Bayesianas, conectando a estrutura probabilística dos grafos com interpretações causais.
Nas seções anteriores, redes Bayesianas foram definidas com base em independência condicional e fatoração de distribuições. No entanto, essas mesmas estruturas também podem ser interpretadas como modelos causais.
Uma rede Bayesiana pode ser interpretada como um modelo causal quando:
- cada aresta representa uma relação de causa e efeito;
- a estrutura do grafo reflete mecanismos causais do sistema;
- as dependências probabilísticas correspondem a relações causais reais.
Assim, uma rede Bayesiana pode ser vista como um modelo que descreve como os dados foram gerados.
Para que uma rede Bayesiana tenha interpretação causal, alguns pressupostos são necessários.
Condição de Markov Causal
Uma variável é independente de seus não-efeitos (não descendentes), dado seus pais. Isso implica que:
- os pais de um nó representam suas causas diretas;
- toda a informação relevante flui através dessas relações.
Faithfulness (ou fidelidade)
A estrutura do grafo deve refletir corretamente as dependências presentes nos dados.
Ausência de variáveis latentes
Não deve haver variáveis ocultas (não observadas) influenciando o sistema.
Caso existam variáveis latentes:
- podem surgir dependências espúrias;
- a estrutura aprendida pode não refletir relações causais reais.
OBS: Na prática, esses pressupostos são difíceis de verificar. Por isso, inferência causal a partir de dados observacionais não é muito comum.
# ------------------------------------------------------------
# 15) BLOCO OPCIONAL: continuar o aprendizado com score Bayesiano
# ------------------------------------------------------------
# usando apenas bnlearn
#library(bnlearn)
#data(marks)
# Hill-climbing com score BGe
bn.hc.bge <- hc(marks, score = "bge", iss = 5)
bn.hc.bge##
## Bayesian network learned via Score-based methods
##
## model:
## [MECH][VECT|MECH][ALG|VECT][ANL|ALG][STAT|ALG]
## nodes: 5
## arcs: 4
## undirected arcs: 0
## directed arcs: 4
## average markov blanket size: 1.60
## average neighbourhood size: 1.60
## average branching factor: 0.80
##
## learning algorithm: Hill-Climbing
## score: Bayesian Gaussian (BGe)
## graph prior: Uniform
## imaginary sample size (normal): 1
## imaginary sample size (Wishart): 7
## tests used in the learning procedure: 26
## optimized: TRUE## [1] -1798.872##
## Bayesian network learned via Score-based methods
##
## model:
## [MECH][VECT|MECH][ALG|MECH:VECT][ANL|ALG][STAT|ALG:ANL]
## nodes: 5
## arcs: 6
## undirected arcs: 0
## directed arcs: 6
## average markov blanket size: 2.40
## average neighbourhood size: 2.40
## average branching factor: 1.20
##
## learning algorithm: Hill-Climbing
## score: BIC (Gauss.)
## penalization coefficient: 2.238668
## tests used in the learning procedure: 34
## optimized: TRUE## [1] -1731.4075.7. Aplicação : perfis de expressão gênica (gene expression profile)
Nesta seção, apresentamos uma aplicação prática de redes Bayesianas estáticas a dados de expressão gênica, conforme discutido no capítulo. O objetivo é ilustrar como essas redes podem ser utilizadas na análise de dados biológicos, com ênfase em três aspectos principais: a construção de redes médias por model averaging, a escolha do limiar de significância e o tratamento de dados intervencionais.
A aplicação é baseada no conjunto de dados de Sachs et al. (2005), amplamente utilizado na literatura para estudos de inferência de redes de sinalização celular.
# ============================================================
# Application to Gene Expression Profiles
# Profa. Glaucia Bressan e Profa. Elisangela Lizzi
# Foco em bnlearn
#FONTE: BAYESIAN NETWOKS IN R: with applications in System Biology
#Radhakrishnan Nagarajan; Marco Scutari; Sophie Lèbre
# ============================================================
# ------------------------------------------------------------
# 0) Pacotes
# ------------------------------------------------------------
# Instale, se necessário:
# install.packages("bnlearn")
# install.packages("Rgraphviz") # opcional para gráficos; pode exigir Bioconductor
library(bnlearn)
# Se quiser gráficos:
library(Rgraphviz)
# ------------------------------------------------------------
# 1) Defina o diretório onde estão os arquivos do Sachs
# ------------------------------------------------------------
# Ajuste este caminho no seu computador, se necessário:
setwd("/Users/elisangelalizzi/Documents/01UTFPR- Parcial/2026-01/Pós graduação- PPGAB e PPGBIOINFO/Completo-Métodos Avançados")
# Arquivos esperados:
obs_file <- "sachs.data.txt"
int_file <- "sachs.interventional.txt"
# ------------------------------------------------------------
# 2) Leitura dos dados observacionais
# ------------------------------------------------------------
# dados observacionais de sinalização proteica.
# Aqui lemos o arquivo bruto.
sachs <- read.table(obs_file, header = TRUE)
str(sachs)## 'data.frame': 853 obs. of 11 variables:
## $ Raf : num 26.4 35.9 59.4 73 33.7 18.8 44.9 47.4 104 21.1 ...
## $ Mek : num 13.2 16.5 44.1 82.8 19.8 3.75 36.5 15 61.5 21.5 ...
## $ Plcg: num 8.82 12.3 14.6 23.1 5.19 17.6 10.4 14.6 10.6 1.88 ...
## $ PIP2: num 18.3 16.8 10.2 13.5 9.73 22.1 132 30.5 21.1 205 ...
## $ PIP3: num 58.8 8.13 13 1.29 24.8 10.9 16.3 17.5 41.8 43.7 ...
## $ Erk : num 6.61 18.6 14.9 5.83 21.1 11.9 8.66 20.2 11.5 13.2 ...
## $ Akt : num 17 32.5 32.5 11.8 46.1 25.7 17.9 45.3 23.5 135 ...
## $ PKA : num 414 352 403 528 305 610 835 466 445 213 ...
## $ PKC : num 17 3.37 11.4 13.7 4.66 13.7 15 6.44 29.2 14.6 ...
## $ P38 : num 44.9 16.5 31.9 28.6 25.7 49.1 35.9 24.4 61 26.7 ...
## $ Jnk : num 40 61.5 19.5 23.1 81.3 57.8 18.1 20 25.3 101 ...## Raf Mek Plcg PIP2 PIP3 Erk Akt PKA PKC P38 Jnk
## 1 26.4 13.20 8.82 18.30 58.80 6.61 17.0 414 17.00 44.9 40.0
## 2 35.9 16.50 12.30 16.80 8.13 18.60 32.5 352 3.37 16.5 61.5
## 3 59.4 44.10 14.60 10.20 13.00 14.90 32.5 403 11.40 31.9 19.5
## 4 73.0 82.80 23.10 13.50 1.29 5.83 11.8 528 13.70 28.6 23.1
## 5 33.7 19.80 5.19 9.73 24.80 21.10 46.1 305 4.66 25.7 81.3
## 6 18.8 3.75 17.60 22.10 10.90 11.90 25.7 610 13.70 49.1 57.8# ------------------------------------------------------------
# 3) Discretização dos dados
# ------------------------------------------------------------
# discretização de Hartemink com:
# - 3 níveis finais
# - 60 quebras iniciais
# - discretização inicial por quantis
dsachs <- discretize(
sachs,
method = "hartemink",
breaks = 3,
ibreaks = 60,
idisc = "quantile"
)
str(dsachs)## 'data.frame': 853 obs. of 11 variables:
## $ Raf : Factor w/ 3 levels "[1.61,39.52]",..: 1 1 2 2 1 1 2 2 3 1 ...
## $ Mek : Factor w/ 3 levels "[1,15.1]","(15.1,24.8]",..: 1 2 3 3 2 1 3 1 3 2 ...
## $ Plcg: Factor w/ 3 levels "[1,10.8]","(10.8,24.8]",..: 1 2 2 2 1 2 1 2 1 1 ...
## $ PIP2: Factor w/ 3 levels "[1.11,23.3]",..: 1 1 1 1 1 1 2 2 1 3 ...
## $ PIP3: Factor w/ 3 levels "[1,18.88]","(18.88,45.3]",..: 3 1 1 1 2 1 1 1 2 2 ...
## $ Erk : Factor w/ 3 levels "[1,15.32]","(15.32,25.3]",..: 1 2 1 1 2 1 1 2 1 1 ...
## $ Akt : Factor w/ 3 levels "[1.7,23.5]","(23.5,39.6]",..: 1 2 2 1 3 2 1 3 1 3 ...
## $ PKA : Factor w/ 3 levels "[1.95,192.2]",..: 2 2 2 2 2 3 3 2 2 2 ...
## $ PKC : Factor w/ 3 levels "[1,8.06]","(8.06,16.58]",..: 3 1 2 2 1 2 2 1 3 2 ...
## $ P38 : Factor w/ 3 levels "[1.53,25]","(25,44.34]",..: 3 1 2 2 2 3 2 1 3 2 ...
## $ Jnk : Factor w/ 3 levels "[1,9.82]","(9.82,35.9]",..: 3 3 2 2 3 3 2 2 2 3 ...## Raf Mek Plcg PIP2 PIP3 Erk
## 1 [1.61,39.52] [1,15.1] [1,10.8] [1.11,23.3] (45.3,764] [1,15.32]
## 2 [1.61,39.52] (15.1,24.8] (10.8,24.8] [1.11,23.3] [1,18.88] (15.32,25.3]
## 3 (39.52,74.3] (24.8,389] (10.8,24.8] [1.11,23.3] [1,18.88] [1,15.32]
## 4 (39.52,74.3] (24.8,389] (10.8,24.8] [1.11,23.3] [1,18.88] [1,15.32]
## 5 [1.61,39.52] (15.1,24.8] [1,10.8] [1.11,23.3] (18.88,45.3] (15.32,25.3]
## 6 [1.61,39.52] [1,15.1] (10.8,24.8] [1.11,23.3] [1,18.88] [1,15.32]
## Akt PKA PKC P38 Jnk
## 1 [1.7,23.5] (192.2,547] (16.58,106] (44.34,170] (35.9,343]
## 2 (23.5,39.6] (192.2,547] [1,8.06] [1.53,25] (35.9,343]
## 3 (23.5,39.6] (192.2,547] (8.06,16.58] (25,44.34] (9.82,35.9]
## 4 [1.7,23.5] (192.2,547] (8.06,16.58] (25,44.34] (9.82,35.9]
## 5 (39.6,3555] (192.2,547] [1,8.06] (25,44.34] (35.9,343]
## 6 (23.5,39.6] (547,4491] (8.06,16.58] (44.34,170] (35.9,343]# ------------------------------------------------------------
# 4) Model averaging com bootstrap
#
# ------------------------------------------------------------
# Vamos aprender várias redes por bootstrap usando hill-climbing
# e score BDe.
set.seed(123)
boot <- boot.strength(
data = dsachs,
R = 500,
algorithm = "hc",
algorithm.args = list(score = "bde", iss = 10)
)
# Visualizando as arestas mais fortes
subset(boot, strength > 0.85 & direction > 0.5)## from to strength direction
## 11 Mek Raf 1.000 0.5150000
## 23 Plcg PIP2 1.000 0.5080000
## 24 Plcg PIP3 0.998 0.5240481
## 34 PIP2 PIP3 1.000 0.5090000
## 57 Erk PKA 0.990 0.5262626
## 67 Akt PKA 1.000 0.5330000
## 89 PKC P38 1.000 0.5090000
## 100 P38 Jnk 1.000 0.5030000
## 109 Jnk PKC 1.000 0.5380000# ------------------------------------------------------------
# 5) Construindo a rede média
# ------------------------------------------------------------
avg.boot <- averaged.network(boot, threshold = 0.85)
avg.boot##
## Random/Generated Bayesian network
##
## model:
## [partially directed graph]
## nodes: 11
## arcs: 9
## undirected arcs: 1
## directed arcs: 8
## average markov blanket size: 1.64
## average neighbourhood size: 1.64
## average branching factor: 0.73
##
## generation algorithm: Model Averaging
## significance threshold: 0.85## [1] "[Mek][Plcg][Akt][PKC][Raf|Mek][PIP2|Plcg][Erk|Akt][P38|PKC][PIP3|Plcg:PIP2][PKA|Erk:Akt][Jnk|P38]"#modelstring(avg.boot)
# Plot
library(Rgraphviz)
graphviz.plot(avg.boot, main = "Averaged network (bootstrap, threshold = 0.85)")
# ------------------------------------------------------------
# 6) Comparando com uma busca hill-climbing simples
# ------------------------------------------------------------
# Isso mostra a diferença entre:
# - uma única busca
# - uma rede média baseada em bootstrap
bn.hc <- hc(dsachs, score = "bde", iss = 10)
bn.hc##
## Bayesian network learned via Score-based methods
##
## model:
## [Raf][Erk][PKC][Mek|Raf][Akt|Erk][P38|PKC][Plcg|Akt][PKA|Erk:Akt]
## [Jnk|PKC:P38][PIP2|Plcg][PIP3|Plcg:PIP2]
## nodes: 11
## arcs: 11
## undirected arcs: 0
## directed arcs: 11
## average markov blanket size: 2.00
## average neighbourhood size: 2.00
## average branching factor: 1.00
##
## learning algorithm: Hill-Climbing
## score: Bayesian Dirichlet (BDe)
## graph prior: Uniform
## imaginary sample size: 10
## tests used in the learning procedure: 165
## optimized: TRUE## [1] "[Raf][Erk][PKC][Mek|Raf][Akt|Erk][P38|PKC][Plcg|Akt][PKA|Erk:Akt][Jnk|PKC:P38][PIP2|Plcg][PIP3|Plcg:PIP2]"
# ------------------------------------------------------------
# 7) Comparando classes de equivalência
# ------------------------------------------------------------
# Duas redes podem ter arcos orientados de forma diferente, mas ainda
# pertencer à mesma classe de equivalência de Markov.
all.equal(cpdag(avg.boot), cpdag(bn.hc))## [1] "Different number of directed/undirected arcs"# obter um DAG completo a partir da rede média:
avg.boot.dag <- cextend(cpdag(avg.boot))
avg.boot.dag##
## Random/Generated Bayesian network
##
## model:
## [Mek][PIP3][PKA][Jnk][Raf|Mek][PIP2|PIP3][Akt|PKA][P38|Jnk][Plcg|PIP2:PIP3]
## [Erk|Akt:PKA][PKC|P38]
## nodes: 11
## arcs: 9
## undirected arcs: 0
## directed arcs: 9
## average markov blanket size: 1.64
## average neighbourhood size: 1.64
## average branching factor: 0.82
##
## generation algorithm: Model Averaging
## significance threshold: 0.85## [1] "[Mek][PIP3][PKA][Jnk][Raf|Mek][PIP2|PIP3][Akt|PKA][P38|Jnk][Plcg|PIP2:PIP3][Erk|Akt:PKA][PKC|P38]"## [1] -8513.144# ------------------------------------------------------------
# 8) Escolha do limiar de significância
#
# ------------------------------------------------------------
# a rede média pode ser pouco sensível a alguns
# valores de threshold. Vamos reproduzir isso.
avg.050 <- averaged.network(boot, threshold = 0.50)
avg.070 <- averaged.network(boot, threshold = 0.70)
avg.def <- averaged.network(boot) # threshold estimado automaticamente
avg.050.dag <- cextend(avg.050)
avg.070.dag <- cextend(avg.070)
avg.def.dag <- cextend(avg.def)
modelstring(avg.050.dag)## [1] "[Mek][Plcg][PKC][Raf|Mek][PIP2|Plcg][Akt|Plcg][P38|PKC][PIP3|Plcg:PIP2][Erk|Akt][Jnk|P38][PKA|Erk:Akt]"## [1] "[Mek][Plcg][Akt][PKC][Raf|Mek][PIP2|Plcg][Erk|Akt][P38|PKC][PIP3|Plcg:PIP2][PKA|Erk:Akt][Jnk|P38]"## [1] "[Mek][Plcg][Akt][PKC][Raf|Mek][PIP2|Plcg][Erk|Akt][P38|PKC][PIP3|Plcg:PIP2][PKA|Erk:Akt][Jnk|P38]"## [1] "Different number of directed/undirected arcs"## [1] TRUE
Interpretação dos Resultados da Aplicação em Perfis de Expressão Gênica
Os resultados obtidos nesta aplicação demonstram como redes Bayesianas podem ser utilizadas para inferir relações estruturais entre variáveis biológicas a partir de dados de sinalização celular e expressão gênica.
Após a leitura do conjunto de dados de Sachs et al. (2005), observou-se a presença de múltiplas variáveis contínuas representando níveis de ativação de diferentes proteínas envolvidas em vias de sinalização intracelular. Como alguns algoritmos de aprendizado estrutural assumem variáveis discretas, realizou-se inicialmente a discretização dos dados por meio do método de Hartemink, agrupando os valores contínuos em três categorias ordenadas. Essa etapa preserva relações estatísticas importantes ao mesmo tempo em que adapta os dados ao processo de inferência estrutural discreta.
Em seguida, foi aplicada a técnica de bootstrap combinada com o algoritmo Hill-Climbing utilizando o score BDe, permitindo gerar múltiplas redes Bayesianas a partir de diferentes reamostragens dos dados. O objetivo dessa abordagem é identificar conexões estruturais consistentes ao longo de várias execuções, reduzindo a influência de variações aleatórias específicas da amostra original.
A análise das forças das arestas (boot.strength) permitiu verificar quais relações apresentaram maior estabilidade estatística. Por exemplo, conexões como Mek–Raf, Plcg–PIP2, PIP2–PIP3, PKC–P38 e P38–Jnk apresentaram valores de força próximos de 1, indicando que essas associações foram identificadas de forma extremamente consistente nas redes aprendidas ao longo dos reamostramentos. Isso sugere forte evidência estatística de dependência estrutural entre essas variáveis.
A construção da rede média (averaged network) utilizando threshold de 0.85 resultou em uma estrutura consenso contendo apenas as arestas mais robustas. Essa rede representa uma estimativa mais conservadora e estável da estrutura causal/probabilística subjacente, removendo conexões menos frequentes e potencialmente espúrias.
Quando comparada com a rede obtida por uma única execução do Hill-Climbing tradicional, observou-se que a rede média apresentou menor número de arcos, indicando que o processo de model averaging atua como mecanismo de regularização estrutural, priorizando robustez em detrimento de complexidade. Isso evidencia que nem toda relação encontrada em uma única busca é necessariamente estável ou reproduzível sob pequenas perturbações dos dados.
Além disso, a análise de diferentes limiares de significância mostrou que a estrutura média permaneceu praticamente inalterada entre alguns thresholds testados, sugerindo estabilidade estrutural da rede inferida e robustez na escolha do parâmetro de corte.
De forma geral, os resultados sugerem que as redes Bayesianas foram capazes de recuperar padrões relevantes de associação entre proteínas envolvidas nas vias de sinalização celular, produzindo estruturas compatíveis com mecanismos biológicos conhecidos e demonstrando o potencial dessa abordagem como ferramenta para modelagem de sistemas biológicos complexos, inferência de redes regulatórias e geração de hipóteses em bioinformática e biologia de sistemas.
Referências bibliográficas do módulo 05
Bøttcher, S. G., & Dethlefsen, C. (2003). deal: a package for learning Bayesian networks. Journal of Statistical Software, 8(20), 1–40.
Bouckaert, R. R. (1995). Bayesian belief networks: from construction to inference (Tese de Doutorado). Utrecht University, The Netherlands.
Castillo, E., Gutiérrez, J. M., & Hadi, A. S. (1997). Expert systems and probabilistic network models. New York: Springer.
Chickering, D. M. (1995). A transformational characterization of equivalent Bayesian network structures. In Proceedings of the 11th Conference on Uncertainty in Artificial Intelligence (UAI95) (pp. 87–98).
Cover, T. M., & Thomas, J. A. (2006). Elements of information theory (2nd ed.). Hoboken: Wiley.
Edwards, D. I. (2000). Introduction to graphical modelling (2nd ed.). New York: Springer.
Geiger, D., & Heckerman, D. (1994). Learning Gaussian networks (Technical Report). Microsoft Research.
Hausser, J., & Strimmer, K. (2009). Entropy inference and the James-Stein estimator, with application to nonlinear gene association networks. Journal of Machine Learning Research, 10, 1469–1484.
Heckerman, D., Geiger, D., & Chickering, D. M. (1995). Learning Bayesian networks: the combination of knowledge and statistical data. Machine Learning, 20(3), 197–243.
Jensen, F. V. (2001). Bayesian networks and decision graphs. New York: Springer.
Koller, D., & Friedman, N. (2009). Probabilistic graphical models: principles and techniques. Cambridge: MIT Press.
Korb, K. B., & Nicholson, A. E. (2010). Bayesian artificial intelligence (2nd ed.). Boca Raton: Chapman and Hall.
Kullback, S. (1968). Information theory and statistic. New York: Dover Publications.
Larranaga, P., Sierra, B., Gallego, M. J., Michelena, M. J., & Picaza, J. M. (1997). Learning Bayesian networks by genetic algorithms: a case study in the prediction of survival in malignant skin melanoma. In Proceedings of the 6th Conference on Artificial Intelligence in Medicine in Europe (AIME ’97) (pp. 261–272). Springer.
Legendre, P. (2000). Comparison of permutation methods for the partial correlation and partial mantel tests. Journal of Statistical Computation and Simulation, 67, 37–73.
Margaritis, D. (2003). Learning Bayesian network model structure from data (Tese de Doutorado). Carnegie Mellon University, Pittsburgh, PA.
Neapolitan, R. E. (2003). Learning Bayesian networks. Englewood Cliffs: Prentice Hall.
Pearl, J. (1988). Probabilistic reasoning in intelligent systems: networks of plausible inference.
Rissanen, J. (2007). Information and complexity in statistical models. New York: Springer.
Sachs, K., Perez, O., Pe’er, D., Lauffenburger, D. A., & Nolan, G. P. (2005). Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721), 523–529.
Scutari, M., & Brogini, A. (2012). Bayesian network structure learning with permutation tests. Communications in Statistics - Theory and Methods, 41(16–17), 3233–3243.
Shafer, J., & Strimmer, K. (2005). A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics. Statistical Applications in Genetics and Molecular Biology, 4, 32.
Spirtes, P., Glymour, C., & Scheines, R. (2001). Causation, prediction, and search (2nd ed.). Cambridge: MIT Press.
Tsamardinos, I., Brown, L. E., & Aliferis, C. F. (2006). The max-min hill-climbing Bayesian network structure learning algorithm. Machine Learning, 65(1), 31–78.
Verma, T. S., & Pearl, J. (1991). Equivalence and synthesis of causal models. Uncertainty in Artificial Intelligence, 6, 255–268.
Whittaker, J. (1990). Graphical models in applied multivariate statistics. Wiley.
Yaramakala, S., & Margaritis, D. (2005). Speculative Markov blanket discovery for optimal feature selection. In Proceedings of the 5th IEEE International Conference on Data Mining (ICDM ’05) (pp. 809–812). IEEE Computer Society.