# One-way ANOVA followed by Tukey test - confidence intervals

I want to test 3 types of treatments (T1, T2, T3) along with a control/placebo treatment (T4). In order to do so, I'm performing One-Way ANOVA to my data and then proceeding to the Tukey test.

     # generating artificial data for this example
y <- c(rnorm(10,21,1),rnorm(10,17,2),rnorm(10,15,2),rnorm(10,18,3))
Type <- rep(c("T1","T2","T3","T4"),each=10)  # Variável explicativa (categórica)
data <- data.frame(y,Type)
data

> data
y    Type
1  22.26349   T1
2  18.95373   T1
3  20.53918   T1
4  19.12539   T1
5  20.21258   T1
6  19.66778   T1
7  20.51515   T1
8  19.66407   T1
9  21.76689   T1
10 21.62950   T1
11 18.59258   T2
12 21.16842   T2
13 20.14588   T2
14 18.36684   T2
15 16.49939   T2
16 14.71304   T2
17 16.53829   T2
18 15.56384   T2
19 17.95939   T2
20 16.62988   T2
21 16.14444   T3
22 14.31957   T3
23 16.42375   T3
24 14.37554   T3
25 14.07891   T3
26 14.48411   T3
27 13.86916   T3
28 11.20102   T3
29 15.65992   T3
30 14.85290   T3
31 21.96175   T4
32 24.26759   T4
33 16.52458   T4
34 17.19564   T4
35 17.63090   T4
36 17.81682   T4
37 19.63599   T4
38 21.06390   T4
39 12.96356   T4
40 20.75986   T4


From the summary of my ANOVA model I can see that the p-value is 3.059e-06 < 0.05, which indicates that we can reject the null hypothesis that all means across all groups are equal, i.e, there is a significant difference among group means:

         anova.model <- lm(y ~ Type, data= data)
> anova(anova.model)
Analysis of Variance Table

Response: y
Df Sum Sq Mean Sq F value    Pr(>F)
Type       3 189.54  63.179  14.117 3.059e-06 ***
Residuals 36 161.12   4.475
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Another way to get this p-value would be:

        > summary(anova.model)

Call:
lm(formula = y ~ Type, data = data)

Residuals:
Min      1Q  Median      3Q     Max
-6.0185 -1.1301 -0.1111  1.2301  5.2855

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  20.4338     0.6690  30.544  < 2e-16 ***
TypeT2       -2.8160     0.9461  -2.976  0.00519 **
TypeT3       -5.8928     0.9461  -6.229 3.44e-07 ***
TypeT4       -1.4517     0.9461  -1.534  0.13367
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.116 on 36 degrees of freedom
Multiple R-squared:  0.5405, Adjusted R-squared:  0.5022
F-statistic: 14.12 on 3 and 36 DF,  p-value: 3.059e-06


Since this only indicates that there is a difference among the groups without telling us which ones are significant, I decided to run a Tukey test:

       > TukeyHSD(aov(y ~ Type, data= data))
Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = y ~ Type, data = data)

