# Comparing two multinomial distributions

Background: Imagine a pizza cut into 8 slices.

On each straight edge of the slice, I insert a magnet with opposite polarities facing outwards. If I separate these components, prevent flipping and shake them, they should form a full pizza.

Now, if I put an additional slice (same dimensions, 1/8th of a full pizza), a full pizza will not always form. It could form clusters of: 4 & 5, 3 & 6, 2 & 7 and 1 & 8.

A model (given by Hosokawa et al. (1994)) gives me the probability of each cluster forming. To validate the model, I do some physical experiments. I conduct 20 trials for each experimental condition.

My results look like this:

Cluster  Theoretical  Physical
3,6:    6.01961132827  4
1,8:    2.77455224377  5
4,5:    6.62198848501  5
2,7:    4.58384794294  6


The aforementioned data is a multinomial distribution (akin to a distribution obtained when rolling a dice). When there are 9 slices, each trial can end in one of 4 states.

In addition to the 9 slice experiment, I have data for a 40 slice (and a couple others) experiment as well. (Please let me know if you would like me to include it here)

Problem: To test the goodness of fit, I would run a Pearson's Chi-squared test. However, since the means of both distributions are "close", I cannot reject the null hypothesis. However, I cannot accept the null hypothesis either.

Question: Is there a better way show how "close" the model gets to the physical experiments? The multinomial equivalent of "standard deviation", or perhaps a confidence interval? Regression with confidence intervals?

Update: A colleague of mine suggested the following approach for regression analysis in R:

d=read.csv("data.csv")
length(unique(d$abs_state)) nrow(d) d$n_componentsf=as.factor(d$n_components) ncomps=9 dsubs=d[d$n_components==ncomps,]
# using exact multinomial test in EMT (better)
library(EMT)
# using Chi square test statistics
EMT::multinomial.test(dsubs$freq_obs, dsubs$freq_pred/sum(dsubs$freq_pred), useChisq=TRUE,MonteCarlo=FALSE) # using p value calculated via Monte Carlo approach EMT::multinomial.test(dsubs$freq_obs,
dsubs$freq_pred/sum(dsubs$freq_pred),
ntrial = 100000,
MonteCarlo=TRUE)

# other exact multinomial test implementation
library(RVAideMemoire)
RVAideMemoire::multinomial.test(dsubs$freq_obs, p=dsubs$freq_pred/sum(dsubs$freq_pred)) # if you're interested in 95% confidence limits # on your multinomial proportions you can # calculate these like this library(DescTools) MultinomCI(dsubs$freq_obs, conf.level=0.95)

library(ggplot2)
library(lattice)

# if you would like to analyse all data together
# you could also analyse the observed frequency
# counts as having approx Poisson expectation,
# using a Poisson GLM
# (or a binomial GLM if you put cbind(freq_obs,freq_notobs) as your dependent var and use family=binomial)
library(afex)
library(MASS)
library(car)
library(effects)
set_sum_contrasts() # use effect coding, as data is unbalanced

fit=glm(freq_obs ~ freq_pred,
family=quasipoisson, data=d)
summary(fit)
Anova(fit, test="LR", type="III")

plot(allEffects(fit),type="response", rug=F)
plot1=plot(Effect(c("freq_pred"),
mod=fit),type="response",rug=F)
plot2=xyplot(freq_obs ~ freq_pred,
data=d,
type=c("p"),
par.settings=
simpleTheme(col=c("blue","red"),pch=16,cex=1.1))
library(latticeExtra)
plot1+plot2 # nice correspondence between observed and expected frequencies

# sticking model predictions in data frame and
# plotting this in ggplot2 instead
df=data.frame(effect(term="freq_pred",
mod=fit,
xlevels=list(freq_pred=seq(0,15,1)),type="response"))
df

library(ggplot2)
library(ggthemes)
ggplot(data=df) +
geom_ribbon(data = df,
aes(x=freq_pred, ymin=lower, ymax=upper), linetype=2, alpha=0.1) +
geom_line(data = df, aes(x=freq_pred, y=fit)) +
geom_point(data = d, aes(x=freq_pred, y=freq_obs, color=n_componentsf), cex=3) +
ylab("Predicted frequency") +
xlab("Observed frequency") +
theme_few(base_size=16)


Data.csv has:

