Estatística Avançada - Trabalho Final
Emanuelly e Nicolas
2026-07-01
Objetivo
Entender fatores influenciando no aumento do tamanho do genoma.
Carregamento de pacotes
library(ggplot2)
library(corrplot)
library(dplyr)Carregamento dos dados
glm_pangenome <- read.csv("~/dev/anvio-summary/glm_pangenome.csv")
colnames(glm_pangenome)## [1] "genome_name"
## [2] "total_length"
## [3] "gc_content"
## [4] "percent_redundancy"
## [5] "num_genes"
## [6] "avg_gene_length"
## [7] "num_genes_per_kb"
## [8] "singleton_gene_clusters"
## [9] "num_gene_clusters"
## [10] "As"
## [11] "Cd"
## [12] "Co"
## [13] "Cr"
## [14] "Cu"
## [15] "Fe"
## [16] "Metal_genes"
## [17] "Mg"
## [18] "Mn"
## [19] "Mo"
## [20] "Ni"
## [21] "Skip"
## [22] "Te"
## [23] "Zn"
## [24] "Amino_acid_transport_and_metabolism"
## [25] "Carbohydrate_transport_and_metabolism"
## [26] "Cell_cycle_control_cell_division_chromosome_partitioning"
## [27] "Cell_motility"
## [28] "Cell_wall_membrane_envelope_biogenesis"
## [29] "Chromatin_structure_and_dynamics"
## [30] "Coenzyme_transport_and_metabolism"
## [31] "Cytoskeleton"
## [32] "Defense_mechanisms"
## [33] "Energy_production_and_conversion"
## [34] "Extracellular_structures"
## [35] "Function_unknown"
## [36] "General_function_prediction_only"
## [37] "Inorganic_ion_transport_and_metabolism"
## [38] "Intracellular_trafficking_secretion_and_vesicular_transport"
## [39] "Lipid_transport_and_metabolism"
## [40] "Mobilome_prophages_transposons"
## [41] "Nucleotide_transport_and_metabolism"
## [42] "Posttranslational_modification_protein_turnover_chaperones"
## [43] "RNA_processing_and_modification"
## [44] "Replication_recombination_and_repair"
## [45] "Secondary_metabolites_biosynthesis_transport_and_catabolism"
## [46] "Signal_transduction_mechanisms"
## [47] "Transcription"
## [48] "Translation_ribosomal_structure_and_biogenesis"
## [49] "HMA_genes"Gráfico de Correlação
You can also embed plots, for example:
glm_pangenome_filt <- glm_pangenome %>%
select(total_length, gc_content,
Metal_genes, HMA_genes,
Defense_mechanisms,
Mobilome_prophages_transposons,
Replication_recombination_and_repair)
# Compute correlation matrix
cor_mat <- cor(glm_pangenome_filt, use = "pairwise.complete.obs")
short_names <- c("Size", "GC%", "Metal Res.", "HMA (CxxC)", "Defense", "Mobilome", "Rep.")
colnames(cor_mat) <- short_names
rownames(cor_mat) <- short_names
# Plot
corrplot(cor_mat, diag = FALSE)
cor_mat## Size GC% Metal Res. HMA (CxxC) Defense Mobilome
## Size 1.0000000 -0.7747406 0.6987942 0.8538561 0.8633926 0.7992820
## GC% -0.7747406 1.0000000 -0.4424625 -0.5904111 -0.7504806 -0.6742407
## Metal Res. 0.6987942 -0.4424625 1.0000000 0.5674816 0.6959402 0.5847626
## HMA (CxxC) 0.8538561 -0.5904111 0.5674816 1.0000000 0.7523670 0.6818718
## Defense 0.8633926 -0.7504806 0.6959402 0.7523670 1.0000000 0.8026213
## Mobilome 0.7992820 -0.6742407 0.5847626 0.6818718 0.8026213 1.0000000
## Rep. 0.6840843 -0.7230336 0.4908238 0.6353768 0.6699162 0.7281611
## Rep.
## Size 0.6840843
## GC% -0.7230336
## Metal Res. 0.4908238
## HMA (CxxC) 0.6353768
## Defense 0.6699162
## Mobilome 0.7281611
## Rep. 1.0000000Modelagem
\(Y\): Tamanho do genoma
ggplot(glm_pangenome, aes(x = total_length)) +
geom_histogram(bins=10, fill = "skyblue")
mod_1 <- glm(total_length ~ 1, family = poisson, data = glm_pangenome)mod_metres <- glm(total_length ~ Metal_genes, family = poisson, data = glm_pangenome)
summary(mod_metres)##
## Call:
## glm(formula = total_length ~ Metal_genes, family = poisson, data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.476e+01 2.356e-03 6263.7 <2e-16 ***
## Metal_genes 3.536e-03 1.107e-05 319.4 <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: 207370 on 32 degrees of freedom
## Residual deviance: 105675 on 31 degrees of freedom
## AIC: 106251
##
## Number of Fisher Scoring iterations: 3mod_hma <- glm(total_length ~ HMA_genes, family = poisson, data = glm_pangenome)
summary(mod_hma)##
## Call:
## glm(formula = total_length ~ HMA_genes, family = poisson, data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.477e+01 1.899e-03 7776.7 <2e-16 ***
## HMA_genes 1.971e-03 5.049e-06 390.3 <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: 207370 on 32 degrees of freedom
## Residual deviance: 55685 on 31 degrees of freedom
## AIC: 56261
##
## Number of Fisher Scoring iterations: 3mod_gc <- glm(total_length ~ gc_content, family = poisson, data = glm_pangenome)
summary(mod_gc)##
## Call:
## glm(formula = total_length ~ gc_content, family = poisson, data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 26.84964 0.03202 838.4 <2e-16 ***
## gc_content -32.22405 0.09099 -354.2 <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: 207370 on 32 degrees of freedom
## Residual deviance: 82949 on 31 degrees of freedom
## AIC: 83526
##
## Number of Fisher Scoring iterations: 3mod_mobilome <- glm(total_length ~ Mobilome_prophages_transposons, family = poisson, data = glm_pangenome)
summary(mod_mobilome)##
## Call:
## glm(formula = total_length ~ Mobilome_prophages_transposons,
## family = poisson, data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.541e+01 2.721e-04 56640.2 <2e-16 ***
## Mobilome_prophages_transposons 1.478e-03 4.045e-06 365.4 <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: 207370 on 32 degrees of freedom
## Residual deviance: 74233 on 31 degrees of freedom
## AIC: 74810
##
## Number of Fisher Scoring iterations: 3mod_rep <- glm(total_length ~ Replication_recombination_and_repair, family = poisson, data = glm_pangenome)
summary(mod_rep)##
## Call:
## glm(formula = total_length ~ Replication_recombination_and_repair,
## family = poisson, data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.517e+01 1.069e-03 14194.8 <2e-16 ***
## Replication_recombination_and_repair 1.978e-03 6.324e-06 312.7 <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: 207370 on 32 degrees of freedom
## Residual deviance: 110548 on 31 degrees of freedom
## AIC: 111124
##
## Number of Fisher Scoring iterations: 3mod_def <- glm(total_length ~ Defense_mechanisms, family = poisson, data = glm_pangenome)
summary(mod_def)##
## Call:
## glm(formula = total_length ~ Defense_mechanisms, family = poisson,
## data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.511e+01 1.002e-03 15083.3 <2e-16 ***
## Defense_mechanisms 2.380e-03 6.031e-06 394.7 <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: 207370 on 32 degrees of freedom
## Residual deviance: 52654 on 31 degrees of freedom
## AIC: 53230
##
## Number of Fisher Scoring iterations: 3AIC(mod_1,
mod_metres,
mod_hma,
mod_gc,
mod_mobilome,
mod_rep,
mod_def)## df AIC
## mod_1 1 207944.61
## mod_metres 2 106251.28
## mod_hma 2 56260.97
## mod_gc 2 83525.51
## mod_mobilome 2 74809.86
## mod_rep 2 111124.15
## mod_def 2 53230.15glm_pangenome$y_hat <- predict(mod_def, newdata = glm_pangenome, type = "response")
ggplot(glm_pangenome, aes(x = Defense_mechanisms, y = total_length)) +
geom_point() +
geom_line(aes(y = y_hat), color = "red", size = 1)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
glm_pangenome$y_hat <- predict(mod_hma, newdata = glm_pangenome, type = "response")
ggplot(glm_pangenome, aes(x = HMA_genes, y = total_length)) +
geom_point() +
geom_line(aes(y = y_hat), color = "red", size = 1)
glm_pangenome$y_hat <- predict(mod_mobilome, newdata = glm_pangenome, type = "response")
ggplot(glm_pangenome, aes(x = Mobilome_prophages_transposons, y = total_length)) +
geom_point() +
geom_line(aes(y = y_hat), color = "red", size = 1)
glm_pangenome$y_hat <- predict(mod_metres, newdata = glm_pangenome, type = "response")
ggplot(glm_pangenome, aes(x = Metal_genes, y = total_length)) +
geom_point() +
geom_line(aes(y = y_hat), color = "red", size = 1)
Usando defesa como modelo base
Defesa + Mobiloma
mod_def_mob <- glm(total_length ~ Defense_mechanisms + Mobilome_prophages_transposons, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_mob, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + Mobilome_prophages_transposons
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2.2e-16 ***
## 3 30 45863 1 6791 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Defesa + GC
mod_def_gc <- glm(total_length ~ Defense_mechanisms + gc_content, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_gc, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + gc_content
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2.2e-16 ***
## 3 30 45402 1 7251 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Defesa + MetRes
mod_def_metres <- glm(total_length ~ Defense_mechanisms + Metal_genes, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_metres, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + Metal_genes
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2.2e-16 ***
## 3 30 48634 1 4020 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Defesa + HMA
mod_def_hma <- glm(total_length ~ Defense_mechanisms + HMA_genes, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_hma, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + HMA_genes
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2.2e-16 ***
## 3 30 32538 1 20116 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Defesa + Rep
mod_def_rep <- glm(total_length ~ Defense_mechanisms + Replication_recombination_and_repair, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_rep, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + Replication_recombination_and_repair
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2.2e-16 ***
## 3 30 48448 1 4205 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Usando mobiloma como modelo base
Construindo novos modelos com mais de uma var explicativa, utilizando mobiloma como ponto de partida
Mob + GC
# Mobiloma + gc
mod_mob_gc <- glm(total_length ~ Mobilome_prophages_transposons + gc_content, family = poisson, data = glm_pangenome)
anova(mod_1, mod_mobilome, mod_mob_gc, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Mobilome_prophages_transposons
## Model 3: total_length ~ Mobilome_prophages_transposons + gc_content
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 74233 1 133137 < 2.2e-16 ***
## 3 30 53663 1 20570 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mob + MetRes
# Mobiloma + met res
mod_mob_metres <- glm(total_length ~ Mobilome_prophages_transposons + Metal_genes, family = poisson, data = glm_pangenome)
anova(mod_1, mod_mobilome, mod_mob_metres, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Mobilome_prophages_transposons
## Model 3: total_length ~ Mobilome_prophages_transposons + Metal_genes
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 74233 1 133137 < 2.2e-16 ***
## 3 30 57164 1 17070 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mob + HMA
# Mobiloma + hma
mod_mob_hma <- glm(total_length ~ Mobilome_prophages_transposons + HMA_genes, family = poisson, data = glm_pangenome)
anova(mod_1, mod_mobilome, mod_mob_hma, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Mobilome_prophages_transposons
## Model 3: total_length ~ Mobilome_prophages_transposons + HMA_genes
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 74233 1 133137 < 2.2e-16 ***
## 3 30 37323 1 36910 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mob + Rep
# Mobiloma + rep
mod_mob_rep <- glm(total_length ~ Mobilome_prophages_transposons + Replication_recombination_and_repair, family = poisson, data = glm_pangenome)
anova(mod_1, mod_mobilome, mod_mob_rep, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Mobilome_prophages_transposons
## Model 3: total_length ~ Mobilome_prophages_transposons + Replication_recombination_and_repair
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 74233 1 133137 < 2.2e-16 ***
## 3 30 69762 1 4471 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No modelo aditivo sobre mob, hma explicou melhor que as outras variáveis:
Interações
# Investigar interação
# Melhorou, mas não muito
mod_mobXhma <- glm(total_length ~ Mobilome_prophages_transposons * HMA_genes, family = poisson, data = glm_pangenome)
anova(mod_1, mod_mobilome, mod_mob_hma, mod_mobXhma, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Mobilome_prophages_transposons
## Model 3: total_length ~ Mobilome_prophages_transposons + HMA_genes
## Model 4: total_length ~ Mobilome_prophages_transposons * HMA_genes
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 74233 1 133137 < 2.2e-16 ***
## 3 30 37323 1 36910 < 2.2e-16 ***
## 4 29 37231 1 92 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_defXhma <- glm(total_length ~ Defense_mechanisms * HMA_genes, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_hma, mod_defXhma, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + HMA_genes
## Model 4: total_length ~ Defense_mechanisms * HMA_genes
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2e-16 ***
## 3 30 32538 1 20116 < 2e-16 ***
## 4 29 32534 1 4 0.04516 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_defXmob <- glm(total_length ~ Defense_mechanisms * Mobilome_prophages_transposons, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_mob, mod_defXmob, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + Mobilome_prophages_transposons
## Model 4: total_length ~ Defense_mechanisms * Mobilome_prophages_transposons
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2.2e-16 ***
## 3 30 45863 1 6791 < 2.2e-16 ***
## 4 29 44733 1 1130 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_completo <- glm(total_length ~ Defense_mechanisms + Mobilome_prophages_transposons + HMA_genes, family = poisson, data = glm_pangenome)
anova(mod_1, mod_def, mod_def_mob, mod_completo, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: total_length ~ 1
## Model 2: total_length ~ Defense_mechanisms
## Model 3: total_length ~ Defense_mechanisms + Mobilome_prophages_transposons
## Model 4: total_length ~ Defense_mechanisms + Mobilome_prophages_transposons +
## HMA_genes
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 207370
## 2 31 52654 1 154716 < 2.2e-16 ***
## 3 30 45863 1 6791 < 2.2e-16 ***
## 4 29 29465 1 16398 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod_completo)##
## Call:
## glm(formula = total_length ~ Defense_mechanisms + Mobilome_prophages_transposons +
## HMA_genes, family = poisson, data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.494e+01 2.384e-03 6266.81 <2e-16 ***
## Defense_mechanisms 1.016e-03 1.146e-05 88.61 <2e-16 ***
## Mobilome_prophages_transposons 3.839e-04 6.926e-06 55.43 <2e-16 ***
## HMA_genes 1.004e-03 7.836e-06 128.12 <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: 207370 on 32 degrees of freedom
## Residual deviance: 29465 on 29 degrees of freedom
## AIC: 30046
##
## Number of Fisher Scoring iterations: 3Validação
anova(mod_def_hma, test = "Chisq")## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: total_length
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 32 207370
## Defense_mechanisms 1 154716 31 52654 < 2.2e-16 ***
## HMA_genes 1 20116 30 32538 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod_def_hma)##
## Call:
## glm(formula = total_length ~ Defense_mechanisms + HMA_genes,
## family = poisson, data = glm_pangenome)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.487e+01 2.009e-03 7399.6 <2e-16 ***
## Defense_mechanisms 1.397e-03 9.183e-06 152.1 <2e-16 ***
## HMA_genes 1.092e-03 7.686e-06 142.0 <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: 207370 on 32 degrees of freedom
## Residual deviance: 32538 on 30 degrees of freedom
## AIC: 33116
##
## Number of Fisher Scoring iterations: 3VALORES AJUSTADOS E RESÍDUOS
ajustados <- fitted(mod_def_hma)
eta_hat <- predict(mod_def_hma, type = "link")
res_bruto <- residuals(mod_def_hma, type = "response")
res_pearson <- residuals(mod_def_hma, type = "pearson")
res_desvio <- residuals(mod_def_hma, type = "deviance")
res_working <- residuals(mod_def_hma, type = "working")
# Resíduos padronizados / studentizados
res_std_pearson <- rstandard(mod_def_hma, type = "pearson")
res_std_dev <- rstandard(mod_def_hma, type = "deviance")
res_stud <- rstudent(mod_def_hma)
# Medidas de influência
lev <- hatvalues(mod_def_hma)
cook <- cooks.distance(mod_def_hma)
dff <- dffits(mod_def_hma)
# Tabela-resumo
diagnostico <- data.frame(
id = glm_pangenome$genome_name,
y = glm_pangenome$total_length,
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
## 1 B_paramycoides_BD_11 5483749 5375280 108468.841 46.784765 46.628726
## 2 B_paramycoides_BD_12 5774287 5714329 59958.270 25.082250 25.038577
## 3 B_paramycoides_BD_13 5627347 5778538 -151190.862 -62.895064 -63.172364
## 4 B_paramycoides_BD_15 5363035 5295194 67841.175 29.481694 29.419075
## 5 B_paramycoides_BD_31 5364188 5353060 11127.943 4.809654 4.807989
## 6 B_paramycoides_BD_331 5199139 5219656 -20516.802 -8.980256 -8.986149
## 7 B_paramycoides_BD_333 5182738 5240899 -58161.051 -25.405594 -25.452802
## 8 B_paramycoides_BD_335 5363943 5269584 94359.186 41.105130 40.983362
## 9 B_paramycoides_BD_336 5294777 5288483 6293.933 2.736885 2.736343
## 10 B_paramycoides_BD_9 5299588 5388673 -89084.648 -38.376200 -38.482673
## res_std_pearson res_std_dev res_stud leverage cook dffits
## 1 48.484620 48.322913 48.334070 0.06889019 57.9753271 0.40729360
## 2 26.625826 26.579466 26.584690 0.11258497 29.9804801 0.28575562
## 3 -69.391504 -69.697447 -69.642942 0.17847549 348.6979609 -1.05150529
## 4 31.610274 31.543134 31.551880 0.13014204 49.8315661 0.36994593
## 5 4.969755 4.968035 4.968144 0.06339234 0.5572215 0.03860059
## 6 -9.372112 -9.378262 -9.377759 0.08187361 2.6109298 -0.08372128
## 7 -26.360732 -26.409715 -26.406233 0.07115386 17.7438706 -0.22059792
## 8 42.727130 42.600558 42.609998 0.07448257 48.9729477 0.37129377
## 9 2.841120 2.840557 2.840598 0.07203011 0.2088518 0.02362939
## 10 -39.766083 -39.876412 -39.868845 0.06868134 38.8727137 -0.33149077IDENTIFICAÇÃO DE PONTOS SUSPEITOS
# Regras práticas
n = 33
limite_res <- 2
limite_lev <- 2 * length(coef(mod_def_hma)) / n
limite_cook <- 4 / n
limite_dff <- 2 * sqrt(length(coef(mod_def_hma)) / n)
cat("\nLimites práticos:\n")##
## Limites práticos:cat(" |resíduo studentizado| >", limite_res, "\n")## |resíduo studentizado| > 2cat(" leverage >", round(limite_lev, 4), "\n")## leverage > 0.1818cat(" Cook >", round(limite_cook, 4), "\n")## Cook > 0.1212cat(" |DFFITS| >", round(limite_dff, 4), "\n\n")## |DFFITS| > 0.603suspeitos <- 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
## 18 B_paramycoides_DE0103 5782815 5549336 233478.747 99.112136
## 3 B_paramycoides_BD_13 5627347 5778538 -151190.862 -62.895064
## 1 B_paramycoides_BD_11 5483749 5375280 108468.841 46.784765
## 8 B_paramycoides_BD_335 5363943 5269584 94359.186 41.105130
## 21 B_paramycoides_NPDC077662 5450429 5537051 -86622.493 -36.812168
## 10 B_paramycoides_BD_9 5299588 5388673 -89084.648 -38.376200
## 11 B_paramycoides_B_14577 5231975 5308815 -76840.330 -33.349583
## 12 B_paramycoides_B_14578 5227308 5303023 -75715.040 -32.879136
## 13 B_paramycoides_B_14580 5227789 5303023 -75234.040 -32.670263
## 22 B_paramycoides_NPDC093745 5897029 5829974 67054.601 27.771239
## 4 B_paramycoides_BD_15 5363035 5295194 67841.175 29.481694
## 29 B_paramycoides_NPDC094358 5285347 5218064 67283.295 29.454561
## 19 B_paramycoides_LB_RP2 5342545 5411559 -69013.654 -29.667014
## 17 B_paramycoides_CP40 5619544 5685960 -66416.339 -27.853067
## 33 B_paramycoides_h2_1 5718601 5653284 65316.840 27.471018
## 2 B_paramycoides_BD_12 5774287 5714329 59958.270 25.082250
## 7 B_paramycoides_BD_333 5182738 5240899 -58161.051 -25.405594
## 32 B_paramycoides_RZ3MS14 5493110 5438558 54552.000 23.392078
## 25 B_paramycoides_NPDC094098 5286979 5336930 -49950.935 -21.622100
## 16 B_paramycoides_CIP111461T 5436504 5486000 -49496.083 -21.132113
## 20 B_paramycoides_NH24A2 5453107 5492961 -39854.478 -17.004891
## 14 B_paramycoides_B_3435 5660034 5630908 29125.555 12.273968
## 28 B_paramycoides_NPDC094345 5211806 5180835 30970.849 13.606713
## 23 B_paramycoides_NPDC094078 5387444 5413210 -25765.788 -11.074293
## 31 B_paramycoides_NRS_1309 5314253 5290097 24156.351 10.502666
## 24 B_paramycoides_NPDC094088 5409897 5386336 23560.718 10.151770
## 26 B_paramycoides_NPDC094118 5482343 5460692 21650.850 9.265121
## 6 B_paramycoides_BD_331 5199139 5219656 -20516.802 -8.980256
## 30 B_paramycoides_NPDC094368 5585329 5600977 -15648.053 -6.611933
## 15 B_paramycoides_B_615 5560887 5574605 -13717.970 -5.810089
## 5 B_paramycoides_BD_31 5364188 5353060 11127.943 4.809654
## 9 B_paramycoides_BD_336 5294777 5288483 6293.933 2.736885
## res_desvio res_std_pearson res_std_dev res_stud leverage cook
## 18 98.429048 107.371241 106.631230 106.741020 0.14792519 667.1432428
## 3 -63.172364 -69.391504 -69.697447 -69.642942 0.17847549 348.6979609
## 1 46.628726 48.484620 48.322913 48.334070 0.06889019 57.9753271
## 8 40.983362 42.727130 42.600558 42.609998 0.07448257 48.9729477
## 21 -36.908782 -41.155949 -41.263963 -41.242388 0.19994922 141.1062170
## 10 -38.482673 -39.766083 -39.876412 -39.868845 0.06868134 38.8727137
## 11 -33.430523 -34.690258 -34.774452 -34.768078 0.07580046 32.9003018
## 12 -32.957845 -34.340079 -34.422286 -34.415447 0.08327682 35.7081444
## 13 -32.747972 -34.121925 -34.203087 -34.196335 0.08327682 35.2558945
## 22 27.718257 31.624520 31.564186 31.578003 0.22884331 98.9286827
## 4 29.419075 31.610274 31.543134 31.551880 0.13014204 49.8315661
## 29 29.391600 30.702060 30.636432 30.641662 0.07961382 27.1789162
## 19 -29.730409 -30.160597 -30.225046 -30.222956 0.03246246 10.1735414
## 17 -27.907557 -29.310166 -29.367506 -29.361952 0.09695479 30.7450413
## 33 27.418372 28.666589 28.611651 28.616142 0.08167272 24.3618536
## 2 25.038577 26.625826 26.579466 26.584690 0.11258497 29.9804801
## 7 -25.452802 -26.360732 -26.409715 -26.406233 0.07115386 17.7438706
## 32 23.353134 23.769495 23.729923 23.731170 0.03150431 6.1261961
## 25 -21.655961 -22.614068 -22.649482 -22.646446 0.08580601 15.9998348
## 16 -21.164010 -21.660949 -21.693644 -21.692068 0.04823248 7.9257869
## 20 -17.025517 -17.305418 -17.326408 -17.325686 0.03443050 3.5596124
## 14 12.263410 14.438788 14.426367 14.429813 0.27738246 26.6753849
## 28 13.593190 14.341441 14.327188 14.328611 0.09983759 7.6039196
## 23 -11.083096 -11.248737 -11.257678 -11.257403 0.03077520 1.3392528
## 31 10.494688 10.848284 10.840043 10.840560 0.06270348 2.6243117
## 24 10.144383 10.322369 10.314857 10.315104 0.03278099 1.2037455
## 26 9.259009 9.507427 9.501155 9.501471 0.05032244 1.5965786
## 6 -8.986149 -9.372112 -9.378262 -9.377759 0.08187361 2.6109298
## 30 -6.615015 -6.969206 -6.972455 -6.972130 0.09990110 1.7969062
## 15 -5.812474 -6.014282 -6.016751 -6.016586 0.06674988 0.8623801
## 5 4.807989 4.969755 4.968035 4.968144 0.06339234 0.5572215
## 9 2.736343 2.841120 2.840557 2.840598 0.07203011 0.2088518
## dffits
## 18 1.64448595
## 3 -1.05150529
## 1 0.40729360
## 8 0.37129377
## 21 -0.63262820
## 10 -0.33149077
## 11 -0.30299982
## 12 -0.31553162
## 13 -0.31344715
## 22 0.52137513
## 4 0.36994593
## 29 0.27296588
## 19 -0.16765355
## 17 -0.29116244
## 33 0.25800070
## 2 0.28575562
## 7 -0.22059792
## 32 0.12889272
## 25 -0.20881123
## 16 -0.14686068
## 20 -0.09813100
## 14 0.26769523
## 28 0.14289821
## 23 -0.06000528
## 31 0.08385500
## 24 0.05678427
## 26 0.06538490
## 6 -0.08372128
## 30 -0.06939935
## 15 -0.04806567
## 5 0.03860059
## 9 0.02362939Medidas globais de ajuste
cat("\n--- Medidas globais ---\n")##
## --- Medidas globais ---cat("Desvio residual:", deviance(mod_def_hma), "\n")## Desvio residual: 32537.78cat("GL residual:", df.residual(mod_def_hma), "\n")## GL residual: 30cat("p-valor (desvio):", 1 - pchisq(deviance(mod_def_hma), df.residual(mod_def_hma)), "\n")## p-valor (desvio): 0X2 <- sum(res_pearson^2)
cat("Qui-quadrado de Pearson:", X2, "\n")## Qui-quadrado de Pearson: 32636.58cat("p-valor (Pearson):", 1 - pchisq(X2, df.residual(mod_def_hma)), "\n")## p-valor (Pearson): 0Gráficos de diagnóstico
par(mfrow = c(3, 3), mar = c(4, 4, 2, 1))
# (1) Observado vs Ajustado
plot(ajustados, glm_pangenome$total_length,
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(glm_pangenome$Mobilome_prophages_transposons, 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(glm_pangenome$HMA_genes, 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)
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)
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(mod_mob_hma))
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, glm_pangenome$genome_name, ""),
pos = 3, cex = 0.8)
Teste com Redes Bayesianas
library(bnlearn) #pacote para criação das redes
library(Rgraphviz) #pacote para visualização das redesglm_pangenome_filt <- glm_pangenome %>%
select(total_length, gc_content,
Metal_genes, HMA_genes,
Defense_mechanisms,
Mobilome_prophages_transposons,
Replication_recombination_and_repair)# Teste com Redes Bayesianas
glm_pangenome_filt[] <- lapply(glm_pangenome_filt, function(x) {
if (is.integer(x))
as.numeric(x)
else
x
})
dpan <- discretize(
glm_pangenome_filt,
method = "hartemink",
breaks = 3,
ibreaks = 60,
idisc = "quantile"
)
str(dpan)## 'data.frame': 33 obs. of 7 variables:
## $ total_length : Factor w/ 3 levels "[5.18274e+06,5.27823e+06]",..: 3 3 3 2 2 1 1 2 2 2 ...
## $ gc_content : Factor w/ 3 levels "[0.349648,0.351117]",..: 2 1 1 2 2 2 3 2 3 2 ...
## $ Metal_genes : Factor w/ 3 levels "[200,209.067]",..: 2 2 2 1 2 1 2 1 1 1 ...
## $ HMA_genes : Factor w/ 3 levels "[353,365.828]",..: 1 3 3 2 2 1 1 2 2 2 ...
## $ Defense_mechanisms : Factor w/ 3 levels "[149,159.696]",..: 3 3 3 1 1 1 1 1 1 3 ...
## $ Mobilome_prophages_transposons : Factor w/ 3 levels "[36,50.3846]",..: 3 3 3 1 3 1 1 1 1 3 ...
## $ Replication_recombination_and_repair: Factor w/ 3 levels "[150,163.852]",..: 2 3 3 1 3 1 1 1 1 3 ...
## - attr(*, "cutpoints")=List of 7
## ..$ total_length : num [1:4] 5182738 5278231 5425861 5897029
## ..$ gc_content : num [1:4] 0.35 0.351 0.353 0.353
## ..$ Metal_genes : num [1:4] 200 209 216 229
## ..$ HMA_genes : num [1:4] 353 366 375 405
## ..$ Defense_mechanisms : num [1:4] 149 160 163 196
## ..$ Mobilome_prophages_transposons : num [1:4] 36 50.4 69.8 99
## ..$ Replication_recombination_and_repair: num [1:4] 150 164 174 203set.seed(123)
boot <- boot.strength(
data = dpan,
R = 100,
algorithm = "hc",
algorithm.args = list(score = "bde", iss = 10)
)## Warning in check.data(x, allow.missing = TRUE): variable gc_content in the data
## has levels that are not observed in the data.
## Warning in check.data(x, allow.missing = TRUE): variable gc_content in the data
## has levels that are not observed in the data.## Warning in check.data(x, allow.missing = TRUE): variable total_length in the
## data has levels that are not observed in the data.avg.boot <- averaged.network(boot, threshold = 0.7)graphviz.plot(avg.boot, main = "Averaged network (bootstrap, threshold = 0.7)")