Cochran’s and SNK Tests

Let’s start by inputting all of the necessary data.

Remember to set both factors as fixed and that all three are the same length.

solution<-c(rep(1,4),rep(2,4),rep(3,4))
days<-c(seq(1,4),seq(1,4),seq(1,4))
obs<-c(13,22,18,39,
       16,24,17,44,
       5,4,1,22)
solution<-as.fixed(solution)
days<-as.fixed(days)

Now that we have data we can construct our model using each only using each individual factors.

model<-lm(obs~solution)
C.test(model)
## 
##  Cochran test of homogeneity of variances
## 
## data:  model
## C = 0.43732, n = 4, k = 3, p-value = 0.8851
## alternative hypothesis: Group 2 has outlying variance
## sample estimates:
##        1        2        3 
## 127.3333 168.9167  90.0000

From this large p-value (0.8851) we can not reject the Ho that the solution has outlying variance. This is partially because there are so few samples that it becomes hard to reject Ho due to a lack of power in this current test.

model<-lm(obs~days)
C.test(model)
## 
##  Cochran test of homogeneity of variances
## 
## data:  model
## C = 0.58277, n = 3, k = 4, p-value = 0.2905
## alternative hypothesis: Group 3 has outlying variance
## sample estimates:
##        1        2        3        4 
##  20.3333 136.3333 399.0000 129.0000

Performing this same test with the days factor, we still fail to reject although with a lower p-value of 0.2905. Do note that this test suffers from the same lack of power as the first. This is despite the fact that the variance of 3 is so much larger than the rest. It is recommended to perform this test with no less than 15 observations (n).

SNK Test

The snk test utilizes a studentized range test similar to Tukey’s Test but controls experiment-wise error rate which differs from Tukey’s.

To start with SNK we need a model with both factors and no interactions since the days is a block.

Personally, adding the model tab line was necessary to provide the function with the gad data it needed although this differs from other examples.

model<-lm(obs~solution+days)
model.tab <- gad(model)
snk.test(model,term='solution', anova.tab = model.tab)
## Student-Newman-Keuls test for: solution 
## 
## Standard error = 1.4696
## Df = 6 
##               3       1      2    
## Rank order:   1       2      3    
## Ranked means: 8       23     25.25
## Comparisons:                      
## 1             3-1 ***             
## 2             2-1 *** 3-2 ns      
## Summary: 3 < 1 = 2 
## ---
## Signif. codes: ***p < 0.001; **p < 0.01; *p < 0.05; ns, p > 0.05

The results of this test shows that there is a significant difference between rank 3 and 1 indicated by ***. Using this logic we can also see that there is a significant difference between ranks 2 and 1. Exercise caution when reading these outputs as it shows the interaction between the ranks from the Rank Order: line not the header of the table.

Full Code

Here is the full code if you wish to experiment with the code shown in this guide.

library(GAD)
library(tidyr)
library(agricolae)
solution<-c(rep(1,4),rep(2,4),rep(3,4))
days<-c(seq(1,4),seq(1,4),seq(1,4))
obs<-c(13,22,18,39,
       16,24,17,44,
       5,4,1,22)
solution<-as.fixed(solution)
days<-as.fixed(days)
model<-lm(obs~solution+days)
gad(model)

model<-lm(obs~solution)
C.test(model)
model<-lm(obs~days)
C.test(model)


model<-lm(obs~solution+days)
model.tab <- gad(model)
snk.test(model,term='solution', anova.tab = model.tab)