1
$\begingroup$

I have a dataset describing age of wine in 3 classes of quality.

> ages
    class age
1  class1 <10
2  class1 <10
3  class1 <10
4  class1 <10
5  class1 <10
6  class1 <10
7  class1 <10
8  class1 <10
9  class1 mid
10 class1 mid
11 class1 mid
12 class1 mid
13 class1 mid
14 class1 mid
15 class1 >60
16 class1 >60
17 class1 >60
18 class1 >60
19 class1 >60
20 class1 >60
21 class1 >60
22 class1 >60
23 class2 <10
24 class2 <10
25 class2 <10
26 class2 <10
27 class2 mid
28 class2 mid
29 class2 mid
30 class2 mid
31 class2 mid
32 class2 mid
33 class2 mid
34 class2 mid
35 class2 mid
36 class2 mid
37 class2 >60
38 class2 >60
39 class2 >60
40 class2 >60
41 class2 >60
42 class3 <10
43 class3 <10
44 class3 <10
45 class3 mid
46 class3 mid
47 class3 >60
48 class3 >60
49 class3 >60
50 class3 >60

> ages_n
   class age  n
1 class1 <10  8
2 class1 mid  6
3 class1 >60  8
4 class2 <10  4
5 class2 mid 10
6 class2 >60  5
7 class3 <10  3
8 class3 mid  2
9 class3 >60  4

or

> ages_t
     class
age   class1 class2 class3
  <10      8      4      3
  mid      6     10      2
  >60      8      5      4

I want to check, if the age structure differs between the 3 classes.

So I run the chi2 test:

> chisq.test(ages_t, simulate.p.value =FALSE)

    Pearson's Chi-squared test

data:  ages_t
X-squared = 3.8921, df = 4, p-value = 0.4208

Warning message:
In chisq.test(ages_t, simulate.p.value = FALSE) :
  Chi-squared approximation may be incorrect

> chisq.test(ages_t, simulate.p.value = TRUE)

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  ages_t
X-squared = 3.8921, df = NA, p-value = 0.4418

No differences in contingency was found. But what exactly does it mean?

I plotted the data:

enter image description here

and

> corrplot(chisq.test((ages_t))$residuals, is.cor = FALSE)

enter image description here

It seems there is some difference between classes of wine.

So I run a log-linear model:

> loglm(n ~ age + class , data=ages_n)
Call:
loglm(formula = n ~ age + class, data = ages_n)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 3.867862  4 0.424184
Pearson          3.892103  4 0.420805

