Read data

x= read.csv("D:\\OneDrive\\_Researches\\_Thiet bi ho tro ho hap\\Advanced Ambu bag_data\\ambu.csv")
#withdrawal because of pulse
#[1] 22483180 22483180 22483180 22483180
#[5] 22483180 22483180 22483180
#withdrawal because of systolic
#22483180 22483180 22483180 22483180 22483180
#withdrawal because of EtCO2
#22442497 22442497
#two patients who have discontinued the study (id: 22483180, 22442497)
library (dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library (magrittr)
d= x %>% filter(id != 22483180 & id != 22442497)
set.seed(123)

Data description

d$gender = as.factor(d$gender)
d$grady = as.factor (d$grady)
d$alv = as.factor (d$alv)

print(tableone::CreateTableOne(vars = c("age", "gender", "los", "grady", "dia1", "dia2", "sf", "mode", "mv","fio2", "peep", "sc", "alv"), data = d), quote = TRUE, noSpaces = TRUE)
##                                      ""
##  ""                                   "Overall"      
##   "n"                                 "27"           
##   "age (mean (SD))"                   "67.33 (10.61)"
##   "gender = 1 (%)"                    "12 (44.4)"    
##   "los (mean (SD))"                   "9.33 (6.46)"  
##   "grady (%)"                         ""             
##   "   0"                              "14 (51.9)"    
##   "   1"                              "6 (22.2)"     
##   "   2"                              "5 (18.5)"     
##   "   3"                              "2 (7.4)"      
##   "dia1 (%)"                          ""             
##   "   Hemorrhagic stroke"             "10 (37.0)"    
##   "   Ischemic stroke"                "11 (40.7)"    
##   "   Malignant hypertension"         "1 (3.7)"      
##   "   Organophosphate toxicity"       "1 (3.7)"      
##   "   Pneumonia"                      "2 (7.4)"      
##   "   Ruptured intracranial aneurysm" "1 (3.7)"      
##   "   Septic shock"                   "1 (3.7)"      
##   "dia2 (%)"                          ""             
##   "   Hemorrhagic stroke"             "2 (7.4)"      
##   "   Ischemic stroke"                "1 (3.7)"      
##   "   Pleural effusion"               "1 (3.7)"      
##   "   Pneumonia"                      "23 (85.2)"    
##   "sf (mean (SD))"                    "3.25 (0.50)"  
##   "mode (%)"                          ""             
##   "   cpap-psv"                       "10 (37.0)"    
##   "   simv-psv"                       "5 (18.5)"     
##   "   vac"                            "12 (44.4)"    
##   "mv (mean (SD))"                    "7.96 (2.02)"  
##   "fio2 (mean (SD))"                  "30.93 (4.81)" 
##   "peep (mean (SD))"                  "4.96 (0.19)"  
##   "sc (mean (SD))"                    "40.03 (13.22)"
##   "alv (%)"                           ""             
##   "   0"                              "8 (29.6)"     
##   "   1"                              "11 (40.7)"    
##   "   2"                              "7 (25.9)"     
##   "   3"                              "1 (3.7)"

Generate data from pulse

p= dplyr::select(d, c('id', "p13a", "p12a", "p11a", "p10a", "p9a", "p8a", "p7a", "p6a", "p5a", "p4a", "p3a", "p2a", "p1a", "p0a", "p1b", "p2b", "p3b", "p4b", "p5b", "p6b", "p7b", "p8b", "p9b", "p10b", "p11b", "p12b", "p13b", "p14b"))
p1= reshape2::melt(p, id= c("id"), measure.vars= c("p13a", "p12a", "p11a", "p10a", "p9a", "p8a", "p7a", "p6a", "p5a", "p4a", "p3a", "p2a", "p1a", "p0a", "p1b", "p2b", "p3b", "p4b", "p5b", "p6b", "p7b", "p8b", "p9b", "p10b", "p11b", "p12b", "p13b", "p14b"))
names (p1)[names (p1) == "variable"] = "time"
names (p1)[names (p1) == "value"] = "pulse"
p1$time = as.character(p1$time)

p1$time [p1$time == "p13a"] = 1
p1$time [p1$time == "p12a"] = 2
p1$time [p1$time == "p11a"] = 3
p1$time [p1$time == "p10a"] = 4
p1$time [p1$time == "p9a"] = 5
p1$time [p1$time == "p8a"] = 6
p1$time [p1$time == "p7a"] = 7
p1$time [p1$time == "p6a"] = 8
p1$time [p1$time == "p5a"] = 9
p1$time [p1$time == "p4a"] = 10
p1$time [p1$time == "p3a"] = 11
p1$time [p1$time == "p2a"] = 12
p1$time [p1$time == "p1a"] = 13
p1$time [p1$time == "p0a"] = 14
p1$time [p1$time == "p1b"] = 15
p1$time [p1$time == "p2b"] = 16
p1$time [p1$time == "p3b"] = 17
p1$time [p1$time == "p4b"] = 18
p1$time [p1$time == "p5b"] = 19
p1$time [p1$time == "p6b"] = 20
p1$time [p1$time == "p7b"] = 21
p1$time [p1$time == "p8b"] = 22
p1$time [p1$time == "p9b"] = 23
p1$time [p1$time == "p10b"] = 24
p1$time [p1$time == "p11b"] = 25
p1$time [p1$time == "p12b"] = 26
p1$time [p1$time == "p13b"] = 27
p1$time [p1$time == "p14b"] = 28

p1$time = as.numeric(p1$time)

p1$group [p1$time == 1| p1$time == 2|p1$time == 3|p1$time == 4|p1$time == 5|p1$time == 6|p1$time == 7|p1$time == 8|p1$time == 9|p1$time == 10|p1$time == 11|p1$time == 12|p1$time == 13|p1$time == 14] = "Conventional mechanical ventilation"
p1$group [p1$time == 15|p1$time == 16|p1$time == 17|p1$time == 18|p1$time == 19|p1$time == 20|p1$time == 21|p1$time == 22|p1$time == 23|p1$time == 24|p1$time == 25|p1$time == 26|p1$time == 27|p1$time == 28]= "Advanced Ambu bag system"

#Create pulse variable above 120 to find discontinuous cases
p1$p120 [p1$pulse <= 120 & p1$pulse >= 60] = "inclusion"
p1$p120 [p1$pulse > 120 | p1$pulse <60] = "EX"
#Look for cases with pulse > 120 or < 60 for more than 10 minutes
library (dplyr)
library (magrittr)
p1a = p1 %>% count (id, time, group, p120)
p1a$id [p1a$p120 == "EX"]
## [1] 23059529 23059529 23059529 23059529
library (ggplot2)
p1$group= factor(p1$group, levels= c("Conventional mechanical ventilation", "Advanced Ambu bag system"))
plot1= ggplot (p1, aes(x = time, y = pulse, group= id)) + geom_point(aes(x= time, y= pulse, fill= group, col= group), size=1, alpha= 0.2)+ geom_line(aes(x= time, y= pulse, fill= group, col= group)) + stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 2, color = "blue") + xlab (" ") + ylab ("Pulse rate/min") + theme(legend.position = "none")
## Warning in geom_line(aes(x = time, y = pulse, fill = group, col = group)):
## Ignoring unknown aesthetics: fill
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
table1::table1(~ pulse|group, data= p1)
Conventional mechanical ventilation
(N=378)		 Advanced Ambu bag system
(N=378)		 Overall
(N=756)
pulse
Mean (SD) 89.4 (10.7) 92.9 (13.3) 91.1 (12.2)
Median [Min, Max] 92.0 [63.0, 120] 95.0 [54.0, 120] 93.0 [54.0, 120] 
t.test (p1$pulse ~ p1$group)
## 
##  Welch Two Sample t-test
## 
## data:  p1$pulse by p1$group
## t = -4.0879, df = 721.14, p-value = 4.843e-05
## alternative hypothesis: true difference in means between group Conventional mechanical ventilation and group Advanced Ambu bag system is not equal to 0
## 95 percent confidence interval:
##  -5.310150 -1.864453
## sample estimates:
## mean in group Conventional mechanical ventilation 
##                                          89.35185 
##            mean in group Advanced Ambu bag system 
##                                          92.93915
#Impact assessment
library (lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
fit1 = lmer (pulse ~ group + (1|id), data= p1)
summary (fit1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pulse ~ group + (1 | id)
##    Data: p1
## 
## REML criterion at convergence: 5000.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1934 -0.5581 -0.0285  0.5865  3.6116 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 111.88   10.578  
##  Residual              37.52    6.125  
## Number of obs: 756, groups:  id, 27
## 
## Fixed effects:
##                               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                    89.3519     2.0599  26.6189  43.377  < 2e-16 ***
## groupAdvanced Ambu bag system   3.5873     0.4456 728.0000   8.051 3.34e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## grpAdvncAbs -0.108
#icc and 95%CI
performance::icc (fit1); performance::icc (fit1, ci= T)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.749
##   Unadjusted ICC: 0.733
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.749 [0.549, 0.854]
##   Unadjusted ICC: 0.733 [0.525, 0.839]

Generate data from Systolic

sy= dplyr::select(d, c('id', "s13a", "s12a", "s11a", "s10a", "s9a", "s8a", "s7a", "s6a", "s5a", "s4a", "s3a", "s2a", "s1a", "s0a", "s1b", "s2b", "s3b", "s4b", "s5b", "s6b", "s7b", "s8b", "s9b", "s10b", "s11b", "s12b", "s13b", "s14b"))
sy1= reshape2::melt(sy, id= c("id"), measure.vars= c("s13a", "s12a", "s11a", "s10a", "s9a", "s8a", "s7a", "s6a", "s5a", "s4a", "s3a", "s2a", "s1a", "s0a", "s1b", "s2b", "s3b", "s4b", "s5b", "s6b", "s7b", "s8b", "s9b", "s10b", "s11b", "s12b", "s13b", "s14b"))
names (sy1)[names (sy1) == "variable"] = "time"
names (sy1)[names (sy1) == "value"] = "systolic"
sy1$time = as.character(sy1$time)

sy1$time [sy1$time == "s13a"] = 1
sy1$time [sy1$time == "s12a"] = 2
sy1$time [sy1$time == "s11a"] = 3
sy1$time [sy1$time == "s10a"] = 4
sy1$time [sy1$time == "s9a"] = 5
sy1$time [sy1$time == "s8a"] = 6
sy1$time [sy1$time == "s7a"] = 7
sy1$time [sy1$time == "s6a"] = 8
sy1$time [sy1$time == "s5a"] = 9
sy1$time [sy1$time == "s4a"] = 10
sy1$time [sy1$time == "s3a"] = 11
sy1$time [sy1$time == "s2a"] = 12
sy1$time [sy1$time == "s1a"] = 13
sy1$time [sy1$time == "s0a"] = 14
sy1$time [sy1$time == "s1b"] = 15
sy1$time [sy1$time == "s2b"] = 16
sy1$time [sy1$time == "s3b"] = 17
sy1$time [sy1$time == "s4b"] = 18
sy1$time [sy1$time == "s5b"] = 19
sy1$time [sy1$time == "s6b"] = 20
sy1$time [sy1$time == "s7b"] = 21
sy1$time [sy1$time == "s8b"] = 22
sy1$time [sy1$time == "s9b"] = 23
sy1$time [sy1$time == "s10b"] = 24
sy1$time [sy1$time == "s11b"] = 25
sy1$time [sy1$time == "s12b"] = 26
sy1$time [sy1$time == "s13b"] = 27
sy1$time [sy1$time == "s14b"] = 28

sy1$time = as.numeric(sy1$time)

sy1$group [sy1$time == 1| sy1$time == 2|sy1$time == 3|sy1$time == 4|sy1$time == 5|sy1$time == 6|sy1$time == 7|sy1$time == 8|sy1$time == 9|sy1$time == 10|sy1$time == 11|sy1$time == 12|sy1$time == 13|sy1$time == 14] = "Conventional mechanical ventilation"
sy1$group [sy1$time == 15|sy1$time == 16|sy1$time == 17|sy1$time == 18|sy1$time == 19|sy1$time == 20|sy1$time == 21|sy1$time == 22|sy1$time == 23|sy1$time == 24|sy1$time == 25|sy1$time == 26|sy1$time == 27|sy1$time == 28]= "Advanced Ambu bag system"

#Create variable systolic blood pressure > 170 mmHg to find discontinuous cases
sy1$s170 [sy1$systolic <= 170] = "inclusion"
sy1$s170 [sy1$systolic > 170] = "EX"
sy1$id [sy1$s170 == "EX"]
## [1] 22452287 22418722 22452287
sy1a= sy1 %>% count (id, time, group, s170)


library (ggplot2)
sy1$group = factor(sy1$group, levels= c("Conventional mechanical ventilation", "Advanced Ambu bag system"))

plot2= ggplot (sy1, aes(x = time, y = systolic, group= id)) + geom_point(aes(x= time, y= systolic, fill= group, col= group), size=1, alpha= 0.2)+ geom_line(aes(x= time, y= systolic, fill= group, col= group)) + stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 2, color = "blue") + xlab (" ") + ylab ("SBP (mmHg)") + theme(legend.position = "none")
## Warning in geom_line(aes(x = time, y = systolic, fill = group, col = group)):
## Ignoring unknown aesthetics: fill
t.test (sy1$systolic ~ sy1$group)
## 
##  Welch Two Sample t-test
## 
## data:  sy1$systolic by sy1$group
## t = -6.5935, df = 749.54, p-value = 8.105e-11
## alternative hypothesis: true difference in means between group Conventional mechanical ventilation and group Advanced Ambu bag system is not equal to 0
## 95 percent confidence interval:
##  -10.206806  -5.523352
## sample estimates:
## mean in group Conventional mechanical ventilation 
##                                          122.2063 
##            mean in group Advanced Ambu bag system 
##                                          130.0714
table1::table1 (~ sy1$systolic| sy1$group, data= sy1)
Conventional mechanical ventilation
(N=378)		 Advanced Ambu bag system
(N=378)		 Overall
(N=756)
sy1$systolic
Mean (SD) 122 (15.8) 130 (17.0) 126 (16.9)
Median [Min, Max] 121 [84.0, 169] 130 [91.0, 188] 126 [84.0, 188] 
#Impact assessment
library (lme4)
library(lmerTest)
fit2 = lmer (systolic ~ group + (1|id), data= sy1)
summary (fit2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: systolic ~ group + (1 | id)
##    Data: sy1
## 
## REML criterion at convergence: 5793.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7375 -0.6026 -0.0023  0.5906  5.0089 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 164.7    12.83   
##  Residual             109.9    10.48   
## Number of obs: 756, groups:  id, 27
## 
## Fixed effects:
##                               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                   122.2063     2.5278  27.2245   48.34   <2e-16 ***
## groupAdvanced Ambu bag system   7.8651     0.7627 728.0000   10.31   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## grpAdvncAbs -0.151
#icc and 95%CI
performance::icc (fit2); performance::icc (fit2, ci= T)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.600
##   Unadjusted ICC: 0.568
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.600 [0.420, 0.727]
##   Unadjusted ICC: 0.568 [0.389, 0.686]

Generate data from diastolic

di= dplyr::select(d, c('id', "d13a", "d12a", "d11a", "d10a", "d9a", "d8a", "d7a", "d6a", "d5a", "d4a", "d3a", "d2a", "d1a", "d0a", "d1b", "d2b", "d3b", "d4b", "d5b", "d6b", "d7b", "d8b", "d9b", "d10b", "d11b", "d12b", "d13b", "d14b"))
di1= reshape2::melt(di, id= c("id"), measure.vars= c("d13a", "d12a", "d11a", "d10a", "d9a", "d8a", "d7a", "d6a", "d5a", "d4a", "d3a", "d2a", "d1a", "d0a", "d1b", "d2b", "d3b", "d4b", "d5b", "d6b", "d7b", "d8b", "d9b", "d10b", "d11b", "d12b", "d13b", "d14b"))
names (di1)[names (di1) == "variable"] = "time"
names (di1)[names (di1) == "value"] = "diastolic"
di1$time = as.character(di1$time)

di1$time [di1$time == "d13a"] = 1
di1$time [di1$time == "d12a"] = 2
di1$time [di1$time == "d11a"] = 3
di1$time [di1$time == "d10a"] = 4
di1$time [di1$time == "d9a"] = 5
di1$time [di1$time == "d8a"] = 6
di1$time [di1$time == "d7a"] = 7
di1$time [di1$time == "d6a"] = 8
di1$time [di1$time == "d5a"] = 9
di1$time [di1$time == "d4a"] = 10
di1$time [di1$time == "d3a"] = 11
di1$time [di1$time == "d2a"] = 12
di1$time [di1$time == "d1a"] = 13
di1$time [di1$time == "d0a"] = 14
di1$time [di1$time == "d1b"] = 15
di1$time [di1$time == "d2b"] = 16
di1$time [di1$time == "d3b"] = 17
di1$time [di1$time == "d4b"] = 18
di1$time [di1$time == "d5b"] = 19
di1$time [di1$time == "d6b"] = 20
di1$time [di1$time == "d7b"] = 21
di1$time [di1$time == "d8b"] = 22
di1$time [di1$time == "d9b"] = 23
di1$time [di1$time == "d10b"] = 24
di1$time [di1$time == "d11b"] = 25
di1$time [di1$time == "d12b"] = 26
di1$time [di1$time == "d13b"] = 27
di1$time [di1$time == "d14b"] = 28

di1$time = as.numeric(di1$time)

di1$group [di1$time == 1| di1$time == 2|di1$time == 3|di1$time == 4|di1$time == 5|di1$time == 6|di1$time == 7|di1$time == 8|di1$time == 9|di1$time == 10|di1$time == 11|di1$time == 12|di1$time == 13|di1$time == 14] = "Conventional mechanical ventilation"
di1$group [di1$time == 15|di1$time == 16|di1$time == 17|di1$time == 18|di1$time == 19|di1$time == 20|di1$time == 21|di1$time == 22|di1$time == 23|di1$time == 24|di1$time == 25|di1$time == 26|di1$time == 27|di1$time == 28]= "Advanced Ambu bag system"

di1$group = factor(di1$group, levels= c("Conventional mechanical ventilation", "Advanced Ambu bag system"))

library (ggplot2)
plot3= ggplot (di1, aes(x = time, y = diastolic, group= id)) + geom_point(aes(x= time, y= diastolic, fill= group, col= group), size=1, alpha= 0.2)+ geom_line(aes(x= time, y= diastolic, fill= group, col= group)) + stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 2, color = "blue") + xlab (" ") + ylab ("DBP (mmHg)") + theme(legend.position = "none")
## Warning in geom_line(aes(x = time, y = diastolic, fill = group, col = group)):
## Ignoring unknown aesthetics: fill
table1::table1(~ di1$diastolic | di1$group, data= di1)
Conventional mechanical ventilation
(N=378)		 Advanced Ambu bag system
(N=378)		 Overall
(N=756)
di1$diastolic
Mean (SD) 78.6 (14.6) 83.6 (16.9) 81.1 (16.0)
Median [Min, Max] 77.0 [46.0, 120] 84.0 [47.0, 119] 80.0 [46.0, 120] 
t.test (di1$diastolic ~ di1$group)
## 
##  Welch Two Sample t-test
## 
## data:  di1$diastolic by di1$group
## t = -4.3607, df = 737.51, p-value = 1.481e-05
## alternative hypothesis: true difference in means between group Conventional mechanical ventilation and group Advanced Ambu bag system is not equal to 0
## 95 percent confidence interval:
##  -7.270185 -2.756270
## sample estimates:
## mean in group Conventional mechanical ventilation 
##                                          78.55291 
##            mean in group Advanced Ambu bag system 
##                                          83.56614

Generate data from mean of blood pressure

map= merge (di1, sy1, by= c("id", "time", "group"), all.x= T, all.y= T)
map$MAP= ((2*map$diastolic) + map$systolic)/3

library (ggplot2)
plot4= ggplot (map, aes(x = time, y = MAP, group= id)) + geom_point(aes(x= time, y= MAP, fill= group, col= group), size=1, alpha= 0.2)+ geom_line(aes(x= time, y= MAP, fill= group, col= group)) + stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 2, color = "blue") + xlab (" ") + ylab ("MAP (mmHg)") + theme(legend.position = "none")
## Warning in geom_line(aes(x = time, y = MAP, fill = group, col = group)):
## Ignoring unknown aesthetics: fill
table1::table1(~ map$MAP | map$group, data= map)
Conventional mechanical ventilation
(N=378)		 Advanced Ambu bag system
(N=378)		 Overall
(N=756)
map$MAP
Mean (SD) 93.1 (13.3) 99.1 (15.4) 96.1 (14.7)
Median [Min, Max] 90.7 [63.3, 133] 98.8 [69.7, 133] 94.2 [63.3, 133] 
t.test (map$MAP ~ map$group)
## 
##  Welch Two Sample t-test
## 
## data:  map$MAP by map$group
## t = -5.7013, df = 738.99, p-value = 1.719e-08
## alternative hypothesis: true difference in means between group Conventional mechanical ventilation and group Advanced Ambu bag system is not equal to 0
## 95 percent confidence interval:
##  -8.017424 -3.910266
## sample estimates:
## mean in group Conventional mechanical ventilation 
##                                          93.10406 
##            mean in group Advanced Ambu bag system 
##                                          99.06790

