library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.4 v dplyr 1.0.2
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggplot2)
library(ggpubr)
library(vroom)
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
library(survival)
library(ranger)
library(ggplot2)
library(dplyr)
library(ggfortify)
library(survminer)
library(knitr)
db <- vroom("saint.csv")
## Rows: 672
## Columns: 49
## Delimiter: ","
## chr [ 1]: studyno
## dbl [48]: treat, day, age, sex, vsheartrate, vsbloodpresssystolic, vsbloodpressdiastolic, ...
##
## Use `spec()` to retrieve the guessed column specification
## Pass a specification to the `col_types` argument to quiet this message
db1 <- db %>%
select(studyno, treat,day,ct_e,ct_n,vl_n,vl_e) %>% #filter(day > 4) %>%
mutate(treat = factor(treat)) %>%
na.omit() %>% as.data.frame() %>% ungroup() %>%
mutate(ct_e1 = ifelse(ct_e<=30,0,1),
ct_n1 = ifelse(ct_n<=30,0,1),
treat = ifelse(treat==1,"Ivermectina","Placebo"))
db2 <- db %>%
select(studyno, treat,igg) %>% #filter(day > 4) %>%
mutate(treat = factor(treat)) %>%
na.omit() %>% as.data.frame() %>% ungroup() %>%
mutate(treat = ifelse(treat==1,"Ivermectina","Placebo"))
bxp <- db2 %>% ggboxplot(y="igg", x="treat", bxp.errorbar = T, add = "jitter", fill="treat", palette = "Dark2", ylab="IgG", title = "IgG") + theme(legend.position = "right") + labs(fill="Treatment")
## Warning: Ignoring unknown aesthetics: fill
wilcox_test(db2,igg ~ treat, detailed = T)
## # A tibble: 1 x 12
## estimate .y. group1 group2 n1 n2 statistic p conf.low conf.high
## * <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1.47 igg Iverm~ Place~ 12 12 52 0.266 -4.43 1.78
## # ... with 2 more variables: method <chr>, alternative <chr>
stat.test <- db2 %>%
wilcox_test(igg ~ treat) %>%
add_significance()
stat.test
## # A tibble: 1 x 8
## .y. group1 group2 n1 n2 statistic p p.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
## 1 igg Ivermectina Placebo 12 12 52 0.266 ns
db2 %>% wilcox_effsize(igg ~ treat)
## # A tibble: 1 x 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 igg Ivermectina Placebo 0.236 12 12 small
ggdensity(db2, x = "igg", rug = TRUE, fill = "lightgray") +
scale_x_continuous(limits = c(0, 18)) +
stat_central_tendency(type = "median", color = "red", linetype = "dashed") +
geom_vline(xintercept = 10, color = "blue", linetype = "dashed") +
labs(subtitle = get_test_label(stat.test, detailed = TRUE))
stat.test <- stat.test %>% add_xy_position(x = "treat")
bxp +
stat_pvalue_manual(stat.test, tip.length = 0) +
labs(subtitle = get_test_label(stat.test, detailed = TRUE))
db2 %>%
group_by(treat) %>%
get_summary_stats(igg, type = "full") %>% kable
| Ivermectina |
igg |
12 |
2.483 |
15.837 |
4.748 |
3.608 |
8.588 |
4.979 |
2.712 |
6.540 |
4.208 |
1.215 |
2.674 |
| Placebo |
igg |
12 |
2.973 |
17.286 |
7.506 |
4.752 |
9.390 |
4.637 |
3.954 |
7.785 |
3.815 |
1.101 |
2.424 |
library(broom)
fit <- survfit(Surv(day, ct_e1) ~ treat, data = subset(db1))
tidy(fit) %>% arrange(time) %>% filter(time<21) %>% kable()
| 1 |
49 |
0 |
12 |
1.0000000 |
0.0000000 |
1.0000000 |
1.0000000 |
treat=Ivermectina |
| 1 |
49 |
0 |
12 |
1.0000000 |
0.0000000 |
1.0000000 |
1.0000000 |
treat=Placebo |
| 4 |
37 |
4 |
8 |
0.8918919 |
0.0572364 |
0.9977734 |
0.7972463 |
treat=Ivermectina |
| 4 |
37 |
1 |
11 |
0.9729730 |
0.0273998 |
1.0000000 |
0.9220999 |
treat=Placebo |
| 7 |
25 |
8 |
3 |
0.6064865 |
0.1486591 |
0.8116336 |
0.4531920 |
treat=Ivermectina |
| 7 |
25 |
4 |
8 |
0.8172973 |
0.0914866 |
0.9778080 |
0.6831350 |
treat=Placebo |
| 14 |
14 |
7 |
1 |
0.3032432 |
0.3058236 |
0.5522156 |
0.1665227 |
treat=Ivermectina |
| 14 |
13 |
10 |
0 |
0.1886071 |
0.5145678 |
0.5170814 |
0.0687950 |
treat=Placebo |
summary(fit)$table
## records n.max n.start events *rmean *se(rmean) median
## treat=Ivermectina 49 49 49 25 13.04378 1.1632450 14
## treat=Placebo 49 49 49 18 13.96025 0.9482471 14
## 0.95LCL 0.95UCL
## treat=Ivermectina 7 NA
## treat=Placebo 14 NA
res.sum <- surv_summary(fit)
## Warning in .get_data(x, data = data): The `data` argument is not provided. Data
## will be extracted from model fit.
res.sum
## time n.risk n.event n.censor surv std.err upper lower
## 1 1 49 0 12 1.0000000 0.00000000 1.0000000 1.00000000
## 2 4 37 4 8 0.8918919 0.05723638 0.9977734 0.79724631
## 3 7 25 8 3 0.6064865 0.14865912 0.8116336 0.45319203
## 4 14 14 7 1 0.3032432 0.30582365 0.5522156 0.16652275
## 5 21 6 6 0 0.0000000 Inf NA NA
## 6 1 49 0 12 1.0000000 0.00000000 1.0000000 1.00000000
## 7 4 37 1 11 0.9729730 0.02739983 1.0000000 0.92209995
## 8 7 25 4 8 0.8172973 0.09148660 0.9778080 0.68313498
## 9 14 13 10 0 0.1886071 0.51456783 0.5170814 0.06879502
## 10 21 3 3 0 0.0000000 Inf NA NA
## strata treat
## 1 treat=Ivermectina Ivermectina
## 2 treat=Ivermectina Ivermectina
## 3 treat=Ivermectina Ivermectina
## 4 treat=Ivermectina Ivermectina
## 5 treat=Ivermectina Ivermectina
## 6 treat=Placebo Placebo
## 7 treat=Placebo Placebo
## 8 treat=Placebo Placebo
## 9 treat=Placebo Placebo
## 10 treat=Placebo Placebo
attr(res.sum, "table")
## records n.max n.start events *rmean *se(rmean) median
## treat=Ivermectina 49 49 49 25 13.04378 1.1632450 14
## treat=Placebo 49 49 49 18 13.96025 0.9482471 14
## 0.95LCL 0.95UCL
## treat=Ivermectina 7 NA
## treat=Placebo 14 NA
surv_diff <- survdiff(Surv(day, ct_e1) ~ treat, data = subset(db1))
surv_diff
## Call:
## survdiff(formula = Surv(day, ct_e1) ~ treat, data = subset(db1))
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## treat=Ivermectina 49 25 23.3 0.122 0.552
## treat=Placebo 49 18 19.7 0.144 0.552
##
## Chisq= 0.6 on 1 degrees of freedom, p= 0.5
d <- data.frame(time = fit$time,
n.risk = fit$n.risk,
n.event = fit$n.event,
n.censor = fit$n.censor,
surv = fit$surv,
upper = fit$upper,
lower = fit$lower)
a<- ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
xlim = c(0, 21))
ggsurvplot(fit,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c('Dark2'))
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
fit1 <- survfit(Surv(day, ct_n1) ~ treat, data = subset(db1))
fit1
## Call: survfit(formula = Surv(day, ct_n1) ~ treat, data = subset(db1))
##
## n events median 0.95LCL 0.95UCL
## treat=Ivermectina 49 25 14 7 NA
## treat=Placebo 49 18 14 14 NA
summary(fit1)
## Call: survfit(formula = Surv(day, ct_n1) ~ treat, data = subset(db1))
##
## treat=Ivermectina
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 37 4 0.892 0.0510 0.797 0.998
## 7 25 8 0.606 0.0902 0.453 0.812
## 14 14 7 0.303 0.0927 0.167 0.552
## 21 6 6 0.000 NaN NA NA
##
## treat=Placebo
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 37 1 0.973 0.0267 0.922 1.000
## 7 25 5 0.778 0.0807 0.635 0.954
## 14 13 9 0.240 0.1027 0.103 0.555
## 21 3 3 0.000 NaN NA NA
summary(fit1)$table
## records n.max n.start events *rmean *se(rmean) median
## treat=Ivermectina 49 49 49 25 13.04378 1.163245 14
## treat=Placebo 49 49 49 18 14.04407 1.034350 14
## 0.95LCL 0.95UCL
## treat=Ivermectina 7 NA
## treat=Placebo 14 NA
d1 <- data.frame(time = fit1$time,
n.risk = fit1$n.risk,
n.event = fit1$n.event,
n.censor = fit1$n.censor,
surv = fit1$surv,
upper = fit1$upper,
lower = fit1$lower)
b<-ggsurvplot(fit1,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
xlim = c(0, 21))
ggsurvplot(fit1,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c('Dark2'))
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
ggarrange(a$plot,b$plot,common.legend = T, legend = "right")
ggsurvplot(
fit, # survfit object with calculated statistics.
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for
# point estimaes of survival curves.
conf.int.style = "step", # customize style of confidence intervals
xlab = "Time in days", # customize X axis label.
break.time.by = 4, # break X axis in time intervals by 200.
ggtheme = theme_light(), # customize plot and risk table with a theme.
risk.table = "abs_pct", # absolute number and percentage at risk.
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.y.text = FALSE,# show bars instead of names in text annotations
# in legend of risk table.
ncensor.plot = TRUE, # plot the number of censored subjects at time t
surv.median.line = "hv", # add the median survival pointer.
palette =
c("#E7B800", "#2E9FDF") # custom color palettes.
)
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
ggsurvplot(
fit1, # survfit object with calculated statistics.
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for
# point estimaes of survival curves.
conf.int.style = "step", # customize style of confidence intervals
xlab = "Time in days", # customize X axis label.
break.time.by = 4, # break X axis in time intervals by 200.
ggtheme = theme_light(), # customize plot and risk table with a theme.
risk.table = "abs_pct", # absolute number and percentage at risk.
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.y.text = FALSE,# show bars instead of names in text annotations
# in legend of risk table.
ncensor.plot = TRUE, # plot the number of censored subjects at time t
surv.median.line = "hv", # add the median survival pointer.
palette =
c("#E7B800", "#2E9FDF") # custom color palettes.
)
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
xlim = c(0, 21))
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
fun = "event")
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
fun = "cumhaz")
m <- length(fit$strata)
df <- data.frame(time = c(rep(0, m), fit$time),
surv = c(rep(1, m), fit$surv),
group = c(names(fit$strata),
rep(names(fit$strata), fit$strata)))
ggplot(data = df) +
geom_step(aes(x = time, y = surv, colour = as.factor(group))) +
ylim(0, 1) +
labs(colour = "Curvas", y = "Proporção de sobreviventes",
x = "Tempo")
