pigs2

bayes

qm2

qm2-pruefung

regression

exam-22

Published

January 11, 2023

Aufgabe

library(tidyverse)  # Datenjudo
library(rstanarm)  # Bayes-Inferenz
library(easystats)  # Komfort

Laden Sie den Datensatz toothgrowth, z.B. von dieser Quelle. In dem Datensatz sind die Daten eines einfaches Experiments berichtet.

In dem Experiment wird der (kausale) Effekt verschiedener Darreichungsformen und Konzentrationen von Vitamin C auf das Zahnwachstum von Meerschweinchen untersucht.

Forschungsfrage:

Hat die Dosis (dose) einen (kausalen) Effekt auf die AV, Zahnlänge (len)?

Wir gehen mal einfach davon aus, dass der Faktor experimentell (also randomisiert und auf Störeffekte hin kontrollliert) veraeicht wurde. Sonst wäre eine Kausalinterpretation nicht (ohne Weiteres) möglich.

Aufgabe: Berechnen Sie die Breite eines 95%-HDI für den Effekt!

Hinweise

Lösung

Schritt: Modell berechnen

lm2 <- stan_glm(len ~ dose, data = d,
                seed = 42,
                refresh = 0)

Zur Erinnerung: Bei der Regressionsformel gilt immer av ~ uv.

Schritt: Posteriori-Verteilung betrachten

Mit parameters() kriegt man einen guten Überblick über die Modellparameter:

parameters(lm2)

Parameter   | Median |        95% CI |   pd |  Rhat |     ESS |                   Prior
---------------------------------------------------------------------------------------
(Intercept) |   7.40 | [4.84,  9.96] | 100% | 0.999 | 3881.00 | Normal (18.81 +- 19.12)
dose        |   9.76 | [7.81, 11.76] | 100% | 0.999 | 4144.00 |  Normal (0.00 +- 30.41)

Das Modell zeigt einen positiven Effekt für dose:

Pro Einheit von dose steigt die Zahnlänge (len) um ca. 8-12 mm im Schnitt (laut unserem Modell).

Null ist nicht im Intervall enthalten; die Nullhypothese ist demnach auszuschließen (falls das jemanden interessiert). Man sagt, man verwirft die Nullhypothese (oder weist sie zurück).

Das können wir auch plotten:

plot(parameters(lm2))

Man kann sich auch ein HDI direkt ausgeben, ohne die sonstigen Informationen, die parameters() ausgibt:

hdi(lm2, ci = .95)

Highest Density Interval

Parameter   |       95% HDI
---------------------------
(Intercept) | [4.96, 10.06]
dose        | [7.77, 11.70]

Wie wir sehen, wird im Standard ein 95%-Intervall berichtet, wie in der Aufgabenstellung verlangt.

Wieder sehen wir, die Null ist nicht im Intervall enthalten. Null ist also kein plausibler Wert für den gesuchten Effekt (laut unserem Modell).

Schauen wir uns mal zum Vergleich die Stichproben-Daten an:

d %>% 
  ggplot(aes(y = len, x = dose)) +
  geom_point() +
  geom_smooth(method = "lm")

Man sieht deutlich einen positiven Effekt: Die Regressionsgerade steigt.

Die Breite des Schätzintervalls für den Effekt beträgt also:

sol <- 11.70 - 7.77
sol

[1] 3.93

Categories:

bayes
~
regression
exam-22

--- exname: pigs2 extype: num extol: 1 exsolution: r sol categories: - bayes - qm2 - qm2-pruefung - regression - exam-22 date: '2023-01-11' slug: pigs2 title: pigs2 --- ```{r libs, include = FALSE} library(tidyverse) ``` ```{r global-knitr-options, include=FALSE} knitr::opts_chunk$set(fig.pos = 'H', fig.asp = 0.618, fig.width = 4, fig.cap = "", fig.path = "", echo = TRUE, message = FALSE, fig.show = "hold") ``` # Aufgabe ```{r message=FALSE} library(tidyverse) # Datenjudo library(rstanarm) # Bayes-Inferenz library(easystats) # Komfort ``` Laden Sie den Datensatz `toothgrowth`, z.B. von [dieser Quelle](https://vincentarelbundock.github.io/Rdatasets/csv/datasets/ToothGrowth.csv). In dem Datensatz sind die Daten eines einfaches Experiments berichtet. ```{r echo = FALSE} d_path <- "https://vincentarelbundock.github.io/Rdatasets/csv/datasets/ToothGrowth.csv" d <- read_csv(d_path) ``` In dem Experiment wird der (kausale) Effekt verschiedener Darreichungsformen und Konzentrationen von Vitamin C auf das Zahnwachstum von Meerschweinchen untersucht. Forschungsfrage: Hat die Dosis (`dose`) einen (kausalen) Effekt auf die AV, Zahnlänge (`len`)? Wir gehen mal einfach davon aus, dass der Faktor experimentell (also randomisiert und auf Störeffekte hin kontrollliert) veraeicht wurde. Sonst wäre eine Kausalinterpretation nicht (ohne Weiteres) möglich. Aufgabe: *Berechnen Sie die Breite eines 95%-HDI für den Effekt!* [Hinweise](https://start-bayes.netlify.app/pruefung#bearbeitungshinweise) # Lösung 1. Schritt: Modell berechnen ```{r} lm2 <- stan_glm(len ~ dose, data = d, seed = 42, refresh = 0) ``` Zur Erinnerung: Bei der Regressionsformel gilt immer `av ~ uv`. 2. Schritt: Posteriori-Verteilung betrachten Mit `parameters()` kriegt man einen guten Überblick über die Modellparameter: ```{r} parameters(lm2) ``` Das Modell zeigt einen positiven Effekt für `dose`: Pro Einheit von `dose` steigt die Zahnlänge (`len`) um ca. 8-12 mm im Schnitt (laut unserem Modell). *Null* ist nicht im Intervall enthalten; die Nullhypothese ist demnach *auszuschließen* (falls das jemanden interessiert). Man sagt, man *verwirft* die Nullhypothese (oder weist sie zurück). Das können wir auch plotten: ```{r} plot(parameters(lm2)) ``` Man kann sich auch ein HDI direkt ausgeben, ohne die sonstigen Informationen, die `parameters()` ausgibt: ```{r} hdi(lm2, ci = .95) ``` Wie wir sehen, wird im Standard ein 95%-Intervall berichtet, wie in der Aufgabenstellung verlangt. Wieder sehen wir, die Null ist *nicht* im Intervall enthalten. Null ist also *kein* plausibler Wert für den gesuchten Effekt (laut unserem Modell). Schauen wir uns mal zum Vergleich die Stichproben-Daten an: ```{r} d %>% ggplot(aes(y = len, x = dose)) + geom_point() + geom_smooth(method = "lm") ``` Man sieht deutlich einen positiven Effekt: Die Regressionsgerade steigt. Die Breite des Schätzintervalls für den Effekt beträgt also: ```{r} sol <- 11.70 - 7.77 sol ``` --- Categories: - bayes - ~ - regression - exam-22