Analysis of proportions in repeated experiments

Question

What would be the best solutions to analyse proportions of multiple classes.

The classes are mutually exclusive, with each observation falling into one of the classes (i.e. no other class is possible) and are unordered.

For example, having the following data, where observations in the same experiments are dependent (i.e. for the same Exp A+B+C = 100% in each condition)

Percentage	Counts	Class	Group	Exp
80	2475	A	Control	1
81	2906	A	Control	2
76	5102	A	Control	3
11	351	B	Control	1
11	405	B	Control	2
14	938	B	Control	3
9	280	C	Control	1
8	282	C	Control	2
10	639	C	Control	3
75	2380	A	Treated	1
76	2584	A	Treated	2
69	2210	A	Treated	3
15	480	B	Treated	1
15	510	B	Treated	2
16	530	B	Treated	3
10	310	C	Treated	1
9	290	C	Treated	2
15	480	C	Treated	3

prop_data <- data.frame(Percentage = c(80, 81, 76, 
                                       11, 11, 14, 
                                       9, 8, 10,
                                   75, 76, 69, 
                                   15, 15, 16, 
                                   10, 9, 15   
        ),
        Counts = c(2475, 2906, 5102,
                   351, 405, 938,
                   280, 282, 639,
                   2380, 2584, 2210,
                   480, 510, 530,
                   310, 290, 480
        ),
        Class = factor(rep(c(&quot;A&quot;, &quot;B&quot;, &quot;C&quot;), each = 3)),
        Group = factor(rep(c(&quot;Control&quot;, &quot;Treated&quot;), each = 9)),
        Exp = 1:3)

I would like to see whether the proportion of A,B and C is dependent on the Group

I thought of doing something like

model <- glm(Percentage ~ Class * Group, data = prop_data, family = poisson)

I am not, however, completely sure whether this is correct. For example, I am not dealing with counts but rather with percentages (so, not integers) here. Also, should I include the experiments as a random factor (I tried using glmer but I get a singular fit)?

I have also considered using a multinomial model but am not quite sure how to set it up.

EDIT:

I have added raw counts to the R code, as suggested in the comments

Answer 1

With "A," "B," and "C" as mutually exclusive and unordered outcome classes in the trials of each combination of experiment (Exp) and treatment Group, this is a standard application of multinomial logistic regression.

One Class is taken as the reference and regression coefficients then express predictor-associated differences in log-odds of each of the other classes with respect to that reference class. This differs from a set of binomial regressions against the reference class, as the probabilities among classes are constrained to sum to 1. The principles are nicely outlined in this UCLA web page, which recommends the multinom() function in the R nnet package

As your data are already aggregated, you use the Counts as weights for each of the outcomes. Treat Exp as a fixed effect after making sure it's interpreted as a factor rather than numeric. The following allows for the effect of the treatment Group to differ among experiments.

library(net)
prop_data$Exp <- factor(prop_data$Exp)
mnMod <- multinom(Class~Group*Exp,data=prop_data,weights=Counts)
## some output about convergence is omitted
summary(mnMod)$coefficients
##   (Intercept) GroupTreated        Exp2      Exp3
## B   -1.953216    0.3521456 -0.01742925 0.2595691
## C   -2.179205    0.1409320 -0.15352828 0.1017433
##   GroupTreated:Exp2 GroupTreated:Exp3
## B      -0.004179719       -0.08637392
## C       0.004549808        0.40957042

The overall reference is for Class "A" under Group "Control" in Exp 1. Intercepts for "B" and "C" are log-odds differences from "A" at those reference predictor values. Regression coefficients are further log-odds differences associated with Group, Exp, and the interactions between them.

For these sample data you would need to pay attention to systematic differences among the experiments, as some Exp and Group:Exp (interaction) coefficients are quite significant:

summary(mnMod)$coefficients/summary(mnMod)$standard.errors
##   (Intercept) GroupTreated       Exp2     Exp3
## B   -34.24556     4.641310 -0.2237779 3.862921
## C   -34.56243     1.614329 -1.7310543 1.343344
##   GroupTreated:Exp2 GroupTreated:Exp3
## B       -0.04000224        -0.8928776
## C        0.03672621         3.7516681

The coefficient estimates are assumed to have a limiting normal distribution, so these ratios to their standard errors are z-scores. A ratio with absolute value > 1.96 would be considered significant at p < 0.05.

You might wonder what happened to all of the information about Class "A." That's implicit in the constraint that the probabilities among all classes sum to 1.

Post-modeling software like the emmeans package allows you to extract that information and perform appropriate post-hoc tests. For example, if one could ignore the experiment-specific outcomes and average over the experiments, you would get the following report of Class probabilities and standard errors as a function of Group:

library(emmeans)
emmeans(mnMod,~Class|Group)
## Group = Control:
##  Class   prob      SE df lower.CL upper.CL
##  A     0.7898 0.00368 12   0.7818   0.7979
##  B     0.1221 0.00295 12   0.1156   0.1285
##  C     0.0881 0.00257 12   0.0825   0.0937
## 
## Group = Treated:
##  Class   prob      SE df lower.CL upper.CL
##  A     0.7336 0.00446 12   0.7239   0.7433
##  B     0.1556 0.00367 12   0.1476   0.1636
##  C     0.1109 0.00317 12   0.1039   0.1178
## 
## Results are averaged over the levels of: Exp 
## Confidence level used: 0.95 

That suggests that "Treated" leads to an increase in probabilities of "B" and "C" with a corresponding decrease in probability of "A."