5-1 HoRM{datasets}

install.packages(“HoRM”)

library(HoRM)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

HoRM: Supplemental Functions and Datasets for “Handbook of Regression Methods”

data(BAC, package="HoRM")
dta <- reshape(data.frame(BAC, id=1:15), idvar='id', varying=list(1:2), 
                  direction = 'long', timevar="test", v.names="bac")
dta
##      id test   bac
## 1.1   1    1 0.160
## 2.1   2    1 0.170
## 3.1   3    1 0.180
## 4.1   4    1 0.100
## 5.1   5    1 0.170
## 6.1   6    1 0.100
## 7.1   7    1 0.060
## 8.1   8    1 0.100
## 9.1   9    1 0.170
## 10.1 10    1 0.056
## 11.1 11    1 0.111
## 12.1 12    1 0.162
## 13.1 13    1 0.143
## 14.1 14    1 0.079
## 15.1 15    1 0.006
## 1.2   1    2 0.145
## 2.2   2    2 0.156
## 3.2   3    2 0.181
## 4.2   4    2 0.108
## 5.2   5    2 0.180
## 6.2   6    2 0.112
## 7.2   7    2 0.081
## 8.2   8    2 0.104
## 9.2   9    2 0.176
## 10.2 10    2 0.048
## 11.2 11    2 0.092
## 12.2 12    2 0.144
## 13.2 13    2 0.121
## 14.2 14    2 0.065
## 15.2 15    2 0.000

DT函數datatable()–視覺化呈現

DT::datatable(dta[, c(2,3,1)], rownames=FALSE, fillContainer=FALSE, options=list(pageLength=9))

5-2 Summary statistics

aggregate(bac ~ test, data=dta, FUN=mean) |> knitr::kable()
test bac
1 0.1178
2 0.1142 
aggregate(bac ~ test, data=dta, FUN=sd) |> knitr::kable()
test bac
1 0.0524911
2 0.0519810

5-3 Design plot

dplyr 套件中的 mutate() 函數，新增變數的資料框，也可更改規則 bty6種邊框 bty=“u”繪製圖形的左邊框、下邊框和右邊框；

library(dplyr)
## 
## 載入套件：'dplyr'
## 下列物件被遮斷自 'package:stats':
## 
##     filter, lag
## 下列物件被遮斷自 'package:base':
## 
##     intersect, setdiff, setequal, union
dta <- dta |> mutate(id = as.factor(id),
                     test = as.factor(test))

dta |> plot.design (bty='n',
                    ylab="Mean of BAC")
grid()

5-4 Subject by group effect

with(dta, 
     interaction.plot(id, 
                      test, 
                      bac, 
                      bty='n', 
                      ylab="BAC"))
grid()
i <- 1:15
segments(as.numeric(dta$id[i]), 
         dta$bac[i], 
         as.numeric(dta$id[i]), 
         dta$bac[i+15],
         col=1, 
         lty=3, 
         lwd=1.2)

5-5 Order subject by response means

x <- aggregate(bac~ id, data=dta, max)
x[order(x$bac),] |> t() |> knitr::kable()
15 10 14 7 8 4 11 6 13 1 12 2 9 5 3
id 15 10 14 7 8 4 11 6 13 1 12 2 9 5 3
bac 0.006 0.056 0.079 0.081 0.104 0.108 0.111 0.112 0.143 0.160 0.162 0.170 0.176 0.180 0.181 
dta$id <- factor(dta$id, levels(dta$id)[x[order(x$bac),1]])

5-6 Profile (of means) plot

library(lattice)
xyplot(id ~ bac, 
       groups=test, 
       data=dta, 
       lb=dta$bac[dta$test==1],
       ub=dta$bac[dta$test==2],
       panel=function(x, y, lb, ub, ...){
         panel.segments(lb, y, ub, y, col="plum")
         panel.points(x, y, pch=1, col=as.numeric(dta$test)+4)
       },
       ylab="Subject ID",
       xlab="BAC",
       par.settings=list(superpose.symbol=list(pch=1,
                                               col=c("#28E2E5", "#CD0BBC"))),
       auto.key=list(space='top', columns=2))

5-7 Slope plot

pacman::p_load(PairedData)
dta2 <- with(dta, 
             data.frame(id=id[test==1],
                        T1=bac[test==1],
                        T2=bac[test==2]))
paired.plotProfiles(dta2,
                    "T1", "T2",
                    subjects="id")+ 
  labs(x="Test", 
       y="BAC")+
  scale_y_continuous(limits=c(0.00, 0.2),
                     breaks=seq(0.00, 0.2, by=0.01))+
  geom_point()+
  theme_classic()

