Antecedentes: Imagine uma pizza cortada em 8 fatias.
[
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.
Respostas:
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:
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.
fonte