R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

Version history

-Update from the first version (Jan to Oct): –changed data file from CSV to xlsx; –Removed dummy time zero (theoretical 450ml of study fluid) because it has been considered misleading; –Added code to help adding new formulas for estimating gastric emptying; –Added mean and sd of baseline data; –RMANOVA was divided for each pair of groups analysis.

Additional codes used in this analysis

read.xlsx2 = function(bdata){
  nv=length(bdata)
  ischar=rep(FALSE,nv)
  for(i in c(1:nv)){
    if(is.character(bdata[,i][0])){
      bdata[,i]=factor(bdata[,i])
    }
  }
  bdata
}
sm = function (var,group){
  print("Summary statistics for ")
print("Median (min - max)")
print(paste(round(median(var),2),"(",round(min(var),2),"-",round(max(var),2),")"))
print("Mean(sd)")
print(paste(round(mean(var),2),"(",round(sd(var),2),")"))
}
knitr::opts_chunk$set(echo = TRUE)
library(ez)
library(ggplot2)
library(lattice)
library(openxlsx)
library(plyr)
# Data can be downloaded from https://data.mendeley.com/datasets/yj333czkbw/2 
bda=bd=read.xlsx2(read.xlsx("D:\\DataAnalysis\\NCT03280706\\NCT03280706.xlsx"))

Antral volume (in mm²) to residual gastric volume (in ml) formulas

bda$residual.volume.Perlas=27+14.6*bda$antral.area-1.28*bda$age
bda$residual.volume.Perlas.ml.kg=bda$residual.volume.Perlas/bda$weight
bda$fasted.volume.Perlas=27+14.6*bda$fasted.antral.area-1.28*bda$age
bda$fasted.volume.Perlas.ml.kg=bda$fasted.volume.Perlas/bda$weight

Subsets

# restrict to per protocol analysis
bdpa=subset(bda,bda$selection=="included")
# for non-repeated-measures analysis
bdpa2=subset(bdpa,bdpa$t==5)
oj=subset(bdpa2, bdpa2$study.group=="OJ")
m=subset(bdpa2, bdpa2$study.group=="M")
cm=subset(bdpa2, bdpa2$study.group=="CM")

Data summary at time 5

