Rozpoczynam od postawienia hipotez:

H0: Relacje pomiędzy rodzajem przewoźnika a miastem są nieistotne. H1: Relacje pomiędzy rodzajem przewoźnika a miastem są istotne.

Ładuję dane:

dane = data.frame(
  Usluga = c(rep("DHL", 8), rep("Paczkomat", 16),
              rep("PocztaPolska", 12), rep("DPD", 12)),
  Miasto = c(rep(c("Gdansk", "Warszawa", "Szczecin",
                        "Wroclaw", "Krakow", "Poznan"), 8)),
  Czas = c(15.23, 14.32, 14.77, 15.12, 14.05, 17.12, 18.15, 20.05,
           15.48, 14.13, 14.46, 15.62, 14.23, 15.19, 14.67, 14.48, 15.34, 14.22,
           16.66, 16.27, 16.35, 16.93, 15.05, 16.98, 16.43, 15.95, 16.73, 15.62,
           16.53, 16.26, 15.69, 16.97, 15.37, 17.12, 16.65, 15.73, 17.77, 15.52,
           16.15, 16.86, 15.18, 17.96, 15.26, 16.36, 16.44, 14.82, 17.62, 15.04)
)

Testuję normalność rozkładu:

shapiro.test(c(15.23, 14.32, 14.77, 15.12, 14.05, 17.12, 18.15, 20.05,
           15.48, 14.13, 14.46, 15.62, 14.23, 15.19, 14.67, 14.48, 15.34, 14.22,
           16.66, 16.27, 16.35, 16.93, 15.05, 16.98, 16.43, 15.95, 16.73, 15.62,
           16.53, 16.26, 15.69, 16.97, 15.37, 17.12, 16.65, 15.73, 17.77, 15.52,
           16.15, 16.86, 15.18, 17.96, 15.26, 16.36, 16.44, 14.82, 17.62, 15.04))
## 
##  Shapiro-Wilk normality test
## 
## data:  c(15.23, 14.32, 14.77, 15.12, 14.05, 17.12, 18.15, 20.05, 15.48, 14.13, 14.46, 15.62, 14.23, 15.19, 14.67, 14.48, 15.34, 14.22, 16.66, 16.27, 16.35, 16.93, 15.05, 16.98, 16.43, 15.95, 16.73, 15.62, 16.53, 16.26, 15.69, 16.97, 15.37, 17.12, 16.65, 15.73, 17.77, 15.52, 16.15, 16.86, 15.18, 17.96, 15.26, 16.36, 16.44, 14.82, 17.62, 15.04)
## W = 0.95322, p-value = 0.0538

Wartość p-value>alfa (0.0538>0.05), zatem stwierdzam, że jest rozkład normalny.

bartlett.test(Czas~Usluga, dane)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Czas by Usluga
## Bartlett's K-squared = 15.134, df = 3, p-value = 0.001706

Wynik testu pokazuje, że p-value<alfa, czyli występuje niejednorodność wariancji.

Sprawdzam, czy są obserwacje odstające:

dane %>%
  group_by(Usluga,Miasto) %>%
  identify_outliers(Czas)
## [1] Usluga     Miasto     Czas       is.outlier is.extreme
## <0 wierszy> (lub 'row.names' o zerowej długości)

# Nie ma żadnych obserwacji odstających.

Mamy 3 próby i dwie zmienne grupujące, więc wykonujemy test anova dwuczynnikowy:

dane %>%
  anova_test(Czas~Usluga*Miasto)
## Coefficient covariances computed by hccm()
## ANOVA Table (type II tests)
## 
##          Effect DFn DFd     F     p p<.05   ges
## 1        Usluga   3  24 1.079 0.377       0.119
## 2        Miasto   5  24 0.291 0.914       0.057
## 3 Usluga:Miasto  15  24 0.336 0.984       0.174

#P-value > 0.05, więc nie mogę odrzucić hipotezy zerowej, czyli stwierdzam, że interakcje między rodzajem przewoźnika a miastem nie są istotne z prawdopodobieństwem 95%