knitr::opts_chunk$set(echo = TRUE)
require(ggplot2)
require(tigerstats)
Let’s read in the data:
lines <-
"TREAT DAYS
1 8
1 7
1 7
1 10
1 8
2 4
2 6
2 6
2 3
2 5
2 6
3 6
3 4
3 4
3 5
3 4
3 3
4 7
4 4
4 6
4 6
4 3
4 5
5 9
5 3
5 5
5 7
5 7
5 6"
con <- textConnection(lines)
data <- read.csv(con, sep="")
data$TREAT <- relevel(factor(data$TREAT), ref=5)
The data consists of observation on patients who were given five different treatments. For each patient, we have the treatment they were given, and the day it took for the patient to heal.
Here is a boxplot:
bwplot(TREAT~DAYS, data=data)
Let’s fit a linear model:
fit <- lm(DAYS~TREAT, data=data)
summary(fit)
##
## Call:
## lm(formula = DAYS ~ TREAT, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1667 -1.0000 0.0000 0.8333 2.8333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1667 0.5951 10.362 2.44e-10 ***
## TREAT1 1.8333 0.8827 2.077 0.0487 *
## TREAT2 -1.1667 0.8416 -1.386 0.1784
## TREAT3 -1.8333 0.8416 -2.178 0.0394 *
## TREAT4 -1.0000 0.8416 -1.188 0.2464
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.458 on 24 degrees of freedom
## Multiple R-squared: 0.4606, Adjusted R-squared: 0.3707
## F-statistic: 5.124 on 4 and 24 DF, p-value: 0.003959
Here is the output from running ANOVA
anova(fit)
## Analysis of Variance Table
##
## Response: DAYS
## Df Sum Sq Mean Sq F value Pr(>F)
## TREAT 4 43.552 10.888 5.1237 0.003959 **
## Residuals 24 51.000 2.125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Some diagnostic plots:
plot(fit)
Here are the pairwise comparisons:
with(data, pairwise.t.test(DAYS, TREAT, p.adj = "none"))
##
## Pairwise comparisons using t tests with pooled SD
##
## data: DAYS and TREAT
##
## 5 1 2 3
## 1 0.04868 - - -
## 2 0.17843 0.00237 - -
## 3 0.03944 0.00036 0.43605 -
## 4 0.24639 0.00375 0.84469 0.33198
##
## P value adjustment method: none
with(data, pairwise.t.test(DAYS, TREAT, p.adj = "bonf"))
##
## Pairwise comparisons using t tests with pooled SD
##
## data: DAYS and TREAT
##
## 5 1 2 3
## 1 0.4868 - - -
## 2 1.0000 0.0237 - -
## 3 0.3944 0.0036 1.0000 -
## 4 1.0000 0.0375 1.0000 1.0000
##
## P value adjustment method: bonferroni
TukeyHSD(aov(fit), conf.level = 0.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = fit)
##
## $TREAT
## diff lwr upr p adj
## 1-5 1.8333333 -0.7671386 4.4338052 0.2621906
## 2-5 -1.1666667 -3.6461193 1.3127860 0.6419253
## 3-5 -1.8333333 -4.3127860 0.6461193 0.2217409
## 4-5 -1.0000000 -3.4794527 1.4794527 0.7578974
## 2-1 -3.0000000 -5.6004719 -0.3995281 0.0182053
## 3-1 -3.6666667 -6.2671386 -1.0661948 0.0029935
## 4-1 -2.8333333 -5.4338052 -0.2328614 0.0279925
## 3-2 -0.6666667 -3.1461193 1.8127860 0.9304016
## 4-2 0.1666667 -2.3127860 2.6461193 0.9996305
## 4-3 0.8333333 -1.6461193 3.3127860 0.8571322