n_components,abs_state,freq_pred,freq_obs,freq_notobs
9,"3,6",6.019611328,4,16
9,"1,8",2.774552244,5,15
9,"4,5",6.621988485,5,15
9,"2,7",4.583847943,6,14
15,"3,6,6",1.81009773,0,20
15,"4,5,6",3.927622683,7,13
15,"7,8",13.49657189,13,7
15,"5,5,5",0.765707695,0,20
40,"4,5,5,5,6,7,8",0.803987454,0,20
40,"5,6,6,7,8,8",3.376961833,3,17
40,"2,7,7,8,8,8",0.230595232,0,20
40,"5,7,7,7,7,7",0.175319358,0,20
40,"5,6,7,7,7,8",2.709196442,1,19
40,"5,5,5,5,6,7,7",0.170797287,0,20
40,"4,5,5,6,6,7,7",0.847746645,1,19
40,"8,8,8,8,8",0.230119576,0,20
40,"4,6,7,7,8,8",1.795116571,0,20
40,"6,6,6,6,8,8",0.483846181,0,20
40,"5,5,5,5,6,6,8",0.185675465,1,19
40,"4,7,7,7,7,8",0.4072505,0,20
40,"5,5,5,6,6,6,7",0.274814936,1,19
40,"5,5,6,8,8,8",0.708881447,1,19
40,"4,6,6,8,8,8",0.479526825,1,19
40,"4,5,5,5,5,8,8",0.071967085,0,20
40,"5,5,7,7,8,8",1.233981848,1,19
40,"4,5,6,6,6,6,7",0.34756786,1,19
40,"6,6,6,7,7,8",2.208449237,2,18
40,"4,5,5,6,6,6,8",0.494611498,1,19
40,"5,5,5,5,5,7,8",0.061650486,0,20
40,"4,5,5,5,5,5,5,6",0.010322162,0,20
40,"3,6,6,6,6,6,7",0.071274376,0,20
40,"4,6,6,6,6,6,6",0.015244456,0,20
40,"6,6,7,7,7,7",0.656508868,1,19
40,"4,5,5,5,7,7,7",0.155256863,2,18
40,"5,5,6,6,6,6,6",0.04917738,0,20
40,"3,6,7,8,8,8",0.634760319,0,20
40,"3,7,7,7,8,8",0.430171526,1,19
40,"5,5,5,5,5,5,5,5",0.00022644,0,20


The code creates 3 impressive [to me at least :) ] plots. However, I am skeptical about their correctness. From what I understand, the regression curve represents the correspondence between the predictions and observations. This curve is probably the "best fit" given the correspondences. It seems, however, it concerns me that bulk of the corresponding points are made at the start. Is this concern valid?

With references to the regression curve, how is the confidence interval computed? I found the Bauer (1972) paper, but wasn't able to understand it satisfactorily. Help understanding how the confidence interval is computed will be appreciated.

Other (possibly rudimentary) questions: - What is the effect of not observing a state on the regression plot? For instance, the 1st and 4th states of the 15 slice experiment are not observed physically? - In the third plot (which you obtain after executing the R code), is there a particular reason why the axes might be flipped?

Finally, what conclusions can I draw from the plots? Most of the correspondence points on the regression plot are outside the confidence interval. Is this a foolproof indication of the model's (in)ability to make predictions?

Lastly and obviously, the data is insufficient. In order to find the required amount of data, one would have to conduct sufficient number of trials until the confidence interval is reduced to a desirable width (I've come to realize that choosing this width is not trivial). I've tried random replication of predicted frequencies (poor man's bootstrapping?), but once again, I run into the problem of not understanding how the confidence interval is computed from the Pearson's chi-squared statistics.

• This comment might not be entirely helpful, but since the multinomial distribution is multivariate, the analog of the variance would be the (co)variance matrix en.wikipedia.org/wiki/Covariance_matrix I am not familiar enough with multivariate statistics to give recommendations for specific tests to use with the covariance matrix. – Chill2Macht Apr 15 '17 at 18:11
• @Chill2Macht thanks for your super quick response! I did consider the covariance matrix, but a couple of concerns: (1) I'm not sure what in the covariance matrix I would have to use, and (2) "standard deviation" is (var)^0.5. So, what would that be in a covariance matrix? My grip on stats is not good enough to figure this out by myself. – troisquatre Apr 15 '17 at 18:19
• Well a matrix is a covariance matrix of some distribution if and only if it is positive semi-definite, and every positive semi-definite matrix has a "square root" linear.ups.edu/jsmath/0210/fcla-jsmath-2.10li108.html -- as for how to apply this fact in a statistical context, I have no idea. It would be the blind leading the blind. I wasn't sure if you had thought of the covariance matrix, and sometimes knowing terminology can make accidentally stumbling on the right result on google a lot easier. But other than that I don't really know, I'm sorry. – Chill2Macht Apr 16 '17 at 10:30
• @Chill2Macht : Thanks a lot for the input. If you know someone who might be able to shed light on my problem, please do share my questions with that person. – troisquatre Apr 17 '17 at 14:59
• Are you interested in a measure of uncertainty, inference or a descriptive measure? – Todd D Dec 10 '17 at 5:39

## 1 Answer

I believe the best test for your data is the multinomial test, which you already estimated.

For a "measure" of the difference between the expected and observed distributions, consider the Kullback-Leibler Divergence. For your first experiment:

require(entropy)
cv <- d[1:4,]
set.seed(42)
probs <- cv$freq_pred/20 cv$freq_exp <- table(sample(c(1:4), 10000, prob=probs, replace=T))
KL.empirical(cv$freq_exp, cv$freq_obs)


The KL Divergence is 0.072. If 0, the 2nd distribution is similar to the 1st distribution. If 1, the 2nd distribution is less related to the 1st distribution. More here for a finer interpretation based in information theory.