0
$\begingroup$

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.

$\endgroup$
1
+50
$\begingroup$
  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
$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$ – Numbermind Mar 22 at 9:14
  • $\begingroup$ I found about Rmisc package. Can't I just use the CI function to each of the groups to get the confidence interval? $\endgroup$ – Numbermind Mar 22 at 11:39
  • $\begingroup$ 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. $\endgroup$ – R Carnell Mar 22 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.