# 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")
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"))
# 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%) |
# 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] |
# 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
# 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
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
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
)
)