BAC_RueiNing
Ruei-Ning
2021/10/17
BAC replication
The dataset, BAC{HoRM}, is from a study to compare the blood alcohol concentration (BAC) of subjects using two different methods (Breathalyzer Model 5000 or labtest BAC in a laboratory) on 15 subjects. Replicate the analysis by following the sleep example of the class note.
library(HoRM)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoodata(BAC, package="HoRM")
dta <- reshape(data.frame(BAC, id=1:15), idvar='id', varying=list(1:2),direction = 'long', timevar="test", v.names="bac")BAC data
Take a look of the BAC data.
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.000The effect of using two methods (test 1 and test 2) in BAC level on 15 participants. I combined the data by names of rows and show 15 participants in every page.
DT::datatable(dta, rownames=FALSE, fillContainer=FALSE, options=list(pageLength=15))Summary statistics
Mean of BAC in two tests
Calculate the mean and the sd of BAC in two tests. Show the outcomes by table.
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 |
Plot the BAC level by factors
Plot the means of the BAC level by id and by test. Give the dotted line dot on the graph by using the code “grid()”.
library(dplyr)##
## 載入套件：'dplyr'## 下列物件被遮斷自 'package:stats':
##
## filter, lag## 下列物件被遮斷自 'package:base':
##
## intersect, setdiff, setequal, uniondta <- dta |> mutate(id = as.factor(id),
test = as.factor(test))
dta |> plot.design(bty='n',
ylab="Mean of bac level")
grid() #given the dotted line dot 
# Plot the BAC level of two tests by id Plot the BAC level of two tests by id. Give the dotted line dot of the difference between test 1 and teat 2 by each id. We show the differences in red.
with(dta,
interaction.plot(id,
test,
bac,
bty='n',
ylab="BAC level"))
grid()
i <- 1:15
segments(as.numeric(dta$id[i]),
dta$bac[i],
as.numeric(dta$id[i]),
dta$bac[i+15],
col=2,
lty=3,
lwd=2.4)
# 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 |
x <- aggregate(bac ~ id, data=dta, min)
x[order(x$bac),] |> t() |> knitr::kable()| 15 | 10 | 7 | 14 | 11 | 4 | 6 | 8 | 13 | 12 | 1 | 2 | 5 | 9 | 3 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| id | 15 | 10 | 7 | 14 | 11 | 4 | 6 | 8 | 13 | 12 | 1 | 2 | 5 | 9 | 3 |
| bac | 0.000 | 0.048 | 0.060 | 0.065 | 0.092 | 0.100 | 0.100 | 0.100 | 0.121 | 0.144 | 0.145 | 0.156 | 0.170 | 0.170 | 0.180 |
dta$id <- factor(dta$id, levels(dta$id)[x[order(x$bac),1]])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="Level of BAC",
par.settings=list(superpose.symbol=list(pch=1,
col=c("#28E2E5", "#CD0BBC"))),
auto.key=list(space='top', columns=2))
## 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 type",
y="Level of BAC")+
scale_y_continuous(limits=c(0, 0.2),
breaks=seq(0, 0.2, by=0.01))+
geom_point()+
theme_classic()
Data vectors
Do pair-t test.
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.0036another code of pair-t test
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.0036Formula 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.0036aov(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 ' ' 1Linear model
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 ' ' 1Difference Plot
pacman::p_load(multicon)
with(dta,
diffPlot(bac[test==1],
bac[test==2],
paired=TRUE,
grp.names=c("Test 1", "Test 2"),
xlab="",
ylab="Level of BAC",
main="Difference of paired means",
sub="Arms are 95 percent CIs"))
## Differences in BAC
ehs_d <- with(dta, bac[test==2] - bac[test==1])
ehs_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 |
Difference of BAC
stripchart(ehs_d,
frame=F, pch=1,
method="stack",
xlim=c(-0.03, 0.03),
xlab="Differences in two tests of BAC level",
main="Differences of BAC")
abline(v=0, lty=2)
boxplot(ehs_d,
horizontal=TRUE,
frame=F,
add=TRUE,
at=.6,
pars=list(boxwex=0.5,
staplewex=0.25))
knitr::kable(t(sort(ehs_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 |
Comparing the mean differences to null
t.test(ehs_d)##
## One Sample t-test
##
## data: ehs_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.0036Equivalence test (Schuirmann, 1987)
H0: The mean difference of BAC between tests is different beyond the equivalence bound of (0, 1) H1: The mean difference of BAC between tests is the same within (0, 1)
pacman::p_load(TOSTER)
TOSTone(m=mean(ehs_d),
mu=0,
sd=sd(ehs_d),
n=length(ehs_d),
low_eqbound_d=0,
high_eqbound_d=1,
alpha=.05,
plot=TRUE,
verbose=FALSE)
## Equivalence bound revised H0: The mean difference of BAC between tests is different beyond the equivalence bound of (-0.1, 0.1) H1: The mean difference of BAC between tests is the same within (-0.1, 0.1
TOSTone(m=mean(ehs_d),
mu=0,
sd=sd(ehs_d),
n=length(ehs_d),
low_eqbound_d=-0.1,
high_eqbound_d=0.1,
alpha=.05,
plot=TRUE,
verbose=FALSE)
Conclusion
There is no difference in BAC level by using Breathalyzer Model 5000 or labtest BAC in a laboratory (t = 1.0447, df = 14, p-value = 0.3139).