Open behinger opened 6 years ago
weiterer Gedanke (und sorry für die Sprachwechsel): Da wir pro Subject nur wenige Trials haben, ist die Wahrscheinlichkeit hoch, dass wir die Varianz zwischen ihnen überschätzen, richtig? Ich würde daher den Prior für die sn_idx group etwas enger schnüren als mit cauchy(0,1) -- z.B. exp(1). Was denkst du dazu? Allerdings bekomme ich dann 17 divergent transitions. Wie geht man damit um?
bzgl. BayesFactor vs. WAIC: Ich denke bei Modality und Abstraction macht es Sinn, das ganze Ding rein / raus zu nehmen. Dann können wir evtl. die Aussage treffen: "Die Modality / das Abstraction Level zu kennen beeinflusst unsere Predictions nicht nennenswert". Und das wäre doch ne schöne klare Aussage :)
Ach ja, und unten bei den Plots waren noch ein paar Fragen die sich darauf beziehen, ob ich die Plots richtig verstanden habe ;-)
Problem der Varianz Schätzung: das wäre mir neu. Eher das man die einzelnen Subjectvarianzen nicht schätzen kann.
Loo können wir gerne machen. Einfach loo(modelWithA,modelWithoutA) (package loo)
Plot fragen schau ich mir Zuhause an
Sind die Ergebnisse mit cauchy(0,1) denn deiner Meinung nach publizierbar? Oder gibt es da noch Bedarf etwas zu ändern? Ich würde die Auswertung eigentlich gern übers Wochenende abschließen und Montag Gordon präsentieren.
Die Ergebnisse der validation zeigen leicht niedrigere (bessere) Werte für die stärkeren exponential priors, allerdings gibt's da im größten Model (mit modality und abstraction) einige divergent transitions.
Ich habe, wie du im Code siehst, auch eine 10-fold cross-validation gemacht. Leider ist nicht das model ganz ohne modality und abstraction das beste. Schau mal rein und verrat mir, was du meinst, wie sich das auf unsere Aussagen auswirkt. Oder ob wir doch eher mit BayesFactor arbeiten sollten.
Ich hatte übrigens Probleme, einen posterior predictive check durchzuführen. Ich kenne einfach diese ganzen Model-Objekte nicht gut, und weiß daher nicht genau was da passiert und wie ich das am besten inspiziere, was da raus kommt. An dem Punkt könnte ich noch mal Hilfe gebrauchen.
Ja cauchy ist gut (wenn keine divergent transitions dabei waren). Ich würde statt exp(1) entweder ne normalverteilung oder student_t(4) oder so nehmen. Cauchy hat sehr große tails. Exp(1) fällt sehr schnell ab (wo hast du den exp(1) her?). An diesem Punkt muss man aber ehrlich gesagt sagen, dass unser Prior partially von den Daten/Ergebnissen abhängt. Das sollte man im Paper diskutieren
das es ~80 datenpunkte gibt die nicht vom Model gut modeliert werden können, ist etwas schräg. Ich könnte mir vorstellen, dass wen du nen logit von +10 hast, dann entspricht das 99.9999% erwartung einer "1". Wenn dann doch ne 0 kommt, ist das ein Problem. Bin aber nicht tiefgenug in importance sampling drin um das genauer zu wissen (in der doku steht, dass WAIC nicht reliable ist in dem Kontext).
loo_model_weights <-- kannst du noch weights (wie bei WAIC) kriegen Es scheint als unterscheiden sich die Modelle nicht wirklich (sehr große SE bei den differenzen), das liegt vllt wiederum daran, dass wir nahe an separability dran sind und somit schwer zu sagen ist, ob sie sich nicht doch unterscheiden(?). Ich würde das reporten, diskutieren und mich dann auf die posterior-density estimates berufen und zeigen dass quasi alle abstraction/modality um null rum sind.
hier ist schwer für mich zu sagen weil ich kein modelobject da habe. Werde ich versuchen später zu überprüfen.
Die Errorbars sind genau mean+ SE (was nicht superviel sinn macht bei binären daten, aber das ist am schnellsten zu berechnen). Solltest du daten plotten wollen, sollte man bootstrapped CIs berechnen (geom_summary(method="boot_ci")
)
Deine Interpretation der Plots klingt gut. Alleine die Datenplots zu zeigen würde mich schon überzeugen ;-)
Prior exponential(1) habe ich aus statistical rethinking, da nimmt er es um die heavy tails der cauchy distribution zu vermeiden. aber wir können es bei cauchy(0,1) belassen, denke ich.
Model Comparison Meinst du nicht, 10-fold cross-validation ist angebrachter als loo in diesem Fall? Ich würde die vorschlagen, die Werte der vier verschiedenen Models aus der 10-fold cross-validation zu reporten und zusätzlich die Posterior distributions der Parameter, um zu argumentieren, dass sich da nicht viel tut.
Posterior Predictive Hast du schon was? :-)
Accuracy Ich denke es macht vielleicht Sinn, accuracy Werte des Models zu berichten, also einfach mean accuracy in einer 10-fold cross-validation.
Ja sorry, den k10 kannst du auch in die loo_models_weights reinwerfen.
######
##Posterior Predictives
######
pp = posterior_predict(mres_full,nsamples = 1000)
d1.hum.effect$pp = apply(pp,2,mean)
d1.pp = reshape2::melt(d1.hum.effect,measure.vars=c("pp","choice_left"))
ggplot(d1.pp,aes(x=young_diff,y=value,group=variable,color=factor(variable)))+stat_summary(position=position_dodge(width=0.2))+geom_hline(yintercept=0.5)+facet_wrap(~abstraction)
dnew = d1.hum.effect %>% group_by(modality,abstraction,sex_diff,visonset_left,elderly_diff,young_diff)%>%summarise(choice_left=mean(choice_left))
#dnew = expand.grid(abstraction=c(-0.5,0.5),modality=c(-0.5,0.5),sex_diff=c(-0.5,0,0.5),visonset_left=c(-0.5,0,0.5),elderly_diff=c(-0.5,0,0.5),young_diff=c(-0.5,0,0.5))
pp_newsub = posterior_predict(mres_full,newdata = data.frame(dnew),allow_new_levels = TRUE,sample_new_levels="gaussian")
pp_uncert = posterior_predict(mres_full,newdata = data.frame(dnew),allow_new_levels = TRUE,sample_new_levels="uncertainty")
dnew$pp_newsub = apply(pp_newsub,2,mean)
dnew$pp_uncert = apply(pp_uncert,2,mean)
tmp = d1.hum.effect %>% group_by(modality,abstraction,sex_diff,visonset_left,elderly_diff,young_diff,sn_idx)%>%summarise(choice_left=mean(choice_left))
pp_samesub = posterior_predict(mres_full,newdata = data.frame(tmp),allow_new_levels = FALSE,sample_new_levels="gaussian")
tmp$pp_samesub = apply(pp_samesub,2,mean)
dnew$pp_samesub = tmp %>% group_by(modality,abstraction,sex_diff,visonset_left,elderly_diff,young_diff)%>%summarise(pp_samesub=mean(pp_samesub))%>%.$pp_samesub
#dnew$pp_samesub = apply(pp_samesub,2,mean)
dnew.pp = reshape2::melt(dnew,measure.vars=c("pp_samesub","pp_newsub","pp_uncert","choice_left"))
ggplot(dnew.pp,aes(x=young_diff,y=value,group=variable,color=factor(variable)))+stat_summary(position=position_dodge(width=0.2))+geom_hline(yintercept=0.5)+facet_grid(modality~abstraction)+ylim(c(0,1))
Die Farben:
Habe alle parameter einmal durchgecheckt, mir ist nur hier was aufgefallen. Wie du rechts unten siehst, stimmen model & daten nicht unbedingt überein.
Was das bedeuten könnte ist das die Subjects nicht der multivariaten normalverteilung folgen. Konsequenz ist, dass das model den young_diff effekt leicht unterschätzt (wenn ich das alles richtig verstehe - mache das aber auch so zum ersten mal ;-)).
ich glaube das ist nichts worüber wir uns sorgen machen müssen. Reporten, diskutieren, nicht so wichtig
bin mir nicht sicher was du mit accuracy meinst?
Accuracy Damit meine ich, wie viel % er korrekt vorhersagt in einer 10-fold x-validation. Also parameter und trial (dh obstacles, visonset lane, etc.) in die Funktion stecken, und dann die Probability bei 50% thresholden, um damit einen 0/1 output zu erzwingen, das dann mit der Entscheidung des Probanden abgleichen. Einfach damit man sagen kann "dieses Model hätte in X% der Fälle das Verhalten der Probanden korrekt vorher gesagt". Kann die k-fold Funktion so einen Accuracy-Wert irgendwie ausspucken, oder müsste ich das manuell machen? Man müsste dafür natürlich auch irgendwie jedes Sample aus den Daten mit einer Distribution über die Subjects validieren, wenn man nicht das genaue Subject-Parameter-Matching betreiben will. Ein bisschen wie mit den posterior predictives, falls du verstehst was ich meine?
Sonstiges In logreg_analysis_4.R findest du jetzt die Sachen, die ich vorhabe im Paper zu behandeln bzgl. dieses Datensatzes. Ich will das ganze nicht extrem detailliert mit allem drum herum auseinander Pflücken, also sowas wie die pp checks eher nicht erwähnen (außer vllt. im Appendix, wenn du meinst, dass die wichtig sind). Schau mal ob du mit allem d'accord bist.
Accuracy Die deviance ist doch eigentlich genau das, nur das man nicht nen 50% cutoff macht, sondern eine Gewichtung hat (d.h. wenn ich 90% ne 1 vorhersage, ist ne "1" zu sehen besser, als wenn ich 50% vorhersage). Ich sehe den vorteil von accuracy nicht gegenüber likelihood
sonstiges Stimme dir da überein, man könnte ein online notebook erstellen mit allen posterior predictive checks. Die interpretation ist ein bisschen tricky zu beschreiben. Im endeffekt denke ich aber dass beides keinen einfluss hat
Plots
Welche plots willst du ins paper nehmen? Ich würde die Daten zusätzlich noch plotten (nicht nur die Parameter estimates)
Accuracy Ich meinte nicht stattdessen, sondern nur zusätzlich, um ne einfache Zahl zu haben, mit der jeder was anfangen kann.
Plots Arbeite ich später dran, melde mich wenn ich was hab ;-)
Single-Level vs. Multi-Level Das single-level model entspricht ja sozusagen einem multi-level model mit unendlich starken 0-priors für die Varianz der sn_idx parameter, richtig? Je schwächer der Prior im multi-level model wird, desto größere Varianz kriegen wir im den sn_idx-posteriors (bis zu einem gewissen Punkt wo der Prior uninformative wird). Gleichzeitig werden die Population-Level effects größer, um weiterhin die "gleichen" predictions zu machen. Letzteres ist etwas unintuitiv, dass bei größerer Varianz auf dem sn_idx level sich der population mean verschiebt, was aber wohl am logit-link liegt. Varianz auf dem sn_idx-level in Richtung weg von der Null fällt sozusagen nicht mehr ins Gewicht, weil die Predictions eh schon bei 99%+ sind, in die andere Richtung aber fällt sie ins Gewicht -- da bedeutet eine std dann mitunter 20-30% Wahrscheinlichkeit. Insofern macht es sinn, dass sich der mean weg von der Null verschiebt im Multi-Level Model.
Die Schätzungen der Varianzen auf sn_idx-level sind übrigens ziemlich die gleichen im zweiten, größeren Datensatz mit wesentlich mehr Datenpunkten pro Proband. In sofern glaube ich nicht, dass wir hier die Varianz falsch Schätzen.
Dennoch klingen die Population-Level Effects im multi-level model natürlich zunächst extrem groß. Ich denke, was wir hier sehen ist, dass die meisten Probanden kategorisch immer die jüngeren Personen schützen, und ein paar wenige dort eher zufällig entscheiden oder sogar die älteren Schützen. Das könnte man als Erkenntnisgewinn und Ansatz für Überlegungen im Design von Nachfolge-Studien verstehen (und da habe ich ein paar).
Accuracy
# get 1000 random samples from the posterior (already 0&1's)
p = posterior_predict(mres_full,nsamples=1000)
#calculate mean for each data point over posterior (mean prediction accuracy for each data point)
d1.hum.effect$accuracy = apply(apply(p,1,function(x)(x==d1.hum.effect$choice_left)),1,mean)
#plot it subjectwise (that is the accuracy each subject has
ggplot(d1.hum.effect%>%group_by(sn_idx)%>%summarise(acc=mean(accuracy)),aes(x=0,y=acc))+ggbeeswarm::geom_quasirandom() + stat_summary(method="mean_cl_boot",color='red',size=1)+xlim(c(-1,1))+ylim(c(0,1))
Single-Level vs. Multi-Level single-level, genau - man geht davon aus dass die subjects sich nicht unterschieden, d.h. varianz=0 (mega starker prior)
Dass die population level effects größer werden wenn man von single auf multilevel geht ist überhaupt nicht klar. IIch stimme dir zu, liegt vermutlich am logitlink. Bei normalen mixed models sind die "mean-effects" genau gleichgroß.
Ich würde mir die single-subject effect-mean estimates gerne anschauen. Vielleicht können wir so subjects identifizieren die gegen den Strom schwimmen und die großen Varianzen (und somit Effekte) erzeugen.
Ich mache jetzt einfach noch mal eine Feststellung bzgl. der großen Haupteffekte:
Die Haupteffekte sind im single-level model schon groß genug, um die Alterseffekte in den Sättigungsbereich der Sigmoidfunktion zu schieben (also wo sie sehr flach wird). Dadurch Wirkt sich die per-subject Varianz für die Predictions nur noch in eine Richtung aus, nämlich dann, wenn die subjects wenig oder keine Unterschiede zwischen den Altersstufen machen. Das sind unsere gegen-den-Strom-Schwimmer.
Es sieht also so aus, dass die meisten Probanden einfach kategorisch immer die jüngeren Schützen, und daher gemeinsam am einen Extrem der Streuung sitzen. Die Gegenschwimmer sind sozusagen die Outlier, die überhaupt dafür sorgen, dass es keine absolut eindeutigen Entscheidungen sind. Wenn ich das richtig sehe, fittet das Model aufgrund dieser Gegebenheit auf den Haupteffekten ziemlich genau den "Strom", also nicht wirklich die echte Mitte der Population, weil es die Gegenschwimmer in die Varianz "auslagern" kann.
Um diese These zu überprüfen, müssten wir also einen Überblick über die Shapes der Verteilungen bekommen.
Es sieht wirklich so aus wie wir uns das überlegt hatte (der population mean ist noch nicht draufgerechnet).
qplot(ranef(mres_full)$sn_idx[,1,]%>%data.frame()%>%.$elderly_diff)
Jap. Ich wäre daher dafür, eher das single level model als Haupt-Modell zu reporten, und dann die Varianz in der Population in einem eigenen Paragraphen zu bringen, und dabei vielleicht einfach die geschätzte sd für die parameter aus dem multi-level model zu reporten, evtl. mit dem Plot den du grad gemacht hast für alle 3 Parameter als Unterstützung. Was meinst du?
Bayes Factor Das Thema würde ich gern noch mal eben auf machen, bzgl. "Signifikanz" von einzelnen Parametern. Ich werd da direkt noch mal nach googlen, und komme dann evtl. mit Fragen zurück ;)
Ich bin stark gegen das single-level model. Das single-level model ist einfach falsch, weil within-subject Lieber habe ich Verteilungen die auf dem höheren level nicht normalverteilt sind. Ich würde das model so reporten und die nicht-normalität diskutieren. Bei den Effektgrößen macht das nicht viel und zieht den Effekt auch eher noch gegen 0 hin.
Bayes factor tutorial in brms
Ok, können wir im Model irgendwie die normal-assumption ändern? Ich habe immer noch die Befürchtung, dass das Multi-Level Model die Population-Level Effekte überschätzt aufgrund der Thematik, die wir am Board besprochen haben.
siehe neues Ticket, ich glaube langsam dass das ein Single subject model die Effekte unterschätzt.
Ja das ist theorethisch möglich, ich muss mal schauen ob man das direkt in BRMS machen kann
BTW aus dem brms tutorial: "Here’s a short post on..." hust ;-)
oh und bezüglich Intercept-Prios in brms:
The reason is that brms by default uses a little trick in parameterizing the intercept which speeds up the MCMC sampling. In order to specify a prior for the intercept, you’ll have to take the default intercept out (0 +), and use the reserved string intercept to say that you mean the regular intercept.
Ich glaube das würde ich gerne in unseren models auch machen, um einen Prior auf das Intercept setzen zu können.
Bayes Factor Das Tutorial habe ich durchgelesen und versucht anzuwenden, aber ich bekomme nur NAs und weiß nicht warum:
# intercept
h1 <- hypothesis(mres_full, "intercept = 0")
h1
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (intercept) = 0 -0.17 0.16 -0.49 0.15 NA NA
---
'*': The expected value under the hypothesis lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
Wenn ich die Hypothese mit < formuliere bekomme ich Werte, aber die entsprechen ja nicht ner sauberen Nullhypothese.
Ich lese gerade auch viel, dass Leute den BF nicht mögen, weil er so stark von den Priors abhängig ist. Meinst du nicht, wir können anhand der distributions, bzw. MAP/SD ratio bzw. CIs argumentieren? Damit würde ich mich aktuell wohler fühlen.
Habe das "cowplot" package benutzt weli ich grid_arrange nicht kenne :)
# extract all random effects, add the fixed effects and put them to a nice data frame
r = ranef(mres_full)$sn_idx
for (s in seq(dim(r)[1])){
r[s,,] = r[s,,]+t(fixef(mres_full))
}
r = reshape2::melt(r) %>% magrittr::set_colnames(c('subject','column','effect','value'))
r = reshape2::dcast(r,subject+effect~column)
r$type = 'hierarchical'
colnames(r)[c(5,6)] = c("2.5%","97.5%")
theme_Publication <- function(base_size=14, base_family="helvetica") {
library(grid)
library(ggthemes)
(theme_foundation(base_size=base_size, base_family=base_family)
+ theme(plot.title = element_text(face = "bold",
size = rel(1.2), hjust = 0.5),
text = element_text(),
panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour = NA),
axis.title = element_text(face = "bold",size = rel(1)),
axis.title.y = element_text(angle=90,vjust =2),
axis.title.x = element_text(vjust = -0.2),
axis.text = element_text(),
axis.line = element_line(colour="black"),
axis.ticks = element_line(),
panel.grid.major = element_line(colour="#f0f0f0"),
panel.grid.minor = element_blank(),
legend.key = element_rect(colour = NA),
legend.position = "bottom",
legend.direction = "horizontal",
legend.key.size= unit(0.2, "cm"),
legend.margin = unit(0, "cm"),
legend.title = element_text(face="italic"),
plot.margin=unit(c(10,5,5,5),"mm"),
strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
strip.text = element_text(face="bold")
))
}
p_individual = ggplot(r%>%filter(effect%in%c("sex_diff","young_diff","elderly_diff")),aes(x=Estimate))+
geom_histogram(binwidth=0.25)+
scale_x_continuous(labels=exp,breaks = log(10^seq(-9,9,3)),name = "Odds Ratios")+
facet_grid(effect~.,scales = "free_y") +
expand_limits(x = 0)+ # be sure to include the "0", in the way I plot the data now this is not necessary
scale_y_continuous(name="",breaks = NULL)+
geom_vline(xintercept = 0,linetype='dotted')+
theme_Publication()
# plot
cowplot::plot_grid(p1+theme_Publication(),p_individual,nrow=1,rel_widths=c(2,1))
Ich habe noch die lines bei "0" durch dotted lines ersetzt, aber irgendwie sit mein R abgestürzt und plottet nichtsmehr...
Wie findest du den plot?
Yes exactly. They are the mean of the hierarchical subject level
The effect of abstraction (changing from text to image, is
0.38
. That is, to get the parameterestimate of abstraction one would do: Intercept + 0.5abstraction (or intercept - 0.5abstraction).Yes.
questions / to do:
vision*abstract => vision + abstract + vision:abstract
I thought it is more explicit to just add the interaction here, we could add the "abstract" main effect twice (because later we write ...abstract*modality) and the program is smart enough to eliminate one of them. Personal preference.
- is the coding on the correct scale to make predictions as suggested above?
Yes absolutely! We could look into posterior predictives. Check out the
pp_check
function.Bayesfactor is a tool usually used to compare the prior-only model and the posterior-model. One often uses the savage-dicky ratio which logic goes as follows: Define a H_0 (e.g. parameter = 0). BF = How much do I believe the H_0 with just the prior / How much do I believe the H_0 after seeing the data.
It's really just dividing those two posterior values. To do so is super simple in BRMS:
h <- hypothesis(m, "intercept = 0.5")
More info & tutorial here
Yes. (I would prefere "loo" values from the package "loo", but we could also do that.) It would answer a slightliy different question though. If you test a single parameter you ask the question: Does this interaction tell me something. If you remove everything related to "A" (so also all the interactions), you would ask: "does knowing 'A' tell me something".
Not sure which one you prefer :)