Comparando duas distribuições multinomiais

7

Antecedentes: Imagine uma pizza cortada em 8 fatias.

[8 fatias de pizza]

Em cada extremidade reta da fatia, insiro um ímã com polaridades opostas voltadas para fora. Se eu separar esses componentes, evitar sacudi-los e sacudi-los, eles deverão formar uma pizza cheia.

Agora, se eu colocar uma fatia adicional (mesmas dimensões, 1/8 de uma pizza completa), uma pizza completa nem sempre será formada. Poderia formar grupos de: 4 e 5, 3 e 6, 2 e 7 e 1 e 8.

Um modelo (dado por Hosokawa et al. (1994)) me fornece a probabilidade de formação de cada cluster. Para validar o modelo, faço alguns experimentos físicos. Eu conduzo 20 tentativas para cada condição experimental.

Meus resultados são assim:

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

Os dados acima mencionados são uma distribuição multinomial (semelhante a uma distribuição obtida ao rolar um dado). Quando existem 9 fatias, cada avaliação pode terminar em um dos quatro estados.

Além do experimento de 9 fatias, também tenho dados para um experimento de 40 fatias (e algumas outras). (Entre em contato se você quiser que eu inclua aqui)

Problema: Para testar a qualidade do ajuste, eu faria um teste qui-quadrado de Pearson. No entanto, como os meios de ambas as distribuições são "próximos", não posso rejeitar a hipótese nula. No entanto, também não posso aceitar a hipótese nula.

Pergunta: Existe uma maneira melhor de mostrar o quão "próximo" o modelo chega dos experimentos físicos? O equivalente multinomial de "desvio padrão", ou talvez um intervalo de confiança? Regressão com intervalos de confiança?

Atualização: Um colega meu sugeriu a seguinte abordagem para análise de regressão em 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 possui:

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

O código cria três gráficos impressionantes [para mim pelo menos :)]. No entanto, sou cético em relação à correção deles. Pelo que entendi, a curva de regressão representa a correspondência entre as previsões e as observações. Essa curva é provavelmente o "melhor ajuste", dadas as correspondências. Parece, no entanto, que me preocupa que grande parte dos pontos correspondentes seja feita no início. Esta preocupação é válida?

Com referências à curva de regressão, como é calculado o intervalo de confiança? Encontrei o artigo de Bauer (1972), mas não o compreendi satisfatoriamente. Ajudar a entender como o intervalo de confiança é calculado será apreciado.

Outras questões (possivelmente rudimentares): - Qual é o efeito de não observar um estado no gráfico de regressão? Por exemplo, o 1º e o 4º estados do experimento de 15 fatias não são observados fisicamente? - No terceiro gráfico (que você obtém após a execução do código R), existe uma razão específica para que os eixos possam ser invertidos?

Finalmente, que conclusões posso tirar das parcelas? A maioria dos pontos de correspondência no gráfico de regressão está fora do intervalo de confiança. Isso é uma indicação infalível da (in) capacidade do modelo de fazer previsões?

Por fim e obviamente, os dados são insuficientes. Para encontrar a quantidade necessária de dados, seria necessário realizar um número suficiente de tentativas até o intervalo de confiança ser reduzido para uma largura desejável (cheguei à conclusão de que escolher essa largura não é trivial). Eu tentei a replicação aleatória das frequências previstas (bootstrapping do pobre homem?), Mas, mais uma vez, me deparo com o problema de não entender como o intervalo de confiança é calculado a partir das estatísticas qui-quadrado de Pearson.

troisquatre
fonte
Esse comentário pode não ser totalmente útil, mas como a distribuição multinomial é multivariada, o analógico da variância seria a matriz de (co) variância en.wikipedia.org/wiki/Covariance_matrix Não estou familiarizado o suficiente com estatísticas multivariadas para fornecer recomendações para. testes específicos para usar com a matriz de covariância.
precisa saber é o seguinte
11
@ Chill2Macht obrigado pela sua resposta super rápida! Eu considerei a matriz de covariância, mas algumas preocupações: (1) não tenho certeza do que precisaria usar na matriz de covariância e (2) "desvio padrão" é (var) ^ 0,5. Então, o que seria isso em uma matriz de covariância? Meu controle sobre as estatísticas não é bom o suficiente para descobrir isso sozinho.
troisquatre
Bem, uma matriz é uma matriz de covariância de alguma distribuição se, e somente se, é semi-definida positiva, e toda matriz semi-definida positiva tem uma "raiz quadrada" linear.ups.edu/jsmath/0210/fcla-jsmath-2.10li108 .html - não sei como aplicar esse fato em um contexto estatístico. Seria o cego guiando o cego. Eu não tinha certeza se você tinha pensado na matriz de covariância e, às vezes, conhecer a terminologia pode tornar muito mais fácil encontrar acidentalmente o resultado certo no google. Mas fora isso, eu realmente não sei, me desculpe.
precisa saber é o seguinte
@ Chill2Macht: Muito obrigado pela contribuição. Se você conhece alguém que possa esclarecer meu problema, compartilhe minhas perguntas com essa pessoa.
troisquatre
Você está interessado em uma medida de incerteza, inferência ou uma medida descritiva?
Todd D

Respostas:

0

Acredito que o melhor teste para seus dados é o teste multinomial, que você já estimou.

Para uma "medida" da diferença entre as distribuições esperadas e observadas, considere a divergência de Kullback-Leibler. Para sua primeira experiência:

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)

A divergência KL é 0,072. Se 0, a 2ª distribuição é semelhante à 1ª distribuição. Se 1, a 2ª distribuição está menos relacionada à 1ª distribuição. Mais aqui para uma interpretação mais refinada, baseada na teoria da informação.

Todd D
fonte