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
## **Here we show how to perform a Hierarchical Bayesian post-hoc analysis on the data.**
23
+
24
+
This type of analysis is similar to the post-hoc single fit bayesian analysis, however here we fit all our subjects together.
25
+
26
+
This is analogous to the difference between running many separate linear regression and then aggregating the results, compared to running a linear mixed effects model. (This script is what a linear mixed effects model is).
27
+
28
+
29
+
30
+
## *Running this model*
31
+
32
+
For the moment, running this type of hierarchical model, opens the black-box of the sampling algorithm from the single fit analysis.
33
+
34
+
Firstly, we need to specify the model:
35
+
36
+
```{r}
37
+
mod = cmdstanr::cmdstan_model(here::here("docs","source","examples","R","src","Hierarchical Cummulative Normal.stan"))
38
+
```
39
+
40
+
and load the data
41
+
42
+
```{r}
43
+
data = read.csv(here::here("docs","source","examples","R","data","Hierarchical data","Hierarchical_data.csv"))
44
+
```
45
+
46
+
47
+
This data set is simulated and represents a within subject design with two sessions. We can visualize the aggregated (population level) effects for each session:
48
+
49
+
```{r}
50
+
data %>% ggplot(aes(x = Alpha, y = resp, col = as.factor(sessions)))+
51
+
geom_smooth(method = "glm",
52
+
method.args = list(family = "binomial"),
53
+
se = FALSE)+theme_classic()
54
+
```
55
+
56
+
To run the model one needs to run the following:
57
+
58
+
This just aggregates the data for the stan model to run a bit faster. Note that you will need a participant_id column that is a unique identifier for each subjects and a sessions column that is a unique session identifier.
59
+
60
+
```{r}
61
+
data = data %>% mutate(sessions = ifelse(sessions == 1, 0 ,1))
62
+
63
+
data = transform_data_to_stan(data) %>% arrange(participant_id,sessions)
64
+
```
65
+
66
+
67
+
Now one just needs to put it all into one big list as below and run the model!
68
+
69
+
70
+
Note! Here 3 matrices are specified "X_lapse", "X_alpha", "X_beta". These represent parameterizations of these three parameters of the model.
71
+
72
+
What this means is that you can specify you desired constast of interrest by adding a column to this matrix, the model will then give the estimate (as in a linear model), of how much this "covariate" explains. Here we are interrested in the difference in slope (beta) and threshold (alpha) between sessions.
Copy file name to clipboardExpand all lines: docs/source/examples/R/Example scripts/Example_analysis_bayesian.Rmd
+79-13Lines changed: 79 additions & 13 deletions
Original file line number
Diff line number
Diff line change
@@ -6,7 +6,7 @@ output: html_document
6
6
---
7
7
8
8
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!**
9
+
# **Here we show how to perform a Bayesian post-hoc analysis on the data, 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. 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
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:
37
+
The difference between this example and the simple analysis is that here we set Bayesian equal to T (TRUE) and specify a model.
38
+
39
+
The models can be found in the src directory inside the .stan files. These are probabilistic models written in stan.
40
+
41
+
42
+
There are two options at the moment for re-fitting the data using this Bayesian model, there is the standard cumulative normal as well as a cumulative normal that incorporates a lapse rate, that specifies the minimum and maximum of the tails of the psychometric.
43
+
44
+
45
+
A lapse rate of 5% (0.05) means that the psychometric function (the cumulative normal) on the lower end is 5% and on the upper end is 95%.
46
+
47
+
The reason to include a lapse rate is that if responses are made that are attentional slips or miss-clicks, then this is greatly influence the slope of the function if not accounted for.
48
+
49
+
In order to run a Bayesian probabilistic model, one needs to specify priors on the free parameters. Here there are 2 or 3 depending on the model.
50
+
51
+
The priors of the Bayesian model is as follows:
38
52
39
53
alpha ~ normal(0,20);
40
54
beta ~ normal(0,3);
41
55
lambda ~ normal(-4,2);
42
56
43
-
This means that the parameters in the constrained space look like this:
57
+
Note that all the parameters are specified in the unconstrained space, this means that the slope of the psychometric i.e. beta is going to be constrained to be strictly positive and is therefore exponentially transformed.
58
+
The lapse is constrained between 0 and 0.5 meaning its inverse logit transformed and then divided by 2:
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:
69
+
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.
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.
105
+
The addition of the lapse rate is visually obvious when looking at the very steep psychometric functions, which do not asymtote at 0 and 1, but at the lambda value and 1-lambda. This is what is called prior predictive checks.
106
+
107
+
### **Changing the priors**
108
+
109
+
If you want to change the priors of the Bayesian model, this has to be done inside the Stan scripts.
110
+
111
+
Open the .stan file and then navigate to the last couple of lines of code where the syntax is the same as above i.e.
112
+
113
+
alpha ~ normal(0,20);
114
+
beta ~ normal(0,3);
115
+
lambda ~ normal(-4,2);
116
+
117
+
These lines of code are the priors of the model and it is therefore possible to first visualize what the prior distributions for the parameters entail (prior predictive checks) for the shape of the psychometric and then changing them inside the Stan scripts themselves.
87
118
88
119
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**
120
+
**Running the analysis**
121
+
122
+
123
+
Using this bayesian fit invovles the same as for the simple analysis with two addition arguments.
124
+
125
+
Firstly, the flag Bayesian needs to be set to T (TRUE), and a model has to be specified.
126
+
127
+
There are at the moment two different models to choose from, one with the lapse rate and one without. These are called
**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**
148
+
The results list now also contains a new index called bayesian_plot.
149
+
150
+
This is a list of either 1 or 3 plots. There will be 1 if you only have one Modality and 2 if you have two i.e. both Extero and Intero.
106
151
107
152
Lets look at them individually:
108
153
109
154
```{r}
110
155
results$bayesian_plot[[1]]
111
156
```
112
157
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)**
114
158
115
159
```{r}
116
160
results$bayesian_plot[[2]]
117
161
```
118
-
##**Here is the number of mean in both conditions divergences:**
162
+
# **Convergence and trust in the model**
163
+
164
+
NOTE:
165
+
166
+
When running a Bayesian model like this convergence is not a given. It is therefore important to check good model covergence!
167
+
168
+
In-order to check for good model convergence watch the upper plots in the above two plots:
169
+
170
+
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), indicating good convergence.
171
+
172
+
Lastly, one needs to investigate whether there are divergences in the sampling process.
173
+
174
+
This information is 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.
175
+
176
+
Dealing with divergences for single subjects fits, involves changing priors and or the model itself (i.e. leaving out or including the lapse rate).
177
+
178
+
Other approaches for dealing with these divergences exist, but is out of the scope of this tutorial [see](https://discourse.mc-stan.org/t/divergent-transitions-a-primer/17099)
179
+
180
+
181
+
182
+
## **Here is the number of mean in both conditions divergences:**
119
183
```{r}
120
184
results$stats$divergences
121
185
```
122
-
Indicating that there are divergences here so perhaps runnning without the Lapse rate would be preferable, or changing the priors.
186
+
Indicating that there are divergences here so perhaps running without the Lapse rate would be preferable, or changing the priors.
187
+
188
+
123
189
124
190
125
191
**Of cause this can be run through several subjects like the "simple" analysis**
0 commit comments