Prova cristiano_aguzzoli

1. Leitura preparação

# install packages 
packages <- c(
  "tidyverse", "ggpubr", "table1", "nortest", "patchwork", "readxl", "gtable", "ggplot2", "patchwork")

# Load installed packages
lapply(packages, require, character.only = TRUE)
## Loading required package: tidyverse
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'tidyverse'
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 4.3.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.3
## Loading required package: table1
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
## Loading required package: nortest
## Warning: package 'nortest' was built under R version 4.3.3
## Loading required package: patchwork
## Warning: package 'patchwork' was built under R version 4.3.3
## Loading required package: readxl
## Loading required package: gtable
## Warning: package 'gtable' was built under R version 4.3.3
## [[1]]
## [1] FALSE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] TRUE
## 
## [[5]]
## [1] TRUE
## 
## [[6]]
## [1] TRUE
## 
## [[7]]
## [1] TRUE
## 
## [[8]]
## [1] TRUE
## 
## [[9]]
## [1] TRUE
# load df 
df <- read_excel("~/Documents/pucrs_phd/df.xlsx")

2. Recodifique cidade e tabagismo e altere para factor

df$cidade[df$cidade == "0"] <- "Porto Alegre"
df$cidade[df$cidade == "1"] <- "Canoas"
df$cidade[df$cidade == "2"] <- "Guaiba"
df$cidade[df$cidade == "3"] <- "Gravatai"
df$cidade <- factor(df$cidade, levels = c("Porto Alegre", "Canoas", "Guaiba", "Gravatai"))

df$tabagismo[df$tabagismo == "0"] <- "Não Fumante"
df$tabagismo[df$tabagismo == "1"] <- "Ex-Fumante"
df$tabagismo[df$tabagismo == "2"] <- "Fumante"
df$tabagismo <- factor(df$tabagismo, levels = c("Não Fumante", "Ex-Fumante", "Fumante"))

3. Tabela demográfica

# demographic characteristics 
table1(
  ~ idade + altura + peso + imc + tabagismo | cidade,
  data = df,
  overall = "Overall")
Porto Alegre
(N=75)
Canoas
(N=75)
Guaiba
(N=75)
Gravatai
(N=75)
Overall
(N=300)
idade
Mean (SD) 34.0 (11.0) 35.6 (11.2) 36.1 (11.6) 29.8 (10.3) 33.9 (11.3)
Median [Min, Max] 35.5 [18.0, 53.8] 35.7 [18.6, 62.9] 34.8 [18.3, 55.9] 26.1 [18.1, 54.2] 33.0 [18.0, 62.9]
altura
Mean (SD) 173 (7.78) 175 (6.78) 172 (7.86) 175 (8.13) 174 (7.69)
Median [Min, Max] 174 [157, 192] 175 [162, 196] 172 [157, 197] 174 [150, 197] 174 [150, 197]
peso
Mean (SD) 82.0 (15.4) 82.8 (14.2) 80.8 (15.6) 81.1 (13.5) 81.7 (14.6)
Median [Min, Max] 82.0 [54.0, 125] 82.0 [50.0, 132] 80.0 [58.0, 140] 80.0 [54.0, 122] 80.0 [50.0, 140]
imc
Mean (SD) 26.9 (5.70) 27.0 (3.96) 27.2 (4.37) 26.6 (4.32) 26.9 (4.62)
Median [Min, Max] 27.0 [0.00257, 38.7] 26.6 [17.7, 39.0] 26.1 [19.6, 40.9] 25.8 [18.6, 41.7] 26.4 [0.00257, 41.7]
tabagismo
Não Fumante 46 (61.3%) 44 (58.7%) 52 (69.3%) 53 (70.7%) 195 (65.0%)
Ex-Fumante 14 (18.7%) 13 (17.3%) 17 (22.7%) 5 (6.7%) 49 (16.3%)
Fumante 15 (20.0%) 18 (24.0%) 6 (8.0%) 17 (22.7%) 56 (18.7%)

4. Tabela demográfica da função pulmonar por cidade

# funcao pulmonar por cidade 
table1(
  ~ fvcz + fev1z + fev1fvcz + fef2575z | cidade,
  data = df,
  overall = "Overall")