It agreed with the chi2 (Pearson's) test, as expected.

Then I fit the Poisson model with log link, to try a different approach. Out of curiosity I run the analysis of deviance over the log linear model. I know I shouldn't test hypotheses ad hoc, but this is purely exploratory and I'm just curious about the numbers and patterns, just for learning.

First, I set the contrast appropriately:

options(contrasts = c("contr.sum", "contr.poly"))

> car::Anova(glm(n ~ class * age , family=poisson(link="log"), 
                 data=age_n), type=3)
Analysis of Deviance Table (Type III tests)

Response: n
          LR Chisq Df Pr(>Chisq)  
class       5.9499  2    0.05105 .
age         0.2065  2    0.90188  
class:age   3.8679  4    0.42418  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Again, the interaction is non-significant, the p-value matches the LRT test from the log-linear model. So it looks like the Chi2 tests checked if "no ages differ at each levels of the class". But now I can see the class has p-value close to the significance level.

This is a small dataset, so I know such analysis lacks of power. Since I know there seems to be no interaction, I can run a more powerful type 2 ANOVA.

> car::Anova(glm(n ~ class + age , family=poisson(link="log"),  
                 data=age_n), type=2)
Analysis of Deviance Table (Type II tests)

Response: n
      LR Chisq Df Pr(>Chisq)  
class   6.1035  2    0.04728 *
age     0.2831  2    0.86802  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(by the way, the type-1 sequential ANOVA showed the same: class = 0.04728, interaction 0.42418)

Yes, as expected, now the class seems to differentiate the mean log(counts). It doesn't matter whether it's significant or not, here I treat p like a measure of the discrepancy between the data with the null hypothesis and the result is 10x lower the others.

I also plotted the predicted counts:

library(emmeans)
emmip(emmeans(glm(n ~ class * age , family=poisson(link="log"), 
          data=age_n), specs = ~class*age), formula = ~age|class, 
          type="response")

enter image description here

It seems the interaction is "masked" by opposite directions in all groups, but there is a visible discrepancy across classes.

Summary: Both the classic chi2 test and the log-linear analysis tell me, that there is no differences causes by the interaction of age and class.

Both type-3 and type-2 ANOVA tell me instead, that there is difference caused by the class in presence of age.

Now I'm curious, which of the two approaches:

  1. testing for the interaction (chi2, log-linear, type-1/type-3 ANOVA) - neither showed it
  2. testing for the main effect only (all ANOVAs showed it, type-2 had more power)

tests the hypothesis that the class of quality of wine doesn't matter?

Or simpler - does it show the lack of homogeneity in age across the 3 classes?


Please note, it's not about asking "which anova is better", rather which kind of analysis tests the hypothesis, that "is there any heterogeneity in age across classes?" or "are all classes similar in terms of age"? - the interaction (non-existent here) or the main effect (let's pretend it's significant)

The illustration suggested that there is some effect, but I expected rather the interaction to be significant.


EDIT: It seems I found the answer.

The GLM, as it could be expected, predicted conditional mean of the counts. And YES, the ABSOLUTE counts differed across the classes. And the ANOVA showed it. Indeed, class3 has the smallest number of cases.

But my question was about the percentage structure, NOT the counts itself, and the differences in % weren't that big, to make it statistically significant. I realized it looking at the descriptive summary:

> age_n %>% group_by(class, age) %>% summarize(n=n) %>% 
  group_by(class) %>% mutate(mean_n = mean(n), p=n/sum(n)) %>% 
  mutate(p=sprintf("%.1f%%", p*100))
`summarise()` has grouped output by 'class'. You can override 
 using the `.groups` argument.
# A tibble: 9 x 5
# Groups:   class [3]
  class  age       n mean_n p    
  <fct>  <fct> <dbl>  <dbl> <chr>
1 class1 <10       8   7.33 36.4%
2 class1 mid       6   7.33 27.3%
3 class1 >60       8   7.33 36.4%
4 class2 <10       4   6.33 21.1%
5 class2 mid      10   6.33 52.6%
6 class2 >60       5   6.33 26.3%
7 class3 <10       3   3    33.3%
8 class3 mid       2   3    22.2%
9 class3 >60       4   3    44.4%

So it seems, the answer is:

  1. YES, class differentiates the mean count of ages
  2. NO, class does NOT differentiate the % structure of ages.

So, the answer is - the chi2 test + the log-linear analysis (and the interaction term in ANOVA over the Poisson GLM).

$\endgroup$
2
  • 1
    $\begingroup$ Your classes seems to be ordinal, and that could be used in the analysis. Maybe add the tag ordinal-data $\endgroup$ Dec 6, 2022 at 12:51
  • $\begingroup$ I added an answer for treating Age as an ordinal variable. You might post your own answer for treating Age as a nominal variable. $\endgroup$ Dec 6, 2022 at 14:27

1 Answer 1

0
$\begingroup$

Your variable Age is probably best treated as ordinal. Your statement, "I want to check, if the age structure differs between the 3 classes" suggests that you can treat Age as a dependent variable and Class as an independent variable. But I don't know if makes actual sense for your data. Some tests you might look into are the extended Cochran–Armitage test, and Kruskal-Wallis. But for the KW test, you need to Age as the dependent variable, not Count.

R code for some of these tests.

Data = read.table(header=TRUE, stringsAsFactors=TRUE, text="
OBS class age
1  class1 <10
2  class1 <10
3  class1 <10
4  class1 <10
5  class1 <10
6  class1 <10
7  class1 <10
8  class1 <10
9  class1 mid
10 class1 mid
11 class1 mid
12 class1 mid
13 class1 mid
14 class1 mid
15 class1 >60
16 class1 >60
17 class1 >60
18 class1 >60
19 class1 >60
20 class1 >60
21 class1 >60
22 class1 >60
23 class2 <10
24 class2 <10
25 class2 <10
26 class2 <10
27 class2 mid
28 class2 mid
29 class2 mid
30 class2 mid
31 class2 mid
32 class2 mid
33 class2 mid
34 class2 mid
35 class2 mid
36 class2 mid
37 class2 >60
38 class2 >60
39 class2 >60
40 class2 >60
41 class2 >60
42 class3 <10
43 class3 <10
44 class3 <10
45 class3 mid
46 class3 mid
47 class3 >60
48 class3 >60
49 class3 >60
50 class3 >60
")

Data$age = ordered(Data$age, levels=c("<10", "mid", ">60"))

kruskal.test(as.numeric(age) ~ class, data=Data)

   ### Kruskal-Wallis rank sum test
   ###
   ### Kruskal-Wallis chi-squared = 0.12691, df = 2, p-value = 0.9385

plot(as.numeric(age) ~ class, data=Data)

Table = xtabs(~ age + class, data=Data)

Table

   ###      class
   ### age   class1 class2 class3
   ### <10      8      4      3
   ### mid      6     10      2
   ### >60      8      5      4

library(coin)

chisq_test(Table, scores = list("age" = c(-1, 0, 1)))

   ### Asymptotic Generalized Pearson Chi-Squared Test
   ###
   ### data:  class by age (<10 < mid < >60)
   ### chi-squared = 0.13118, df = 2, p-value = 0.9365

library(coin)

spineplot(Table)
$\endgroup$
1
  • $\begingroup$ Thank you! Well, the Kruskal-Wallis does not test what I need. It looks at ranked values, roughly related to the age value itself, while I need to know if the % structure of the 3 classes change. In such analysis it's more important when there is exceeding % of certain category, as their importance is not linearly related to the value of age itself. Just if there are more "oldies" than others. The two hypotheses may result in similar results, but are not equivalent. In each class of quality the % should be same, but mean/median ages may differ. This is why I focused on the %. $\endgroup$
    – CroissSant
    Dec 6, 2022 at 16:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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