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