Porto Alegre
(N=75)
Canoas
(N=75)
Guaiba
(N=75)
Gravatai
(N=75)
Overall
(N=300)
fvcz
Mean (SD) -0.980 (1.02) -0.577 (0.830) -0.921 (0.920) -1.24 (1.02) -0.930 (0.974)
Median [Min, Max] -0.975 [-4.11, 1.61] -0.537 [-2.84, 1.26] -0.872 [-3.82, 1.33] -1.06 [-3.61, 0.700] -0.883 [-4.11, 1.61]
fev1z
Mean (SD) -0.760 (1.05) -0.493 (0.848) -0.700 (0.847) -1.16 (1.02) -0.779 (0.972)
Median [Min, Max] -0.772 [-3.20, 1.34] -0.467 [-2.73, 1.37] -0.589 [-3.39, 1.57] -0.980 [-4.13, 1.23] -0.724 [-4.13, 1.57]
fev1fvcz
Mean (SD) 0.423 (1.18) 0.192 (1.13) 0.444 (1.03) 0.158 (1.30) 0.304 (1.16)
Median [Min, Max] 0.348 [-2.48, 3.55] 0.173 [-2.57, 2.46] 0.470 [-2.51, 3.61] 0.176 [-3.65, 3.74] 0.319 [-3.65, 3.74]
fef2575z
Mean (SD) -0.184 (0.954) -0.137 (0.899) -0.0772 (0.782) -0.491 (1.09) -0.222 (0.946)
Median [Min, Max] -0.126 [-2.36, 1.93] -0.178 [-2.57, 1.59] -0.0775 [-2.53, 1.42] -0.440 [-4.38, 1.68] -0.175 [-4.38, 1.93]

5. Normalidade das variaveis de funcao pulmonar com Lilliefors (Kolmogorov-Smirnov) Test

# normalidade das variaveis de funcao pulmonar com  Lilliefors (Kolmogorov-Smirnov) Test 
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

6. Testar diferenças entre funcao pulmonar entre duas cidades com Kruskal.test

# testar diferenças entre funcao pulmonar entre duas cidades com Kruskal.test
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
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
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

Teste post-hoc comparando a função pulmonar entre as cidades com pairwise.wilcox.test

# teste post-hoc comparando a função pulmonar entre as cidades com pairwise.wilcox.test
pairwise.wilcox.test(df$fvcz, df$cidade, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$fvcz and df$cidade 
## 
##          Porto Alegre Canoas Guaiba
## Canoas   0.0540       -      -     
## Guaiba   1.0000       0.2041 -     
## Gravatai 1.0000       0.0013 0.4600
## 
## P value adjustment method: bonferroni
pairwise.wilcox.test(df$fev1z, df$cidade, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$fev1z and df$cidade 
## 
##          Porto Alegre Canoas  Guaiba 
## Canoas   0.45626      -       -      
## Guaiba   1.00000      0.97223 -      
## Gravatai 0.21576      0.00026 0.01835
## 
## P value adjustment method: bonferroni
pairwise.wilcox.test(df$fef2575z, df$cidade, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$fef2575z and df$cidade 
## 
##          Porto Alegre Canoas Guaiba
## Canoas   1.000        -      -     
## Guaiba   1.000        1.000  -     
## Gravatai 0.427        0.250  0.035 
## 
## P value adjustment method: bonferroni

Gráfico combinado de cada variável de função pulmonar por cidade

# gráfico combinado de cada variável de função pulmonar por cidade 

plot1 <- ggplot(data = df, aes(x = cidade, y = fvcz, fill = cidade)) +
  geom_boxplot(alpha = 0.5) +
  stat_compare_means(method = "kruskal.test") +
  scale_fill_viridis_d(option = "D") +
  theme_classic() +
  theme(legend.position = "none") +
  labs(x = "Município", y = "FVC (z score)")

plot2 <- ggplot(data = df, aes(x = cidade, y = fev1z, fill = cidade)) +
  geom_boxplot(alpha = 0.5) +
  stat_compare_means(method = "kruskal.test") +
  scale_fill_viridis_d(option = "D") +
  theme_classic() +
  theme(legend.position = "none") +
  labs(x = "Município", y = "FEV1 (z score)")

plot3 <- ggplot(data = df, aes(x = cidade, y = fev1fvcz, fill = cidade)) +
  geom_boxplot(alpha = 0.5) +
  stat_compare_means(method = "kruskal.test") +
  scale_fill_viridis_d(option = "D") +
  theme_classic() +
  theme(legend.position = "none") +
  labs(x = "Município", y = "FEV1/FVC (z score)")

plot4 <- ggplot(data = df, aes(x = cidade, y = fef2575z, fill = cidade)) +
  geom_boxplot(alpha = 0.5) +
  stat_compare_means(method = "kruskal.test") +
  scale_fill_viridis_d(option = "D") +
  theme_classic() +
  theme(legend.position = "none") +
  labs(x = "Município", y = "FEF25-75 (z score)")


  (plot1 + plot2) /
  (plot3 + plot4) +
  plot_annotation(
    title = "Função pulmonar por cidade",
    tag_levels = "A"
  ) &
  theme(
    plot.margin = margin(5, 5, 5, 5),
    plot.tag = element_text(size = 14, face = "bold"),
    plot.title = element_text(
      size = 14,
      face = "bold",
      hjust = 0.5
    )
  )