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.