5-8 Data vectors

t.test(x=subset(dta, test=='1')$bac,
       y=subset(dta, test=='2')$bac, 
       paired=TRUE)
## 
##  Paired t-test
## 
## data:  subset(dta, test == "1")$bac and subset(dta, test == "2")$bac
## t = 1.0447, df = 14, p-value = 0.3139
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.003790739  0.010990739
## sample estimates:
## mean of the differences 
##                  0.0036
with(dta,
     t.test(bac[test==1],
            bac[test==2], 
            pair=T))
## 
##  Paired t-test
## 
## data:  bac[test == 1] and bac[test == 2]
## t = 1.0447, df = 14, p-value = 0.3139
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.003790739  0.010990739
## sample estimates:
## mean of the differences 
##                  0.0036

5-9 Formula approach

t.test(bac ~ test, pair=T, data=dta)
## 
##  Paired t-test
## 
## data:  bac by test
## t = 1.0447, df = 14, p-value = 0.3139
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.003790739  0.010990739
## sample estimates:
## mean of the differences 
##                  0.0036

5-10 ANOVA

aov(bac ~ id + test, data=dta) |> anova()
## Analysis of Variance Table
## 
## Response: bac
##           Df   Sum Sq   Mean Sq F value    Pr(>F)    
## id        14 0.075156 0.0053683 60.2791 4.851e-10 ***
## test       1 0.000097 0.0000972  1.0914    0.3139    
## Residuals 14 0.001247 0.0000891                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5-11 Linear models

lm(bac ~ id + test, data=dta) |> anova()
## Analysis of Variance Table
## 
## Response: bac
##           Df   Sum Sq   Mean Sq F value    Pr(>F)    
## id        14 0.075156 0.0053683 60.2791 4.851e-10 ***
## test       1 0.000097 0.0000972  1.0914    0.3139    
## Residuals 14 0.001247 0.0000891                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5-12 Difference plot

pacman::p_load(multicon)
with(dta,
  diffPlot(bac[test==1],
           bac[test==2],
           paired=TRUE,
           grp.names=c("1", "2"), 
           xlab="", 
           ylab="BAC", 
           main="Difference of paired means", 
           sub="Arms are 95 percent CIs"))

5-13 Differences in extra hours of sleep

bac_d <- with(dta, bac[test==2] - bac[test==1])
bac_d |> knitr::kable()
x
-0.015
-0.014
0.001
0.008
0.010
0.012
0.021
0.004
0.006
-0.008
-0.019
-0.018
-0.022
-0.014
-0.006

5-14 Prolongation of BAC

stripchart(bac_d, 
           frame=F, pch=1, 
           method="stack", 
           xlim=c(-0.05, 0.05),
           xlab="Differences in BAC(test 2-1)",
           main="BAC")
abline(v=0, lty=2)
boxplot(bac_d, 
        horizontal=TRUE, 
        frame=F,
        add=TRUE,
        at=.6, 
        pars=list(boxwex=0.5, 
                  staplewex=0.25))

knitr::kable(t(sort(bac_d)))
-0.022 -0.019 -0.018 -0.015 -0.014 -0.014 -0.008 -0.006 0.001 0.004 0.006 0.008 0.01 0.012 0.021

5-15 Comparing the mean differences to null

t.test(bac_d)
## 
##  One Sample t-test
## 
## data:  bac_d
## t = -1.0447, df = 14, p-value = 0.3139
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.010990739  0.003790739
## sample estimates:
## mean of x 
##   -0.0036

5-16 Equivalence test (Schuirmann, 1987)

Two One-Sided Tests (TOST) Equivalence Testing

H0: The mean of BAC between groups is different beyond the equivalence bound of (0, 1)

H1: The mean of BAC is the same within (0, 1)

pacman::p_load(TOSTER)
TOSTone(m=mean(bac_d),
        mu=0,
        sd=sd(bac_d),
        n=length(bac_d),
        low_eqbound_d = 0, 
        high_eqbound_d = 0.05, 
        alpha=.05,
        plot=TRUE, 
        verbose=FALSE)

5-17 Equivalence bound revised

H0: The mean of BAC between groups is different beyond the equivalence bound of (0, 2)

H1: The mean of BAC is the same within (0, 1)

TOSTone(m=mean(bac_d),
        mu=0,
        sd=sd(bac_d),
        n=length(bac_d),
        low_eqbound_d=0, 
        high_eqbound_d=1, 
        alpha=.05,
        plot=TRUE, 
        verbose=FALSE)

Conclusions

No difference before and after BAC measurement (P=0.3139).