1. Leitura e preparação do banco

knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE
)

pacman::p_load(readxl, tidyverse, ggpubr, table1, nortest, patchwork, janitor, gtsummary, dplyr)

df = read_xlsx("df.xlsx") %>%
  clean_names()

df$cidade <- factor(df$cidade,
                    levels = c(0, 1, 2, 3),
                    labels = c("Porto Alegre", "Canoas", "Guaiba", "Gravatai"))
table(df$cidade)
## 
## Porto Alegre       Canoas       Guaiba     Gravatai 
##           75           75           75           75
df$tabagismo <- factor(df$tabagismo,
                       levels = c(0, 1, 2),
                       labels = c("Não Fumante", "Ex-Fumante", "Fumante"))
table(df$tabagismo)
## 
## Não Fumante  Ex-Fumante     Fumante 
##         195          49          56

2. Tabela demográfica

tabela_demo <- df %>%
  select(idade, altura, peso, imc, tabagismo, cidade) %>%
  tbl_summary(
    by = cidade,
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 1
  )

tabela_demo
Characteristic Porto Alegre
N = 751		 Canoas
N = 751		 Guaiba
N = 751		 Gravatai
N = 751
idade 34.0 (11.0) 35.6 (11.2) 36.1 (11.6) 29.8 (10.3)
altura 173.5 (7.8) 175.0 (6.8) 172.2 (7.9) 174.5 (8.1)
peso 82.0 (15.4) 82.8 (14.2) 80.8 (15.6) 81.1 (13.5)
imc 26.9 (5.7) 27.0 (4.0) 27.2 (4.4) 26.6 (4.3)
tabagismo



    Não Fumante 46 (61%) 44 (59%) 52 (69%) 53 (71%)
    Ex-Fumante 14 (19%) 13 (17%) 17 (23%) 5 (6.7%)
    Fumante 15 (20%) 18 (24%) 6 (8.0%) 17 (23%)
1 Mean (SD); n (%)

3. Tabela função pulmonar

tabela_pulmonar <- df %>%
  select(fev1z, fvcz, fev1fvcz, fef2575z, cidade) %>%
  tbl_summary(
    by = cidade,
    statistic = list(
      all_continuous() ~ "{mean} ({sd})"
    ),
    digits = all_continuous() ~ 1
  )

tabela_pulmonar
Characteristic Porto Alegre
N = 751		 Canoas
N = 751		 Guaiba
N = 751		 Gravatai
N = 751
fev1z -0.8 (1.0) -0.5 (0.8) -0.7 (0.8) -1.2 (1.0)
fvcz -1.0 (1.0) -0.6 (0.8) -0.9 (0.9) -1.2 (1.0)
fev1fvcz 0.4 (1.2) 0.2 (1.1) 0.4 (1.0) 0.2 (1.3)
fef2575z -0.2 (1.0) -0.1 (0.9) -0.1 (0.8) -0.5 (1.1)
1 Mean (SD)

4. Teste de normalidade

lillie.test(df$fvcz)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  df$fvcz
## D = 0.058471, p-value = 0.01504
lillie.test(df$fev1z)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  df$fev1z
## D = 0.050336, p-value = 0.06404
lillie.test(df$fev1fvcz)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  df$fev1fvcz
## D = 0.050541, p-value = 0.06195
lillie.test(df$fef2575z)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  df$fef2575z
## D = 0.056845, p-value = 0.02054

5. Diferença de função pulmonar entre as cidades

kruskal.test(fev1z ~ cidade, data = df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  fev1z by cidade
## Kruskal-Wallis chi-squared = 17.793, df = 3, p-value = 0.0004852
pairwise.wilcox.test(df$fev1z, df$cidade, p.adjust.method = "hochberg")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$fev1z and df$cidade 
## 
##          Porto Alegre Canoas  Guaiba 
## Canoas   0.22813      -       -      
## Guaiba   0.50107      0.32408 -      
## Gravatai 0.14384      0.00026 0.01529
## 
## P value adjustment method: hochberg
kruskal.test(fvcz ~ cidade, data = df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  fvcz by cidade
## Kruskal-Wallis chi-squared = 15.091, df = 3, p-value = 0.00174
pairwise.wilcox.test(df$fvcz, df$cidade, p.adjust.method = "hochberg")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$fvcz and df$cidade 
## 
##          Porto Alegre Canoas Guaiba
## Canoas   0.0450       -      -     
## Guaiba   0.5011       0.1360 -     
## Gravatai 0.4297       0.0013 0.2300
## 
## P value adjustment method: hochberg
kruskal.test(fev1fvcz ~ cidade, data = df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  fev1fvcz by cidade
## Kruskal-Wallis chi-squared = 4.733, df = 3, p-value = 0.1924
kruskal.test(fef2575z ~ cidade, data = df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  fef2575z by cidade
## Kruskal-Wallis chi-squared = 7.9169, df = 3, p-value = 0.04776
pairwise.wilcox.test(df$fef2575z, df$cidade, p.adjust.method = "hochberg")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$fef2575z and df$cidade 
## 
##          Porto Alegre Canoas Guaiba
## Canoas   0.781        -      -     
## Guaiba   0.781        0.781  -     
## Gravatai 0.285        0.208  0.035 
## 
## P value adjustment method: hochberg

6. Gráfico

df_long <- df %>%
  select(cidade, fvcz, fev1z, fev1fvcz, fef2575z) %>%
  pivot_longer(
    cols = -cidade,
    names_to = "variavel",
    values_to = "valor"
  )

labels <- c(
  fvcz = "FVC (z score)",
  fev1z = "FEV1 (z score)",
  fev1fvcz = "FEV1/FVC (z score)",
  fef2575z = "FEF25–75 (z score)"
)

ggplot(df_long, aes(x = cidade, y = valor, fill = cidade)) +
  geom_boxplot(color = "black") +
  facet_wrap(~variavel, scales = "free_y", labeller = as_labeller(labels)) +
  labs(
    x = "Cidade",
    y = "",
    fill = "Cidade"
  ) +
  scale_fill_manual(values = c(
    "Porto Alegre" = "#1b9e77",
    "Canoas"       = "#d95f02",
    "Guaiba"       = "#7570b3",
    "Gravatai"     = "#e7298a"
  )) +
  theme_bw() +
theme(
  legend.position = "right",
  axis.text.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.x = element_blank()
)