summary(bdpa2)
##        id           selection 
##  Min.   : 1.00   excluded: 0  
##  1st Qu.:14.25   included:54  
##  Median :27.50                
##  Mean   :29.13                
##  3rd Qu.:40.75                
##  Max.   :62.00                
##                               
##                              selection.reason study.group      age       
##  cesarean indicated before completion: 0      CM:18       Min.   :18.00  
##  early birth                         : 0      M :18       1st Qu.:21.25  
##  ingested another solution           : 0      OJ:18       Median :25.50  
##  vomited                             : 0                  Mean   :26.56  
##  NA's                                :54                  3rd Qu.:30.00  
##                                                           Max.   :42.00  
##                                                                          
##      weight          height      gestational.age    parityG     
##  Min.   :53.30   Min.   :1.480   Min.   :24.42   Min.   :1.000  
##  1st Qu.:65.33   1st Qu.:1.573   1st Qu.:38.03   1st Qu.:1.000  
##  Median :71.15   Median :1.600   Median :38.92   Median :2.000  
##  Mean   :72.16   Mean   :1.605   Mean   :38.52   Mean   :2.093  
##  3rd Qu.:76.97   3rd Qu.:1.637   3rd Qu.:39.85   3rd Qu.:2.750  
##  Max.   :99.00   Max.   :1.720   Max.   :41.57   Max.   :7.000  
##                                                                 
##     parityP          parityC          parityA       labor.induction labor 
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   N:38            S:54  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   S:16                  
##  Median :0.0000   Median :0.0000   Median :0.0000                         
##  Mean   :0.7037   Mean   :0.1296   Mean   :0.2593                         
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000                         
##  Max.   :6.0000   Max.   :1.0000   Max.   :2.0000                         
##                                                                           
##   labor.hours                    morbidity  labor.analgesia
##  Min.   : 1.000   ArterialHypertension: 4   N:51           
##  1st Qu.: 5.000   N                   :50   S: 3           
##  Median : 6.000                                            
##  Mean   : 6.375                                            
##  3rd Qu.: 9.000                                            
##  Max.   :12.000                                            
##  NA's   :46                                                
##  advese.outcomes       t      antral.area     fasted.antral.area
##  N     :53       Min.   :5   Min.   : 8.629   Min.   :0.8973    
##  nausea: 1       1st Qu.:5   1st Qu.:12.399   1st Qu.:1.2966    
##                  Median :5   Median :19.234   Median :1.5068    
##                  Mean   :5   Mean   :18.667   Mean   :1.6165    
##                  3rd Qu.:5   3rd Qu.:24.070   3rd Qu.:1.8260    
##                  Max.   :5   Max.   :29.297   Max.   :3.8863    
##                                                                 
##   time.5.area      time.30.area     time.60.area      time.90.area   
##  Min.   : 8.629   Min.   : 2.629   Min.   : 0.8274   Min.   :0.7877  
##  1st Qu.:12.399   1st Qu.: 4.260   1st Qu.: 1.9685   1st Qu.:1.3452  
##  Median :19.234   Median :10.453   Median : 5.3240   Median :2.5897  
##  Mean   :18.667   Mean   : 9.873   Mean   : 4.8970   Mean   :2.5161  
##  3rd Qu.:24.070   3rd Qu.:13.670   3rd Qu.: 5.8671   3rd Qu.:3.3510  
##  Max.   :29.297   Max.   :18.845   Max.   :12.8384   Max.   :5.3452  
##                                                                      
##  time.120.area    residual.volume.Perlas residual.volume.Perlas.ml.kg
##  Min.   :0.7685   Min.   :119.7          Min.   :1.826               
##  1st Qu.:0.9548   1st Qu.:179.3          1st Qu.:2.568               
##  Median :1.5363   Median :272.9          Median :3.693               
##  Mean   :1.6167   Mean   :265.5          Mean   :3.748               
##  3rd Qu.:2.0969   3rd Qu.:347.9          3rd Qu.:4.766               
##  Max.   :4.2493   Max.   :431.7          Max.   :7.574               
##                                                                      
##  fasted.volume.Perlas fasted.volume.Perlas.ml.kg
##  Min.   : 0.00        Min.   :0.0000            
##  1st Qu.: 9.99        1st Qu.:0.1292            
##  Median :15.35        Median :0.2217            
##  Mean   :16.61        Mean   :0.2400            
##  3rd Qu.:23.02        3rd Qu.:0.3116            
##  Max.   :60.70        Max.   :1.0649            
## 

Data distribution between groups

