20211004 Homework
Chang-Chun Chen
2021-10-14
xfun::pkg_load2(c("htmltools", "mime"))1 Homework 1
2 Homework 2
3 Homework 3
3.1 Homework 3-1
5 Homework 5
5.1 Data Description
{HoRM} BAC dataset
This dataset is from a study to compare the blood alcohol concentration (BAC) of subjects using two different methods.
- Format
This data frame consists of 2 variables measured on 15 subjects:
Column 1: breath BAC obtained using the Breathalyzer Model 5000.
Column 2: labtest BAC based on a breath estimate in a laboratory.
# install.packages("HoRM")
library(HoRM)Input data and create id & test
test= 1 is the breath BAC value
test= 2 is the labtest BAC value
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")The BAC value of two test from 15 subjects.
DT::datatable(dta[,], rownames=FALSE, fillContainer=FALSE, options=list(pageLength=15))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 |
The mean and SD value of BAC from Breathalyzer and breath estimate in a laboratory are similar.
5.3 Design plot
library(dplyr)dta <- dta |> mutate(id = as.factor(id),
test = as.factor(test))
dta |> plot.design (bty='n',
ylab="Mean BAC value")
grid()
5.4 Subject by test 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 value",
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 value")+
scale_y_continuous(limits=c(-0.01, 0.2),
breaks=seq(-0.01, 0.2, by=0.01))+
geom_point()+
theme_classic()
5.8 Slope plot with ID
ggplot(dta) +
geom_boxplot(aes(x = as.numeric(test), y = bac, group = test), width=0.1)+
geom_point(aes(x = as.numeric(test), y = bac)) +
geom_line(aes(x = as.numeric(test), y = bac, group = id)) +
geom_text(aes(x = as.numeric(test), y = bac, label = id), nudge_x = 0.1)+
scale_x_continuous(breaks = c(1,2), labels = c("T1", "T2"))+
scale_y_continuous(limits = c(-0.01,0.2))+
xlab("Test")+
theme_minimal()
5.9 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.00365.10 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.00365.11 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 ' ' 15.12 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 ' ' 15.13 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 value",
main="Difference of paired means",
sub="Arms are 95 percent CIs"))
5.14 Differences in BAC
(Lab-Breathalyzer)
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.15 BAC
stripchart(bac_d,
frame=F, pch=1,
method="stack",
xlim=c(-0.05, 0.05),
xlab="Differences in BAC",
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 |
誤差範圍 +- 0.02 ?
5.16 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.0036No difference between lab and the equipment
5.17 Equivalence test
Two One-Sided Tests (TOST) Equivalence Testing.
Null hypothesis significance testing (NHST)
H0: The mean BAC between groups is different beyond the equivalence bound of (0, 0.05)
H1: The mean BAC is the same within (0, 0.05)
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.18 Equivalence bound revised
H0: The mean BAC between groups is different beyond the equivalence bound of (-1, 0.5)
H1: The mean BAC is the same within (-1, 0.5)
TOSTone(m=mean(bac_d),
mu=0,
sd=sd(bac_d),
n=length(bac_d),
low_eqbound_d=-1,
high_eqbound_d=0.5,
alpha=.05,
plot=TRUE,
verbose=FALSE)
5.19 Conclusion
-The BAC value measured from equipment is same as the estimated from laboratory. (Equipment can be an alternative to measure BAC?)
5.20 Reference
- Source
Krishnamoorthy, K., Kulkarni, P. M., and Mathew, T. (2001), Multiple Use One-Sided Hypotheses Testing in Univariate Linear Calibration, Journal of Statistical Planning and Inference, 93, 211– 223.