library(tidyverse)
library(ggplot2)
O dataset Pima Indians Diabetes Database, disponibilizado na plataforma Kaggle, é um dos conjuntos de dados mais clássicos e utilizados no cenário de Machine Learning para tarefas de classificação. Contudo, no presente trabalho, adotamos uma abordagem estatística ao aplicar o teste inferencial de Wilcoxon não pareado para uma amostra. O objetivo desta análise é verificar quais parâmetros clínicos da população Pima saudável se distanciam dos valores de referência universalmente aceitos na literatura médica, evidenciando características basais desta população geneticamente predisposta à obesidade e resistência insulínica.
A população Pima dos Estados Unidos é amplamente reconhecida na literatura epidemiológica por apresentar uma das maiores prevalências mundiais de diabetes tipo 2 e obesidade, fenótipos atribuídos a uma combinação de fatores genéticos, históricos e ambientais. Estudos comparativos demonstram que mesmo indivíduos Pima com tolerância normal à glicose apresentam índice de massa corporal médio de aproximadamente 34,1 ± 8,0 kg/m² e insulina sérica de 2 horas pós-OGTT de 64,1 µU/mL - valores substancialmente superiores aos padrões universais de referência (Valencia et al., 2010). Essa divergência é central para a interpretação dos resultados: a rejeição da hipótese nula para determinados exames não indica inadequação do marcador clínico, mas sim que a população estudada possui um baseline metabólico e antropométrico distinto do adulto médio global.
Definimos a H₀ como: a mediana dos resultados de exames dos pacientes saudáveis é estatisticamente igual ao respectivo valor de referência médico universal, e a H₁ como: a mediana difere significativamente do valor esperado. Com o teste sobre as hipóteses mostradas, podemos identificar quais parâmetros já apresentam desvio sistêmico no grupo saudável, fornecendo subsídios sobre a predisposição fisiológica desta população. Adicionalmente, para mensurar a força dessa distorção, utilizamos o coeficiente descritivo de tamanho de efeito r de Rosenthal para indicar a magnitude da diferença, enquanto a direção do afastamento (para mais ou para menos) é determinada pela análise comparativa direta entre a mediana observada e o valor de referência.
A escolha por esta base de dados para a aplicação do teste mencionado justifica-se, primeiramente, pelo fato de as variáveis analisadas não seguirem uma distribuição normal, como demonstramos pelo teste de Shapiro-Wilk e por meio de visualizações diagnósticas. Sendo assim, os dados compreendem um cenário ideal para a aplicação de um teste não paramétrico. Estruturalmente, o conjunto de dados é composto por resultados de exames clínicos, como Glicose e Pressão Arterial, além de fatores demográficos e históricos, como a idade da paciente e o número de gestações. Por fim, a variável alvo Outcome (Resultado) atua como o indicador binário que define a presença ou ausência do diagnóstico de diabetes.
Para nossa pesquisa selecionamos apenas os pacientes saudáveis (Outcome = 0) e descartamos as variáveis de idade e gestações, visto que para elas não faria sentido comparar seus valores com uma referência clínica universal estabelecida. O mesmo se aplica para a variável DiabetesPedigreeFunction, que não possui valor de referencia por se tratar de uma métrica original dos pesquisadores que montaram o dataset. Também não segmentamos os pacientes por idade, sob a visão de que o “teste para uma amostra” não deveria englobar subgrupos, preservando a unidade amostral. Especificamente comparamos os valores observados com valores estáticos de referência, pois entendemos que o objetivo exigido do trabalho (teste para uma amostra) não compreendia a comparação entre grupos, requerida explicitamente para outra equipe. Os valores de referência adotados são provenientes de diretrizes internacionais (OMS, ADA) e estudos epidemiológicos em populações gerais, sendo intencionalmente conservadores para evidenciar as divergências basais da população Pima.
| Variável | Descrição | Valor de referência | Faixa normal | Justificativa | Referência |
|---|---|---|---|---|---|
| Glucose | Concentração de glicose plasmática após teste oral de tolerância à glicose. | 100 mg/dL | 60-140 | Ponto médio da faixa normal pós-OGTT (ADA/OMS). | https://www.ncbi.nlm.nih.gov/books/NBK532915/ |
| BloodPressure | Pressão arterial diastólica. | 70 mmHg | 60–80 | Média da faixa diastólica normal em adultos. | https://jamanetwork.com/journals/jamanetworkopen/fullarticle/2776530 |
| SkinThickness | Espessura da dobra cutânea do tríceps. | 22 mm | 16–28 | Ponto médio para adultos saudáveis. | https://pmc.ncbi.nlm.nih.gov/articles/PMC9127233/ |
| Insulin | Nível de insulina sérica em 2 horas. | 40 µU/mL | 16–166 | Média de insulina 2h (41,6 ± 30,8). | https://pmc.ncbi.nlm.nih.gov/articles/PMC5831567/ |
| BMI | Índice de Massa Corporal (IMC). | 21,7 kg/m² | 18,5–24,9 | Média exata da faixa OMS de peso normal. | https://www.who.int/news-room/fact-sheets/detail/obesity-and-overweight |
valores_referencia <- c(
Glucose = 100,
BloodPressure = 70,
SkinThickness = 22,
Insulin = 40,
BMI = 21.7
)
df <- read.csv('diabetes.csv')
head(df)
## Pregnancies Glucose BloodPressure SkinThickness Insulin BMI
## 1 6 148 72 35 0 33.6
## 2 1 85 66 29 0 26.6
## 3 8 183 64 0 0 23.3
## 4 1 89 66 23 94 28.1
## 5 0 137 40 35 168 43.1
## 6 5 116 74 0 0 25.6
## DiabetesPedigreeFunction Age Outcome
## 1 0.627 50 1
## 2 0.351 31 0
## 3 0.672 32 1
## 4 0.167 21 0
## 5 2.288 33 1
## 6 0.201 30 0
summary(df)
## Pregnancies Glucose BloodPressure SkinThickness
## Min. : 0.000 Min. : 0.0 Min. : 0.00 Min. : 0.00
## 1st Qu.: 1.000 1st Qu.: 99.0 1st Qu.: 62.00 1st Qu.: 0.00
## Median : 3.000 Median :117.0 Median : 72.00 Median :23.00
## Mean : 3.845 Mean :120.9 Mean : 69.11 Mean :20.54
## 3rd Qu.: 6.000 3rd Qu.:140.2 3rd Qu.: 80.00 3rd Qu.:32.00
## Max. :17.000 Max. :199.0 Max. :122.00 Max. :99.00
## Insulin BMI DiabetesPedigreeFunction Age
## Min. : 0.0 Min. : 0.00 Min. :0.0780 Min. :21.00
## 1st Qu.: 0.0 1st Qu.:27.30 1st Qu.:0.2437 1st Qu.:24.00
## Median : 30.5 Median :32.00 Median :0.3725 Median :29.00
## Mean : 79.8 Mean :31.99 Mean :0.4719 Mean :33.24
## 3rd Qu.:127.2 3rd Qu.:36.60 3rd Qu.:0.6262 3rd Qu.:41.00
## Max. :846.0 Max. :67.10 Max. :2.4200 Max. :81.00
## Outcome
## Min. :0.000
## 1st Qu.:0.000
## Median :0.000
## Mean :0.349
## 3rd Qu.:1.000
## Max. :1.000
str(df)
## 'data.frame': 768 obs. of 9 variables:
## $ Pregnancies : int 6 1 8 1 0 5 3 10 2 8 ...
## $ Glucose : int 148 85 183 89 137 116 78 115 197 125 ...
## $ BloodPressure : int 72 66 64 66 40 74 50 0 70 96 ...
## $ SkinThickness : int 35 29 0 23 35 0 32 0 45 0 ...
## $ Insulin : int 0 0 0 94 168 0 88 0 543 0 ...
## $ BMI : num 33.6 26.6 23.3 28.1 43.1 25.6 31 35.3 30.5 0 ...
## $ DiabetesPedigreeFunction: num 0.627 0.351 0.672 0.167 2.288 ...
## $ Age : int 50 31 32 21 33 30 26 29 53 54 ...
## $ Outcome : int 1 0 1 0 1 0 1 0 1 1 ...
sum(is.na(df))
## [1] 0
sum(duplicated(df))
## [1] 0
Nenhum valor nulo explícito ou duplicado
cols_com_zero_invalido <- c("Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI")
n_invalidas <- df |>
filter(if_any(all_of(cols_com_zero_invalido), ~ . == 0)) |>
nrow()
cat("Linhas com valores inválidos:", n_invalidas, "\n")
## Linhas com valores inválidos: 376
cat("Linhas restantes:", nrow(df) - n_invalidas, "\n")
## Linhas restantes: 392
df <- df |> filter(!if_any(all_of(cols_com_zero_invalido), ~ . == 0))
Muitos valores para exames são iguais a zero, porém isso é biologicamente impossível, logo consideramos como dados faltantes e os removemos
table(df$Outcome)
##
## 0 1
## 262 130
nrow(df[df$Age < 21, ])
## [1] 0
Nenhum paciente menor de 21 anos foi encontrado no dataset.
Temos 262 saudáveis, os demais (130 diabéticos) são descartados da análise
df <- df |> filter(Outcome == 0) |> select(-Outcome)
nrow(df)
## [1] 262
summary(df)
## Pregnancies Glucose BloodPressure SkinThickness
## Min. : 0.000 Min. : 56.0 Min. : 24.00 Min. : 7.00
## 1st Qu.: 1.000 1st Qu.: 94.0 1st Qu.: 60.00 1st Qu.:18.25
## Median : 2.000 Median :107.5 Median : 70.00 Median :27.00
## Mean : 2.721 Mean :111.4 Mean : 68.97 Mean :27.25
## 3rd Qu.: 4.000 3rd Qu.:126.0 3rd Qu.: 76.00 3rd Qu.:34.00
## Max. :13.000 Max. :197.0 Max. :106.00 Max. :60.00
## Insulin BMI DiabetesPedigreeFunction Age
## Min. : 15.0 Min. :18.20 Min. :0.0850 Min. :21.00
## 1st Qu.: 66.0 1st Qu.:26.12 1st Qu.:0.2610 1st Qu.:22.00
## Median :105.0 Median :31.25 Median :0.4135 Median :25.00
## Mean :130.9 Mean :31.75 Mean :0.4722 Mean :28.35
## 3rd Qu.:163.8 3rd Qu.:36.10 3rd Qu.:0.6242 3rd Qu.:30.00
## Max. :744.0 Max. :57.30 Max. :2.3290 Max. :81.00
É possível visualizar algumas distâncias dos valores esperados para saudáveis como no caso da glicose e insulina, mas apenas com o teste inferencial podemos definir se esses valores podem ou não ter surgido ao acaso.
As colunas mencionadas também são removidas
df <- df |> select(-Pregnancies, -Age, -DiabetesPedigreeFunction)
num_cols <- names(df)
df_long <- df |>
pivot_longer(cols = everything(), names_to = "Variavel", values_to = "Valor")
ggplot(df_long, aes(x = Valor)) +
geom_histogram(
aes(y = after_stat(density)),
bins = 30,
fill = "steelblue",
alpha = 0.5,
color = "white"
) +
geom_density(color = "steelblue", linewidth = 0.8) +
facet_wrap(~Variavel, scales = "free", ncol = 3) +
labs(
title = "Distribuição das variáveis (grupo saudável)",
x = NULL,
y = "Densidade"
) +
theme_minimal(base_size = 12)
A assimetria e caudas longas presentes nos gráficos já demonstram uma possível não normalidade de algumas distribuições, porém a análise visual ainda não é confiável para estabelecer uma decisão.
ggplot(df_long, aes(x = factor(""), y = Valor)) +
geom_boxplot(fill = "steelblue", alpha = 0.7, outlier.alpha = 0.4) +
facet_wrap(~Variavel, scales = "free_y", ncol = 3) +
labs(
title = "Boxplots das variáveis (grupo saudável)",
x = NULL,
y = NULL
) +
theme_minimal(base_size = 12) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
Notavelmente existem valores muito discrepantes dos intervalos interquartílicos aceitáveis, principalmente para a variável Insulin
resultados_shapiro <- map_dfr(num_cols, function(col) {
dados <- df[[col]]
teste <- shapiro.test(dados)
tibble(
Variavel = col,
`Estatística` = round(teste$statistic, 4),
`p-valor` = teste$p.value,
`Normal? (α=0.05)` = ifelse(teste$p.value > 0.05, "Sim", "Não"),
`Normal? (α=0.01)` = ifelse(teste$p.value > 0.01, "Sim", "Não")
)
})
print(resultados_shapiro)
## # A tibble: 5 × 5
## Variavel Estatística `p-valor` `Normal? (α=0.05)` `Normal? (α=0.01)`
## <chr> <dbl> <dbl> <chr> <chr>
## 1 Glucose 0.966 6.46e- 6 Não Não
## 2 BloodPressure 0.990 7.29e- 2 Sim Sim
## 3 SkinThickness 0.978 4.91e- 4 Não Não
## 4 Insulin 0.771 7.45e-19 Não Não
## 5 BMI 0.984 4.44e- 3 Não Não
Apenas BloodPressure demonstrou normalidade perante o teste, logo a descartamos já que não faria sentido aplicar wilcoxon para dados normais.
df <- df |> select(-BloodPressure)
Usamos a função wilcox.test do R base. Infelizmente uma
limitação da função é não retornar a estatística Z associada ao teste, e
como queríamos expô-la, foi preciso criar uma função que realiza o teste
de Wilcoxon manualmente, mas não em totalidade, apenas o suficiente para
gerarmos a estatística Z. Também nos atentamos a manter essa nossa
função (z_stat) alinhada com o funcionamento da função
wilcox.test utilizada, para garantir o mesmo valor de Z
calculado internamente pelo método.
Ademais utilizamos o valor de Z para calcularmos também o R de rosenthal, que nos permite definir o tamanho de efeito da variação, se muito alto, alto, médio ou baixo. Na tabela apresentada o sinal de rosenthal (+) ou (-) é calculado a partir da diferença da mediana observada com o valor referencia, já que R não mostra o sentido da variação, apenas seu nível.
x — Vetor numérico de dados (a amostra
observada).
y = NULL — Segundo vetor numérico de
dados, opcional. Se fornecido, realiza o teste de Wilcoxon de duas
amostras (Mann-Whitney U). Se NULL, realiza o teste de uma
amostra ou pareado, dependendo de paired.
alternative = "two.sided" — Hipótese
alternativa: "two.sided" detecta diferença em qualquer
direção; "less" testa se a locação é menor que
mu; "greater" se é maior.
mu = 0 — Parâmetro de locação sob H₀.
No teste de uma amostra ou pareado é o valor de referência; no de duas
amostras é a diferença de locação esperada.
paired = FALSE — Se TRUE,
realiza o teste de postos com sinal para amostras pareadas — requer
x e y com o mesmo comprimento. Se
FALSE, realiza o teste de soma de postos (Mann-Whitney U)
quando y é fornecido.
exact = NULL — Controla o cálculo do
p-valor: TRUE usa a distribuição exata dos postos;
FALSE usa a aproximação normal; NULL escolhe
automaticamente — exato se n < 50 e sem empates, aproximação caso
contrário.
correct = TRUE — Aplica correção de
continuidade na aproximação normal. Só tem efeito quando
exact = FALSE.
conf.int = FALSE — Se
TRUE, retorna intervalo de confiança e estimativa pontual
da locação (pseudomediana para uma amostra; diferença de Hodges-Lehmann
para duas amostras).
conf.level = 0.95 — Nível de confiança
para o intervalo retornado por conf.int. Padrão: 0,95.
formula = NULL — Interface de fórmula
do tipo y ~ grupo, onde grupo é um fator com
dois níveis. Alternativa à especificação direta de x e
y.
data = NULL — Data frame onde as
variáveis da fórmula serão buscadas. Só tem efeito com
formula.
subset = NULL — Vetor lógico ou
numérico para selecionar um subconjunto de observações quando
formula é usada.
na.action = NULL — Função que define o
tratamento de NAs. O padrão segue
getOption("na.action"); por exemplo, na.omit
remove as observações ausentes silenciosamente.
Para as demais configurações mantemos os valores padrão, por se tratarem de configurações estruturais para se adequarem ao estilo de código implementado pelo cientista.
rosenthal_r <- function(z, n) {
z / sqrt(n)
}
rosenthal_label <- function(r) {
abs_r <- abs(r)
stars <- dplyr::case_when(
abs_r < 0.10 ~ "",
abs_r < 0.30 ~ "*",
abs_r < 0.50 ~ "**",
TRUE ~ "***"
)
sinal <- ifelse(r >= 0, "+", "-")
paste0("(", sinal, ")", stars)
}
z_stat <- function(diferencas) {
diferencas_nao_zero <- diferencas[diferencas != 0]
n <- length(diferencas_nao_zero)
if (n == 0) return(list(z = NA_real_, n = n))
abs_ranks <- rank(abs(diferencas_nao_zero))
T_pos <- sum(abs_ranks[diferencas_nao_zero > 0])
media_w <- n * (n + 1) / 4
sd_w <- sqrt(n * (n + 1) * (2 * n + 1) / 24)
z <- ((T_pos - media_w) - 0.5) / sd_w
list(z = z, n = n)
}
num_cols <- intersect(names(valores_referencia), names(df))
resultados_wilcoxon <- map_dfr(num_cols, function(col) {
dados <- df[[col]]
referencia <- valores_referencia[col]
diferencas <- dados - referencia
teste <- wilcox.test(
dados,
mu = referencia,
alternative = "two.sided",
exact = FALSE,
correct = TRUE
)
zn <- z_stat(diferencas)
z <- zn$z
n <- zn$n
r <- rosenthal_r(z, n)
tibble(
Variavel = col,
`Mediana Observada` = round(median(dados, na.rm = TRUE), 4),
`Valor Esperado` = referencia,
n = n,
W = round(teste$statistic, 4),
Z = round(z, 4),
r = rosenthal_label(r),
`p-valor` = round(teste$p.value, 6),
`α=0.01` = ifelse(teste$p.value < 0.01, "Rejeita H₀", "Não Rejeita H₀"),
`α=0.05` = ifelse(teste$p.value < 0.05, "Rejeita H₀", "Não Rejeita H₀"),
`α=0.10` = ifelse(teste$p.value < 0.10, "Rejeita H₀", "Não Rejeita H₀")
)
})
resultados_wilcoxon <- column_to_rownames(resultados_wilcoxon, var = "Variavel")
print(resultados_wilcoxon)
## Mediana Observada Valor Esperado n W Z r
## Glucose 107.50 100.0 251 23194.5 6.4106 (+)**
## SkinThickness 27.00 22.0 256 24688.0 6.9480 (+)**
## Insulin 105.00 40.0 260 33531.0 13.6485 (+)***
## BMI 31.25 21.7 262 34181.5 13.8096 (+)***
## p-valor α=0.01 α=0.05 α=0.10
## Glucose 0 Rejeita H₀ Rejeita H₀ Rejeita H₀
## SkinThickness 0 Rejeita H₀ Rejeita H₀ Rejeita H₀
## Insulin 0 Rejeita H₀ Rejeita H₀ Rejeita H₀
## BMI 0 Rejeita H₀ Rejeita H₀ Rejeita H₀
Rejeitamos as hipóteses nulas para todos os testes com graus médios/altos de magnitude para a diferença entre os valores do grupo e o valor esperado. Em todos os casos os pacientes de Pima superaram os valores mundialmente tidos como normais para os atributos estudados, para os 3 níveis de significância clássicos, revelando como a população, mesmo saudável, é muito distoante do padrão mundial de saúde.
resultados_ttest <- map_dfr(num_cols, function(col) {
dados <- na.omit(df[[col]])
referencia <- valores_referencia[col]
teste <- t.test(
dados,
mu = referencia,
alternative = "two.sided"
)
n <- length(dados)
tibble(
Variavel = col,
`Média Observada` = round(mean(dados), 4),
`Valor Esperado` = referencia,
t = round(teste$statistic, 4),
gl = n - 1,
`p-valor` = round(teste$p.value, 6),
`α=0.01` = ifelse(teste$p.value < 0.01, "Rejeita H₀", "Não Rejeita H₀"),
`α=0.05` = ifelse(teste$p.value < 0.05, "Rejeita H₀", "Não Rejeita H₀"),
`α=0.10` = ifelse(teste$p.value < 0.10, "Rejeita H₀", "Não Rejeita H₀")
)
})
resultados_ttest <- column_to_rownames(resultados_ttest, var = "Variavel")
print(resultados_ttest)
## Média Observada Valor Esperado t gl p-valor α=0.01
## Glucose 111.4313 100.0 7.5088 261 0 Rejeita H₀
## SkinThickness 27.2519 22.0 8.1473 261 0 Rejeita H₀
## Insulin 130.8550 40.0 14.3298 261 0 Rejeita H₀
## BMI 31.7508 21.7 23.9421 261 0 Rejeita H₀
## α=0.05 α=0.10
## Glucose Rejeita H₀ Rejeita H₀
## SkinThickness Rejeita H₀ Rejeita H₀
## Insulin Rejeita H₀ Rejeita H₀
## BMI Rejeita H₀ Rejeita H₀
O relativo mesmo resultado, aqui também rejeitamos as hipóteses nulas com p-valores muito pequenos. Não foi possível observar uma diferença significativa entre os testes
numiqo - https://www.youtube.com/watch?v=NZsL2eDQiDQ&pp=ygUOdGVzdGUgd2lsY294b24%3D
Valencia, M. E., et al. (2010). Differences in Insulin Resistance in Mexican and U.S. Pima Indians with Normal Glucose Tolerance. Diabetes Care, 33(2), 285–289. https://pmc.ncbi.nlm.nih.gov/articles/PMC2968731/
World Health Organization. (2021). Obesity and overweight. WHO Fact Sheets. https://www.who.int/news-room/fact-sheets/detail/obesity-and-overweight
American Diabetes Association. (2024). Standards of Care in Diabetes—2024. Diabetes Care, 47(Suppl. 1). https://www.ncbi.nlm.nih.gov/books/NBK532915/