Generate data from spo2

s= dplyr::select(d, c('id', "sp13a", "sp12a", "sp11a", "sp10a", "sp9a", "sp8a", "sp7a", "sp6a", "sp5a", "sp4a", "sp3a", "sp2a", "sp1a", "sp0a", "sp1b", "sp2b", "sp3b", "sp4b", "sp5b", "sp6b", "sp7b", "sp8b", "sp9b", "sp10b", "sp11b", "sp12b", "sp13b", "sp14b"))
s1= reshape2::melt(s, id= c("id"), measure.vars= c("sp13a", "sp12a", "sp11a", "sp10a", "sp9a", "sp8a", "sp7a", "sp6a", "sp5a", "sp4a", "sp3a", "sp2a", "sp1a", "sp0a", "sp1b", "sp2b", "sp3b", "sp4b", "sp5b", "sp6b", "sp7b", "sp8b", "sp9b", "sp10b", "sp11b", "sp12b", "sp13b", "sp14b"))
names (s1)[names (s1) == "variable"] = "time"
names (s1)[names (s1) == "value"] = "spo2"

s1$time = as.character(s1$time)

s1$time [s1$time == "sp13a"] = 1
s1$time [s1$time == "sp12a"] = 2
s1$time [s1$time == "sp11a"] = 3
s1$time [s1$time == "sp10a"] = 4
s1$time [s1$time == "sp9a"] = 5
s1$time [s1$time == "sp8a"] = 6
s1$time [s1$time == "sp7a"] = 7
s1$time [s1$time == "sp6a"] = 8
s1$time [s1$time == "sp5a"] = 9
s1$time [s1$time == "sp4a"] = 10
s1$time [s1$time == "sp3a"] = 11
s1$time [s1$time == "sp2a"] = 12
s1$time [s1$time == "sp1a"] = 13
s1$time [s1$time == "sp0a"] = 14
s1$time [s1$time == "sp1b"] = 15
s1$time [s1$time == "sp2b"] = 16
s1$time [s1$time == "sp3b"] = 17
s1$time [s1$time == "sp4b"] = 18
s1$time [s1$time == "sp5b"] = 19
s1$time [s1$time == "sp6b"] = 20
s1$time [s1$time == "sp7b"] = 21
s1$time [s1$time == "sp8b"] = 22
s1$time [s1$time == "sp9b"] = 23
s1$time [s1$time == "sp10b"] = 24
s1$time [s1$time == "sp11b"] = 25
s1$time [s1$time == "sp12b"] = 26
s1$time [s1$time == "sp13b"] = 27
s1$time [s1$time == "sp14b"] = 28