aggregate(age~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group  age.Min. age.1st Qu. age.Median  age.Mean age.3rd Qu.
## 1          CM 18.000000   21.000000  24.500000 26.333333   31.500000
## 2           M 18.000000   22.500000  27.000000 27.333333   30.000000
## 3          OJ 19.000000   21.250000  24.000000 26.000000   30.750000
##    age.Max.    age.sd
## 1 41.000000  7.129062
## 2 42.000000  6.543969
## 3 35.000000  5.335784
aggregate(weight~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group weight.Min. weight.1st Qu. weight.Median weight.Mean
## 1          CM   57.000000      64.500000     70.025000   72.000000
## 2           M   53.300000      66.250000     71.650000   72.250000
## 3          OJ   59.300000      68.000000     72.000000   72.216667
##   weight.3rd Qu. weight.Max. weight.sd
## 1      77.450000   99.000000 11.183036
## 2      76.400000   98.000000 10.100393
## 3      76.550000   89.000000  8.455506
aggregate(height~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group height.Min. height.1st Qu. height.Median height.Mean
## 1          CM  1.49000000     1.58000000    1.59500000  1.60500000
## 2           M  1.48000000     1.57250000    1.60000000  1.60388889
## 3          OJ  1.50000000     1.56250000    1.62000000  1.60555556
##   height.3rd Qu. height.Max.  height.sd
## 1     1.62000000  1.72000000 0.05933455
## 2     1.63500000  1.72000000 0.05637468
## 3     1.64500000  1.69000000 0.05627895
aggregate(gestational.age~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group gestational.age.Min. gestational.age.1st Qu.
## 1          CM            37.000000               38.035000
## 2           M            24.420000               38.070000
## 3          OJ            25.420000               38.210000
##   gestational.age.Median gestational.age.Mean gestational.age.3rd Qu.
## 1              38.565000            38.790000               39.532500
## 2              38.850000            38.266111               40.030000
## 3              39.280000            38.494444               39.962500
##   gestational.age.Max. gestational.age.sd
## 1            40.870000           1.120688
## 2            41.570000           3.721232
## 3            40.850000           3.442585
aggregate(parityG~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group parityG.Min. parityG.1st Qu. parityG.Median parityG.Mean
## 1          CM    1.0000000       1.0000000      2.0000000    1.8333333
## 2           M    1.0000000       1.0000000      2.0000000    2.1666667
## 3          OJ    1.0000000       1.0000000      2.0000000    2.2777778
##   parityG.3rd Qu. parityG.Max. parityG.sd
## 1       2.0000000    4.0000000  0.9235481
## 2       2.7500000    7.0000000  1.5811388
## 3       2.7500000    7.0000000  1.7083034
aggregate(parityP~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group parityP.Min. parityP.1st Qu. parityP.Median parityP.Mean
## 1          CM    0.0000000       0.0000000      0.0000000    0.5555556
## 2           M    0.0000000       0.0000000      0.0000000    0.8888889
## 3          OJ    0.0000000       0.0000000      0.0000000    0.6666667
##   parityP.3rd Qu. parityP.Max. parityP.sd
## 1       1.0000000    2.0000000  0.7838234
## 2       1.0000000    6.0000000  1.5296631
## 3       1.0000000    4.0000000  1.0846523
aggregate(parityC~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group parityC.Min. parityC.1st Qu. parityC.Median parityC.Mean
## 1          CM    0.0000000       0.0000000      0.0000000    0.1111111
## 2           M    0.0000000       0.0000000      0.0000000    0.1666667
## 3          OJ    0.0000000       0.0000000      0.0000000    0.1111111
##   parityC.3rd Qu. parityC.Max. parityC.sd
## 1       0.0000000    1.0000000  0.3233808
## 2       0.0000000    1.0000000  0.3834825
## 3       0.0000000    1.0000000  0.3233808
aggregate(parityA~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group parityA.Min. parityA.1st Qu. parityA.Median parityA.Mean
## 1          CM    0.0000000       0.0000000      0.0000000    0.1666667
## 2           M    0.0000000       0.0000000      0.0000000    0.1111111
## 3          OJ    0.0000000       0.0000000      0.0000000    0.5000000
##   parityA.3rd Qu. parityA.Max. parityA.sd
## 1       0.0000000    2.0000000  0.5144958
## 2       0.0000000    1.0000000  0.3233808
## 3       1.0000000    2.0000000  0.7859052
aggregate(labor.hours~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group labor.hours.Min. labor.hours.1st Qu. labor.hours.Median
## 1          CM         6.000000            6.000000           6.000000
## 2           M         1.000000            1.000000           1.000000
## 3          OJ         2.000000            5.500000           9.000000
##   labor.hours.Mean labor.hours.3rd Qu. labor.hours.Max. labor.hours.sd
## 1         6.750000            6.750000         9.000000       1.500000
## 2         1.000000            1.000000         1.000000             NA
## 3         7.666667           10.500000        12.000000       5.131601
aggregate(fasted.antral.area~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group fasted.antral.area.Min. fasted.antral.area.1st Qu.
## 1          CM               0.8972603                  1.2924658
## 2           M               0.9479452                  1.2116438
## 3          OJ               0.9849315                  1.3034247
##   fasted.antral.area.Median fasted.antral.area.Mean
## 1                 1.6191781               1.7021309
## 2                 1.4760274               1.4846271
## 3                 1.5520548               1.6627093
##   fasted.antral.area.3rd Qu. fasted.antral.area.Max. fasted.antral.area.sd
## 1                  2.0061644               3.8863014             0.7181185
## 2                  1.6904110               2.2986301             0.3565778
## 3                  1.8315068               2.8972603             0.4800739
aggregate(time.5.area~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group time.5.area.Min. time.5.area.1st Qu. time.5.area.Median
## 1          CM        12.298630           18.204452          19.897945
## 2           M         8.628767           10.749315          11.731507
## 3          OJ        18.900000           22.745205          24.384247
##   time.5.area.Mean time.5.area.3rd Qu. time.5.area.Max. time.5.area.sd
## 1        20.465677           21.073973        29.297260       4.322567
## 2        11.695129           12.572603        14.273973       1.452768
## 3        23.838965           25.650000        27.698630       2.726833
aggregate(time.30.area~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group time.30.area.Min. time.30.area.1st Qu. time.30.area.Median
## 1          CM          6.668493            13.123973           13.819863
## 2           M          2.628767             2.934932            3.226027
## 3          OJ          7.689041            10.397603           11.698630
##   time.30.area.Mean time.30.area.3rd Qu. time.30.area.Max. time.30.area.sd
## 1         13.466591            14.190068         18.845205        2.909226
## 2          3.809209             3.910616          7.846575        1.416304
## 3         12.343798            13.470890         17.302740        2.906101
aggregate(time.60.area~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group time.60.area.Min. time.60.area.1st Qu. time.60.area.Median
## 1          CM         5.0945205            5.6390411           5.8561644
## 2           M         0.8273973            1.3736301           1.4952055
## 3          OJ         3.5191781            5.2832192           5.3815068
##   time.60.area.Mean time.60.area.3rd Qu. time.60.area.Max. time.60.area.sd
## 1         6.6705479            6.7335616        12.8383562       2.1660518
## 2         1.7369102            1.9325342         4.0068493       0.7179637
## 3         6.2836377            5.9770548        10.6246575       2.1479291
aggregate(time.90.area~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group time.90.area.Min. time.90.area.1st Qu. time.90.area.Median
## 1          CM         2.7465753            3.1739726           3.5424658
## 2           M         0.7876712            0.8972603           1.0479452
## 3          OJ         1.7986301            2.2140411           2.5897260
##   time.90.area.Mean time.90.area.3rd Qu. time.90.area.Max. time.90.area.sd
## 1         3.6579909            4.1027397         5.3452055       0.6501418
## 2         1.1471081            1.2869863         1.8328767       0.3394416
## 3         2.7431507            3.0027397         4.3684932       0.6727453
aggregate(time.120.area~study.group,FUN = function(i) c(summary(i, na.rm=TRUE),sd=(sd(i,na.rm=TRUE))),data=bdpa2)
##   study.group time.120.area.Min. time.120.area.1st Qu.
## 1          CM          1.2876712             2.0962329
## 2           M          0.7684932             0.8493151
## 3          OJ          0.9410959             1.2167808
##   time.120.area.Median time.120.area.Mean time.120.area.3rd Qu.
## 1            2.1801370          2.3063927             2.4996575
## 2            0.8726027          1.0310502             0.9739726
## 3            1.5363014          1.5127854             1.7465753
##   time.120.area.Max. time.120.area.sd
## 1          4.2493151        0.6369628
## 2          1.8328767        0.3433957
## 3          2.0972603        0.3638995

Hypothesis test

Main outcome

Subseting for data analysis

cmoj=droplevels(subset(bdpa,bdpa$study.group!="M"))
cmm=droplevels(subset(bdpa,bdpa$study.group!="OJ"))
ojm=droplevels(subset(bdpa,bdpa$study.group!="CM"))

Maltodextrin and Pulp-Free Orange Juice comparison

ezANOVA(data=ojm,dv=antral.area,wid=factor(id),between=study.group,within=t)
## Warning: Converting "id" to factor for ANOVA.
## Warning: "t" will be treated as numeric.
## Warning: There is at least one numeric within variable, therefore aov()
## will be used for computation and no assumption checks will be obtained.
## $ANOVA
##          Effect DFn DFd         F            p p<.05       ges
## 1   study.group   1  34  199.0169 9.010976e-16     * 0.7965497
## 2             t   1  34 1591.3384 3.858672e-30     * 0.9393870
## 3 study.group:t   1  34  244.8965 4.189734e-17     * 0.7045840

Maltodextrin and Coffee with Milk comparison

ezANOVA(data=cmm,dv=antral.area,wid=factor(id),between=study.group,within=t)
## Warning: Converting "id" to factor for ANOVA.
## Warning: "t" will be treated as numeric.
## Warning: There is at least one numeric within variable, therefore aov()
## will be used for computation and no assumption checks will be obtained.
## $ANOVA
##          Effect DFn DFd         F            p p<.05       ges
## 1   study.group   1  34 138.82363 1.491273e-13     * 0.7204511
## 2             t   1  34 761.38156 7.367411e-25     * 0.8919959
## 3 study.group:t   1  34  79.24876 2.099244e-10     * 0.4622593

Coffee with Milk and Pulp-Free Orange juice comparison

ezANOVA(data=cmoj,dv=antral.area,wid=factor(id),between=study.group,within=t)
## Warning: Converting "id" to factor for ANOVA.
## Warning: "t" will be treated as numeric.
## Warning: There is at least one numeric within variable, therefore aov()
## will be used for computation and no assumption checks will be obtained.
## $ANOVA
##          Effect DFn DFd            F            p p<.05          ges
## 1   study.group   1  34 8.106373e-04 9.774524e-01       1.510557e-05
## 2             t   1  34 1.095666e+03 1.880455e-27     * 9.219257e-01
## 3 study.group:t   1  34 6.971124e+00 1.241844e-02     * 6.987981e-02

Figures (intention to treat)

Main outcome

bdse=ddply(bdpa,.(study.group,t),summarize, sd=sd(antral.area), antral.area=mean(antral.area))

p=ggplot(bdse,aes(x=t,y=antral.area,group=study.group,color=study.group))+geom_line()+geom_point(size=5)+geom_errorbar(aes(ymin=antral.area-sd, ymax=antral.area+sd),width=0.2,position=position_dodge(0.05))
p= p + labs(x = "Time in minutes", y = "Antral area in mm²", colour = "Study group")
print(p)

bdse=ddply(bdpa,.(study.group,t),summarize, min=min(residual.volume.Perlas.ml.kg),max=max(residual.volume.Perlas.ml.kg),sd=sd(residual.volume.Perlas.ml.kg,na.rm = TRUE), residual.volume.Perlas.ml.kg=mean(residual.volume.Perlas.ml.kg))

p=ggplot(bdse,aes(x=t,y=bdse$residual.volume.Perlas.ml.kg,group=bdse$study.group,color=bdse$study.group))+geom_line()+geom_point(size=5)+geom_errorbar(aes(ymin=residual.volume.Perlas.ml.kg-sd, ymax=residual.volume.Perlas.ml.kg+sd),width=0.2,position=position_dodge(0.05))
p= p + labs(x = "Time in minutes", y = "Estimated residual gastric volume in ml/kg(Perlas)", colour = "Study group")
p=p+ geom_hline(yintercept = 1.5) + annotate("text", min(0), 1.5, vjust = -1, label = "    1.5ml/kg")
p=p+ geom_hline(yintercept = 0.8) + 
    annotate("text", min(0), 0.8, vjust = -1, label = "    0.8ml/kg")

print(p)

Fitted curves

p=ggplot(data = bdpa, aes(x = bdpa$t, y = bdpa$antral.area, group = bdpa$id, color=bdpa$study.group)) + labs(
    x = "Time in minutes",
    y = "Antral area in mm2",
    colour = "Study group"
   )
p=p + stat_smooth(aes(group = bdpa$study.group))
p=p+ theme(legend.background = element_rect(fill=rgb(1,1,1), size=0.5, linetype=2));
print(p)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

p=ggplot(data = bdpa, aes(x = bdpa$t, y = bdpa$residual.volume.Perlas.ml.kg, group = bdpa$id, color=bdpa$study.group)) + labs(
    x = "Time in minutes",
    y = "Estimated antral volume in ml/kg (Perlas formula)",
    colour = "Study group"
   )
p=p + stat_smooth(aes(group = bdpa$study.group))
p=p+ geom_hline(yintercept = 1.5) + annotate("text", min(0), 1.5, vjust = -1, label = "    1.5ml/kg")
p=p+ geom_hline(yintercept = 0.8) + 
    annotate("text", min(0), 0.8, vjust = -1, label = "    0.8ml/kg")
p=p+ theme(legend.background = element_rect(fill=rgb(1,1,1), size=0.5, linetype=2));
print(p)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'