Hintergrund: Stellen Sie sich eine in 8 Scheiben geschnittene Pizza vor.
[
An jeder geraden Kante der Scheibe füge ich einen Magneten mit entgegengesetzten Polaritäten nach außen ein. Wenn ich diese Komponenten trenne, ein Umdrehen verhindere und sie schüttle, sollten sie eine volle Pizza bilden.
Wenn ich jetzt ein zusätzliches Stück (gleiche Abmessungen, 1/8 einer vollen Pizza) lege, bildet sich nicht immer eine volle Pizza. Es könnte Cluster bilden aus: 4 & 5, 3 & 6, 2 & 7 und 1 & 8.
Ein Modell (gegeben von Hosokawa et al. (1994)) gibt mir die Wahrscheinlichkeit, dass sich jeder Cluster bildet. Um das Modell zu validieren, mache ich einige physikalische Experimente. Ich führe 20 Versuche für jede experimentelle Bedingung durch.
Meine Ergebnisse sehen folgendermaßen aus:
Cluster Theoretical Physical
3,6: 6.01961132827 4
1,8: 2.77455224377 5
4,5: 6.62198848501 5
2,7: 4.58384794294 6 Die oben genannten Daten sind eine Multinomialverteilung (ähnlich einer Verteilung, die beim Würfeln erhalten wird). Wenn 9 Scheiben vorhanden sind, kann jeder Versuch in einem von 4 Zuständen enden.
Zusätzlich zum 9-Slice-Experiment habe ich auch Daten für ein 40-Slice-Experiment (und einige andere). (Bitte lassen Sie mich wissen, ob ich es hier aufnehmen soll.)
Problem: Um die Passgenauigkeit zu testen, würde ich einen Pearson-Chi-Quadrat-Test durchführen. Da die Mittelwerte beider Verteilungen jedoch "nahe" liegen, kann ich die Nullhypothese nicht ablehnen. Die Nullhypothese kann ich jedoch auch nicht akzeptieren.
Frage: Gibt es eine bessere Möglichkeit zu zeigen, wie nahe das Modell den physikalischen Experimenten kommt? Das multinomiale Äquivalent der "Standardabweichung" oder vielleicht ein Konfidenzintervall? Regression mit Konfidenzintervallen?
Update: Ein Kollege von mir schlug den folgenden Ansatz für die Regressionsanalyse in R vor:
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 hat:
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,20Der Code erstellt 3 beeindruckende [zumindest für mich :)] Diagramme. Ich bin jedoch skeptisch hinsichtlich ihrer Richtigkeit. Soweit ich weiß, repräsentiert die Regressionskurve die Entsprechung zwischen den Vorhersagen und Beobachtungen. Diese Kurve ist angesichts der Entsprechungen wahrscheinlich die "beste Anpassung". Es scheint mir jedoch ein Anliegen zu sein, dass ein Großteil der entsprechenden Punkte zu Beginn gemacht wird. Ist dieses Anliegen berechtigt?
Wie wird das Konfidenzintervall unter Bezugnahme auf die Regressionskurve berechnet? Ich fand das Papier von Bauer (1972), konnte es aber nicht zufriedenstellend verstehen. Es wird hilfreich sein, zu verstehen, wie das Konfidenzintervall berechnet wird.
Andere (möglicherweise rudimentäre) Fragen: - Wie wirkt es sich aus, wenn ein Zustand auf dem Regressionsplot nicht beobachtet wird? Zum Beispiel werden der 1. und 4. Zustand des 15-Schicht-Experiments physikalisch nicht beobachtet? - Gibt es im dritten Diagramm (das Sie nach Ausführung des R-Codes erhalten) einen bestimmten Grund, warum die Achsen möglicherweise umgedreht werden?
Welche Schlussfolgerungen kann ich schließlich aus den Darstellungen ziehen? Die meisten Korrespondenzpunkte im Regressionsdiagramm liegen außerhalb des Konfidenzintervalls. Ist dies ein narrensicherer Hinweis auf die (Un-) Fähigkeit des Modells, Vorhersagen zu treffen?
Schließlich sind die Daten offensichtlich unzureichend. Um die erforderliche Datenmenge zu finden, müsste eine ausreichende Anzahl von Versuchen durchgeführt werden, bis das Konfidenzintervall auf eine gewünschte Breite reduziert ist (ich habe festgestellt, dass die Auswahl dieser Breite nicht trivial ist). Ich habe versucht, vorhergesagte Frequenzen zufällig zu replizieren (Bootstrapping des armen Mannes?), Aber ich stoße erneut auf das Problem, nicht zu verstehen, wie das Konfidenzintervall aus den Chi-Quadrat-Statistiken von Pearson berechnet wird.
Antworten:
Ich glaube, der beste Test für Ihre Daten ist der Multinomialtest, den Sie bereits geschätzt haben.
Betrachten Sie für ein "Maß" der Differenz zwischen der erwarteten und der beobachteten Verteilung die Kullback-Leibler-Divergenz. Für Ihr erstes Experiment:
Die KL-Divergenz beträgt 0,072. Bei 0 ähnelt die 2. Verteilung der 1. Verteilung. Wenn 1, ist die 2. Verteilung weniger mit der 1. Verteilung verbunden. Mehr hier für eine feinere Interpretation basierend auf Informationstheorie.
quelle