s1$time = as.numeric(s1$time)

s1$group [s1$time == 1| s1$time == 2|s1$time == 3|s1$time == 4|s1$time == 5|s1$time == 6|s1$time == 7|s1$time == 8|s1$time == 9|s1$time == 10|s1$time == 11|s1$time == 12|s1$time == 13|s1$time == 14] = "Conventional mechanical ventilation"
s1$group [s1$time == 15|s1$time == 16|s1$time == 17|s1$time == 18|s1$time == 19|s1$time == 20|s1$time == 21|s1$time == 22|s1$time == 23|s1$time == 24|s1$time == 25|s1$time == 26|s1$time == 27|s1$time == 28]= "Advanced Ambu bag system"

s1$group = factor(s1$group, levels= c("Conventional mechanical ventilation", "Advanced Ambu bag system"))
library (ggplot2)
plot5= ggplot (s1, aes(x = time, y = spo2, group= id)) + geom_point(aes(x= time, y= spo2, fill= group, col= group), size=1, alpha= 0.2)+ geom_line(aes(x= time, y= spo2, fill= group, col= group)) + stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 2, color = "blue") + xlab ("Times") + ylab ("SpO2 (%)") + theme(legend.position = "none")
## Warning in geom_line(aes(x = time, y = spo2, fill = group, col = group)):
## Ignoring unknown aesthetics: fill
table1::table1(~ s1$spo2 | s1$group, data= s1)
Conventional mechanical ventilation
(N=378)		 Advanced Ambu bag system
(N=378)		 Overall
(N=756)
s1$spo2
Mean (SD) 98.2 (1.82) 98.8 (1.41) 98.5 (1.65)
Median [Min, Max] 99.0 [88.0, 100] 99.0 [93.0, 100] 99.0 [88.0, 100] 
t.test (s1$spo2 ~ s1$group)
## 
##  Welch Two Sample t-test
## 
## data:  s1$spo2 by s1$group
## t = -4.9372, df = 708.47, p-value = 9.894e-07
## alternative hypothesis: true difference in means between group Conventional mechanical ventilation and group Advanced Ambu bag system is not equal to 0
## 95 percent confidence interval:
##  -0.8171516 -0.3521605
## sample estimates:
## mean in group Conventional mechanical ventilation 
##                                          98.20899 
##            mean in group Advanced Ambu bag system 
##                                          98.79365
#impact assessment
library (lme4)
library(lmerTest)
options (scipen = 999)
fit4 = lmer (spo2 ~ group + (1|id), data= s1)
summary (fit4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: spo2 ~ group + (1 | id)
##    Data: s1
## 
## REML criterion at convergence: 2482.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.3564 -0.4107  0.0903  0.5438  2.6360 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 1.312    1.145   
##  Residual             1.384    1.176   
## Number of obs: 756, groups:  id, 27
## 
## Fixed effects:
##                                Estimate Std. Error        df t value
## (Intercept)                    98.20899    0.22860  27.92013 429.618
## groupAdvanced Ambu bag system   0.58466    0.08556 728.00000   6.833
##                                           Pr(>|t|)    
## (Intercept)                   < 0.0000000000000002 ***
## groupAdvanced Ambu bag system      0.0000000000175 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## grpAdvncAbs -0.187
#icc and 95%CI
set.seed(123)
performance::icc (fit4); performance::icc (fit4, ci= T)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.487
##   Unadjusted ICC: 0.472
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.487 [0.256, 0.660]
##   Unadjusted ICC: 0.472 [0.242, 0.638]
#Find cases of prolonged hypoxia.
s1$s90 [s1$spo2 < 90] = "EX"
s1$s90 [s1$spo2 >= 90] = "inclusion"
s1$id [s1$s90 == "EX"]
## [1] 23000205 22446528
s1a= s1 %>% count (id, time, group, s90)

Generate data from etco2

e= dplyr::select(d, c("id","e13a", "e12a", "e11a", "e10a", "e9a", "e8a", "e7a", "e6a", "e5a", "e4a", "e3a", "e2a", "e1a", "e0a", "e1b", "e2b", "e3b", "e4b", "e5b", "e6b", "e7b", "e8b", "e9b", "e10b", "e11b", "e12b", "e13b", "e14b"))
e1= reshape2::melt(e, id= c("id"), measure.vars= c("e13a", "e12a", "e11a", "e10a", "e9a", "e8a", "e7a", "e6a", "e5a", "e4a", "e3a", "e2a", "e1a", "e0a", "e1b", "e2b", "e3b", "e4b", "e5b", "e6b", "e7b", "e8b", "e9b", "e10b", "e11b", "e12b", "e13b", "e14b"))
names (e1) [names(e1) == "variable"] = "time"
names (e1) [names(e1) == "value"] = "etco2"

e1$time = as.character(e1$time)

e1$time [e1$time == "e13a"] = 1
e1$time [e1$time == "e12a"] = 2
e1$time [e1$time == "e11a"] = 3
e1$time [e1$time == "e10a"] = 4
e1$time [e1$time == "e9a"] = 5
e1$time [e1$time == "e8a"] = 6
e1$time [e1$time == "e7a"] = 7
e1$time [e1$time == "e6a"] = 8
e1$time [e1$time == "e5a"] = 9
e1$time [e1$time == "e4a"] = 10
e1$time [e1$time == "e3a"] = 11
e1$time [e1$time == "e2a"] = 12
e1$time [e1$time == "e1a"] = 13
e1$time [e1$time == "e0a"] = 14
e1$time [e1$time == "e1b"] = 15
e1$time [e1$time == "e2b"] = 16
e1$time [e1$time == "e3b"] = 17
e1$time [e1$time == "e4b"] = 18
e1$time [e1$time == "e5b"] = 19
e1$time [e1$time == "e6b"] = 20
e1$time [e1$time == "e7b"] = 21
e1$time [e1$time == "e8b"] = 22
e1$time [e1$time == "e9b"] = 23
e1$time [e1$time == "e10b"] = 24
e1$time [e1$time == "e11b"] = 25
e1$time [e1$time == "e12b"] = 26
e1$time [e1$time == "e13b"] = 27
e1$time [e1$time == "e14b"] = 28

e1$time = as.numeric(e1$time)

e1$Ventilation_therapy [e1$time == 1|e1$time == 2|e1$time == 3|e1$time == 4|e1$time == 5|e1$time == 6|e1$time == 7|e1$time == 8|e1$time == 9|e1$time == 10|e1$time == 11|e1$time == 12|e1$time == 13|e1$time == 14]= "Conventional mechanical ventilation"
e1$Ventilation_therapy [e1$time == 15|e1$time == 16|e1$time == 17|e1$time == 18|e1$time == 19|e1$time == 20|e1$time == 21|e1$time == 22|e1$time == 23|e1$time == 24|e1$time == 25|e1$time == 26|e1$time == 27|e1$time == 28] = "Advanced Ambu bag system"

e1$Ventilation_therapy = factor(e1$Ventilation_therapy, levels= c("Conventional mechanical ventilation", "Advanced Ambu bag system"))

library(ggplot2)
plot6= ggplot (e1, aes(x = time, y = etco2, group= id)) + geom_point(aes(x= time, y= etco2, fill= Ventilation_therapy, col= Ventilation_therapy), size=1, alpha= 0.2)+ geom_line(aes(x= time, y= etco2, fill= Ventilation_therapy, col= Ventilation_therapy)) + stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 2, color = "blue") + xlab("Times") + ylab ("EtCO2 (mmHg)") + theme(legend.position = "none")
## Warning in geom_line(aes(x = time, y = etco2, fill = Ventilation_therapy, :
## Ignoring unknown aesthetics: fill
plot7= ggplot (e1, aes(x = time, y = etco2, group= id)) + geom_point(aes(x= time, y= etco2, fill= Ventilation_therapy, col= Ventilation_therapy), size=1, alpha= 0.2)+ geom_line(aes(x= time, y= etco2, fill= Ventilation_therapy, col= Ventilation_therapy)) + stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 2, color = "blue") + xlab("Times") + ylab ("EtCO2 (mmHg)") + theme(legend.position = "bottom")
## Warning in geom_line(aes(x = time, y = etco2, fill = Ventilation_therapy, :
## Ignoring unknown aesthetics: fill
#Impact assessment
library (lme4)
library(lmerTest)
fit5 = lmer (etco2 ~ Ventilation_therapy + (1|id), data= e1)
summary (fit5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: etco2 ~ Ventilation_therapy + (1 | id)
##    Data: e1
## 
## REML criterion at convergence: 4282.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0140 -0.6014  0.1050  0.6368  4.9734 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 10.46    3.235   
##  Residual             15.22    3.901   
## Number of obs: 756, groups:  id, 27
## 
## Fixed effects:
##                                             Estimate Std. Error       df
## (Intercept)                                  31.9153     0.6540  28.6290
## Ventilation_therapyAdvanced Ambu bag system  -4.9233     0.2838 728.0000
##                                             t value            Pr(>|t|)    
## (Intercept)                                   48.80 <0.0000000000000002 ***
## Ventilation_therapyAdvanced Ambu bag system  -17.35 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Vntltn_AAbs -0.217
#icc and 95%CI
set.seed(123)
performance::icc (fit5); performance::icc (fit5, ci= T)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.407
##   Unadjusted ICC: 0.330
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.407 [0.192, 0.588]
##   Unadjusted ICC: 0.330 [0.143, 0.506]
table1::table1 (~ etco2| Ventilation_therapy, data= e1)
Conventional mechanical ventilation
(N=378)		 Advanced Ambu bag system
(N=378)		 Overall
(N=756)
etco2
Mean (SD) 31.9 (5.07) 27.0 (4.99) 29.5 (5.60)
Median [Min, Max] 31.0 [15.0, 51.0] 28.0 [9.00, 48.0] 30.0 [9.00, 51.0] 
t.test (e1$etco2 ~ e1$Ventilation_therapy)
## 
##  Welch Two Sample t-test
## 
## data:  e1$etco2 by e1$Ventilation_therapy
## t = 13.451, df = 753.84, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group Conventional mechanical ventilation and group Advanced Ambu bag system is not equal to 0
## 95 percent confidence interval:
##  4.204729 5.641832
## sample estimates:
## mean in group Conventional mechanical ventilation 
##                                          31.91534 
##            mean in group Advanced Ambu bag system 
##                                          26.99206
e1$etex [e1$etco2 > 45 | e1$etco2 < 15] = "EX"
e1$etex [e1$etco2 <= 45 & e1$etco2 >= 15]= "inclusion"
e1$id [e1$etex == "exclusion"]
## integer(0)
e1a= e1 %>% count (id, time, Ventilation_therapy, etex)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(cowplot)
grid_arrange <- plot_grid(plot1, plot2, plot3, plot4,plot5, plot6, ncol = 2)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
final_plot <- plot_grid(
  grid_arrange,
  get_legend(plot7), 
  ncol = 1,  
  rel_heights = c(1, 0.2) 
)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#final_plot
print(final_plot)