You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: docs/source/examples/R/Example scripts/Example_analysis_bayesian.Rmd
+69-11Lines changed: 69 additions & 11 deletions
Original file line number
Diff line number
Diff line change
@@ -9,8 +9,12 @@ output: html_document
9
9
# **Here we show how to perform a bayesian analysis on the data after it has been collected which is very similar to the "simple" analysis. Please see the "simple analysis before this!**
The only difference here is that we set bayesian equal to T (TRUE) and specify the model. Here the model is a predefined Stan model that is inside the src folder called first_model.stan
35
+
The only difference here is that we set bayesian equal to T (TRUE) and specify the model. The models can be found in the src folder inside the .stan files. These are probabilitic models written in stan that are compiled and does the sampling. There are two options at the moment for re-fitting there is the standard cummulative normal aswell as a cummulative normal that incorporates a lapse rate, that sepcifies the minimum and maximum of the tails of the psychometric. This means that a lapse rate of 5% (0.05) means that the psychometric on the lower end is 5% and on the upper end is 95%. The reason to include a lapse rate is that if responses are made that are attentional slips or misclicks in the far end of the stimulus spectrum (high or low) then this is greatly influence the slope of the psychometric if modelled without the lapse rate.
36
+
37
+
The priors of the bayesian model is as follows in the unconstrained space (beta is constained to be positive so is exponentially transformed, and the lapse is constrained between 0 and 0.5 meaning its inv_logit transformed and then devided by 2:
38
+
39
+
alpha ~ normal(0,20);
40
+
beta ~ normal(0,3);
41
+
lambda ~ normal(-4,2);
42
+
43
+
This means that the parameters in the constrained space look like this:
Below there is a visualization of what this extra lapse rate does as well as what the priors of the model means when looking at the psychometric function itself:
y = list(psychometric_nolapse(seq(-80,80,0.1), alpha, beta))
76
+
) %>%
77
+
ungroup %>%
78
+
mutate(id = 1:n()) %>%
79
+
unnest(cols = c(x, y)) %>% mutate(lapse = F) %>%
80
+
ggplot(aes(x = x, y = y, group = id))+
81
+
geom_line(alpha = 0.5)+theme_classic()+ggtitle("Without Lapse rate")
82
+
83
+
84
+
```
85
+
86
+
If you want to change the priors of the Bayesian model, this has to be done inside the Stan scripts. By opening the .stan File and then changing the last couple of lines of code where the syntax is the same as above, it is therefore possible to visualize what the prior distributions for the parameters and also see what they entail (prior predictive checks) for the shape of the psychometric here in the markdown script and then changing them inside the Stan scripts themselves.
87
+
32
88
33
-
**Doing the same as for the simple analysis with bayesian = T**
89
+
**Running the analysis using this bayesian fit invovles the same as for the simple analysis with two addition arguments firstly the flage bayesian needs to be set to T (TRUE), and a model has to be specified. There are at the moment two different models to choose from, one with the lapse rate and one without**
**The results list now also contains a new index called bayesian_plot. This is a list of either 1 or 3 plots. There'll be 1 if you only have one Morality and 3 if you have two (Extero and Intero). Here there is 3 plots**
105
+
**The results list now also contains a new index called bayesian_plot. This is a list of either 1 or 3 plots. There will be 1 if you only have one Morality and 2 if you have two (Extero and Intero). Here there is 3 plots**
47
106
48
107
Lets look at them individually:
49
108
50
109
```{r}
51
110
results$bayesian_plot[[1]]
52
111
```
53
112
54
-
**NOTE: The Import thing to look at for good model convergence is the upper plots: Here we see that all the 4 chains (to the left) seem to capture the same posterior distribution. It is also clear from the trace-plots to the upper right that the chains mix well (hairy catterpillars), meaning good convergence**
113
+
**NOTE: The Import thing to look at for good model convergence is the upper plots: Here we see that all the 4 chains (to the left) seem to capture the same posterior distribution. It is also clear from the trace-plots to the upper right that the chains mix well (hairy catterpillars), meaning good convergence. Lastly looking into whether there are divergences in the sampling process is pivotal, these are stored in the stats file under divergences, if this column is not 0, then trusting the estimates even with good looking chains is not advised. Dealing with divergences for single subjects fits like here involves changing priors and or the model itself (i.e. leaving out or including the lapse rate)**
55
114
56
115
```{r}
57
116
results$bayesian_plot[[2]]
58
117
```
59
-
60
-
61
-
And the combined plot can be found in the last index
62
-
```{r, fig.height=8,fig.width=14}
63
-
results$bayesian_plot[[3]]
118
+
##**Here is the number of mean in both conditions divergences:**
119
+
```{r}
120
+
results$stats$divergences
64
121
```
122
+
Indicating that there are divergences here so perhaps runnning without the Lapse rate would be preferable, or changing the priors.
65
123
66
124
67
125
**Of cause this can be run through several subjects like the "simple" analysis**
results = single_sub_analysis(df, #The raw dataframe
63
-
interoPost = interoPost, #numpy array for the intero (NA if not avaliable)
64
-
exteroPost = exteroPost, #numpy array for the extero (NA if not avaliable)
65
-
bayesian = F, #Bayesian Analysis (TRUE/FALSE)
66
-
model = NA, #Bayesian model here a stan script (NA if Bayesian is FALSE)
67
-
out = here::here("docs","source","examples","R")) #Output directory for results
62
+
results = single_sub_analysis(df, #The raw dataframe
63
+
interoPost = interoPost, #numpy array for the intero (NA if not avaliable)
64
+
exteroPost = exteroPost, #numpy array for the extero (NA if not avaliable)
65
+
bayesian = F, #Bayesian Analysis (TRUE/FALSE)
66
+
model = NA, #Bayesian model here a stan script (NA if Bayesian is FALSE)
67
+
out = here::here("docs","source","examples","R")) #Output directory for results
68
68
```
69
69
70
70
**Note that these analyses can also be run with only one "Modality", important is that either the interopost or exteropost then gets set to NA i.e. the modality you do not have access to!**
0 commit comments