$Type diff lwr upr p adj T2-T1 -2.816021 -5.364065 -0.2679777 0.0255089 T3-T1 -5.892842 -8.440886 -3.3447987 0.0000020 T4-T1 -1.451716 -3.999759 1.0963275 0.4281661 T3-T2 -3.076821 -5.624865 -0.5287776 0.0127180 T4-T2 1.364305 -1.183738 3.9123487 0.4822404 T4-T3 4.441126 1.893083 6.9891697 0.0002152  From the p-adj got in the Tukey test, one concludes that there is a significant statistical difference between T1 and T2, T1 and T3, T3 and T2, and between T4 and T3. However, I got some questions regarding this procedure. 1. What do the Pr(>|t|) values mean when using the summary command? 2. Are the confidence intervals obtained in the Tukey test referring to the sample or population? For example, for the confidence interval between T3 and T4, I got $$]1.893083,6.9891697[$$. What does it literally mean? Does this mean that $$95\%$$ of the mean differences between T3 and T4 in the sample are in that range of values, or are they saying that for the population? 3. How can I get a confidence interval for each one of the treatments? For example, a $$95\%$$ CI for T1, i.e, a CI that says that people in treatment 1 will have their values in that interval of values. ## 1 Answer 1. What do the Pr(>|t|) values mean when using the summary command? $$P(\beta_{Tj} \neq 0)$$ for a two sided test of $$H_o: E(T_1) = E(T_j)$$ Using this model: $$y = \beta_0 + \beta_{T2} x_{T2} + \beta_{T3} x_{T3} + \beta_{T4} x_{T4} + \epsilon$$ In words, the Pr(>|t|) is the probability that the coefficient on $$T_j$$ is non-zero when compared to $$T_1$$, or it is the probability that $$T_j$$ is significantly difference from $$T_1$$ since $$T_1$$ is the baseline condition. 1. Are the confidence intervals obtained in the Tukey test referring to the sample or population? These are confidence intervals calculated from the sample, but they are estimates of the true parameters (entire population). In other words, (lwr, upr) is a 95% confidence interval for the true difference $$T_2-T_1$$ 1. How can I get a confidence interval for each one of the treatments? For example, a 95% CI for T1, i.e, a CI that says that people in treatment 1 will have their values in that interval of values. See below: set.seed(1973783) # generating artificial data for this example dat <- data.frame(y = c(rnorm(10, 21, 1), rnorm(10, 17, 2), rnorm(10, 15, 2), rnorm(10, 18, 3)), Type = rep(paste0("T", 1:4), each = 10)) lm1 <- lm(y ~ Type, data = dat) anova(lm1) #> Analysis of Variance Table #> #> Response: y #> Df Sum Sq Mean Sq F value Pr(>F) #> Type 3 196.05 65.350 14.792 1.95e-06 *** #> Residuals 36 159.04 4.418 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 TukeyHSD(aov(y ~ Type, data = dat)) #> Tukey multiple comparisons of means #> 95% family-wise confidence level #> #> Fit: aov(formula = y ~ Type, data = dat) #> #>$Type
#>             diff        lwr        upr     p adj
#> T2-T1 -3.2010636 -5.7326501 -0.6694770 0.0085056
#> T3-T1 -6.2567418 -8.7883283 -3.7251552 0.0000005
#> T4-T1 -2.9592189 -5.4908055 -0.4276324 0.0166078
#> T3-T2 -3.0556782 -5.5872647 -0.5240916 0.0127626
#> T4-T2  0.2418447 -2.2897419  2.7734312 0.9939214
#> T4-T3  3.2975228  0.7659363  5.8291094 0.0064627

# confidence interval for the mean of T1 - 95% chance that the true mean is in this interval
with(subset(dat, Type == "T1"), mean(y) + qt(c(0.025, 0.975), length(y) - 1) * sd(y) / sqrt(length(y)))
#> [1] 20.38922 21.54059
Rmisc::CI(dat$$y[dat$$Type == "T1"])
#>   upper     mean    lower
#> 21.54059 20.96491 20.38922
by(dat$$y, dat$$Type, Rmisc::CI)
#> dat$$Type: T1 #> upper mean lower #> 21.54059 20.96491 20.38922 #> -------------------------------------------------------------------------------------------------- #> dat$$Type: T2
#>    upper     mean    lower
#> 18.94730 17.76384 16.58038
#> --------------------------------------------------------------------------------------------------
#> dat$$Type: T3 #> upper mean lower #> 16.19860 14.70817 13.21773 #> -------------------------------------------------------------------------------------------------- #> dat$$Type: T4
#>    upper     mean    lower
#> 20.26172 18.00569 15.74966

# interval for the population - 95% of the observations fall in this interval
with(subset(dat, Type == "T1"), quantile(y, probs = c(0.025, 0.975)))
#>     2.5%    97.5%
#> 19.50665 21.94073

# prediction interval - 95% chance that the next observation is in this interval
with(subset(dat, Type == "T1"), mean(y) + qt(c(0.025, 0.975), length(y) - 1) * sd(y) * sqrt(1 + 1 / length(y)))
#> [1] 19.05557 22.87424

• Thank you, but can you explain better (by words) my first question? The answer was kinda confusing. Also, you presented 3 different IC's. What is the standard/usual IC? I'll give you the bounty – Numbermind Mar 22 at 9:14
• I found about Rmisc package. Can't I just use the CI function to each of the groups to get the confidence interval? – Numbermind Mar 22 at 11:39
• Yes, you can use the Rmisc package. I added it to my answer. It produces a confidence interval for the mean. I also added more words to my answer to your first question. The "usual" confidence interval is the confidence interval for the mean because it is how you normallly compare groups. I also added the other confidence intervals because you asked about "people in treatment 1 will have values in that interval of values". That is a confidence interval for the sample or the population. That is why I added the other two types. – R Carnell Mar 22 at 15:32