data management example: BAC
Jenny Hsu
2021/10/11
data information
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.
#install.packages("HoRM")
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")
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.000reshape是將資料變型的函數，direction設定轉為long form varying是指欄位1:2要進行transform的欄位
The effect of two methods on 15 patients.
DT::datatable(dta, rownames=FALSE, fillContainer=FALSE, options=list(pageLength=10))Summary statistics
caculate mean and sd of bac by test
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 |
Design plot
library(dplyr)dta <- dta |> mutate(id = as.factor(id),
test = as.factor(test))
dta |> plot.design(bty='n',
ylim = c(0.005,0.2),
ylab=" BAC value")
grid()
Subject by test effect
with(dta,
interaction.plot(id,
test,
bac,
bty='n',
ylab=" BAC value"))
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)
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]])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="ID",
xlab="BAC values",
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",
y="BAC values")+
scale_y_continuous(limits=c(-0.01, 0.2),
breaks=seq(-0.01, 0.2, by=0.01))+
geom_point()+
theme_classic()
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.0036with(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.0036ANOVA
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 ' ' 1Linear 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 ' ' 1Difference plot
pacman::p_load(multicon)with(dta,
diffPlot(bac[test==1],
bac[test==2],
paired=TRUE,
grp.names=c("1", "2"),
xlab="",
ylab="BAC values",
main="Difference of paired means",
sub="Arms are 95 percent CIs"))
Differences in bac
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 |
boxplot of BAC
# one dimensional scatter charts
stripchart(bac_d,
frame=F, pch=1,
method="stack",
xlim=c(-0.05, 0.05),
xlab="Differences in BAC value",
main="BAC")
abline(v=0, lty=2)
boxplot(bac_d,
horizontal=TRUE,
frame=F,
add=T,
at=.6,
pars=list(boxwex=0.5,
staplewex=0.25))
stripchart, When dealing with small sample sizes (few data points)
abline = Add a vertical line lty = line types
add = T, means combine stripchart and boxplot
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 |
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.0036Equivalence test (Schuirmann, 1987)
H0: The mean of BAC between tests is different beyond the equivalence bound of (0, 0.5)
H1: The mean of BAC is the same within (0, 0.5)
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.5,
alpha=.05,
plot=TRUE,
verbose=FALSE)
Equivalence bound revised
H0: The mean of BAC between tests is different beyond the equivalence bound of (-0.05, 0.05)
H1: The mean of BAC is the same within (-0.05, 0.05)
TOSTone(m=mean(bac_d),
mu=0,
sd=sd(bac_d),
n=length(bac_d),
low_eqbound_d=-0.05,
high_eqbound_d=0.05,
alpha=.05,
plot=TRUE,
verbose=FALSE)
Conclusions
blood alcohol concentration (BAC) of subjects in two different methods (Breathalyzer Model 5000 or labtest BAC in a laboratory) did not significance different on 15 subjects.