Permalink
Browse files

more work on vignettes

  • Loading branch information...
1 parent 8632b52 commit 3eef5f6d13b24a330cb96940b7606cdc12867104 @paul-buerkner committed Dec 29, 2016
@@ -7,10 +7,10 @@ opts_chunk$set(
comment = NA,
message = FALSE,
warning = FALSE,
- eval = FALSE, # params$EVAL,
+ eval = params$EVAL,
dev = "png",
dpi = 150,
- fig.asp = 0.618,
+ fig.asp = 0.8,
fig.width = 5,
out.width = "60%",
fig.align = "center"
@@ -22,11 +22,11 @@ opts_chunk$set(
# dat1 <- data.frame(group, symptom_post)
# head(dat1)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit1 <- brm(bf(symptom_post ~ group, sigma ~ group),
# data = dat1, family = gaussian())
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# summary(fit1)
# plot(fit1)
# plot(marginal_effects(fit1), points = TRUE)
@@ -45,7 +45,7 @@ opts_chunk$set(
# zinb <- read.csv("http://www.ats.ucla.edu/stat/data/fish.csv")
# head(zinb)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit_zinb1 <- brm(count ~ persons + child + camper,
# data = zinb, family = zero_inflated_poisson())
@@ -54,7 +54,7 @@ opts_chunk$set(
# plot(fit_zinb1)
# plot(marginal_effects(fit_zinb1), ask = FALSE)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit_zinb2 <- brm(bf(count ~ persons + child + camper, zi ~ child),
# data = zinb, family = zero_inflated_poisson())
@@ -68,7 +68,7 @@ opts_chunk$set(
# dat_smooth <- gamSim(eg = 6, n = 200, scale = 2)
# head(dat_smooth)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit_smooth1 <- brm(bf(y ~ s(x1) + s(x2) + (1|fac), sigma ~ s(x0) + (1|fac)),
# data = dat_smooth, family = gaussian(),
# chains = 2, control = list(adapt_delta = 0.95))
@@ -6,7 +6,8 @@ output:
rmarkdown::html_vignette:
toc: yes
params:
- EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true")
+ # EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true")
+ EVAL: !r FALSE
---
<!--
@@ -46,14 +47,14 @@ head(dat1)
The following model estimates the effect of `group` on both the mean and the residual standard deviation of the normal response distribution.
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
fit1 <- brm(bf(symptom_post ~ group, sigma ~ group),
data = dat1, family = gaussian())
```
Useful summary statistics and plots can be obtained via
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
summary(fit1)
plot(fit1)
plot(marginal_effects(fit1), points = TRUE)
@@ -89,7 +90,7 @@ head(zinb)
As predictors we choose the number of people per group, the number of children, as well as whether the group consists of campers. Many groups may not catch any fish and so we directly fit a zero-inflated Poisson model. For now, we assume a constant zero-inflation probability across observations.
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
fit_zinb1 <- brm(count ~ persons + child + camper,
data = zinb, family = zero_inflated_poisson())
```
@@ -106,7 +107,7 @@ According to the parameter estimates, larger groups catch more fish, campers cat
Now, we try to additionally predict the zero-inflation probability by the number of children. But why should we predict the zero-inflation probability at all? One line of reasong might go as follows. We assume two processes behind our response, the first being the basic Poisson process and the second being the zero-inflation process. For the smoker example, it is easy to interprete the zero-inflation process as the decision between smoking and non-smoking in general, which is definitely worth investigating using regression techniques. For other applications of zero-inflated models, for instance in the fishing example, it may not necessarily be that easy to interprete the zero-inflation process (or maybe it is, I have no idea of fishing), but the statistical idea remains the same: Each of the two processes comes with its own implied distribution and the zero-inflated Poisson distribution is a mixture of the two. From this perspective, predicting both parts of the model is natural and reasonable to fully understand the data.
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
fit_zinb2 <- brm(bf(count ~ persons + child + camper, zi ~ child),
data = zinb, family = zero_inflated_poisson())
```
@@ -137,7 +138,7 @@ head(dat_smooth)
The data contains the predictors `x0` to `x3` as well as the grouping factor `fac` indicating the nested structure of the data. We predict the response variable `y` using smooth terms of `x1` and `x2` and a varying intercept of `fac`. In addition, we assume the residual standard deviation `sigma` to vary by a smoothing term of `x0` and a varying intercept of `fac`.
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
fit_smooth1 <- brm(bf(y ~ s(x1) + s(x2) + (1|fac), sigma ~ s(x0) + (1|fac)),
data = dat_smooth, family = gaussian(),
chains = 2, control = list(adapt_delta = 0.95))
@@ -12,7 +12,7 @@
<meta name="author" content="Paul Bürkner" />
-<meta name="date" content="2016-12-28" />
+<meta name="date" content="2016-12-29" />
<title>Estimating Distributional Models with brms</title>
@@ -70,7 +70,7 @@
<h1 class="title toc-ignore">Estimating Distributional Models with brms</h1>
<h4 class="author"><em>Paul Bürkner</em></h4>
-<h4 class="date"><em>2016-12-28</em></h4>
+<h4 class="date"><em>2016-12-29</em></h4>
<div id="TOC">
@@ -12,7 +12,7 @@
<meta name="author" content="Paul Bürkner" />
-<meta name="date" content="2016-12-28" />
+<meta name="date" content="2016-12-29" />
<title>Parameterization of Response Distributions in brms</title>
@@ -32,7 +32,7 @@
<h1 class="title toc-ignore">Parameterization of Response Distributions in brms</h1>
<h4 class="author"><em>Paul Bürkner</em></h4>
-<h4 class="date"><em>2016-12-28</em></h4>
+<h4 class="date"><em>2016-12-29</em></h4>
@@ -7,10 +7,10 @@ opts_chunk$set(
comment = NA,
message = FALSE,
warning = FALSE,
- eval = FALSE, # params$EVAL,
+ eval = params$EVAL,
dev = "png",
dpi = 150,
- fig.asp = 0.618,
+ fig.asp = 0.8,
fig.width = 5,
out.width = "60%",
fig.align = "center"
@@ -24,7 +24,7 @@ opts_chunk$set(
# ls <- mean_ls[income] + rnorm(100, sd = 7)
# dat <- data.frame(income, ls)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# library(brms)
# fit1 <- brm(ls ~ monotonic(income), data = dat)
@@ -33,14 +33,14 @@ opts_chunk$set(
# plot(fit1, pars = "simplex")
# plot(marginal_effects(fit1))
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# dat$income_num <- as.numeric(dat$income)
# fit2 <- brm(ls ~ income_num, data = dat)
## ------------------------------------------------------------------------
# summary(fit2)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# contrasts(dat$income) <- contr.treatment(4)
# fit3 <- brm(ls ~ income, data = dat)
@@ -50,7 +50,7 @@ opts_chunk$set(
## ------------------------------------------------------------------------
# LOO(fit1, fit2, fit3)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# prior4 <- prior(dirichlet(c(2, 1, 1)), class = "simplex", coef = "income")
# fit4 <- brm(ls ~ monotonic(income), data = dat,
# prior = prior4, sample_prior = TRUE)
@@ -66,7 +66,7 @@ opts_chunk$set(
# var_city <- rnorm(10, sd = 10)
# dat$ls <- dat$ls + var_city[dat$city]
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit5 <- brm(ls ~ mo(income) + (1 | city) + (mo(income) | city), data = dat)
## ------------------------------------------------------------------------
@@ -6,7 +6,8 @@ output:
rmarkdown::html_vignette:
toc: yes
params:
- EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true")
+ # EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true")
+ EVAL: !r FALSE
---
<!--
@@ -43,7 +44,7 @@ dat <- data.frame(income, ls)
We now proceed with analyzing the data modeling `income` as a monotonic effect.
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
library(brms)
fit1 <- brm(ls ~ monotonic(income), data = dat)
```
@@ -60,7 +61,7 @@ The distributions of the simplex parameter of `income`, as shown in the `plot` m
Now, let's compare of monotonic model with two common alternative models. (a) Assume `income` to be continuous:
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
dat$income_num <- as.numeric(dat$income)
fit2 <- brm(ls ~ income_num, data = dat)
```
@@ -72,7 +73,7 @@ summary(fit2)
or (b) Assume `income` to be an unordered factor:
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
contrasts(dat$income) <- contr.treatment(4)
fit3 <- brm(ls ~ income, data = dat)
```
@@ -93,7 +94,7 @@ The monotonic model fits better than the continuous model, which is not surprisi
In the previous monotonic model, we have implicitly assumed that all differences between adjacent categories were a-priori the same, or formulated correctly, had the same prior distribution. In the following, we want to show how to change this assumption. The canonical prior distribution of a simplex parameter is the Dirchlet distribution, a multivariate generalization of the beta distribution. It is non-zero for all valid simplexes (i.e., $\zeta_i \in [0,1]$ and $\sum_{i = 1}^K \zeta_i = 1$) and zero otherwise. The Dirichlet prior has a single parameter $\alpha$ of the same length as $\zeta$. The higher $\alpha_i$ the higher the a-priori probability of higher values of $\zeta_i$. Suppose that, before looking at the data, we expected that the same amount of additional money matters more for people who generally have less money. This translates into a higher a-priori values of $\zeta_1$ (difference between 'below_20' and '20_to_40') and hence into higher values of $\alpha_1$. We choose $\alpha_1 = 2$ and $\alpha_2 = \alpha_3 = 1$, the latter being the default value of $\alpha$. To fit the model we write:
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
prior4 <- prior(dirichlet(c(2, 1, 1)), class = "simplex", coef = "income")
fit4 <- brm(ls ~ monotonic(income), data = dat,
prior = prior4, sample_prior = TRUE)
@@ -123,7 +124,7 @@ dat$ls <- dat$ls + var_city[dat$city]
With the following code, we fit a multilevel model assuming the intercept and the effect of `income` to vary by city:
-```{r, results='hide', message=FALSE, warning = FALSE}
+```{r, results='hide'}
fit5 <- brm(ls ~ mo(income) + (1 | city) + (mo(income) | city), data = dat)
```
For technical reasons, monotonic group-level effects have to be specified in separate terms in the model formula. Further, we have used the abbrevation `mo` for `monotonic`, which helps in shortening the formula. The summary output
@@ -12,7 +12,7 @@
<meta name="author" content="Paul Bürkner" />
-<meta name="date" content="2016-12-28" />
+<meta name="date" content="2016-12-29" />
<title>Estimating Monotonic Effects with brms</title>
@@ -70,7 +70,7 @@
<h1 class="title toc-ignore">Estimating Monotonic Effects with brms</h1>
<h4 class="author"><em>Paul Bürkner</em></h4>
-<h4 class="date"><em>2016-12-28</em></h4>
+<h4 class="date"><em>2016-12-29</em></h4>
<div id="TOC">
@@ -7,10 +7,10 @@ opts_chunk$set(
comment = NA,
message = FALSE,
warning = FALSE,
- eval = FALSE, # params$EVAL,
+ eval = params$EVAL,
dev = "png",
dpi = 150,
- fig.asp = 0.618,
+ fig.asp = 0.8,
fig.width = 5,
out.width = "60%",
fig.align = "center"
@@ -22,7 +22,7 @@ opts_chunk$set(
# y <- rnorm(100, mean = b[1] * exp(b[2] * x))
# dat1 <- data.frame(x, y)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# library(brms)
# prior1 <- c(prior(normal(1, 2), nlpar = "b1"),
# prior(normal(0, 2), nlpar = "b2"))
@@ -34,7 +34,7 @@ opts_chunk$set(
# plot(fit1)
# plot(marginal_effects(fit1), points = TRUE)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit2 <- brm(y ~ x, data = dat1)
## ------------------------------------------------------------------------
@@ -52,7 +52,7 @@ opts_chunk$set(
# loss <- read.csv(url)
# head(loss)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit_loss <- brm(bf(cum ~ ult * (1 - exp(-(dev/theta)^omega)),
# ult ~ 1 + (1|AY), omega ~ 1, theta ~ 1,
# nl = TRUE),
@@ -81,14 +81,14 @@ opts_chunk$set(
# answer <- ifelse(runif(300, 0, 1) < p, 1, 0)
# dat_ir <- data.frame(ability, answer)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit_ir1 <- brm(answer ~ ability, data = dat_ir, family = bernoulli())
## ------------------------------------------------------------------------
# summary(fit_ir1)
# plot(marginal_effects(fit_ir1), points = TRUE)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit_ir2 <- brm(bf(answer ~ 0.33 + 0.67 * inv_logit(eta),
# eta ~ ability, nl = TRUE),
# data = dat_ir, family = bernoulli("identity"),
@@ -101,7 +101,7 @@ opts_chunk$set(
## ------------------------------------------------------------------------
# LOO(fit_ir1, fit_ir2)
-## ---- results='hide', message=FALSE, warning = FALSE---------------------
+## ---- results='hide'-----------------------------------------------------
# fit_ir3 <- brm(bf(answer ~ guess + (1 - guess) * inv_logit(eta),
# eta ~ 0 + ability, guess ~ 1, nl = TRUE),
# data = dat_ir, family = bernoulli("identity"),
Oops, something went wrong.

0 comments on commit 3eef5f6

Please sign in to comment.