class: center, middle, inverse, title-slide # Bayesianische Statistik ## Teil 4 <br/> Hierarchische Modelle mit brms ### Andrew Ellis ### Methodenkurs Neurowissenschaft im Computerlab ### 2021-04-14 --- layout: true <!-- Home icon --> <div class="my-footer"> <span> <a href="https://kogpsy.github.io/neuroscicomplab" target="_blank"><svg aria-hidden="true" role="img" viewBox="0 0 576 512" style="height:1em;width:1.12em;vertical-align:-0.125em;margin-right:0.2em;font-size:inherit;fill:#0F4C81;overflow:visible;position:relative;"><path d="M280.37 148.26L96 300.11V464a16 16 0 0 0 16 16l112.06-.29a16 16 0 0 0 15.92-16V368a16 16 0 0 1 16-16h64a16 16 0 0 1 16 16v95.64a16 16 0 0 0 16 16.05L464 480a16 16 0 0 0 16-16V300L295.67 148.26a12.19 12.19 0 0 0-15.3 0zM571.6 251.47L488 182.56V44.05a12 12 0 0 0-12-12h-56a12 12 0 0 0-12 12v72.61L318.47 43a48 48 0 0 0-61 0L4.34 251.47a12 12 0 0 0-1.6 16.9l25.5 31A12 12 0 0 0 45.15 301l235.22-193.74a12.19 12.19 0 0 1 15.3 0L530.9 301a12 12 0 0 0 16.9-1.6l25.5-31a12 12 0 0 0-1.7-16.93z"/></svg></a> Methodenkurs Neurowissenschaft im Computerlab </span> </div> <!-- Name (left) --> <!-- <div class="my-footer"> --> <!-- <span> --> <!-- Andrew Ellis - <a href="https://kogpsy.github.io/neuroscicomplab" target="_blank">kogpsy.github.io/neuroscicomplab</a> --> <!-- </span> --> <!-- </div> --> <!-- slide separator (for xaringan) --> --- ## Hierarchische Modelle <!-- library(lmerTest) --> <!-- options(contrasts=c('contr.sum', 'contr.poly')) --> <!-- anova(lmer(y ~ x + (1|subj),data=myDat),ddf="Kenward-Roger") --> <!-- library(afex) --> <!-- mixed(y ~ x + (1|subj), type=3,method="KR",data=myDat) --> In dieser Sitzung behandeln wir messwiederholte Daten. Dabei stellen wir uns vor allem die Frage: Wie können wir Messwiederholung als lineares Modell formulieren? -- Wir fangen mit einem simplen Beispiel an, und untersuchen im zweiten Teil ein etwas komplexeres experimentelles Design. --- class: middle .pull-left-narrow[ .huge-blue-number[1] ] .pull-right-wide[ .larger[ Mittelwerte schätzen ] ] --- ## Wieder mal ein IQ Beispiel Wir haben dieses Mal 3 Personen, und von jeder Person 3 Messungen. Das Beispiel kommt dieses Mal aus dem Buch von Lee & Wagenmakers. .footnote[ Lee, Michael D., and Eric-Jan Wagenmakers. Bayesian Cognitive Modeling: A Practical Course. Cambridge: Cambridge University Press, 2014. https://doi.org/10.1017/CBO9781139087759. ] .panelset[ .panel[.panel-name[Create data] ```r IQwide <- tribble( ~A, ~B, ~C, 110, 105, 115, 105, 112, 108, 102, 113, 130 ) ``` ] .panel[.panel-name[Pivot longer] ```r IQdata <- IQwide %>% pivot_longer(everything(), names_to = "Person", values_to = "IQ") %>% mutate(Person = as_factor(Person)) %>% arrange(Person) ``` ] .panel[.panel-name[Data] ```r library(kableExtra) IQdata %>% kbl() %>% scroll_box(width = "500px", height = "200px") ``` <div style="border: 1px solid #ddd; padding: 0px; overflow-y: scroll; height:200px; overflow-x: scroll; width:500px; "><table> <thead> <tr> <th style="text-align:left;position: sticky; top:0; background-color: #FFFFFF;"> Person </th> <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> IQ </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> A </td> <td style="text-align:right;"> 110 </td> </tr> <tr> <td style="text-align:left;"> A </td> <td style="text-align:right;"> 105 </td> </tr> <tr> <td style="text-align:left;"> A </td> <td style="text-align:right;"> 102 </td> </tr> <tr> <td style="text-align:left;"> B </td> <td style="text-align:right;"> 105 </td> </tr> <tr> <td style="text-align:left;"> B </td> <td style="text-align:right;"> 112 </td> </tr> <tr> <td style="text-align:left;"> B </td> <td style="text-align:right;"> 113 </td> </tr> <tr> <td style="text-align:left;"> C </td> <td style="text-align:right;"> 115 </td> </tr> <tr> <td style="text-align:left;"> C </td> <td style="text-align:right;"> 108 </td> </tr> <tr> <td style="text-align:left;"> C </td> <td style="text-align:right;"> 130 </td> </tr> </tbody> </table></div> ] ] --- ## Punktschätzung ```r se <- function(x) sd(x)/sqrt(length(x)) IQdata %>% group_by(Person) %>% summarise(mean = mean(IQ), sd = sd(IQ), se = se(IQ)) ``` ``` ## # A tibble: 3 x 4 ## Person mean sd se ## <fct> <dbl> <dbl> <dbl> ## 1 A 106. 4.04 2.33 ## 2 B 110 4.36 2.52 ## 3 C 118. 11.2 6.49 ``` --- ```r IQdata %>% ggplot(aes(Person, IQ)) + geom_point() ``` <img src="04-bayesian-stats-slides_files/figure-html/unnamed-chunk-6-1.png" width="100%" /> --- .discussion[ Wie würden Sie diese Daten analysieren? - Die Daten sind messwiederholt. Für jede Person haben wir 3 Messungen. - Welche Frage könnte hier interessant sein? - Welche Methode(n) würden Sie hier anwenden? - Bei Person 3 gibt einen _Ausreisser_. Was würden Sie tun? ] <div class="countdown" id="timer_6076afb7" style="right:0;bottom:0;" data-warnwhen="0"> <code class="countdown-time"><span class="countdown-digits minutes">03</span><span class="countdown-digits colon">:</span><span class="countdown-digits seconds">00</span></code> </div> --- ## Parameterschätzung Das Ziel ist es nun, sowohl die 3 latenten IQs zu schätzen, als auch den Gruppenmittelwert Wir haben folgende Möglichkeiten: - __No pooling:__ Wir tun so, als ob die Personen alle nichts miteinander zu tun haben. - __Complete pooling:__ Wir tun so, als ob alle Messungen von derselben Person kommen. - __Partial pooling:__ Wir berücksichtigen, dass die Messungen von 3 verschiedenen Personen kommen. Wir glauben, dass Information, welche wir über eine Person haben, auch die Schätzung der anderen Personen beeinflussen sollte. --- ## No pooling ```r library(brms) get_prior(IQ ~ 0 + Person, data = IQdata) ``` ``` ## prior class coef group resp dpar nlpar bound source ## (flat) b default ## (flat) b PersonA (vectorized) ## (flat) b PersonB (vectorized) ## (flat) b PersonC (vectorized) ## student_t(3, 0, 7.4) sigma default ``` --- ## No pooling ```r m_no_pool <- brm(IQ ~ 0 + Person, data = IQdata, file = "models/m_no_pool") ``` --- ## Complete pooling ```r get_prior(IQ ~ 1, data = IQdata) ``` ``` ## prior class coef group resp dpar nlpar bound source ## student_t(3, 110, 7.4) Intercept default ## student_t(3, 0, 7.4) sigma default ``` --- ## Complete pooling ```r m_comp_pool <- brm(IQ ~ 1, data = IQdata, file = "models/m_comp_pool") ``` --- ## Partial pooling ```r get_prior(IQ ~ 1 + (1 | Person), data = IQdata) ``` ``` ## prior class coef group resp dpar nlpar bound ## student_t(3, 110, 7.4) Intercept ## student_t(3, 0, 7.4) sd ## student_t(3, 0, 7.4) sd Person ## student_t(3, 0, 7.4) sd Intercept Person ## student_t(3, 0, 7.4) sigma ## source ## default ## default ## (vectorized) ## (vectorized) ## default ``` --- ## Partial pooling <img src="04-bayesian-stats-slides_files/figure-html/unnamed-chunk-13-1.png" width="100%" /> --- ## Partial pooling ```r m_part_pool <- brm(IQ ~ 1 + (1 | Person), data = IQdata, file = "models/m_part_pool") ``` --- .panelset[ .panel[.panel-name[No Pooling] ```r m_no_pool ``` ``` ## Family: gaussian ## Links: mu = identity; sigma = identity ## Formula: IQ ~ 0 + Person ## Data: IQdata (Number of observations: 9) ## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; ## total post-warmup samples = 4000 ## ## Population-Level Effects: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## PersonA 105.67 4.90 95.63 115.37 1.00 2708 1960 ## PersonB 110.03 4.96 100.05 119.93 1.00 3531 2408 ## PersonC 117.66 4.92 107.87 127.13 1.00 3645 2285 ## ## Family Specific Parameters: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sigma 8.17 2.55 4.76 14.42 1.00 2152 2189 ## ## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1). ``` ] .panel[.panel-name[Complete Pooling] ```r m_comp_pool ``` ``` ## Family: gaussian ## Links: mu = identity; sigma = identity ## Formula: IQ ~ 1 ## Data: IQdata (Number of observations: 9) ## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; ## total post-warmup samples = 4000 ## ## Population-Level Effects: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## Intercept 110.93 2.71 105.68 116.47 1.00 2390 2126 ## ## Family Specific Parameters: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sigma 8.78 2.22 5.64 14.04 1.00 2397 2354 ## ## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1). ``` ] .panel[.panel-name[Partial Pooling] ```r m_part_pool ``` ``` ## Family: gaussian ## Links: mu = identity; sigma = identity ## Formula: IQ ~ 1 + (1 | Person) ## Data: IQdata (Number of observations: 9) ## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; ## total post-warmup samples = 4000 ## ## Group-Level Effects: ## ~Person (Number of levels: 3) ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sd(Intercept) 5.00 3.90 0.18 15.33 1.00 1187 1793 ## ## Population-Level Effects: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## Intercept 111.03 3.52 104.03 118.26 1.00 1537 1491 ## ## Family Specific Parameters: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sigma 8.09 2.29 4.88 13.45 1.00 1733 2515 ## ## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1). ``` ] ] --- ## Partial Pooling ```r f <- fixef(m_part_pool, summary = FALSE) r <- ranef(m_part_pool, summary = FALSE) ``` ```r library(tidybayes) get_variables(m_part_pool) ``` ``` ## [1] "b_Intercept" "sd_Person__Intercept" "sigma" ## [4] "r_Person[A,Intercept]" "r_Person[B,Intercept]" "r_Person[C,Intercept]" ## [7] "lp__" "accept_stat__" "stepsize__" ## [10] "treedepth__" "n_leapfrog__" "divergent__" ## [13] "energy__" ``` ```r person_effects <- m_part_pool %>% spread_draws(b_Intercept, r_Person[Person, Intercept]) %>% # add the grand mean to the person-specific deviations mutate(mu = b_Intercept + r_Person) ``` --- ```r person_effects %>% median_qi(mu) ``` ``` ## # A tibble: 3 x 8 ## # Groups: Person [3] ## Person Intercept mu .lower .upper .width .point .interval ## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr> ## 1 A Intercept 109. 100. 116. 0.95 median qi ## 2 B Intercept 111. 104. 118. 0.95 median qi ## 3 C Intercept 114. 106. 122. 0.95 median qi ``` ```r fixef(m_no_pool) ``` ``` ## Estimate Est.Error Q2.5 Q97.5 ## PersonA 105.6677 4.903999 95.63449 115.3670 ## PersonB 110.0306 4.962230 100.04958 119.9328 ## PersonC 117.6587 4.917794 107.87280 127.1344 ``` --- .panelset[ .panel[.panel-name[Plot Code] ```r person_effects %>% ggplot(aes(mu, Person)) + stat_halfeye(fill = "#A2BF8A") + geom_vline(xintercept = fixef(m_part_pool)[1, 1], color = "#839496", size = 1) + geom_vline(xintercept = fixef(m_part_pool)[1, 3:4], color = "#839496", linetype = 2) + labs(x = expression("Personen-spezifische Mittelwerte"), y = "Personen") + theme(panel.grid = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_text(hjust = 0)) ``` ] .panel[.panel-name[Plot] <img src="04-bayesian-stats-slides_files/figure-html/unnamed-chunk-23-1.png" width="100%" /> ] ] --- ## Vergleich mit No Pooling: Shrinkage .panelset[ .panel[.panel-name[Plot Code] ```r col <- viridis::viridis(3, begin = 0.2, end = 0.8) person_effects %>% ggplot(aes(mu, Person, fill = Person)) + stat_halfeye(alpha = 0.6) + geom_vline(xintercept = fixef(m_part_pool)[1, 1], color = "#ECCC87", size = 1) + geom_vline(xintercept = fixef(m_no_pool)[1, 1], color = col[1], size = 1) + geom_vline(xintercept = fixef(m_no_pool)[2, 1], color = col[2], size = 1) + geom_vline(xintercept = fixef(m_no_pool)[3, 1], color = col[3], size = 1) + scale_fill_viridis_d(begin = 0.2, end = 0.8) + labs(x = expression("Personen-spezifische Mittelwerte"), y = "Personen") + theme(panel.grid = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_text(hjust = 0)) ``` ] .panel[.panel-name[Plot] <img src="04-bayesian-stats-slides_files/figure-html/unnamed-chunk-25-1.png" width="100%" /> ] ] --- class: middle .pull-left-narrow[ .huge-blue-number[2] ] .pull-right-wide[ .larger[ Mixed Models ] ] --- .panelset[ .panel[.panel-name[Create dataset] ```r library(tidyverse) intervention <- rep(c('treat', 'control'), each = 5) pre <- c(20, 10, 60, 20, 10, 50, 10, 40, 20, 10) post <- c(70, 50, 90, 60, 50, 20, 10, 30, 50, 10) ``` ```r dwide <- tibble(id = factor(1:10), intervention, pre, post) %>% mutate(diff = post - pre, id = as_factor(id), intervention = factor(intervention, levels = c("control", "treat"))) ``` ] .panel[.panel-name[Dataframe] ```r dwide %>% paged_table() ``` <div data-pagedtable="false"> <script data-pagedtable-source type="application/json"> {"columns":[{"label":["id"],"name":[1],"type":["fct"],"align":["left"]},{"label":["intervention"],"name":[2],"type":["fct"],"align":["left"]},{"label":["pre"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["post"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["diff"],"name":[5],"type":["dbl"],"align":["right"]}],"data":[{"1":"1","2":"treat","3":"20","4":"70","5":"50"},{"1":"2","2":"treat","3":"10","4":"50","5":"40"},{"1":"3","2":"treat","3":"60","4":"90","5":"30"},{"1":"4","2":"treat","3":"20","4":"60","5":"40"},{"1":"5","2":"treat","3":"10","4":"50","5":"40"},{"1":"6","2":"control","3":"50","4":"20","5":"-30"},{"1":"7","2":"control","3":"10","4":"10","5":"0"},{"1":"8","2":"control","3":"40","4":"30","5":"-10"},{"1":"9","2":"control","3":"20","4":"50","5":"30"},{"1":"10","2":"control","3":"10","4":"10","5":"0"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}} </script> </div> ]] --- ```r d <- dwide %>% select(-diff) %>% pivot_longer(cols = pre:post, names_to = "time", values_to = "score") %>% mutate(time = as_factor(time)) d %>% paged_table() ``` <div data-pagedtable="false"> <script data-pagedtable-source type="application/json"> {"columns":[{"label":["id"],"name":[1],"type":["fct"],"align":["left"]},{"label":["intervention"],"name":[2],"type":["fct"],"align":["left"]},{"label":["time"],"name":[3],"type":["fct"],"align":["left"]},{"label":["score"],"name":[4],"type":["dbl"],"align":["right"]}],"data":[{"1":"1","2":"treat","3":"pre","4":"20"},{"1":"1","2":"treat","3":"post","4":"70"},{"1":"2","2":"treat","3":"pre","4":"10"},{"1":"2","2":"treat","3":"post","4":"50"},{"1":"3","2":"treat","3":"pre","4":"60"},{"1":"3","2":"treat","3":"post","4":"90"},{"1":"4","2":"treat","3":"pre","4":"20"},{"1":"4","2":"treat","3":"post","4":"60"},{"1":"5","2":"treat","3":"pre","4":"10"},{"1":"5","2":"treat","3":"post","4":"50"},{"1":"6","2":"control","3":"pre","4":"50"},{"1":"6","2":"control","3":"post","4":"20"},{"1":"7","2":"control","3":"pre","4":"10"},{"1":"7","2":"control","3":"post","4":"10"},{"1":"8","2":"control","3":"pre","4":"40"},{"1":"8","2":"control","3":"post","4":"30"},{"1":"9","2":"control","3":"pre","4":"20"},{"1":"9","2":"control","3":"post","4":"50"},{"1":"10","2":"control","3":"pre","4":"10"},{"1":"10","2":"control","3":"post","4":"10"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}} </script> </div> --- Wir haben 10 Personen .panelset[ .panel[.panel-name[Subjects] ```r d %>% summarize(id = n_distinct(id)) ``` ``` ## # A tibble: 1 x 1 ## id ## <int> ## 1 10 ``` ] .panel[.panel-name[Trials per subject] mit 2 Messungen pro Person. ```r d %>% group_by(id, intervention) %>% count() %>% rmarkdown::paged_table() ``` <div data-pagedtable="false"> <script data-pagedtable-source type="application/json"> {"columns":[{"label":["id"],"name":[1],"type":["fct"],"align":["left"]},{"label":["intervention"],"name":[2],"type":["fct"],"align":["left"]},{"label":["n"],"name":[3],"type":["int"],"align":["right"]}],"data":[{"1":"1","2":"treat","3":"2"},{"1":"2","2":"treat","3":"2"},{"1":"3","2":"treat","3":"2"},{"1":"4","2":"treat","3":"2"},{"1":"5","2":"treat","3":"2"},{"1":"6","2":"control","3":"2"},{"1":"7","2":"control","3":"2"},{"1":"8","2":"control","3":"2"},{"1":"9","2":"control","3":"2"},{"1":"10","2":"control","3":"2"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}} </script> </div> ] ] --- .discussion[ Wie würden Sie diese Daten analysieren? - Die Daten sind messwiederholt. Jede Person gibt zu 2 Zeitpunkten eine Antwort. - Wir haben einen `between-subjects` Faktor: `intervention` - Welche Frage könnte hier interessant sein? - Welche Methode(n) würden Sie hier anwenden? ] <div class="countdown" id="timer_6076aee6" style="right:0;bottom:0;" data-warnwhen="0"> <code class="countdown-time"><span class="countdown-digits minutes">03</span><span class="countdown-digits colon">:</span><span class="countdown-digits seconds">00</span></code> </div> --- ## Daten zusammenfassen .panelset[ .panel[.panel-name[Standard error function] ```r se <- function(x) sd(x)/sqrt(length(x)) ``` ] .panel[.panel-name[Summarize (long)] ```r d %>% group_by(intervention, time) %>% summarise(mean = mean(score), sd = sd(score), se = se(score)) ``` ``` ## # A tibble: 4 x 5 ## # Groups: intervention [2] ## intervention time mean sd se ## <fct> <fct> <dbl> <dbl> <dbl> ## 1 control pre 26 18.2 8.12 ## 2 control post 24 16.7 7.48 ## 3 treat pre 24 20.7 9.27 ## 4 treat post 64 16.7 7.48 ``` ] .panel[.panel-name[Summarize (wide)] ```r dwide %>% group_by(intervention) %>% summarise(mean = mean(diff), sd = sd(diff), se = se(diff)) ``` ``` ## # A tibble: 2 x 4 ## intervention mean sd se ## <fct> <dbl> <dbl> <dbl> ## 1 control -2 21.7 9.70 ## 2 treat 40 7.07 3.16 ``` ] ] --- .panelset[ .panel[.panel-name[Plot Code] ```r d %>% ggplot(aes(time, score, color = intervention)) + geom_line(aes(group = id), linetype = 1, size = 1) + geom_point(size = 4) + scale_color_viridis_d(end = 0.8) + theme_bw() ``` ] .panel[.panel-name[Plot] <img src="04-bayesian-stats-slides_files/figure-html/unnamed-chunk-37-1.png" width="100%" /> ] ] --- ## t-Test ```r t.test(diff ~ intervention, data = dwide, var.equal = TRUE) ``` ``` ## ## Two Sample t-test ## ## data: diff by intervention ## t = -4.1184, df = 8, p-value = 0.003351 ## alternative hypothesis: true difference in means between group control and group treat is not equal to 0 ## 95 percent confidence interval: ## -65.51672 -18.48328 ## sample estimates: ## mean in group control mean in group treat ## -2 40 ``` --- ## Mixed Model .panelset[ .panel[.panel-name[Model Code] ```r library(lme4) lme_model <- lmer(score ~ intervention * time + (1|id), data = d) ``` ] .panel[.panel-name[Model summary] ```r lme_model %>% sjPlot::tab_model() ``` <table style="border-collapse:collapse; border:none;"> <tr> <th style="border-top: double; text-align:center; font-style:normal; font-weight:bold; padding:0.2cm; text-align:left; ">&nbsp;</th> <th colspan="3" style="border-top: double; text-align:center; font-style:normal; font-weight:bold; padding:0.2cm; ">score</th> </tr> <tr> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; text-align:left; ">Predictors</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Estimates</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">CI</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">p</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">(Intercept)</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">26.00</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">10.08&nbsp;&ndash;&nbsp;41.92</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>0.001</strong></td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">intervention [treat]</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">&#45;2.00</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">&#45;24.52&nbsp;&ndash;&nbsp;20.52</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.862</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">time [post]</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">&#45;2.00</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">&#45;16.13&nbsp;&ndash;&nbsp;12.13</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.782</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">intervention [treat] *<br>time [post]</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">42.00</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">22.01&nbsp;&ndash;&nbsp;61.99</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>&lt;0.001</td> </tr> <tr> <td colspan="4" style="font-weight:bold; text-align:left; padding-top:.8em;">Random Effects</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">&sigma;²</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">130.00</td> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">&tau;₀₀ _id</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">200.00</td> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">ICC</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">0.61</td> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">N _id</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">10</td> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm; border-top:1px solid;">Observations</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left; border-top:1px solid;" colspan="3">20</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">Marginal R² / Conditional R²</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">0.481 / 0.796</td> </tr> </table> ]] --- ```r library(afex) mixed(score ~ intervention * time + (1|id), type = 3, method = "KR", data = d) ``` ``` ## Fitting one lmer() model. [DONE] ## Calculating p-values. [DONE] ``` ``` ## Mixed Model Anova Table (Type 3 tests, KR-method) ## ## Model: score ~ intervention * time + (1 | id) ## Data: d ## Effect df F p.value ## 1 intervention 1, 8 3.41 .102 ## 2 time 1, 8 13.88 ** .006 ## 3 intervention:time 1, 8 16.96 ** .003 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1 ``` --- ## Was geschieht hier genau? ```r mm <- model.matrix(~ intervention * time, data = d) head(mm) ``` ``` ## (Intercept) interventiontreat timepost interventiontreat:timepost ## 1 1 1 0 0 ## 2 1 1 1 1 ## 3 1 1 0 0 ## 4 1 1 1 1 ## 5 1 1 0 0 ## 6 1 1 1 1 ``` --- ## Erklärung Eigentlich ganz einfach: - Wir sagen den Mittelwert der abhängingen Variable `score` mit einem linearen Modell vorher. - Im linearen Modell zerlegen wir den Mittelwert in + fixed effects: `intervention * time` + random effects (eine Abweichung für jede Person): `(1 | id)` --- ## Mit `brms` .panelset[ .panel[.panel-name[get_prior] ```r library(brms) priors <- get_prior(score ~ intervention*time, data = d) ``` ] .panel[.panel-name[Priors] ```r priors ``` ``` ## prior class coef group resp dpar ## (flat) b ## (flat) b interventiontreat ## (flat) b interventiontreat:timepost ## (flat) b timepost ## student_t(3, 25, 22.2) Intercept ## student_t(3, 0, 22.2) sigma ## nlpar bound source ## default ## (vectorized) ## (vectorized) ## (vectorized) ## default ## default ``` ] ] --- ```r m2 <- brm(score ~ intervention*time + (1 | id), data = d, file = "models/04-treat-time") ``` --- ```r summary(m2) ``` ``` ## Family: gaussian ## Links: mu = identity; sigma = identity ## Formula: score ~ treat * time + (1 | id) ## Data: d (Number of observations: 20) ## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; ## total post-warmup samples = 4000 ## ## Group-Level Effects: ## ~id (Number of levels: 10) ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sd(Intercept) 13.91 6.42 1.61 28.05 1.01 541 723 ## ## Population-Level Effects: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## Intercept 25.29 9.25 6.15 42.92 1.00 1465 1161 ## treattreat -2.18 13.19 -26.89 24.74 1.00 1581 1587 ## timepost -1.99 9.07 -21.02 16.25 1.00 2575 2477 ## treattreat:timepost 42.10 12.97 16.21 68.76 1.00 2479 2398 ## ## Family Specific Parameters: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sigma 13.62 3.86 8.07 22.97 1.00 686 1507 ## ## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1). ``` --- ```r m2 %>% mcmc_plot("b_") ``` <img src="04-bayesian-stats-slides_files/figure-html/unnamed-chunk-47-1.png" width="100%" />