library(tidyverse)
library(ggplot2)

Introdução

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.

Descrição das variáveis e valores de referência selecionados

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)

Configuração do teste de Wilcoxon

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.

Todas as configurações possíveis

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.

Escolha dos parâmetros

  1. x: Os dados observados de cada variável
  2. mu: O respectivo valor de referência clínica
  3. alternative: “two.sided”, pois queremos detectar diferença para mais ou para menos
  4. exact: FALSE, pois usamos a aproximação normal para um teste clássico
  5. correct: TRUE correção de continuidade, já que usa aproximação normal

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.

Comparação com teste T

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

Referências

numiqo - https://www.youtube.com/watch?v=NZsL2eDQiDQ&pp=ygUOdGVzdGUgd2lsY294b24%3D

numiqo - https://www.youtube.com/watch?v=2AqoK8itEFQ&t=223s&pp=ugMGCgJwdBABugUEEgJwdMoFDnRlc3RlIHdpbGNveG9u0gcJCT8LAYcqIYzv2AcB

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/