P-value for Sheffe-test from pooled summary data (means and standard deviation)

Question

I'm stuck, how do I calculate p-values for the Scheffé test from pooled data (means and standard deviations) like "sheffeTest" from library(DescTools)? I didn't find any documentation. The contrasts I'm considering are all possible pairwise comparisons.

I got the following data called Data_sum:

Location  n       mean       sd
<fctr>    <dbl>   <dbl>      <dbl>
TM    5   0.08016    0.0145622
NP    5   0.07550    0.0091643
PB    5   0.10434    0.0195800

and my contrast table looks like this (which is the same as using emmeans(model, pairwise ~ Location, adjust="scheffe"):

 Contrast.Name     TM   NP    PB
 TMvsNP           -1    1     0
 TMvsPB           -1    0     1
 NPvsPB            0   -1     1

I calculated an Anova using anovaMean() from library(HH). I can calculate Sheffe by hand given the values from Anova, but im not able to calulate p-values. Any clues?

Edit: solved, see code below

1. Calculate Anova from summary statistics (means, standard deviations, sample sizes)
2. Build contrasts
3. Do the math as advised by @PBulls

note: since pooled_n = mean(n_precalcs), results from inbalanced design may be inaccurate

library(HH)
library(magrittr)
library(rstatix)

scheffe.test <- function(iv,n,ybar,sd, alpha = 0.05){
  if (sd(n) != 0) warning("Inbalanced Design, results may be inaccurate!\n")
  iv = factor(iv,levels = unique(iv)) # iv to levels
  anova_sum = anovaMean(iv,n,ybar,sd) # anova from mean and sd
  
  # calculate contrasts
  calculate_contrasts <- function(lvl1,lvl2,iv,obs,ybar){
    m <- function(lvl){ybar[iv == lvl] %>% mean()}
    n <- function(lvl){obs[iv == lvl] %>% sum()}
    
    return(list(contrast   = paste(lvl1, "-", lvl2), # generate contrasts
                mean_diffs = m(lvl1) - m(lvl2),      # mean diff between contrasts
                n_precalc  = 1/n(lvl1) + 1/n(lvl2)))}# part of t-statistic
    
  # build contrasts
  contr      = character()
  mean_diffs = numeric()
  n_precalcs = numeric()
    
  for (i in 1:(nlevels(iv)-1)) {
    for (j in (i+1):nlevels(iv)) {
      lvl1 = levels(iv)[i]
      lvl2 = levels(iv)[j]
      res  = calculate_contrasts(lvl1,lvl2,iv,n,ybar)
      contr      = append(contr, res$contrast)
      mean_diffs = append(mean_diffs, res$mean_diffs)
      n_precalcs = append(n_precalcs, res$n_precalc)}}
  
  print(contr)

  #scheffe statistics
  n    = sum(n)      # total number of observations (sum of sample sizes)
  k    = nlevels(iv) # total number of groups/ levels of iv variable
  rank = k-1         # rank of the scheffe-test
  mse      = anova_sum$`Mean Sq`[2] # mean squared error from anovaMean
  pooled_n = mean(n_precalcs)       # if balanced all pooled_n are the same
  se       = sqrt(mse*pooled_n)     # standard error (of the means)
  t.ratio  = mean_diffs/se          # t.ratio
  p.values = pf(t.ratio**2/rank, rank, n-k, lower.tail = F)
  t.crit   = sqrt(rank*qf(1-alpha, rank, n-k)) # critical t value
  lwr.ci = mean_diffs - t.crit * se #upper conf int boundary
  upr.ci = mean_diffs + t.crit * se #lower conf int boundary
  
  return(data.frame(contrast  = contr,
                    estimates = mean_diffs,
                    SE        = se,
                    df        = n-k,
                    t.ratio   = t.ratio,
                    p.value   = p.values,
                    lwr.ci    = lwr.ci,
                    upr.ci    = upr.ci))}

The input for the scheffe.test() function is the same as for anovaMean() from library(HH). The standard alpha level is set to 0.05 and can be changed by passing the alpha level additionally in the scheffe.test() function e.g.:

scheffe.test(iv    = Data_sum$Location, # level names
             n     = Data_sum$n,        # sample size for each level
             ybar  = Data_sum$mean,     # sample mean for each level
             sd    = Data_sum$sd,       # sample standard deviation for each level
             alpha = 0.05) %>%          # change 0.05 alpha level if necessary
  add_significance("p.value") #add significance stars

The result is a dataframe in style of "sheffeTest" from library(DescTools) with confidence intervals, but using summary statistics instead. (note: To generate exact same contrasts as DescTools::scheffeTest(), run factor(levels = unique()) on indipendent variable from raw data.)

Answer 1

The Scheffé statistic follows an F distribution with degrees of freedom of the 'rank' (the number of groups under comparison minus one) and the residuals, and the statistic itself is the squared t-statistic divided by the rank.

You can calculate all of this by hand from the info you have provided plus the ANOVA MSE, which you can also calculate from those data:

n <- 15       ## Number of observations
k <- 3        ## Number of groups
rank <- k-1   ## Rank of the test

## Statistics derived from your data
MSE <- 0.00022647282
diffs <- c(-0.08016+0.07550, -0.08016+0.10434, -0.07550+0.10434)

## Calculate the t-statistics - assuming balanced samples (OK here)
tstat <- diffs / sqrt(MSE * (1/5 + 1/5))

## P-values
pf(tstat**2 / rank, rank, 15-k, lower.tail = FALSE)

#> 0.888 0.0756 0.033

Just to prove that this will match e.g. DescTools::ScheffeTest or emmeans with adjust="scheffe":

f <- summary(a <- aov(breaks ~ tension, data=warpbreaks))  ## MSE = 141.1481
DescTools::ScheffeTest(a)

#>           diff  pval    
#> M-L -10.000000  0.0496 
#> H-L -14.722222  0.0022
#> H-M  -4.722222  0.4960

wMSE <- 141.1481
## N=18 per group (balanced)
tstat <- c(-10, -14.722, -4.722) / sqrt(wMSE * (1/18 + 1/18))
pf(tstat**2 / 2, 2, 54-3, lower.tail = FALSE)

#> 0.04959 0.00221 0.49601

Edit: calculations for the confidence interval were also requested. For this you have to recall that the Scheffé critical value for an alpha-level test on k groups is the F distribution quantile sqrt((k-1)*qf(1-alpha, k-1, n-k)) (this would look much nicer in LaTeX). Combining that with our standard error from before:

alpha <- .05
crit <- sqrt(rank*qf(1-alpha, rank, n-k))   ## 2.787577

vapply(diffs, \(d)
  d + c(-1,1) * crit * sqrt(MSE * (1/5 + 1/5)), numeric(2))
#>          [,1]    [,2]   [,3]
#> lower -0.0312 -0.0024 0.0023
#> upper  0.0219  0.0507 0.0554

Doing the same for the warpbreaks example to confirm:

wcrit <- sqrt(2*qf(1-alpha, 2, 54-3))
vapply(c(-10, -14.722, -4.722),  \(d)
  d + c(-1,1) * wcrit * sqrt(wMSE * (1/18 + 1/18)), numeric(2))
#>          [,1]   [,2]   [,3]
#> lower -19.985 -24.71 -14.71
#> upper  -0.015  -4.73   5.26

## Compare again vs. DescTools
DescTools::ScheffeTest(a)
#>      lwr.ci  upr.ci
#> M-L -19.985 -0.0146
#> H-L -24.707 -4.7368
#> H-M -14.707  5.2631