commit e2b548c00cd850d1fa26c53d2a4c9ba5df0ab5c4
Author: nwickel <n.wickelmaier@iwm-tuebingen.de>
Date:   Fri Jan 9 16:55:02 2026 +0100

    Initialize repository

diff --git a/01_intro/exercises.md b/01_intro/exercises.md
new file mode 100644
index 0000000..8b4fd1c
--- /dev/null
+++ b/01_intro/exercises.md
@@ -0,0 +1,51 @@
+Exercise: Binomial test power simulation
+================
+
+## Birth rates
+
+Kanazawa (2007) claims that beautiful parents have more daughters
+
+- Plan a study and calculate the sample size necessary to
+  - detect a deviation from the global 106:100 male-female sex ratio
+  - with about 80% power
+- Wanted: Substance-matter knowledge
+  - What would be a minimum relevant deviation (effect)?
+  - Considering the literature on birth rates, what would be a realistic
+    deviation?
+- Some background
+  - <https://en.wikipedia.org/wiki/Human_sex_ratio>
+  - Literature cited there (e.g., Davis, Gottlieb, and Stampnitzky 1998;
+    Mathews and Hamilton 2005)
+
+### References
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-DavisGottlieb98" class="csl-entry">
+
+Davis, D. L., M. B. Gottlieb, and J. R. Stampnitzky. 1998. “Reduced
+Ratio of Male to Female Births in Several Industrial Countries: A
+Sentinel Health Indicator?” *Journal of the American Medical
+Association* 279 (13): 1018–23.
+<https://doi.org/10.1001/jama.279.13.1018>.
+
+</div>
+
+<div id="ref-Kanazawa07" class="csl-entry">
+
+Kanazawa, S. 2007. “Beautiful Parents Have More Daughters: A Further
+Implication of the Generalized Trivers–Willard Hypothesis (gTWH).”
+*Journal of Theoretical Biology* 244 (1): 133–40.
+<https://doi.org/10.1016/j.jtbi.2006.07.017>.
+
+</div>
+
+<div id="ref-MathewsHamilton05" class="csl-entry">
+
+Mathews, T. J., and B. E. Hamilton. 2005. “Trend Analysis of the Sex
+Ratio at Birth in the United States.” *National Vital Statistics
+Reports* 53 (20): 1–20.
+
+</div>
+
+</div>
diff --git a/01_intro/powersim-intro.pdf b/01_intro/powersim-intro.pdf
new file mode 100644
index 0000000..9531ec0
Binary files /dev/null and b/01_intro/powersim-intro.pdf differ
diff --git a/02_ttest/exercises-ttest1.md b/02_ttest/exercises-ttest1.md
new file mode 100644
index 0000000..520eb83
--- /dev/null
+++ b/02_ttest/exercises-ttest1.md
@@ -0,0 +1,59 @@
+Exercise: Power simulation and power curves for t-test
+================
+
+## Temporal value asymmetry
+
+“participants … were asked to imagine that they had agreed to spend 5 hr
+entering data into a computer and to indicate how much money it would be
+fair for them to receive. Some participants imagined that they had
+completed the work 1 month previously, and others imagined that they
+would complete the work 1 month in the future . . . Participants
+believed that they should receive 101% more money for work they would do
+1 month later ($M = \$125.04$) than for identical work that they had
+done 1 month previously ($M = \$62.20$), $t(119) = 2.22$, $p = .03$,
+$d = 0.41$” (Caruso, Gilbert, and Wilson 2008, 797)
+
+### Plan a direct replication
+
+1.  What is a plausible standard deviation? Hint: $d = (M1 − M2)/SD$
+2.  What is an interesting minimal effect size (in \$, Euro, or min)?
+3.  Simulate responses for 120 participants in both the *past* and the
+    *future* condition, assuming normal distributions with the same
+    variance. Use the standard deviation and the minimal effect size
+    from 1. and 2.
+4.  Parameter recovery: Repeat the simulation from 3. 2000 times to
+    re-estimate the parameters ($\mu_1, \mu_2, \sigma$) from the
+    simulated responses. Visualize the recovered parameters in box
+    plots.  
+    Hint: $SE = 2/\sqrt{n} \cdot SD$, where $n$ is the total sample
+    size.
+
+``` r
+t <- t.test(x, y, mu = 0, var.equal = TRUE)
+c(t$estimate, sd.pool = sqrt(n) / 2 * t$stderr)
+```
+
+5.  Power simulation: Increase the total sample size to find out the `n`
+    necessary for 80% power for the t-test.
+6.  Power curves:
+    - Write an `R` function that takes sample size `n`, minimal effect
+      `d`, standard deviation `sd`, and number of replications `nrep` as
+      arguments. It should return the simulated power.
+    - Use this function to simulate the power for each combination of 4
+      different standard deviations and 4 sample sizes.
+    - Visualize these power curves in a single plot.
+
+### References
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-CarusoGilbert08" class="csl-entry">
+
+Caruso, E. M., D. T. Gilbert, and T. D. Wilson. 2008. “A Wrinkle in
+Time: Asymmetric Valuation of Past and Future Events.” *Psychological
+Science* 19 (8): 796–801.
+<https://doi.org/10.1111/j.1467-9280.2008.02159.x>.
+
+</div>
+
+</div>
diff --git a/02_ttest/exercises-ttest2.md b/02_ttest/exercises-ttest2.md
new file mode 100644
index 0000000..e7aa9cb
--- /dev/null
+++ b/02_ttest/exercises-ttest2.md
@@ -0,0 +1,31 @@
+Exercise: Power simulation and power curves for t-test
+================
+
+## Directed reading activities
+
+Plan a replication of the “directed reading activities” study on
+<https://jasp-stats.github.io/jasp-data-library/myChapters/chapter_2.html>
+
+1.  Which parameter value represents the minimum relevant effect that
+    the test should be able to detect? How large should this effect be?
+    Briefly say what it means in the context of the study.
+
+2.  What would be a plausible value for the standard deviation?
+
+3.  Simulate and draw power curves that depend on the effect size and
+    $n$. Use at least five different effect sizes and at least four
+    different sample sizes. The minimum relevant effect size should be
+    among the simulation conditions.
+
+4.  From inspecting the power curves, what $n$ in each group is required
+    to detect the minimum relevant effect with at least 80% power?
+
+5.  Create a renderable R script or an R Markdown file that includes
+
+    - a header with title, author, date
+    - at least one section head line
+    - the homework questions and your answers
+    - the simulation code and output
+    - the plot of the power curves.
+
+    Render the R or Rmd file to HTML.
diff --git a/02_ttest/powersim-ttest.R b/02_ttest/powersim-ttest.R
new file mode 100644
index 0000000..8344d51
--- /dev/null
+++ b/02_ttest/powersim-ttest.R
@@ -0,0 +1,90 @@
+#' ---
+#' title: "Independent t-test: Power simulation and power curves"
+#' author: ""
+#' date: "Last modified: 2026-01-09"
+#' ---
+
+
+#' ## Application context
+
+#'
+#' Listening experiment
+#'
+#' - Task of each participant is to repeatedly adjust the frequency of a
+#'   comparison tone to sound equal in pitch to a 1000-Hz standard tone
+#' - Mean adjustment estimates the point of subjective equality $\mu$
+#' - Two participants will take part in the experiment providing adjustments
+#'   $X$ and $Y$
+#' - Goal is to detect a difference between their points of subjective
+#'   equality $\mu_x$ and $\mu_y$ of 4 Hz
+#'
+
+#' ## Model
+
+#'
+#' Assumptions
+#' 
+#' - $X_1, \ldots, X_n \sim N(\mu_x, \sigma_x^2)$ i.i.d.
+#' - $Y_1, \ldots, Y_m \sim N(\mu_y, \sigma_y^2)$ i.i.d.
+#' - both samples independent
+#' - $\sigma_x^2 = \sigma_y^2$ but unknown
+#' 
+#' Hypothesis
+#' 
+#' - H$_0\colon~ \mu_x - \mu_y = \delta = 0$
+#'
+
+#' ## Power simulation
+#+ cache = TRUE
+
+n <- 110; m <- 110
+pval <- replicate(2000, {
+  x <- rnorm(n, mean = 1000 + 4, sd = 10)         # Participant 1 responses
+  y <- rnorm(m, mean = 1000,     sd = 10)         # Participant 2 responses
+  t.test(x, y, mu = 0, var.equal = TRUE)$p.value
+})
+mean(pval < 0.05)
+
+#' ## Power curves
+#'
+#' Turn into a function of n and effect size
+
+pwrFun <- function(n = 30, d = 4, sd = 10, nrep = 50) {
+  n <- n; m <- n
+  pval <- replicate(nrep, {
+    x <- rnorm(n, mean = 1000 + d, sd = sd)
+    y <- rnorm(m, mean = 1000,     sd = sd)
+    t.test(x, y, mu = 0, var.equal = TRUE)$p.value
+  })
+  mean(pval < 0.05)
+}
+
+#'
+#' Set up conditions and call power function
+#+ cache = TRUE
+
+cond <- expand.grid(d = 0:5,
+                    n = c(50, 75, 100, 125))
+cond$pwr <- mapply(pwrFun, n = cond$n, d = cond$d, MoreArgs = list(nrep = 500))
+
+## Plot results
+lattice::xyplot(pwr ~ d, cond, groups = n, type = c("g", "b"),
+                auto.key = list(corner = c(0, 1)))
+
+#' ## Parameter recovery
+
+n <- 100
+
+out <- replicate(2000, {
+  x <- rnorm(n, mean = 1000 + 4, sd = 10)
+  y <- rnorm(n, mean = 1000,     sd = 10)
+  t <- t.test(x, y, mu = 0, var.equal = TRUE)
+  c(
+    effect = -as.numeric(diff(t$estimate)),
+    sd.pool = t$stderr * sqrt(2*n) / 2
+  )
+})
+
+rowMeans(out)
+boxplot(t(out))
+
diff --git a/02_ttest/powersim-ttest.md b/02_ttest/powersim-ttest.md
new file mode 100644
index 0000000..de3c29c
--- /dev/null
+++ b/02_ttest/powersim-ttest.md
@@ -0,0 +1,99 @@
+Independent t-test: Power simulation and power curves
+================
+Last modified: 2026-01-09
+
+## Application context
+
+Listening experiment
+
+- Task of each participant is to repeatedly adjust the frequency of a
+  comparison tone to sound equal in pitch to a 1000-Hz standard tone
+- Mean adjustment estimates the point of subjective equality $\mu$
+- Two participants will take part in the experiment providing
+  adjustments $X$ and $Y$
+- Goal is to detect a difference between their points of subjective
+  equality $\mu_x$ and $\mu_y$ of 4 Hz
+
+## Model
+
+Assumptions
+
+- $X_1, \ldots, X_n \sim N(\mu_x, \sigma_x^2)$ i.i.d.
+- $Y_1, \ldots, Y_m \sim N(\mu_y, \sigma_y^2)$ i.i.d.
+- both samples independent
+- $\sigma_x^2 = \sigma_y^2$ but unknown
+
+Hypothesis
+
+- H$_0\colon~ \mu_x - \mu_y = \delta = 0$
+
+## Power simulation
+
+``` r
+n <- 110; m <- 110
+pval <- replicate(2000, {
+  x <- rnorm(n, mean = 1000 + 4, sd = 10)         # Participant 1 responses
+  y <- rnorm(m, mean = 1000,     sd = 10)         # Participant 2 responses
+  t.test(x, y, mu = 0, var.equal = TRUE)$p.value
+})
+mean(pval < 0.05)
+```
+
+    ## [1] 0.831
+
+## Power curves
+
+Turn into a function of n and effect size
+
+``` r
+pwrFun <- function(n = 30, d = 4, sd = 10, nrep = 50) {
+  n <- n; m <- n
+  pval <- replicate(nrep, {
+    x <- rnorm(n, mean = 1000 + d, sd = sd)
+    y <- rnorm(m, mean = 1000,     sd = sd)
+    t.test(x, y, mu = 0, var.equal = TRUE)$p.value
+  })
+  mean(pval < 0.05)
+}
+```
+
+Set up conditions and call power function
+
+``` r
+cond <- expand.grid(d = 0:5,
+                    n = c(50, 75, 100, 125))
+cond$pwr <- mapply(pwrFun, n = cond$n, d = cond$d, MoreArgs = list(nrep = 500))
+
+## Plot results
+lattice::xyplot(pwr ~ d, cond, groups = n, type = c("g", "b"),
+                auto.key = list(corner = c(0, 1)))
+```
+
+![](powersim-ttest_files/figure-gfm/unnamed-chunk-3-1.png)<!-- -->
+
+## Parameter recovery
+
+``` r
+n <- 100
+
+out <- replicate(2000, {
+  x <- rnorm(n, mean = 1000 + 4, sd = 10)
+  y <- rnorm(n, mean = 1000,     sd = 10)
+  t <- t.test(x, y, mu = 0, var.equal = TRUE)
+  c(
+    effect = -as.numeric(diff(t$estimate)),
+    sd.pool = t$stderr * sqrt(2*n) / 2
+  )
+})
+
+rowMeans(out)
+```
+
+    ##   effect  sd.pool 
+    ## 3.978400 9.993142
+
+``` r
+boxplot(t(out))
+```
+
+![](powersim-ttest_files/figure-gfm/unnamed-chunk-4-1.png)<!-- -->
diff --git a/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-3-1.png b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-3-1.png
new file mode 100644
index 0000000..bfc8331
Binary files /dev/null and b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-3-1.png differ
diff --git a/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-4-1.png b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-4-1.png
new file mode 100644
index 0000000..62867da
Binary files /dev/null and b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-4-1.png differ
diff --git a/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-5-1.png b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-5-1.png
new file mode 100644
index 0000000..8101b10
Binary files /dev/null and b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-5-1.png differ
diff --git a/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-5-2.png b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-5-2.png
new file mode 100644
index 0000000..a419880
Binary files /dev/null and b/02_ttest/powersim-ttest_files/figure-gfm/unnamed-chunk-5-2.png differ
diff --git a/03_anova/exercises-anova.md b/03_anova/exercises-anova.md
new file mode 100644
index 0000000..b476449
--- /dev/null
+++ b/03_anova/exercises-anova.md
@@ -0,0 +1,62 @@
+Exercise: Testing the interaction in two-by-two ANOVA
+================
+
+## Anchoring and adjustment
+
+Items and anchor values (Jacowitz and Kahneman 1995)
+
+- How tall is the largest coast redwood in the world? \[20, 168m\]
+- How many member states belong to the United Nations? \[14, 127
+  members\]
+- How much km/h is the maximum speed of a house cat? \[11, 48km/h\]
+
+### Research question
+
+- Does time pressure (respond within 7s) increase the anchor effect?
+
+### Suggest a minimum relevant effect
+
+- Go to <https://apps.mathpsy.uni-tuebingen.de/fw/pars2eta/>
+- Fix the parameters of the ANOVA model
+
+### Some background
+
+- Open anchoring quest (Röseler et al. 2022), <https://osf.io/ygnvb/>
+
+### Plan the study
+
+- Pick one of the three items
+- Parameter recovery
+  - Make a data frame for the two-by-two design
+  - With the parameter values determined before, simulate responses
+  - Re-estimate the parameters
+- Power simulation
+  - Calculate the sample size necessary to detect the time-pressure
+    effect
+
+### Bonus task: Verify the plausibility of your model
+
+- Download the raw data from the open anchoring quest project
+- Estimate $\sigma$ and compare it to your value
+
+### References
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-JacowitzKahneman95" class="csl-entry">
+
+Jacowitz, K. E., and D. Kahneman. 1995. “Measures of Anchoring in
+Estimation Tasks.” *Personality and Social Psychology Bulletin* 21 (11):
+1161–66. <https://doi.org/10.1177/01461672952111004>.
+
+</div>
+
+<div id="ref-RoeselerWeber22" class="csl-entry">
+
+Röseler, L., L. Weber, K. A. C. Helgerth, E. Stich, M. Günther, P.
+Tegethoff, F. S. Wagner, et al. 2022. “OpAQ: Open Anchoring Quest,
+Version 1.1.50.97.” <https://doi.org/10.17605/OSF.IO/YGNVB>.
+
+</div>
+
+</div>
diff --git a/03_anova/powersim-anova.R b/03_anova/powersim-anova.R
new file mode 100644
index 0000000..7f4ef02
--- /dev/null
+++ b/03_anova/powersim-anova.R
@@ -0,0 +1,93 @@
+#' ---
+#' title: "Two-by-two ANOVA: Power simulation of the interaction test"
+#' author: ""
+#' date: "Last modified: 2026-01-09"
+#' ---
+
+#' ## Application context
+
+#'
+#' Effect of fertilizers
+#'
+#' In an experiment, two fertilizers (A and B, each either low or high dose)
+#' will be combined and the yield of peas (Y) in kg be observed. Goal is to
+#' detect an increase of the Fertilizer-A effect by an additional 12 kg when
+#' combined with a high dose of Fertilizer B (interaction effect).
+#'
+
+#+ echo = FALSE
+dat <- data.frame(
+  A = rep(1:2, each = 2),
+  B = rep(1:2, times = 2),
+  y = c(30, 30 + 5, 30 + 30, 30 + 30 + 5 + 12)
+)
+par(mai = c(.6, .6, .1, .1), mgp = c(2, .7, 0))
+plot(y ~ A, dat, type = "n", xlim = c(0.8, 2.2), ylim = c(20, 80),
+     xlab = "Fertilizer A", ylab = "Yield (kg)", xaxt = "n")
+lines(y ~ A, dat[dat$B == 1, ], col = "darkblue")
+lines(y ~ A, dat[dat$B == 2, ], col = "darkblue")
+lines(1:2, c(30 + 5, 30 + 30 + 5), lty = 2, col = "darkblue")
+axis(1, 1:2, c("low", "high"))
+arrows(2, 30 + 30 + 6, 2, 30 + 30 + 5 + 11, code = 3, length = 0.1,
+       col = "darkgray")
+text(c(1, 2, 1, 2, 2.07), c(27, 55, 40, 80, 65 + 6),
+     c(expression(mu), expression(mu + alpha[2]),
+       expression(mu + beta[2]),
+       expression(mu + alpha[2] + beta[2] + (alpha * beta)[22]),
+       "12 kg")
+)
+text(1.5, 27 + 30/2, "Fertilizer B: low", srt = 21, col = "darkgray")
+text(1.5, 38 + (30 + 12)/2, "Fertilizer B: high", srt = 32, col = "darkgray")
+
+#' ## Model
+
+#'
+#' Assumptions
+#' 
+#' - $Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} +
+#'              \varepsilon_{ijk}$
+#' - $\varepsilon_{ijk} \sim N(0, \sigma^2) \text{ i.i.d.}$
+#' - $i = 1, \dots, I$; $j = 1, \dots, J$; $k = 1, \dots, K$
+#' - $\alpha_1 = \beta_1 := 0$
+#' 
+#' Hypothesis
+#' 
+#' - H$_0^{AB}\colon~ (\alpha\beta)_{ij} = 0 \text{ for all } i,j$
+#' 
+
+#' ## Setup
+
+set.seed(1704)
+n <- 96
+dat <- data.frame(
+  A = factor(rep(1:2, each = n/2), labels = c("low", "high")),
+  B = factor(rep(rep(1:2, each = n/4), 2), labels = c("low", "high"))
+)
+X <- model.matrix(~ A*B, dat)
+unique(X)
+beta <- c(mu = 30, a2 = 30, b2 = 5, ab22 = 12)
+means <- X %*% beta
+
+lattice::xyplot(I(means + rnorm(n, sd = 10)) ~ A, dat, groups = B,
+                type = c("g", "p", "a"), auto.key = TRUE, ylab = "Yield (kg)")
+
+#' ## Parameter recovery
+#+ cache = TRUE
+
+out <- replicate(2000, {
+  y <- means + rnorm(n, sd = 10)   # y = mu + a + b + ab + e
+  m <- aov(y ~ A * B, dat)
+  c(coef(m), sigma = sigma(m))
+})
+boxplot(t(out))
+
+#' ## Power simulation
+#+ cache = TRUE
+
+pval <- replicate(2000, {
+  y <- means + rnorm(n, sd = 10)
+  m <- aov(y ~ A*B, dat)
+  summary(m)[[1]]$"Pr(>F)"[3]      # test of interaction
+})
+mean(pval < 0.05)
+
diff --git a/03_anova/powersim-anova.md b/03_anova/powersim-anova.md
new file mode 100644
index 0000000..5e5761d
--- /dev/null
+++ b/03_anova/powersim-anova.md
@@ -0,0 +1,82 @@
+Two-by-two ANOVA: Power simulation of the interaction test
+================
+Last modified: 2026-01-09
+
+## Application context
+
+Effect of fertilizers
+
+In an experiment, two fertilizers (A and B, each either low or high
+dose) will be combined and the yield of peas (Y) in kg be observed. Goal
+is to detect an increase of the Fertilizer-A effect by an additional 12
+kg when combined with a high dose of Fertilizer B (interaction effect).
+
+![](powersim-anova_files/figure-gfm/unnamed-chunk-1-1.png)<!-- -->
+
+## Model
+
+Assumptions
+
+- $Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} +  \varepsilon_{ijk}$
+- $\varepsilon_{ijk} \sim N(0, \sigma^2) \text{ i.i.d.}$
+- $i = 1, \dots, I$; $j = 1, \dots, J$; $k = 1, \dots, K$
+- $\alpha_1 = \beta_1 := 0$
+
+Hypothesis
+
+- H$_0^{AB}\colon~ (\alpha\beta)_{ij} = 0 \text{ for all } i,j$
+
+## Setup
+
+``` r
+set.seed(1704)
+n <- 96
+dat <- data.frame(
+  A = factor(rep(1:2, each = n/2), labels = c("low", "high")),
+  B = factor(rep(rep(1:2, each = n/4), 2), labels = c("low", "high"))
+)
+X <- model.matrix(~ A*B, dat)
+unique(X)
+```
+
+    ##    (Intercept) Ahigh Bhigh Ahigh:Bhigh
+    ## 1            1     0     0           0
+    ## 25           1     0     1           0
+    ## 49           1     1     0           0
+    ## 73           1     1     1           1
+
+``` r
+beta <- c(mu = 30, a2 = 30, b2 = 5, ab22 = 12)
+means <- X %*% beta
+
+lattice::xyplot(I(means + rnorm(n, sd = 10)) ~ A, dat, groups = B,
+                type = c("g", "p", "a"), auto.key = TRUE, ylab = "Yield (kg)")
+```
+
+![](powersim-anova_files/figure-gfm/unnamed-chunk-2-1.png)<!-- -->
+
+## Parameter recovery
+
+``` r
+out <- replicate(2000, {
+  y <- means + rnorm(n, sd = 10)   # y = mu + a + b + ab + e
+  m <- aov(y ~ A * B, dat)
+  c(coef(m), sigma = sigma(m))
+})
+boxplot(t(out))
+```
+
+![](powersim-anova_files/figure-gfm/unnamed-chunk-3-1.png)<!-- -->
+
+## Power simulation
+
+``` r
+pval <- replicate(2000, {
+  y <- means + rnorm(n, sd = 10)
+  m <- aov(y ~ A*B, dat)
+  summary(m)[[1]]$"Pr(>F)"[3]      # test of interaction
+})
+mean(pval < 0.05)
+```
+
+    ## [1] 0.817
diff --git a/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-1-1.png b/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-1-1.png
new file mode 100644
index 0000000..81c32ff
Binary files /dev/null and b/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-1-1.png differ
diff --git a/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-2-1.png b/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-2-1.png
new file mode 100644
index 0000000..b7cb901
Binary files /dev/null and b/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-2-1.png differ
diff --git a/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-3-1.png b/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-3-1.png
new file mode 100644
index 0000000..d0c2c28
Binary files /dev/null and b/03_anova/powersim-anova_files/figure-gfm/unnamed-chunk-3-1.png differ
diff --git a/04_ancova/exercises-ancova1.md b/04_ancova/exercises-ancova1.md
new file mode 100644
index 0000000..cc63a15
--- /dev/null
+++ b/04_ancova/exercises-ancova1.md
@@ -0,0 +1,33 @@
+Exercise: Analysis and power simulation for baseline/follow-up
+measurements
+================
+
+## Shoulder pain and acupuncture
+
+1.  Reanalyze the original data
+    - Re-estimate the ANCOVA model for the Kleinhenz et al. (1999)
+      [data](../data/kleinhenz.txt)
+2.  Run a power simulation for a replication study
+    1.  Draw plausible pre-CMS values
+    2.  Specify the minimum relevant average treatment effect (ATE)
+    3.  Set the remaining parameters to plausible values
+    4.  What is the sample size required for the test to detect the
+        effect with 80% power?
+    5.  How robust is the power simulation when you repeat it with a new
+        set of pre-CMS values? Try it!
+    6.  Recover the parameters of the ANCOVA model
+
+### References
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-KleinhenzStreitberger99" class="csl-entry">
+
+Kleinhenz, J., K. Streitberger, J. Windeler, A. Güßbacher, G. Mavridis,
+and E. Martin. 1999. “Randomised Clinical Trial Comparing the Effects of
+Acupuncture and a Newly Designed Placebo Needle in Rotator Cuff
+Tendinitis.” *Pain* 83 (2): 235–41.
+
+</div>
+
+</div>
diff --git a/04_ancova/exercises-ancova2.md b/04_ancova/exercises-ancova2.md
new file mode 100644
index 0000000..8a26156
--- /dev/null
+++ b/04_ancova/exercises-ancova2.md
@@ -0,0 +1,45 @@
+Exercise: Analysis and power simulation for baseline/follow-up
+measurements
+================
+
+## MASS anorexia data
+
+1.  Analyze the original data:
+    - In R, see ?MASS::anorexia
+    - Data preparation
+
+``` r
+data(anorexia, package = "MASS")
+dat <-
+  subset(anorexia, Treat != "Cont") |>  # exclude control group
+  droplevels()                          # drop empty factor levels
+lbs2kg <- 0.4535924
+dat$Prewt  <- lbs2kg * dat$Prewt        # to kg
+dat$Postwt <- lbs2kg * dat$Postwt
+lattice::xyplot(Postwt ~ Prewt, dat, groups = Treat,
+                type = c("g", "r", "p"), auto.key = TRUE)
+```
+
+- Estimate the average treatment effect (ATE) for FT relative to CBT.
+- What is the 95% CI for the ATE?
+- What are the pre- and post-weight means for the two groups?
+- What are the baseline-adjusted means for the two groups?
+
+2.  Run a power simulation for a replication study:
+    - Draw plausible pre-weights.
+    - Specify the minimum relevant effect.
+    - Set the remaining parameters to plausible values.
+    - What is the sample size required for the test to detect the effect
+      with 80% power?
+    - How robust is the power simulation when you repeat it with a new
+      set of pre-weights? Try it!
+    - Recover the parameters of the ANCOVA model.
+3.  Create a renderable R script or an R Markdown file that includes
+    - a header with title, author, date
+    - at least one section head line
+    - the homework questions and your answers
+    - the R code, output, and plots (if any)
+
+    Render the R or Rmd file to HTML.
+
+### Reference
diff --git a/04_ancova/intro-ancova.pdf b/04_ancova/intro-ancova.pdf
new file mode 100644
index 0000000..c9da59d
Binary files /dev/null and b/04_ancova/intro-ancova.pdf differ
diff --git a/05_mixed1/exercises-pwrlmm.md b/05_mixed1/exercises-pwrlmm.md
new file mode 100644
index 0000000..b1f1b1d
--- /dev/null
+++ b/05_mixed1/exercises-pwrlmm.md
@@ -0,0 +1,50 @@
+Exercise: Power simulation for longitudinal data
+================
+
+## Risperidone vs. haloperidol and schizophrenia
+
+``` r
+dat       <- read.table("../data/moeller.csv", header = TRUE, sep = ",")
+dat$id    <- factor(dat$id)
+dat$treat <- factor(dat$treat, levels = c("risp", "halo"))
+
+lattice::xyplot(pans ~ week, data = dat, groups = treat, type = c("g", "p", "a"), auto.key = TRUE)
+```
+
+1)  Analyze the original data from [moeller.csv](../data/moeller.csv):
+    - `pans`: Positive and Negative Symptom Scale for schizophrenia
+
+    - `treat`: medication group
+
+      - `risp`: atypical neuroleptic risperidone
+      - `halo`: conventional neuroleptic haloperidol
+
+    - What is the sample size in each treatment group?
+
+    - Estimate the by-group random-slope model.
+
+    - What are the estimates for the fixed effects and variance
+      components?
+
+    - Interpret the interaction effect.
+
+    - Test the interaction effect.
+2)  Run a power simulation for a replication study:
+    - Set up a data frame containing the study design and sample size.
+
+    - Specify the minimum relevant effect.
+
+    - Set the fixed effects and variance components to plausible values.
+
+    - How many participants are required for the test of the interaction
+      to detect the specified effect with a power of 80%?
+
+    - Recover the parameters of the by-group random-slope model for one
+      simulated data set.
+3)  Create a renderable R script or an R Markdown file that includes
+    - a header with title, author, date
+    - at least one section head line
+    - the questions from above and your answers
+    - the R code, output, and plots (if any)
+
+    Render the R or Rmd file to HTML.
diff --git a/05_mixed1/intro-lmm.pdf b/05_mixed1/intro-lmm.pdf
new file mode 100644
index 0000000..f1e4113
Binary files /dev/null and b/05_mixed1/intro-lmm.pdf differ
diff --git a/05_mixed1/powersim-lmm.R b/05_mixed1/powersim-lmm.R
new file mode 100644
index 0000000..fee09f4
--- /dev/null
+++ b/05_mixed1/powersim-lmm.R
@@ -0,0 +1,185 @@
+#' ---
+#' title: "Power for mixed-effects models"
+#' author: ""
+#' date: "Last modified: 2026-01-09"
+#' bibliography: ../lit.bib
+#' ---
+
+library(lattice)
+library(lme4)
+
+#' # Reanalysis
+
+#' ## Application context: Depression and type of diagnosis
+#'
+#' - @ReisbyGram77 studied the effect of Imipramin on 66 inpatients
+#'   treated for depression
+#' - Depression was measured with the Hamilton depression rating scale
+#' - Patients were classified into endogenous and non-endogenous depressed
+#' - Depression was measured weekly for 6 time points
+#'
+#' Data: [reisby.txt](../data/reisby.txt)
+
+dat      <- read.table("../data/reisby.txt", header = TRUE)
+dat$id   <- factor(dat$id)
+dat$diag <- factor(dat$diag, levels = c("nonen", "endog"))
+dat      <- na.omit(dat)     # drop missing values
+head(dat, n = 13)
+
+xyplot(hamd ~ week | id, data = dat, type=c("g", "r", "p"),
+  pch = 16, layout = c(11, 6), ylab = "HDRS score", xlab = "Time (week)")
+
+
+#' ## Random-intercept model
+#'
+#' $$
+#' \begin{aligned}
+#' Y_{ij} &= \beta_0 + \beta_1 \, \mathtt{week}_{ij}
+#'         + \upsilon_{0i}
+#'         + \varepsilon_{ij} \\
+#' \upsilon_{0i} &\sim N(0, \sigma^2_{\upsilon_0}) \text{ i.i.d.} \\
+#' \mathbf{\varepsilon}_i &\sim N(0, \, \sigma^2)  \text{ i.i.d.} \\
+#' i &= 1, \ldots, I, \quad j = 1, \ldots n_i
+#' \end{aligned}
+#' $$
+
+m1 <- lmer(hamd ~ week + (1 | id), data = dat, REML = FALSE)
+summary(m1)
+
+#' ## Random-slope model
+#'
+#' $$
+#' \begin{aligned}
+#' Y_{ij} &= \beta_0 + \beta_1 \, \mathtt{week}_{ij}
+#'         + \upsilon_{0i} + \upsilon_{1i}\, \mathtt{week}_{ij}
+#'         + \varepsilon_{ij} \\
+#' \begin{pmatrix} \upsilon_{0i}\\ \upsilon_{1i} \end{pmatrix} &\sim
+#'    N \left(\begin{pmatrix} 0\\ 0 \end{pmatrix}, \,
+#'            \mathbf{\Sigma}_\upsilon =
+#' \begin{pmatrix}
+#'           \sigma^2_{\upsilon_0} & \sigma_{\upsilon_0 \upsilon_1} \\
+#'           \sigma_{\upsilon_0 \upsilon_1} & \sigma^2_{\upsilon_1} \\
+#' \end{pmatrix} \right)
+#'       \text{ i.i.d.} \\
+#' \mathbf{\varepsilon}_i &\sim N(\mathbf{0}, \, \sigma^2 \mathbf{I}_{n_i})
+#'       \text{ i.i.d.} \\
+#' i &= 1, \ldots, I, \quad j = 1, \ldots n_i
+#' \end{aligned}
+#' $$
+
+m2 <- lmer(hamd ~ week + (week | id), data = dat, REML = FALSE)
+summary(m2)
+
+#' ## Partial pooling
+
+#+ fig.height = 10, fig.width = 6.5, fig.aling = "center"
+indiv <- unlist(
+  sapply(unique(dat$id),
+         function(i) predict(lm(hamd ~ week, dat[dat$id == i, ])))
+)
+
+xyplot(hamd + predict(m2, re.form = ~ 0) + predict(m2) + indiv ~ week | id,
+  data = dat, type = c("p", "l", "l", "l"), pch = 16, grid = TRUE,
+  distribute.type = TRUE, layout = c(11, 6), ylab = "HDRS score",
+  xlab = "Time (week)",
+  # customize colors
+  col = c("#434F4F", "#3CB4DC", "#FF6900", "#78004B"),
+  # add legend
+  key = list(space = "top", columns = 3,
+             text = list(c("Population", "Mixed model", "Within-subject")),
+             lines = list(col = c("#3CB4DC", "#FF6900", "#78004B")))
+  )
+
+#' ## By-group random-slope model
+
+m3 <- lmer(hamd ~ week + diag + (week | id), data = dat, REML = FALSE)
+m4 <- lmer(hamd ~ week * diag + (week | id), data = dat, REML = FALSE)
+anova(m3, m4)
+
+#' ## Means and predicted HDRS score by group
+
+dat2 <- aggregate(hamd ~ week + diag, dat, mean)
+dat2$m4 <- predict(m4, newdata = dat2, re.form = ~ 0)
+
+plot(m4 ~ week, dat2[dat2$diag == "endog", ], type = "l",
+     ylim=c(0, 28), xlab="Week", ylab = "HDRS score")
+lines(m4 ~ week, dat2[dat2$diag == "nonen", ], lty = 2)
+points(hamd ~ week, dat2[dat2$diag == "endog", ], pch = 16)
+points(hamd ~ week, dat2[dat2$diag == "nonen", ], pch = 21, bg = "white")
+legend("topright", c("Endogenous", "Non endogenous"),
+       lty = 1:2, pch = c(16, 21), pt.bg = "white", bty = "n")
+
+#' ## Gory details
+
+fixef(m4)
+getME(m4, "theta")
+t(chol(VarCorr(m4)$id))[lower.tri(diag(2), diag = TRUE)] / sigma(m4)
+
+#' # Power simulation
+
+#' ## Setup
+
+## Study design and sample sizes
+n_week <- 6
+n_subj <- 80
+n <- n_week * n_subj
+
+dat <- data.frame(
+     id = factor(rep(seq_len(n_subj), each = n_week)),
+   week = rep(0:(n_week - 1), times = n_subj),
+  treat = factor(rep(0:1, each = n/2), labels = c("ctr", "trt"))
+)
+
+## Fixed effects and variance components
+beta <- c(22, -2, 0, -1)
+Su <- matrix(c(12, -1.5, -1.5, 2), nrow = 2)
+se <- 3.5
+
+#' ## Power
+
+#+ cache = TRUE, warning = FALSE
+pval <- replicate(200, {
+
+  # Data generation
+  means <- model.matrix( ~ week * treat, dat) %*% beta
+  ranu <- MASS::mvrnorm(n_subj, mu = c(0, 0), Sigma = Su)
+  e <- rnorm(n_subj * n_week, mean = 0, sd = se)
+  
+  y <- means + ranu[dat$id, 1] + ranu[dat$id, 2] * dat$week + e
+  
+  # Fitting model to test H0
+  m0 <- lmer(y ~ week + treat + (1 + week | id), data = dat, REML = FALSE)
+  m1 <- lmer(y ~ week * treat + (1 + week | id), data = dat, REML = FALSE)
+  anova(m0, m1)["m1", "Pr(>Chisq)"]
+  }
+)
+
+mean(pval < 0.05)
+
+#' ## Parameter recovery
+
+#+ cache = TRUE, warning = FALSE
+par <- replicate(200, {
+  
+  means <- model.matrix( ~ week * treat, dat) %*% beta
+  ranu <- MASS::mvrnorm(n_subj, mu = c(0, 0), Sigma = Su)
+  e <- rnorm(n_subj * n_week, mean = 0, sd = se)
+  
+  y <- means + ranu[dat$id, 1] + ranu[dat$id, 2] * dat$week + e
+  
+  m1 <- lmer(y ~ week * treat + (1 + week | id), data = dat, REML = FALSE)
+  list(fixef = fixef(m1), theta = getME(m1, "theta"), sigma = sigma(m1))
+  }, simplify = FALSE
+)
+
+rowMeans(sapply(par, function(x) x$fixef))
+rowMeans(sapply(par, function(x) x$theta))
+mean(sapply(par, function(x) x$sigma))
+
+beta
+Lt <- chol(Su)
+t(Lt)[lower.tri(Lt, diag = TRUE)] / se
+se
+
+#' ### References
+
diff --git a/05_mixed1/powersim-lmm.md b/05_mixed1/powersim-lmm.md
new file mode 100644
index 0000000..a0fc3b9
--- /dev/null
+++ b/05_mixed1/powersim-lmm.md
@@ -0,0 +1,340 @@
+Power for mixed-effects models
+================
+Last modified: 2026-01-09
+
+``` r
+library(lattice)
+library(lme4)
+```
+
+# Reanalysis
+
+## Application context: Depression and type of diagnosis
+
+- Reisby et al. (1977) studied the effect of Imipramin on 66 inpatients
+  treated for depression
+- Depression was measured with the Hamilton depression rating scale
+- Patients were classified into endogenous and non-endogenous depressed
+- Depression was measured weekly for 6 time points
+
+Data: [reisby.txt](../data/reisby.txt)
+
+``` r
+dat      <- read.table("../data/reisby.txt", header = TRUE)
+dat$id   <- factor(dat$id)
+dat$diag <- factor(dat$diag, levels = c("nonen", "endog"))
+dat      <- na.omit(dat)     # drop missing values
+head(dat, n = 13)
+```
+
+    ##     id hamd week  diag endweek
+    ## 1  101   26    0 nonen       0
+    ## 2  101   22    1 nonen       0
+    ## 3  101   18    2 nonen       0
+    ## 4  101    7    3 nonen       0
+    ## 5  101    4    4 nonen       0
+    ## 6  101    3    5 nonen       0
+    ## 7  103   33    0 nonen       0
+    ## 8  103   24    1 nonen       0
+    ## 9  103   15    2 nonen       0
+    ## 10 103   24    3 nonen       0
+    ## 11 103   15    4 nonen       0
+    ## 12 103   13    5 nonen       0
+    ## 13 104   29    0 endog       0
+
+``` r
+xyplot(hamd ~ week | id, data = dat, type=c("g", "r", "p"),
+  pch = 16, layout = c(11, 6), ylab = "HDRS score", xlab = "Time (week)")
+```
+
+![](powersim-lmm_files/figure-gfm/unnamed-chunk-2-1.png)<!-- -->
+
+## Random-intercept model
+
+$$
+\begin{aligned}
+Y_{ij} &= \beta_0 + \beta_1 \, \mathtt{week}_{ij}
+        + \upsilon_{0i}
+        + \varepsilon_{ij} \\
+\upsilon_{0i} &\sim N(0, \sigma^2_{\upsilon_0}) \text{ i.i.d.} \\
+\mathbf{\varepsilon}_i &\sim N(0, \, \sigma^2)  \text{ i.i.d.} \\
+i &= 1, \ldots, I, \quad j = 1, \ldots n_i
+\end{aligned}
+$$
+
+``` r
+m1 <- lmer(hamd ~ week + (1 | id), data = dat, REML = FALSE)
+summary(m1)
+```
+
+    ## Linear mixed model fit by maximum likelihood  ['lmerMod']
+    ## Formula: hamd ~ week + (1 | id)
+    ##    Data: dat
+    ## 
+    ##       AIC       BIC    logLik -2*log(L)  df.resid 
+    ##    2293.2    2308.9   -1142.6    2285.2       371 
+    ## 
+    ## Scaled residuals: 
+    ##     Min      1Q  Median      3Q     Max 
+    ## -3.1739 -0.5876 -0.0342  0.5465  3.5297 
+    ## 
+    ## Random effects:
+    ##  Groups   Name        Variance Std.Dev.
+    ##  id       (Intercept) 16.16    4.019   
+    ##  Residual             19.04    4.363   
+    ## Number of obs: 375, groups:  id, 66
+    ## 
+    ## Fixed effects:
+    ##             Estimate Std. Error t value
+    ## (Intercept)  23.5518     0.6385   36.88
+    ## week         -2.3757     0.1350  -17.60
+    ## 
+    ## Correlation of Fixed Effects:
+    ##      (Intr)
+    ## week -0.524
+
+## Random-slope model
+
+$$
+\begin{aligned}
+Y_{ij} &= \beta_0 + \beta_1 \, \mathtt{week}_{ij}
+        + \upsilon_{0i} + \upsilon_{1i}\, \mathtt{week}_{ij}
+        + \varepsilon_{ij} \\
+\begin{pmatrix} \upsilon_{0i}\\ \upsilon_{1i} \end{pmatrix} &\sim
+   N \left(\begin{pmatrix} 0\\ 0 \end{pmatrix}, \,
+           \mathbf{\Sigma}_\upsilon =
+\begin{pmatrix}
+          \sigma^2_{\upsilon_0} & \sigma_{\upsilon_0 \upsilon_1} \\
+          \sigma_{\upsilon_0 \upsilon_1} & \sigma^2_{\upsilon_1} \\
+\end{pmatrix} \right)
+      \text{ i.i.d.} \\
+\mathbf{\varepsilon}_i &\sim N(\mathbf{0}, \, \sigma^2 \mathbf{I}_{n_i})
+      \text{ i.i.d.} \\
+i &= 1, \ldots, I, \quad j = 1, \ldots n_i
+\end{aligned}
+$$
+
+``` r
+m2 <- lmer(hamd ~ week + (week | id), data = dat, REML = FALSE)
+summary(m2)
+```
+
+    ## Linear mixed model fit by maximum likelihood  ['lmerMod']
+    ## Formula: hamd ~ week + (week | id)
+    ##    Data: dat
+    ## 
+    ##       AIC       BIC    logLik -2*log(L)  df.resid 
+    ##    2231.0    2254.6   -1109.5    2219.0       369 
+    ## 
+    ## Scaled residuals: 
+    ##     Min      1Q  Median      3Q     Max 
+    ## -2.7460 -0.5016  0.0332  0.5177  3.6834 
+    ## 
+    ## Random effects:
+    ##  Groups   Name        Variance Std.Dev. Corr 
+    ##  id       (Intercept) 12.631   3.554         
+    ##           week         2.079   1.442    -0.28
+    ##  Residual             12.216   3.495         
+    ## Number of obs: 375, groups:  id, 66
+    ## 
+    ## Fixed effects:
+    ##             Estimate Std. Error t value
+    ## (Intercept)  23.5769     0.5456   43.22
+    ## week         -2.3771     0.2087  -11.39
+    ## 
+    ## Correlation of Fixed Effects:
+    ##      (Intr)
+    ## week -0.449
+
+## Partial pooling
+
+``` r
+indiv <- unlist(
+  sapply(unique(dat$id),
+         function(i) predict(lm(hamd ~ week, dat[dat$id == i, ])))
+)
+
+xyplot(hamd + predict(m2, re.form = ~ 0) + predict(m2) + indiv ~ week | id,
+  data = dat, type = c("p", "l", "l", "l"), pch = 16, grid = TRUE,
+  distribute.type = TRUE, layout = c(11, 6), ylab = "HDRS score",
+  xlab = "Time (week)",
+  # customize colors
+  col = c("#434F4F", "#3CB4DC", "#FF6900", "#78004B"),
+  # add legend
+  key = list(space = "top", columns = 3,
+             text = list(c("Population", "Mixed model", "Within-subject")),
+             lines = list(col = c("#3CB4DC", "#FF6900", "#78004B")))
+  )
+```
+
+![](powersim-lmm_files/figure-gfm/unnamed-chunk-5-1.png)<!-- -->
+
+## By-group random-slope model
+
+``` r
+m3 <- lmer(hamd ~ week + diag + (week | id), data = dat, REML = FALSE)
+m4 <- lmer(hamd ~ week * diag + (week | id), data = dat, REML = FALSE)
+anova(m3, m4)
+```
+
+    ## Data: dat
+    ## Models:
+    ## m3: hamd ~ week + diag + (week | id)
+    ## m4: hamd ~ week * diag + (week | id)
+    ##    npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
+    ## m3    7 2228.9 2256.4 -1107.5    2214.9                     
+    ## m4    8 2230.9 2262.3 -1107.5    2214.9 0.0042  1     0.9486
+
+## Means and predicted HDRS score by group
+
+``` r
+dat2 <- aggregate(hamd ~ week + diag, dat, mean)
+dat2$m4 <- predict(m4, newdata = dat2, re.form = ~ 0)
+
+plot(m4 ~ week, dat2[dat2$diag == "endog", ], type = "l",
+     ylim=c(0, 28), xlab="Week", ylab = "HDRS score")
+lines(m4 ~ week, dat2[dat2$diag == "nonen", ], lty = 2)
+points(hamd ~ week, dat2[dat2$diag == "endog", ], pch = 16)
+points(hamd ~ week, dat2[dat2$diag == "nonen", ], pch = 21, bg = "white")
+legend("topright", c("Endogenous", "Non endogenous"),
+       lty = 1:2, pch = c(16, 21), pt.bg = "white", bty = "n")
+```
+
+![](powersim-lmm_files/figure-gfm/unnamed-chunk-7-1.png)<!-- -->
+
+## Gory details
+
+``` r
+fixef(m4)
+```
+
+    ##    (Intercept)           week      diagendog week:diagendog 
+    ##    22.47626332    -2.36568746     1.98802087    -0.02705576
+
+``` r
+getME(m4, "theta")
+```
+
+    ##      id.(Intercept) id.week.(Intercept)             id.week 
+    ##           0.9760823          -0.1175194           0.3951992
+
+``` r
+t(chol(VarCorr(m4)$id))[lower.tri(diag(2), diag = TRUE)] / sigma(m4)
+```
+
+    ## [1]  0.9760823 -0.1175194  0.3951992
+
+# Power simulation
+
+## Setup
+
+``` r
+## Study design and sample sizes
+n_week <- 6
+n_subj <- 80
+n <- n_week * n_subj
+
+dat <- data.frame(
+     id = factor(rep(seq_len(n_subj), each = n_week)),
+   week = rep(0:(n_week - 1), times = n_subj),
+  treat = factor(rep(0:1, each = n/2), labels = c("ctr", "trt"))
+)
+
+## Fixed effects and variance components
+beta <- c(22, -2, 0, -1)
+Su <- matrix(c(12, -1.5, -1.5, 2), nrow = 2)
+se <- 3.5
+```
+
+## Power
+
+``` r
+pval <- replicate(200, {
+
+  # Data generation
+  means <- model.matrix( ~ week * treat, dat) %*% beta
+  ranu <- MASS::mvrnorm(n_subj, mu = c(0, 0), Sigma = Su)
+  e <- rnorm(n_subj * n_week, mean = 0, sd = se)
+  
+  y <- means + ranu[dat$id, 1] + ranu[dat$id, 2] * dat$week + e
+  
+  # Fitting model to test H0
+  m0 <- lmer(y ~ week + treat + (1 + week | id), data = dat, REML = FALSE)
+  m1 <- lmer(y ~ week * treat + (1 + week | id), data = dat, REML = FALSE)
+  anova(m0, m1)["m1", "Pr(>Chisq)"]
+  }
+)
+
+mean(pval < 0.05)
+```
+
+    ## [1] 0.805
+
+## Parameter recovery
+
+``` r
+par <- replicate(200, {
+  
+  means <- model.matrix( ~ week * treat, dat) %*% beta
+  ranu <- MASS::mvrnorm(n_subj, mu = c(0, 0), Sigma = Su)
+  e <- rnorm(n_subj * n_week, mean = 0, sd = se)
+  
+  y <- means + ranu[dat$id, 1] + ranu[dat$id, 2] * dat$week + e
+  
+  m1 <- lmer(y ~ week * treat + (1 + week | id), data = dat, REML = FALSE)
+  list(fixef = fixef(m1), theta = getME(m1, "theta"), sigma = sigma(m1))
+  }, simplify = FALSE
+)
+
+rowMeans(sapply(par, function(x) x$fixef))
+```
+
+    ##   (Intercept)          week      treattrt week:treattrt 
+    ##   22.00069056   -1.98276116    0.03001846   -1.03298775
+
+``` r
+rowMeans(sapply(par, function(x) x$theta))
+```
+
+    ##      id.(Intercept) id.week.(Intercept)             id.week 
+    ##           0.9767601          -0.1180799           0.3763657
+
+``` r
+mean(sapply(par, function(x) x$sigma))
+```
+
+    ## [1] 3.479607
+
+``` r
+beta
+```
+
+    ## [1] 22 -2  0 -1
+
+``` r
+Lt <- chol(Su)
+t(Lt)[lower.tri(Lt, diag = TRUE)] / se
+```
+
+    ## [1]  0.9897433 -0.1237179  0.3846546
+
+``` r
+se
+```
+
+    ## [1] 3.5
+
+### References
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-ReisbyGram77" class="csl-entry">
+
+Reisby, N., L. F. Gram, P. Bech, A. Nagy, G. O. Petersen, J. Ortmann, I.
+Ibsen, et al. 1977. “Imipramine: Clinical Effects and Pharmacokinetic
+Variability.” *Psychopharmacology* 54: 263–72.
+<https://doi.org/10.1007/BF00426574>.
+
+</div>
+
+</div>
diff --git a/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-2-1.png b/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-2-1.png
new file mode 100644
index 0000000..17185ff
Binary files /dev/null and b/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-2-1.png differ
diff --git a/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-5-1.png b/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-5-1.png
new file mode 100644
index 0000000..526522a
Binary files /dev/null and b/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-5-1.png differ
diff --git a/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-7-1.png b/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-7-1.png
new file mode 100644
index 0000000..b0867c5
Binary files /dev/null and b/05_mixed1/powersim-lmm_files/figure-gfm/unnamed-chunk-7-1.png differ
diff --git a/06_mixed2/exercises-datsimlmm.md b/06_mixed2/exercises-datsimlmm.md
new file mode 100644
index 0000000..aa653f4
--- /dev/null
+++ b/06_mixed2/exercises-datsimlmm.md
@@ -0,0 +1,81 @@
+Exercises: Data simulation for crossed random-effects models
+================
+
+- Change the data simulation by Baayen, Davidson, and Bates (2008) for
+  $N = 30$ subjects instead of only 3
+- You can use the following script and adjust it accordingly
+- You can choose if you want to use model matrices or create the vectors
+  “manually”
+
+``` r
+library(lattice)
+library(lme4)
+
+#--------------- (1) Create data frame ----------------------------------------
+datsim <- expand.grid(subject = factor(c("s1" , "s2" , "s3" )),
+                      item    = factor(c("w1" , "w2" , "w3" )),
+                      soa     = factor(c("long" , "short" ))) |>
+  sort_by(~ subject)
+
+#--------------- (2) Define parameters ----------------------------------------
+beta0 <- 522.22
+beta1 <- -19
+
+sw  <- 21
+sy0 <- 24
+sy1 <- 7
+ry  <- -0.7
+se  <- 9
+
+#--------------- (3) Create vectors and simulate data -------------------------
+# Fixed effects
+b0 <- rep(beta0, 18)
+b1 <- rep(rep(c(0, beta1), each = 3), 3)
+
+# Draw random effects
+w <- rep(rnorm(3, mean = 0, sd = sw), 6)
+e <- rnorm(18, mean = 0, sd = se)
+
+# Bivariate normal distribution
+sig <- matrix(c(sy0^2, ry * sy0 * sy1, ry * sy0 * sy1, sy1^2), 2, 2)
+y01 <- MASS::mvrnorm(3, mu = c(0, 0), Sigma = sig)
+y0 <- rep(y01[,1], each = 6)
+y1 <- rep(c(0, y01[1,2],
+            0, y01[2,2],
+            0, y01[3,2]), each = 3)
+
+datsim$rt <- b0 + b1 + w + y0 + y1 + e
+
+#--------------- (4) Simulate data using model matrices -----------------------
+X <- model.matrix( ~ soa, datsim)
+Z <- model.matrix( ~ 0 + item + subject + subject:soa, datsim,
+                  contrasts.arg = list(subject = contrasts(datsim$subject,
+                                                        contrasts = FALSE)))
+
+# Fixed effects
+beta <- c(beta0, beta1)
+# Random effects
+u <- c(w = unique(w),
+       y0 = y01[,1],
+       y1 = y01[,2])
+
+datsim$rt2 <- X %*% beta + Z %*% u + e
+
+#--------------- (5) Visualize simulated data ---------------------------------
+xyplot(rt ~ soa | subject, datsim, group = item, type = "b", layout = c(3, 1))
+```
+
+### Reference
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-Baayen08" class="csl-entry">
+
+Baayen, R. H., D. J. Davidson, and D. M. Bates. 2008. “Mixed-Effects
+Modeling with Crossed Random Effects for Subjects and Items.” *Journal
+of Memory and Language* 59 (4): 390–412.
+<https://doi.org/10.1016/j.jml.2007.12.005>.
+
+</div>
+
+</div>
diff --git a/06_mixed2/exercises.md b/06_mixed2/exercises.md
new file mode 100644
index 0000000..d3fd939
--- /dev/null
+++ b/06_mixed2/exercises.md
@@ -0,0 +1,67 @@
+Exercise: Power simulation for LMM
+================
+
+## Physical healing as a function of perceived time
+
+Aungle and Langer (2023) investigate how perceived time influences
+physical healing
+
+- They used cupping to induce bruises on 33 subjects, then took a
+  picture, waited for 28 min and took another picture
+- Subjective time was manipulated to feel like 14, 28, or 56 min
+- The pre and post pictures were presented to 25 raters who rated the
+  amount of healing on a 10-point-scale with 0 = not at all healed, 5 =
+  somewhat healed, 10 = completely healed
+- Subjects participated in all three conditions over a two week period
+
+Data: [healing.RData](../data/healing.RData)
+
+``` r
+load("../data/healing.RData")
+
+str(dat)
+
+# Subject ID
+dat$Subject <- factor(dat$Subject)
+# Rater ID
+dat$ResponseId <- factor(dat$ResponseId)
+```
+
+1.  Visualize the data.
+
+    - Aggregate the data over Raters and plot the data for each subject
+      using `lattice::xyplot()`
+    - Aggregate the data over Subjects and plot one panel for each rater
+    - How would you choose the random effects for a model testing
+      healing over the three conditions
+
+2.  Fit the model you think fits the experimental design best
+
+3.  Test the effects of condition
+
+4.  Run a power simulation for a replication study:
+
+    - Set up a data frame containing the study design and sample size.
+
+    - Specify the minimum relevant effects.
+
+    - Set the fixed effects and variance components to plausible values.
+
+    - How many participants are required to detect the specified effect
+      with a power of 80%?
+
+    - Recover the parameters of the model for one simulated data set.
+
+### Reference
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-Aungle23" class="csl-entry">
+
+Aungle, P., and E. Langer. 2023. “Physical Healing as a Function of
+Perceived Time.” *Scientific Reports* 13 (1): 22432.
+<https://doi.org/10.1038/s41598-023-50009-3>.
+
+</div>
+
+</div>
diff --git a/06_mixed2/intro-datsimlmm.pdf b/06_mixed2/intro-datsimlmm.pdf
new file mode 100644
index 0000000..c7aef9c
Binary files /dev/null and b/06_mixed2/intro-datsimlmm.pdf differ
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..93b92a2
--- /dev/null
+++ b/README.md
@@ -0,0 +1,51 @@
+Power simulations
+================
+Nora Wickelmaier
+January, 21-22, 2026
+
+## Schedule
+
+| Day | Time        | Topic                                                                                                              | Exercises                                                                                                                         |
+|:----|:------------|:-------------------------------------------------------------------------------------------------------------------|:----------------------------------------------------------------------------------------------------------------------------------|
+| Wed | 9:00–10:30  | [Simulation-based power analysis](01_intro/powersim-intro.pdf)                                                     | [Binomial test power simulation](01_intro/exercises.md)                                                                           |
+|     | 10:45–12:15 | [Power (curves) for t-tests](02_ttest/powersim-ttest.md)                                                           | [Temporal value asymmetry](02_ttest/exercises-ttest1.md), [Directed reading activities](02_ttest/exercises-ttest2.md)             |
+|     | 13:00–14:00 | [Power for ANOVAs](03_anova/powersim-anova.md)                                                                     | [Testing the interaction in two-by-two ANOVA](03_anova/exercises-anova.md)                                                        |
+| Thu | 9:00–10:30  | [Power for baseline/follow-up measurements](04_ancova/intro-ancova.pdf)                                            | [Shoulder pain and acupuncture](04_ancova/exercises-ancova1.md), [MASS anorexia data](04_ancova/exercises-ancova2.md)             |
+|     | 10:45–12:15 | [Introduction to LMMs](05_mixed1/intro-lmm.pdf), [Power for longitudinal data analysis](05_mixed1/powersim-lmm.md) | [Risperidone vs. haloperidol and schizophrenia](05_mixed1/exercises-pwrlmm.md)                                                    |
+|     | 13:00–14:00 | [Data simulation for crossed random-effects models](06_mixed2/intro-datsimlmm.pdf)                                 | [Data simulation for crossed random-effects models](06_mixed2/exercises-datsimlmm.md), [Physical healing](06_mixed2/exercises.md) |
+
+## Content
+
+- Power and the significance filter
+- Simulation-based power analysis with R
+- Drawing power curves
+- Power for t-tests, ANOVA, ANCOVA, and mixed-effects models
+
+## Prerequisites
+
+- Introductory knowledge of statistics
+- Basic knowledge of R
+
+## Software
+
+Participants will need to have installed:
+
+- A current R version (<https://CRAN.R-project.org/>)
+- An IDE for R (like RStudio or VSCode) or a text editor with syntax
+  highlighting (like Vim or Notepad++)
+- Additional R package:
+  - lme4: <https://CRAN.R-project.org/package=lme4>
+
+## Literature
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-Wickelmaier22" class="csl-entry">
+
+Wickelmaier, F. 2022. “Simulating the Power of Statistical Tests: A
+Collection of R Examples.” *ArXiv*.
+<https://doi.org/10.48550/arXiv.2110.09836>.
+
+</div>
+
+</div>
diff --git a/data/healing.RData b/data/healing.RData
new file mode 100644
index 0000000..76af9de
Binary files /dev/null and b/data/healing.RData differ
diff --git a/data/kleinhenz.txt b/data/kleinhenz.txt
new file mode 100644
index 0000000..31c97f4
--- /dev/null
+++ b/data/kleinhenz.txt
@@ -0,0 +1,57 @@
+# Reference: Kleinhenz et al., Pain 83 (1999) 235-241
+# Source: Vickers, A. J., & Altman, D. G., BMJ 323 (2001) 1123-1124 
+
+pre post  grp
+ 35   35 plac
+ 54   37 plac
+ 35   40 plac
+ 30   42 plac
+ 44   45 plac
+ 49   47 plac
+ 38   51 plac
+ 52   52 plac
+ 52   52 plac
+ 39   53 plac
+ 44   53 plac
+ 53   53 plac
+ 73   53 plac
+ 48   54 plac
+ 58   54 plac
+ 55   57 plac
+ 73   74 plac
+ 65   76 plac
+ 78   78 plac
+ 80   80 plac
+ 59   81 plac
+ 52   81 plac
+ 46   83 plac
+ 57   85 plac
+ 44   85 plac
+ 63   86 plac
+ 80   95 plac
+ 43   74  acu
+ 75   99  acu
+ 66   88  acu
+ 65   85  acu
+ 74  100  acu
+ 58   58  acu
+ 62   77  acu
+ 64  100  acu
+ 59   61  acu
+ 45   83  acu
+ 70   78  acu
+ 66  100  acu
+ 59  100  acu
+ 68   94  acu
+ 73   60  acu
+ 60   92  acu
+ 41   41  acu
+ 53   53  acu
+ 62   80  acu
+ 74   83  acu
+ 54   57  acu
+ 40   67  acu
+ 32   82  acu
+ 78   78  acu
+ 69  100  acu
+
diff --git a/data/moeller.csv b/data/moeller.csv
new file mode 100644
index 0000000..99a5695
--- /dev/null
+++ b/data/moeller.csv
@@ -0,0 +1,2010 @@
+# Data from Moeller et al., International Journal of
+#   Neuropsychopharmacology (2008), 11, 985-997
+id,treat,week,pans
+1,risp,0,85
+1,risp,1,91
+1,risp,2,86
+1,risp,3,77
+1,risp,4,69
+1,risp,5,59
+1,risp,6,69
+1,risp,7,60
+1,risp,8,64
+2,risp,0,77
+2,risp,1,69
+2,risp,2,59
+2,risp,3,52
+2,risp,4,67
+2,risp,5,59
+2,risp,6,56
+2,risp,7,56
+2,risp,8,61
+3,halo,0,66
+3,halo,1,59
+3,halo,2,66
+3,halo,3,55
+3,halo,4,66
+3,halo,5,61
+3,halo,6,48
+3,halo,7,57
+3,halo,8,57
+4,risp,0,99
+4,risp,1,86
+4,risp,2,66
+4,risp,3,58
+4,risp,4,57
+4,risp,5,57
+4,risp,6,57
+4,risp,7,57
+4,risp,8,57
+5,halo,0,64
+5,halo,1,57
+5,halo,2,46
+5,halo,3,50
+5,halo,4,52
+5,halo,5,52
+5,halo,6,52
+5,halo,7,52
+5,halo,8,52
+6,risp,0,88
+6,risp,1,81
+6,risp,2,83
+6,risp,3,71
+6,risp,4,49
+6,risp,5,72
+6,risp,6,57
+6,risp,7,74
+6,risp,8,74
+7,halo,0,112
+7,halo,1,79
+7,halo,2,93
+7,halo,3,93
+7,halo,4,93
+7,halo,5,93
+7,halo,6,93
+7,halo,7,93
+7,halo,8,93
+8,halo,0,68
+8,halo,1,30
+8,halo,2,82
+8,halo,3,74
+8,halo,4,103
+8,halo,5,103
+8,halo,6,103
+8,halo,7,103
+8,halo,8,103
+9,risp,0,105
+9,risp,1,126
+9,risp,2,126
+9,risp,3,126
+9,risp,4,126
+9,risp,5,126
+9,risp,6,126
+9,risp,7,126
+9,risp,8,126
+10,halo,0,84
+10,halo,1,56
+10,halo,2,64
+10,halo,3,55
+10,halo,4,53
+10,halo,5,44
+10,halo,6,44
+10,halo,7,44
+10,halo,8,44
+11,risp,0,88
+11,risp,1,88
+11,risp,2,93
+11,risp,3,94
+11,risp,4,97
+11,risp,5,99
+11,risp,6,99
+11,risp,7,99
+11,risp,8,99
+12,halo,0,121
+12,halo,1,109
+12,halo,2,105
+12,halo,3,82
+12,halo,4,75
+12,halo,5,75
+12,halo,6,75
+12,halo,7,73
+12,halo,8,73
+13,risp,0,111
+13,risp,2,52
+13,risp,3,45
+13,risp,4,45
+13,risp,5,45
+13,risp,6,45
+13,risp,7,45
+13,risp,8,57
+14,risp,0,118
+14,risp,1,106
+14,risp,2,106
+14,risp,3,106
+14,risp,4,106
+14,risp,5,106
+14,risp,6,106
+14,risp,7,106
+14,risp,8,106
+15,halo,0,102
+15,halo,1,94
+15,halo,2,90
+15,halo,3,89
+15,halo,4,96
+15,halo,5,96
+15,halo,6,96
+15,halo,7,96
+15,halo,8,96
+16,halo,0,96
+16,halo,1,63
+16,halo,2,50
+16,halo,3,73
+16,halo,4,81
+16,halo,5,92
+16,halo,6,92
+16,halo,7,92
+16,halo,8,92
+17,halo,0,73
+17,halo,1,41
+17,halo,2,41
+17,halo,3,41
+17,halo,4,41
+17,halo,5,41
+17,halo,6,41
+17,halo,7,41
+17,halo,8,41
+18,halo,0,84
+18,halo,1,40
+18,halo,2,45
+18,halo,3,35
+18,halo,4,34
+18,halo,5,37
+18,halo,6,32
+18,halo,7,32
+18,halo,8,32
+19,halo,0,74
+19,halo,1,65
+19,halo,2,38
+19,halo,3,31
+19,halo,4,33
+19,halo,5,33
+19,halo,6,33
+19,halo,7,33
+19,halo,8,33
+20,risp,0,64
+20,risp,1,30
+20,risp,2,36
+20,risp,3,34
+20,risp,4,39
+20,risp,5,45
+20,risp,6,41
+20,risp,7,45
+20,risp,8,41
+21,halo,0,52
+21,halo,1,38
+21,halo,2,42
+21,halo,3,41
+21,halo,4,36
+21,halo,5,36
+21,halo,6,43
+21,halo,7,41
+21,halo,8,41
+22,risp,0,94
+22,risp,1,72
+22,risp,2,73
+22,risp,3,72
+22,risp,4,69
+22,risp,5,74
+22,risp,6,81
+22,risp,7,89
+22,risp,8,87
+23,risp,0,89
+23,risp,1,35
+23,risp,2,34
+23,risp,3,52
+23,risp,4,58
+23,risp,5,54
+23,risp,6,54
+23,risp,7,54
+23,risp,8,54
+24,halo,0,72
+24,halo,1,64
+24,halo,2,63
+24,halo,3,65
+24,halo,4,71
+24,halo,5,73
+24,halo,6,62
+24,halo,7,62
+24,halo,8,62
+25,halo,0,87
+25,halo,1,95
+25,halo,2,77
+25,halo,3,88
+25,halo,4,70
+25,halo,5,57
+25,halo,6,70
+25,halo,7,71
+25,halo,8,72
+26,risp,0,119
+26,risp,1,132
+26,risp,2,114
+26,risp,3,88
+26,risp,4,83
+26,risp,5,72
+26,risp,6,65
+26,risp,7,65
+26,risp,8,69
+27,halo,0,109
+27,halo,1,93
+27,halo,2,84
+27,halo,3,90
+27,halo,4,82
+27,halo,5,91
+27,halo,6,92
+27,halo,7,92
+27,halo,8,92
+28,halo,0,106
+28,halo,1,107
+28,halo,2,94
+28,halo,3,70
+28,halo,4,70
+28,halo,5,70
+28,halo,6,70
+28,halo,7,70
+28,halo,8,70
+29,halo,0,92
+29,halo,1,90
+29,halo,2,96
+29,halo,3,96
+29,halo,4,96
+29,halo,5,96
+29,halo,6,96
+29,halo,7,96
+29,halo,8,96
+30,halo,0,77
+30,halo,1,88
+30,halo,2,63
+30,halo,3,67
+30,halo,4,67
+30,halo,5,67
+30,halo,6,67
+30,halo,7,67
+30,halo,8,67
+31,risp,0,62
+31,risp,1,47
+31,risp,2,47
+31,risp,3,47
+31,risp,4,39
+31,risp,5,41
+31,risp,6,41
+31,risp,7,41
+31,risp,8,41
+32,risp,0,59
+32,risp,1,52
+32,risp,2,51
+32,risp,3,49
+32,risp,4,43
+32,risp,5,43
+32,risp,6,39
+32,risp,7,36
+32,risp,8,36
+33,halo,0,100
+33,halo,1,99
+33,halo,2,90
+33,halo,3,68
+33,halo,4,74
+33,halo,5,64
+33,halo,6,55
+33,halo,7,45
+33,halo,8,49
+34,risp,0,99
+34,risp,2,115
+34,risp,3,72
+34,risp,4,68
+34,risp,5,68
+34,risp,6,68
+34,risp,7,56
+34,risp,8,55
+35,risp,0,97
+35,risp,1,66
+35,risp,2,64
+35,risp,3,68
+35,risp,4,53
+35,risp,5,56
+35,risp,6,56
+35,risp,7,53
+35,risp,8,54
+36,halo,0,93
+36,halo,1,71
+36,halo,2,71
+36,halo,3,69
+36,halo,4,64
+36,halo,5,65
+36,halo,6,72
+36,halo,7,72
+36,halo,8,72
+37,risp,0,98
+37,risp,1,91
+37,risp,2,97
+37,risp,3,97
+37,risp,4,97
+37,risp,5,97
+37,risp,6,97
+37,risp,7,97
+37,risp,8,97
+38,halo,0,113
+38,halo,1,105
+38,halo,2,55
+38,halo,3,42
+38,halo,4,37
+38,halo,5,38
+38,halo,6,35
+38,halo,7,35
+38,halo,8,35
+39,halo,0,74
+39,halo,1,58
+39,halo,2,65
+39,halo,3,54
+39,halo,4,49
+39,halo,5,49
+39,halo,6,49
+39,halo,7,49
+39,halo,8,49
+40,risp,0,117
+40,risp,1,113
+40,risp,2,109
+40,risp,3,96
+40,risp,4,84
+40,risp,5,83
+40,risp,6,77
+40,risp,7,70
+40,risp,8,69
+41,halo,0,90
+41,halo,1,51
+41,halo,2,49
+41,halo,3,49
+41,halo,4,74
+41,halo,5,74
+41,halo,6,74
+41,halo,7,74
+41,halo,8,74
+42,halo,0,79
+42,halo,1,89
+42,halo,2,69
+42,halo,3,66
+42,halo,4,57
+42,halo,5,57
+42,halo,6,55
+42,halo,7,59
+42,halo,8,40
+43,risp,0,94
+43,risp,1,76
+43,risp,2,72
+43,risp,3,62
+43,risp,4,50
+43,risp,5,39
+43,risp,6,37
+43,risp,7,32
+43,risp,8,32
+44,risp,0,88
+44,risp,1,85
+44,risp,2,92
+44,risp,3,78
+44,risp,4,75
+44,risp,5,64
+44,risp,6,68
+44,risp,7,70
+44,risp,8,68
+45,risp,0,74
+45,risp,1,49
+45,risp,2,48
+45,risp,3,48
+45,risp,4,40
+45,risp,5,43
+45,risp,6,49
+45,risp,7,53
+45,risp,8,47
+46,halo,0,80
+46,halo,1,73
+46,halo,2,74
+46,halo,3,60
+46,halo,4,65
+46,halo,5,49
+46,halo,6,46
+46,halo,7,49
+46,halo,8,46
+47,halo,0,82
+47,halo,1,53
+47,halo,2,48
+47,halo,3,45
+47,halo,4,56
+47,halo,5,43
+47,halo,6,39
+47,halo,7,39
+47,halo,8,39
+48,risp,0,98
+48,risp,1,95
+48,risp,2,48
+48,risp,3,40
+48,risp,4,53
+48,risp,5,50
+48,risp,6,50
+48,risp,7,56
+48,risp,8,59
+49,halo,0,78
+49,halo,1,87
+49,halo,2,77
+49,halo,3,58
+49,halo,4,44
+49,halo,5,48
+49,halo,6,45
+49,halo,7,39
+49,halo,8,41
+50,halo,0,54
+50,halo,1,44
+50,halo,2,43
+50,halo,3,54
+50,halo,4,47
+50,halo,5,45
+50,halo,6,40
+50,halo,7,50
+50,halo,8,43
+51,risp,0,93
+51,risp,1,83
+51,risp,2,79
+51,risp,3,77
+51,risp,4,55
+51,risp,5,61
+51,risp,6,67
+51,risp,7,56
+51,risp,8,61
+52,risp,0,60
+52,risp,1,63
+52,risp,2,51
+52,risp,3,48
+52,risp,4,41
+52,risp,5,38
+52,risp,6,36
+52,risp,7,37
+52,risp,8,40
+53,halo,0,89
+53,halo,1,89
+53,halo,2,62
+53,halo,3,62
+53,halo,4,62
+53,halo,5,62
+53,halo,6,75
+53,halo,7,75
+53,halo,8,75
+54,risp,0,68
+54,risp,1,53
+54,risp,2,52
+54,risp,3,38
+54,risp,4,35
+54,risp,5,37
+54,risp,6,35
+54,risp,7,37
+54,risp,8,34
+55,halo,0,88
+55,halo,1,110
+55,halo,2,119
+55,halo,3,110
+55,halo,4,115
+55,halo,5,115
+55,halo,6,100
+55,halo,7,109
+55,halo,8,109
+56,risp,0,42
+56,risp,1,45
+56,risp,2,43
+56,risp,3,36
+56,risp,4,34
+56,risp,5,35
+56,risp,6,37
+56,risp,7,39
+56,risp,8,39
+57,halo,0,58
+57,halo,1,39
+57,halo,2,30
+57,halo,3,30
+57,halo,4,30
+57,halo,5,30
+57,halo,6,30
+57,halo,7,30
+57,halo,8,30
+58,risp,0,32
+58,risp,1,30
+58,risp,2,35
+58,risp,3,35
+58,risp,4,35
+58,risp,5,35
+58,risp,6,35
+58,risp,7,35
+58,risp,8,35
+59,halo,0,41
+59,halo,1,38
+59,halo,2,32
+59,halo,3,34
+59,halo,4,32
+59,halo,5,30
+59,halo,6,30
+59,halo,7,30
+59,halo,8,30
+60,risp,0,40
+60,risp,1,69
+60,risp,2,32
+60,risp,3,30
+60,risp,4,30
+60,risp,5,31
+60,risp,6,31
+60,risp,7,31
+60,risp,8,31
+61,halo,0,73
+61,halo,1,64
+61,halo,2,62
+61,halo,3,63
+61,halo,4,61
+61,halo,5,63
+61,halo,6,62
+61,halo,7,62
+61,halo,8,62
+62,halo,0,85
+62,halo,1,72
+62,halo,2,61
+62,halo,3,40
+62,halo,4,91
+62,halo,5,99
+62,halo,6,99
+62,halo,7,101
+62,halo,8,101
+63,risp,0,55
+63,risp,1,42
+63,risp,2,44
+63,risp,3,53
+63,risp,4,53
+63,risp,5,58
+63,risp,6,51
+63,risp,7,55
+63,risp,8,73
+64,risp,0,59
+64,risp,1,56
+64,risp,2,78
+64,risp,3,78
+64,risp,4,78
+64,risp,5,78
+64,risp,6,78
+64,risp,7,78
+64,risp,8,78
+65,risp,0,59
+65,risp,1,56
+65,risp,2,58
+65,risp,3,52
+65,risp,4,66
+65,risp,5,64
+65,risp,6,66
+65,risp,7,59
+65,risp,8,57
+66,halo,0,63
+66,halo,1,62
+66,halo,2,69
+66,halo,3,60
+66,halo,4,60
+66,halo,5,71
+66,halo,6,71
+66,halo,7,71
+66,halo,8,71
+67,risp,0,92
+67,risp,1,76
+67,risp,2,86
+67,risp,3,86
+67,risp,4,80
+67,risp,5,68
+67,risp,6,75
+67,risp,7,73
+67,risp,8,73
+68,halo,0,98
+68,halo,1,77
+68,halo,2,85
+68,halo,3,78
+68,halo,4,79
+68,halo,5,79
+68,halo,6,83
+68,halo,7,84
+68,halo,8,84
+69,risp,0,58
+69,risp,1,32
+69,risp,2,49
+69,risp,3,49
+69,risp,4,49
+69,risp,5,49
+69,risp,6,49
+69,risp,7,49
+69,risp,8,49
+70,risp,0,65
+70,risp,1,54
+70,risp,2,54
+70,risp,3,54
+70,risp,4,54
+70,risp,5,54
+70,risp,6,54
+70,risp,7,54
+70,risp,8,54
+71,halo,0,64
+71,halo,1,51
+71,halo,2,48
+71,halo,3,38
+71,halo,4,38
+71,halo,5,38
+71,halo,6,38
+71,halo,7,33
+71,halo,8,35
+72,halo,0,51
+72,halo,1,48
+72,halo,2,48
+72,halo,3,46
+72,halo,4,40
+72,halo,5,40
+72,halo,6,73
+72,halo,7,73
+72,halo,8,73
+73,risp,0,41
+73,risp,1,36
+73,risp,2,37
+73,risp,3,49
+73,risp,4,47
+73,risp,5,47
+73,risp,6,42
+73,risp,7,42
+73,risp,8,42
+74,risp,0,38
+74,risp,1,32
+74,risp,2,32
+74,risp,3,33
+74,risp,4,33
+74,risp,5,33
+74,risp,6,33
+74,risp,7,33
+74,risp,8,33
+75,risp,0,61
+75,risp,3,49
+75,risp,4,42
+75,risp,5,42
+75,risp,6,42
+75,risp,7,57
+75,risp,8,39
+76,risp,0,87
+76,risp,1,81
+76,risp,2,78
+76,risp,3,71
+76,risp,4,60
+76,risp,5,57
+76,risp,6,56
+76,risp,7,56
+76,risp,8,46
+77,risp,0,82
+77,risp,1,52
+77,risp,2,50
+77,risp,3,40
+77,risp,4,45
+77,risp,5,40
+77,risp,6,34
+77,risp,7,34
+77,risp,8,41
+78,halo,0,101
+78,halo,1,98
+78,halo,2,93
+78,halo,3,78
+78,halo,4,78
+78,halo,5,75
+78,halo,6,45
+78,halo,7,36
+78,halo,8,47
+79,risp,0,74
+79,risp,1,61
+79,risp,2,42
+79,risp,3,43
+79,risp,4,43
+79,risp,5,38
+79,risp,6,38
+79,risp,7,38
+79,risp,8,38
+80,halo,0,82
+80,halo,1,58
+80,halo,2,46
+80,halo,3,40
+80,halo,4,43
+80,halo,5,42
+80,halo,6,42
+80,halo,7,39
+80,halo,8,34
+81,halo,0,47
+81,halo,1,45
+81,halo,2,46
+81,halo,3,50
+81,halo,4,43
+81,halo,5,58
+81,halo,6,47
+81,halo,7,45
+81,halo,8,39
+82,risp,0,72
+82,risp,1,43
+82,risp,2,42
+82,risp,3,48
+82,risp,4,48
+82,risp,5,41
+82,risp,6,45
+82,risp,7,44
+82,risp,8,44
+83,halo,0,67
+83,halo,1,52
+83,halo,2,51
+83,halo,3,42
+83,halo,4,37
+83,halo,5,55
+83,halo,6,56
+83,halo,7,56
+83,halo,8,56
+84,risp,0,82
+84,risp,1,58
+84,risp,2,53
+84,risp,3,42
+84,risp,4,36
+84,risp,5,39
+84,risp,6,39
+84,risp,7,39
+84,risp,8,39
+85,halo,0,60
+85,halo,1,43
+85,halo,2,48
+85,halo,3,46
+85,halo,4,47
+85,halo,5,42
+85,halo,6,42
+85,halo,7,42
+85,halo,8,42
+86,risp,0,54
+86,risp,1,47
+86,risp,2,38
+86,risp,3,40
+86,risp,4,36
+86,risp,5,34
+86,risp,6,34
+86,risp,7,34
+86,risp,8,34
+87,halo,0,57
+87,halo,1,52
+87,halo,2,42
+87,halo,3,58
+87,halo,4,58
+87,halo,5,58
+87,halo,6,58
+87,halo,7,58
+87,halo,8,58
+88,risp,0,69
+88,risp,1,51
+88,risp,2,42
+88,risp,3,36
+88,risp,4,38
+88,risp,5,40
+88,risp,6,40
+88,risp,7,40
+88,risp,8,40
+89,risp,0,64
+89,risp,1,61
+89,risp,2,67
+89,risp,3,74
+89,risp,4,70
+89,risp,5,77
+89,risp,6,70
+89,risp,7,78
+89,risp,8,71
+90,risp,0,92
+90,risp,2,82
+90,risp,3,82
+90,risp,4,92
+90,risp,5,83
+90,risp,6,83
+90,risp,7,83
+90,risp,8,83
+91,halo,0,84
+91,halo,1,66
+91,halo,2,83
+91,halo,3,69
+91,halo,4,69
+91,halo,5,70
+91,halo,6,51
+91,halo,7,54
+91,halo,8,55
+92,risp,0,63
+92,risp,1,40
+92,risp,2,36
+92,risp,3,36
+92,risp,4,37
+92,risp,5,37
+92,risp,6,37
+92,risp,7,37
+92,risp,8,37
+93,risp,0,48
+93,risp,1,46
+93,risp,2,50
+93,risp,3,40
+93,risp,4,40
+93,risp,5,40
+93,risp,6,39
+93,risp,7,38
+93,risp,8,42
+94,halo,0,67
+94,halo,1,65
+94,halo,2,76
+94,halo,3,63
+94,halo,4,53
+94,halo,5,51
+94,halo,6,51
+94,halo,7,51
+94,halo,8,51
+95,halo,0,66
+95,halo,1,44
+95,halo,2,46
+95,halo,3,50
+95,halo,4,58
+95,halo,5,60
+95,halo,6,60
+95,halo,7,60
+95,halo,8,60
+96,risp,0,45
+96,risp,1,47
+96,risp,2,43
+96,risp,3,41
+96,risp,4,52
+96,risp,5,42
+96,risp,6,38
+96,risp,7,38
+96,risp,8,43
+97,risp,0,68
+97,risp,1,63
+97,risp,2,54
+97,risp,3,54
+97,risp,4,65
+97,risp,5,68
+97,risp,6,49
+97,risp,7,48
+97,risp,8,48
+98,halo,0,92
+98,halo,1,86
+98,halo,2,71
+98,halo,3,53
+98,halo,4,60
+98,halo,5,66
+98,halo,6,66
+98,halo,7,66
+98,halo,8,66
+99,halo,0,52
+99,halo,1,49
+99,halo,2,42
+99,halo,3,52
+99,halo,4,48
+99,halo,5,57
+99,halo,6,58
+99,halo,7,58
+99,halo,8,58
+100,risp,0,91
+100,risp,1,74
+100,risp,2,60
+100,risp,3,50
+100,risp,4,55
+100,risp,5,49
+100,risp,6,52
+100,risp,7,53
+100,risp,8,50
+101,halo,0,50
+101,halo,1,44
+101,halo,2,44
+101,halo,3,41
+101,halo,4,36
+101,halo,5,32
+101,halo,6,32
+101,halo,7,42
+101,halo,8,34
+102,halo,0,70
+102,halo,1,44
+102,halo,2,57
+102,halo,3,57
+102,halo,4,57
+102,halo,5,62
+102,halo,6,60
+102,halo,7,59
+102,halo,8,54
+103,halo,0,56
+103,halo,1,56
+103,halo,2,54
+103,halo,3,57
+103,halo,4,53
+103,halo,5,47
+103,halo,6,50
+103,halo,7,50
+103,halo,8,50
+104,halo,0,78
+104,halo,1,40
+104,halo,2,41
+104,halo,3,37
+104,halo,4,33
+104,halo,5,39
+104,halo,6,44
+104,halo,7,42
+104,halo,8,45
+105,risp,0,66
+105,risp,1,44
+105,risp,2,43
+105,risp,3,33
+105,risp,4,30
+105,risp,5,30
+105,risp,6,30
+105,risp,7,30
+105,risp,8,30
+106,risp,0,70
+106,risp,1,53
+106,risp,2,56
+106,risp,3,50
+106,risp,4,37
+106,risp,5,37
+106,risp,6,40
+106,risp,7,38
+106,risp,8,43
+107,halo,0,67
+107,halo,1,44
+107,halo,2,59
+107,halo,3,58
+107,halo,4,43
+107,halo,5,44
+107,halo,6,37
+107,halo,7,34
+107,halo,8,36
+108,risp,0,71
+108,risp,1,80
+108,risp,2,66
+108,risp,3,52
+108,risp,4,54
+108,risp,5,54
+108,risp,6,54
+108,risp,7,54
+108,risp,8,54
+109,halo,0,79
+109,halo,1,40
+109,halo,2,34
+109,halo,3,36
+109,halo,4,32
+109,halo,5,32
+109,halo,6,32
+109,halo,7,32
+109,halo,8,30
+110,halo,0,54
+110,halo,1,45
+110,halo,2,42
+110,halo,3,43
+110,halo,4,44
+110,halo,5,42
+110,halo,6,43
+110,halo,7,41
+110,halo,8,41
+111,risp,0,103
+111,risp,1,68
+111,risp,2,67
+111,risp,3,61
+111,risp,4,49
+111,risp,5,56
+111,risp,6,51
+111,risp,7,51
+111,risp,8,49
+112,halo,0,84
+112,halo,1,71
+112,halo,2,59
+112,halo,3,58
+112,halo,4,58
+112,halo,5,55
+112,halo,6,46
+112,halo,7,42
+112,halo,8,46
+113,risp,0,68
+113,risp,1,56
+113,risp,2,40
+113,risp,3,45
+113,risp,4,44
+113,risp,5,44
+113,risp,6,42
+113,risp,7,42
+113,risp,8,44
+114,halo,0,47
+114,halo,1,58
+114,halo,2,47
+114,halo,3,43
+114,halo,4,43
+114,halo,5,46
+114,halo,6,42
+114,halo,7,37
+114,halo,8,40
+115,risp,0,67
+115,risp,1,68
+115,risp,2,60
+115,risp,3,59
+115,risp,4,47
+115,risp,5,50
+115,risp,6,49
+115,risp,7,57
+115,risp,8,53
+116,risp,0,68
+116,risp,1,52
+116,risp,2,49
+116,risp,3,46
+116,risp,4,40
+116,risp,5,40
+116,risp,6,40
+116,risp,7,40
+116,risp,8,32
+117,halo,0,44
+117,halo,1,42
+117,halo,2,43
+117,halo,3,54
+117,halo,4,48
+117,halo,5,42
+117,halo,6,42
+117,halo,7,40
+117,halo,8,38
+118,risp,0,76
+118,risp,1,37
+118,risp,2,36
+118,risp,3,34
+118,risp,4,34
+118,risp,5,32
+118,risp,6,32
+118,risp,7,32
+118,risp,8,30
+119,risp,0,81
+119,risp,1,74
+119,risp,2,69
+119,risp,3,69
+119,risp,4,69
+119,risp,5,69
+119,risp,6,69
+119,risp,7,69
+119,risp,8,69
+120,halo,0,48
+120,halo,1,34
+120,halo,2,30
+120,halo,3,30
+120,halo,4,30
+120,halo,5,30
+120,halo,6,30
+120,halo,7,30
+120,halo,8,30
+121,halo,0,56
+121,halo,1,47
+121,halo,2,36
+121,halo,3,32
+121,halo,4,38
+121,halo,5,31
+121,halo,6,31
+121,halo,7,33
+121,halo,8,32
+122,risp,0,76
+122,risp,1,54
+122,risp,2,71
+122,risp,3,51
+122,risp,4,37
+122,risp,5,36
+122,risp,6,41
+122,risp,7,40
+122,risp,8,40
+123,halo,0,67
+123,halo,1,72
+123,halo,2,46
+123,halo,3,39
+123,halo,4,36
+123,halo,5,32
+123,halo,6,34
+123,halo,7,35
+123,halo,8,37
+124,halo,0,59
+124,halo,1,42
+124,halo,2,38
+124,halo,3,38
+124,halo,4,34
+124,halo,5,32
+124,halo,6,32
+124,halo,7,30
+124,halo,8,30
+125,risp,0,62
+125,risp,1,56
+125,risp,2,53
+125,risp,3,42
+125,risp,4,42
+125,risp,5,41
+125,risp,6,50
+125,risp,7,44
+125,risp,8,44
+126,halo,0,70
+126,halo,1,43
+126,halo,2,54
+126,halo,3,36
+126,halo,4,40
+126,halo,5,42
+126,halo,6,41
+126,halo,7,41
+126,halo,8,56
+127,risp,0,60
+127,risp,1,54
+127,risp,2,51
+127,risp,3,59
+127,risp,4,59
+127,risp,5,59
+127,risp,6,59
+127,risp,7,59
+127,risp,8,59
+128,halo,0,47
+128,halo,1,40
+128,halo,2,36
+128,halo,3,37
+128,halo,4,34
+128,halo,5,34
+128,halo,6,34
+128,halo,7,34
+128,halo,8,34
+129,risp,0,63
+129,risp,1,56
+129,risp,2,43
+129,risp,3,44
+129,risp,4,44
+129,risp,5,44
+129,risp,6,44
+129,risp,7,44
+129,risp,8,44
+130,risp,0,56
+130,risp,1,47
+130,risp,2,40
+130,risp,3,48
+130,risp,4,48
+130,risp,5,44
+130,risp,6,52
+130,risp,7,56
+130,risp,8,43
+131,risp,0,67
+131,risp,1,66
+131,risp,2,58
+131,risp,3,45
+131,risp,4,52
+131,risp,5,47
+131,risp,6,46
+131,risp,7,38
+131,risp,8,41
+132,risp,0,50
+132,risp,1,50
+132,risp,2,46
+132,risp,3,65
+132,risp,4,61
+132,risp,5,55
+132,risp,6,62
+132,risp,7,64
+132,risp,8,63
+133,halo,0,76
+133,halo,1,59
+133,halo,2,54
+133,halo,3,44
+133,halo,4,44
+133,halo,5,44
+133,halo,6,44
+133,halo,7,44
+133,halo,8,44
+134,risp,0,57
+134,risp,1,46
+134,risp,2,49
+134,risp,3,56
+134,risp,4,68
+134,risp,5,61
+134,risp,6,53
+134,risp,7,48
+134,risp,8,50
+135,risp,0,65
+135,risp,1,67
+135,risp,2,44
+135,risp,3,55
+135,risp,4,44
+135,risp,5,44
+135,risp,6,30
+135,risp,7,50
+135,risp,8,43
+136,halo,0,66
+136,halo,1,54
+136,halo,2,44
+136,halo,3,42
+136,halo,4,42
+136,halo,5,42
+136,halo,6,60
+136,halo,7,54
+136,halo,8,50
+137,halo,0,45
+137,halo,1,38
+137,halo,2,35
+137,halo,3,31
+137,halo,4,31
+137,halo,5,38
+137,halo,6,38
+137,halo,7,38
+137,halo,8,38
+138,risp,0,59
+138,risp,1,50
+138,risp,2,44
+138,risp,3,50
+138,risp,4,66
+138,risp,5,66
+138,risp,6,66
+138,risp,7,66
+138,risp,8,66
+139,halo,0,82
+139,halo,1,82
+139,halo,2,80
+139,halo,3,68
+139,halo,4,64
+139,halo,5,52
+139,halo,6,46
+139,halo,7,46
+139,halo,8,46
+140,risp,0,59
+140,risp,1,44
+140,risp,2,35
+140,risp,3,41
+140,risp,4,53
+140,risp,5,63
+140,risp,6,67
+140,risp,7,66
+140,risp,8,66
+141,halo,0,75
+141,halo,1,59
+141,halo,2,59
+141,halo,3,61
+141,halo,4,54
+141,halo,5,40
+141,halo,6,40
+141,halo,7,35
+141,halo,8,43
+142,risp,0,44
+142,risp,1,45
+142,risp,2,45
+142,risp,3,45
+142,risp,4,45
+142,risp,5,45
+142,risp,6,45
+142,risp,7,45
+142,risp,8,45
+143,halo,0,62
+143,halo,1,56
+143,halo,2,60
+143,halo,3,57
+143,halo,4,56
+143,halo,5,58
+143,halo,6,56
+143,halo,7,51
+143,halo,8,58
+144,risp,0,63
+144,risp,1,56
+144,risp,2,53
+144,risp,3,58
+144,risp,4,52
+144,risp,5,51
+144,risp,6,69
+144,risp,7,49
+144,risp,8,55
+145,halo,0,48
+145,halo,1,46
+145,halo,2,42
+145,halo,3,46
+145,halo,4,44
+145,halo,5,44
+145,halo,6,44
+145,halo,7,44
+145,halo,8,44
+146,risp,0,45
+146,risp,1,45
+146,risp,2,31
+146,risp,3,30
+146,risp,4,30
+146,risp,5,30
+146,risp,6,30
+146,risp,7,30
+146,risp,8,30
+147,halo,0,72
+147,halo,1,33
+147,halo,2,33
+147,halo,3,33
+147,halo,4,33
+147,halo,5,33
+147,halo,6,33
+147,halo,7,33
+147,halo,8,33
+148,risp,0,40
+148,risp,1,43
+148,risp,2,42
+148,risp,3,40
+148,risp,4,39
+148,risp,5,39
+148,risp,6,39
+148,risp,7,39
+148,risp,8,39
+149,halo,0,62
+149,halo,1,54
+149,halo,2,53
+149,halo,3,51
+149,halo,4,46
+149,halo,5,46
+149,halo,6,46
+149,halo,7,46
+149,halo,8,46
+150,halo,0,58
+150,halo,1,30
+150,halo,2,40
+150,halo,3,40
+150,halo,4,40
+150,halo,5,40
+150,halo,6,40
+150,halo,7,40
+150,halo,8,40
+151,risp,0,70
+151,risp,1,70
+151,risp,2,70
+151,risp,3,70
+151,risp,4,70
+151,risp,5,70
+151,risp,6,38
+151,risp,7,38
+151,risp,8,41
+152,halo,0,46
+152,halo,1,43
+152,halo,2,42
+152,halo,3,50
+152,halo,4,50
+152,halo,5,45
+152,halo,6,45
+152,halo,7,45
+152,halo,8,45
+153,risp,0,69
+153,risp,1,41
+153,risp,2,38
+153,risp,3,33
+153,risp,4,45
+153,risp,5,38
+153,risp,6,36
+153,risp,7,35
+153,risp,8,37
+154,halo,0,46
+154,halo,1,30
+154,halo,2,34
+154,halo,3,30
+154,halo,4,30
+154,halo,5,34
+154,halo,6,35
+154,halo,7,39
+154,halo,8,44
+155,halo,0,69
+155,halo,1,34
+155,halo,2,34
+155,halo,3,33
+155,halo,4,30
+155,halo,5,30
+155,halo,6,31
+155,halo,7,34
+155,halo,8,34
+156,risp,0,60
+156,risp,1,56
+156,risp,2,51
+156,risp,3,42
+156,risp,4,63
+156,risp,5,48
+156,risp,6,48
+156,risp,7,35
+156,risp,8,37
+157,halo,0,82
+157,halo,1,45
+157,halo,2,45
+157,halo,3,47
+157,halo,4,41
+157,halo,5,43
+157,halo,6,43
+157,halo,7,58
+157,halo,8,55
+158,risp,0,64
+158,risp,1,41
+158,risp,2,48
+158,risp,3,46
+158,risp,4,46
+158,risp,5,47
+158,risp,6,47
+158,risp,7,47
+158,risp,8,47
+159,risp,0,40
+159,risp,1,38
+159,risp,2,43
+159,risp,3,43
+159,risp,4,43
+159,risp,5,43
+159,risp,6,43
+159,risp,7,43
+159,risp,8,43
+160,risp,0,77
+160,risp,1,57
+160,risp,2,45
+160,risp,3,45
+160,risp,4,45
+160,risp,5,45
+160,risp,6,45
+160,risp,7,45
+160,risp,8,45
+161,halo,0,37
+161,halo,1,40
+161,halo,2,35
+161,halo,3,35
+161,halo,4,35
+161,halo,5,35
+161,halo,6,35
+161,halo,7,35
+161,halo,8,35
+162,risp,0,57
+162,risp,3,50
+162,risp,4,50
+162,risp,5,50
+162,risp,6,56
+162,risp,7,53
+162,risp,8,53
+163,risp,0,38
+163,risp,2,31
+163,risp,3,31
+163,risp,4,31
+163,risp,5,31
+163,risp,6,31
+163,risp,7,31
+163,risp,8,31
+164,risp,0,51
+164,risp,1,42
+164,risp,2,35
+164,risp,3,37
+164,risp,4,34
+164,risp,5,34
+164,risp,6,34
+164,risp,7,34
+164,risp,8,34
+165,risp,0,84
+165,risp,1,55
+165,risp,2,46
+165,risp,3,46
+165,risp,4,37
+165,risp,5,37
+165,risp,6,37
+165,risp,7,37
+165,risp,8,37
+166,halo,0,83
+166,halo,1,50
+166,halo,2,41
+166,halo,3,47
+166,halo,4,37
+166,halo,5,37
+166,halo,6,37
+166,halo,7,37
+166,halo,8,37
+167,risp,0,72
+167,risp,1,62
+167,risp,2,59
+167,risp,3,58
+167,risp,4,58
+167,risp,5,58
+167,risp,6,58
+167,risp,7,58
+167,risp,8,58
+168,halo,0,59
+168,halo,1,67
+168,halo,2,67
+168,halo,3,67
+168,halo,4,67
+168,halo,5,67
+168,halo,6,67
+168,halo,7,67
+168,halo,8,67
+169,halo,0,62
+169,halo,1,57
+169,halo,2,59
+169,halo,3,60
+169,halo,4,63
+169,halo,5,63
+169,halo,6,63
+169,halo,7,63
+169,halo,8,63
+170,halo,0,57
+170,halo,1,50
+170,halo,2,36
+170,halo,3,31
+170,halo,4,31
+170,halo,5,37
+170,halo,6,32
+170,halo,7,35
+170,halo,8,31
+171,risp,0,62
+171,risp,1,81
+171,risp,2,58
+171,risp,3,58
+171,risp,4,59
+171,risp,5,68
+171,risp,6,62
+171,risp,7,42
+171,risp,8,45
+172,halo,0,85
+172,halo,1,41
+172,halo,2,57
+172,halo,3,53
+172,halo,4,45
+172,halo,5,45
+172,halo,6,34
+172,halo,7,34
+172,halo,8,33
+173,risp,0,60
+173,risp,1,49
+173,risp,2,42
+173,risp,3,36
+173,risp,4,32
+173,risp,5,54
+173,risp,6,45
+173,risp,7,50
+173,risp,8,50
+174,halo,0,85
+174,halo,1,81
+174,halo,2,60
+174,halo,3,59
+174,halo,4,51
+174,halo,5,52
+174,halo,6,76
+174,halo,7,77
+174,halo,8,80
+175,halo,0,50
+175,halo,1,84
+175,halo,2,71
+175,halo,3,69
+175,halo,4,83
+175,halo,5,78
+175,halo,6,80
+175,halo,7,80
+175,halo,8,80
+176,risp,0,49
+176,risp,1,51
+176,risp,2,35
+176,risp,3,35
+176,risp,4,54
+176,risp,5,52
+176,risp,6,60
+176,risp,7,64
+176,risp,8,58
+177,risp,0,92
+177,risp,1,78
+177,risp,2,78
+177,risp,3,85
+177,risp,4,92
+177,risp,5,86
+177,risp,6,86
+177,risp,7,86
+177,risp,8,86
+178,halo,0,119
+178,halo,1,115
+178,halo,2,94
+178,halo,3,122
+178,halo,4,129
+178,halo,5,129
+178,halo,6,129
+178,halo,7,129
+178,halo,8,129
+179,risp,0,104
+179,risp,1,113
+179,risp,2,96
+179,risp,3,84
+179,risp,4,77
+179,risp,5,74
+179,risp,6,67
+179,risp,7,72
+179,risp,8,70
+180,risp,0,75
+180,risp,1,72
+180,risp,2,72
+180,risp,3,72
+180,risp,4,72
+180,risp,5,72
+180,risp,6,72
+180,risp,7,72
+180,risp,8,72
+181,halo,0,56
+181,halo,1,38
+181,halo,2,36
+181,halo,3,33
+181,halo,4,44
+181,halo,5,47
+181,halo,6,48
+181,halo,7,45
+181,halo,8,45
+182,risp,0,61
+182,risp,1,38
+182,risp,2,39
+182,risp,3,39
+182,risp,4,39
+182,risp,5,39
+182,risp,6,39
+182,risp,7,39
+182,risp,8,39
+183,risp,0,81
+183,risp,1,31
+183,risp,2,43
+183,risp,3,30
+183,risp,4,30
+183,risp,5,30
+183,risp,6,32
+183,risp,7,31
+183,risp,8,44
+184,risp,0,106
+184,risp,1,58
+184,risp,2,55
+184,risp,3,38
+184,risp,4,58
+184,risp,5,52
+184,risp,6,96
+184,risp,7,112
+184,risp,8,96
+185,halo,0,100
+185,halo,1,82
+185,halo,2,86
+185,halo,3,86
+185,halo,4,86
+185,halo,5,86
+185,halo,6,86
+185,halo,7,86
+185,halo,8,86
+186,risp,0,123
+186,risp,1,101
+186,risp,2,99
+186,risp,3,77
+186,risp,4,76
+186,risp,5,76
+186,risp,6,75
+186,risp,7,93
+186,risp,8,70
+187,risp,0,115
+187,risp,1,87
+187,risp,2,91
+187,risp,3,91
+187,risp,4,99
+187,risp,5,109
+187,risp,6,61
+187,risp,7,68
+187,risp,8,69
+188,risp,0,113
+188,risp,1,89
+188,risp,2,82
+188,risp,3,70
+188,risp,4,75
+188,risp,5,67
+188,risp,6,67
+188,risp,7,67
+188,risp,8,67
+189,halo,0,120
+189,halo,1,101
+189,halo,2,85
+189,halo,3,76
+189,halo,4,77
+189,halo,5,53
+189,halo,6,66
+189,halo,7,67
+189,halo,8,70
+190,halo,0,131
+190,halo,1,102
+190,halo,2,83
+190,halo,3,86
+190,halo,4,49
+190,halo,5,53
+190,halo,6,57
+190,halo,7,77
+190,halo,8,77
+191,halo,0,133
+191,halo,1,67
+191,halo,2,68
+191,halo,3,68
+191,halo,4,56
+191,halo,5,47
+191,halo,6,46
+191,halo,7,46
+191,halo,8,39
+192,risp,0,116
+192,risp,1,95
+192,risp,2,94
+192,risp,3,94
+192,risp,4,66
+192,risp,5,85
+192,risp,6,94
+192,risp,7,94
+192,risp,8,92
+193,halo,0,152
+193,halo,2,122
+193,halo,3,120
+193,halo,4,99
+193,halo,5,99
+193,halo,6,82
+193,halo,7,104
+193,halo,8,104
+194,halo,0,75
+194,halo,1,68
+194,halo,2,48
+194,halo,3,47
+194,halo,4,46
+194,halo,5,46
+194,halo,6,46
+194,halo,7,46
+194,halo,8,46
+195,risp,0,72
+195,risp,1,64
+195,risp,2,56
+195,risp,3,49
+195,risp,4,53
+195,risp,5,52
+195,risp,6,53
+195,risp,7,54
+195,risp,8,49
+196,halo,0,88
+196,halo,1,50
+196,halo,2,50
+196,halo,3,50
+196,halo,4,50
+196,halo,5,50
+196,halo,6,50
+196,halo,7,50
+196,halo,8,50
+197,halo,0,87
+197,halo,1,35
+197,halo,2,35
+197,halo,3,35
+197,halo,4,78
+197,halo,5,62
+197,halo,6,47
+197,halo,7,47
+197,halo,8,47
+198,risp,0,125
+198,risp,1,125
+198,risp,2,125
+198,risp,3,125
+198,risp,4,125
+198,risp,5,125
+198,risp,6,125
+198,risp,7,125
+198,risp,8,125
+199,halo,0,102
+199,halo,1,106
+199,halo,2,90
+199,halo,3,69
+199,halo,4,63
+199,halo,5,46
+199,halo,6,46
+199,halo,7,46
+199,halo,8,46
+200,risp,0,66
+200,risp,1,64
+200,risp,2,74
+200,risp,3,74
+200,risp,4,81
+200,risp,5,64
+200,risp,6,44
+200,risp,7,46
+200,risp,8,54
+201,halo,0,73
+201,halo,1,74
+201,halo,2,76
+201,halo,3,71
+201,halo,4,63
+201,halo,5,56
+201,halo,6,56
+201,halo,7,56
+201,halo,8,56
+202,halo,0,108
+202,halo,1,102
+202,halo,2,112
+202,halo,3,118
+202,halo,4,85
+202,halo,5,107
+202,halo,6,107
+202,halo,7,107
+202,halo,8,107
+203,halo,0,86
+203,halo,1,62
+203,halo,2,52
+203,halo,3,49
+203,halo,4,41
+203,halo,5,40
+203,halo,6,31
+203,halo,7,41
+203,halo,8,38
+204,halo,0,85
+204,halo,1,83
+204,halo,2,93
+204,halo,3,89
+204,halo,4,81
+204,halo,5,56
+204,halo,6,37
+204,halo,7,40
+204,halo,8,45
+205,halo,0,66
+205,halo,1,49
+205,halo,2,46
+205,halo,3,44
+205,halo,4,51
+205,halo,5,57
+205,halo,6,70
+205,halo,7,70
+205,halo,8,54
+206,halo,0,97
+206,halo,1,70
+206,halo,2,73
+206,halo,3,76
+206,halo,4,49
+206,halo,5,58
+206,halo,6,62
+206,halo,7,62
+206,halo,8,39
+207,risp,0,125
+207,risp,1,121
+207,risp,2,111
+207,risp,3,94
+207,risp,4,49
+207,risp,5,43
+207,risp,6,43
+207,risp,7,53
+207,risp,8,53
+208,risp,0,75
+208,risp,1,69
+208,risp,2,48
+208,risp,3,45
+208,risp,4,52
+208,risp,5,52
+208,risp,6,50
+208,risp,7,58
+208,risp,8,58
+209,risp,0,69
+209,risp,1,47
+209,risp,2,49
+209,risp,3,42
+209,risp,4,37
+209,risp,5,36
+209,risp,6,34
+209,risp,7,31
+209,risp,8,32
+210,halo,0,131
+210,halo,1,108
+210,halo,2,106
+210,halo,3,99
+210,halo,4,96
+210,halo,5,91
+210,halo,6,89
+210,halo,7,94
+210,halo,8,101
+211,risp,0,105
+211,risp,1,82
+211,risp,2,74
+211,risp,3,66
+211,risp,4,66
+211,risp,5,66
+211,risp,6,66
+211,risp,7,66
+211,risp,8,66
+212,halo,0,108
+212,halo,1,75
+212,halo,2,56
+212,halo,3,46
+212,halo,4,36
+212,halo,5,36
+212,halo,6,36
+212,halo,7,36
+212,halo,8,36
+213,risp,0,99
+213,risp,1,92
+213,risp,2,75
+213,risp,3,60
+213,risp,4,63
+213,risp,5,63
+213,risp,6,63
+213,risp,7,63
+213,risp,8,63
+214,halo,0,105
+214,halo,1,107
+214,halo,2,107
+214,halo,3,107
+214,halo,4,107
+214,halo,5,107
+214,halo,6,107
+214,halo,7,107
+214,halo,8,107
+215,risp,0,69
+215,risp,1,55
+215,risp,2,61
+215,risp,3,41
+215,risp,4,39
+215,risp,5,58
+215,risp,6,58
+215,risp,7,58
+215,risp,8,58
+216,risp,0,140
+216,risp,1,103
+216,risp,2,87
+216,risp,3,74
+216,risp,4,74
+216,risp,5,74
+216,risp,6,74
+216,risp,7,74
+216,risp,8,74
+217,risp,0,140
+217,risp,1,117
+217,risp,2,97
+217,risp,3,80
+217,risp,4,72
+217,risp,5,72
+217,risp,6,72
+217,risp,7,72
+217,risp,8,72
+218,halo,0,125
+218,halo,1,97
+218,halo,2,89
+218,halo,3,89
+218,halo,4,89
+218,halo,5,89
+218,halo,6,89
+218,halo,7,89
+218,halo,8,89
+219,halo,0,98
+219,halo,1,77
+219,halo,2,65
+219,halo,3,66
+219,halo,4,66
+219,halo,5,66
+219,halo,6,66
+219,halo,7,66
+219,halo,8,66
+220,halo,0,119
+220,halo,1,91
+220,halo,2,91
+220,halo,3,91
+220,halo,4,91
+220,halo,5,91
+220,halo,6,91
+220,halo,7,91
+220,halo,8,91
+221,halo,0,105
+221,halo,1,94
+221,halo,2,81
+221,halo,3,61
+221,halo,4,46
+221,halo,5,44
+221,halo,6,35
+221,halo,7,32
+221,halo,8,32
+222,risp,0,65
+222,risp,1,62
+222,risp,2,60
+222,risp,3,60
+222,risp,4,56
+222,risp,5,56
+222,risp,6,59
+222,risp,7,59
+222,risp,8,59
+223,halo,0,69
+223,halo,1,46
+223,halo,2,46
+223,halo,3,67
+223,halo,4,46
+223,halo,5,53
+223,halo,6,53
+223,halo,7,53
+223,halo,8,53
+224,risp,0,96
+224,risp,1,68
+224,risp,2,68
+224,risp,3,68
+224,risp,4,68
+224,risp,5,68
+224,risp,6,68
+224,risp,7,68
+224,risp,8,68
diff --git a/data/reisby.txt b/data/reisby.txt
new file mode 100644
index 0000000..60f81d1
--- /dev/null
+++ b/data/reisby.txt
@@ -0,0 +1,399 @@
+## Data from Reisby et al., Psychopharmacology 54 (1977) 263-272
+ id hamd week  diag endweek
+101   26    0 nonen       0
+101   22    1 nonen       0
+101   18    2 nonen       0
+101    7    3 nonen       0
+101    4    4 nonen       0
+101    3    5 nonen       0
+103   33    0 nonen       0
+103   24    1 nonen       0
+103   15    2 nonen       0
+103   24    3 nonen       0
+103   15    4 nonen       0
+103   13    5 nonen       0
+104   29    0 endog       0
+104   22    1 endog       1
+104   18    2 endog       2
+104   13    3 endog       3
+104   19    4 endog       4
+104    0    5 endog       5
+105   22    0 nonen       0
+105   12    1 nonen       0
+105   16    2 nonen       0
+105   16    3 nonen       0
+105   13    4 nonen       0
+105    9    5 nonen       0
+106   21    0 endog       0
+106   25    1 endog       1
+106   23    2 endog       2
+106   18    3 endog       3
+106   20    4 endog       4
+106   NA    5 endog       5
+107   21    0 endog       0
+107   21    1 endog       1
+107   16    2 endog       2
+107   19    3 endog       3
+107   NA    4 endog       4
+107    6    5 endog       5
+108   21    0 endog       0
+108   22    1 endog       1
+108   11    2 endog       2
+108    9    3 endog       3
+108    9    4 endog       4
+108    7    5 endog       5
+113   21    0 nonen       0
+113   23    1 nonen       0
+113   19    2 nonen       0
+113   23    3 nonen       0
+113   23    4 nonen       0
+113   NA    5 nonen       0
+114   NA    0 nonen       0
+114   17    1 nonen       0
+114   11    2 nonen       0
+114   13    3 nonen       0
+114    7    4 nonen       0
+114    7    5 nonen       0
+115   NA    0 endog       0
+115   16    1 endog       1
+115   16    2 endog       2
+115   16    3 endog       3
+115   16    4 endog       4
+115   11    5 endog       5
+117   19    0 endog       0
+117   16    1 endog       1
+117   13    2 endog       2
+117   12    3 endog       3
+117    7    4 endog       4
+117    6    5 endog       5
+118   NA    0 endog       0
+118   26    1 endog       1
+118   18    2 endog       2
+118   18    3 endog       3
+118   14    4 endog       4
+118   11    5 endog       5
+120   20    0 nonen       0
+120   19    1 nonen       0
+120   17    2 nonen       0
+120   18    3 nonen       0
+120   16    4 nonen       0
+120   17    5 nonen       0
+121   20    0 nonen       0
+121   22    1 nonen       0
+121   19    2 nonen       0
+121   19    3 nonen       0
+121   12    4 nonen       0
+121   14    5 nonen       0
+123   15    0 nonen       0
+123   15    1 nonen       0
+123   15    2 nonen       0
+123   13    3 nonen       0
+123    5    4 nonen       0
+123    5    5 nonen       0
+501   29    0 endog       0
+501   30    1 endog       1
+501   26    2 endog       2
+501   22    3 endog       3
+501   19    4 endog       4
+501   24    5 endog       5
+502   21    0 endog       0
+502   22    1 endog       1
+502   13    2 endog       2
+502   11    3 endog       3
+502    2    4 endog       4
+502    1    5 endog       5
+504   19    0 nonen       0
+504   17    1 nonen       0
+504   15    2 nonen       0
+504   16    3 nonen       0
+504   12    4 nonen       0
+504   12    5 nonen       0
+505   21    0 nonen       0
+505   11    1 nonen       0
+505   18    2 nonen       0
+505    0    3 nonen       0
+505    0    4 nonen       0
+505    4    5 nonen       0
+507   27    0 endog       0
+507   26    1 endog       1
+507   26    2 endog       2
+507   25    3 endog       3
+507   24    4 endog       4
+507   19    5 endog       5
+603   28    0 nonen       0
+603   22    1 nonen       0
+603   18    2 nonen       0
+603   20    3 nonen       0
+603   11    4 nonen       0
+603   13    5 nonen       0
+604   27    0 nonen       0
+604   27    1 nonen       0
+604   13    2 nonen       0
+604    5    3 nonen       0
+604    7    4 nonen       0
+604   NA    5 nonen       0
+606   19    0 endog       0
+606   33    1 endog       1
+606   12    2 endog       2
+606   12    3 endog       3
+606    3    4 endog       4
+606    1    5 endog       5
+607   30    0 endog       0
+607   39    1 endog       1
+607   30    2 endog       2
+607   27    3 endog       3
+607   20    4 endog       4
+607    4    5 endog       5
+608   24    0 nonen       0
+608   19    1 nonen       0
+608   14    2 nonen       0
+608   12    3 nonen       0
+608    3    4 nonen       0
+608    4    5 nonen       0
+609   NA    0 endog       1
+609   25    1 endog       1
+609   22    2 endog       2
+609   14    3 endog       3
+609   15    4 endog       4
+609    2    5 endog       5
+610   34    0 endog       0
+610   NA    1 endog       1
+610   33    2 endog       2
+610   23    3 endog       3
+610   NA    4 endog       4
+610   11    5 endog       5
+302   18    0 endog       0
+302   22    1 endog       1
+302   16    2 endog       2
+302    8    3 endog       3
+302    9    4 endog       4
+302   12    5 endog       5
+303   21    0 nonen       0
+303   21    1 nonen       0
+303   13    2 nonen       0
+303   14    3 nonen       0
+303   10    4 nonen       0
+303    5    5 nonen       0
+304   21    0 endog       0
+304   27    1 endog       1
+304   29    2 endog       2
+304   NA    3 endog       3
+304   12    4 endog       4
+304   24    5 endog       5
+305   19    0 nonen       0
+305   17    1 nonen       0
+305   15    2 nonen       0
+305   11    3 nonen       0
+305    5    4 nonen       0
+305    1    5 nonen       0
+308   22    0 nonen       0
+308   21    1 nonen       0
+308   18    2 nonen       0
+308   17    3 nonen       0
+308   12    4 nonen       0
+308   11    5 nonen       0
+309   22    0 nonen       0
+309   22    1 nonen       0
+309   16    2 nonen       0
+309   19    3 nonen       0
+309   20    4 nonen       0
+309   11    5 nonen       0
+310   24    0 endog       0
+310   19    1 endog       1
+310   11    2 endog       2
+310    7    3 endog       3
+310    6    4 endog       4
+310   NA    5 endog       5
+311   20    0 endog       0
+311   16    1 endog       1
+311   21    2 endog       2
+311   17    3 endog       3
+311   NA    4 endog       4
+311   15    5 endog       5
+312   17    0 endog       0
+312   NA    1 endog       1
+312   18    2 endog       2
+312   17    3 endog       3
+312   17    4 endog       4
+312    6    5 endog       5
+313   21    0 nonen       0
+313   19    1 nonen       0
+313   10    2 nonen       0
+313   11    3 nonen       0
+313   11    4 nonen       0
+313    8    5 nonen       0
+315   27    0 endog       0
+315   21    1 endog       1
+315   17    2 endog       2
+315   13    3 endog       3
+315    5    4 endog       4
+315   NA    5 endog       5
+316   32    0 endog       0
+316   26    1 endog       1
+316   23    2 endog       2
+316   26    3 endog       3
+316   23    4 endog       4
+316   24    5 endog       5
+318   17    0 endog       0
+318   18    1 endog       1
+318   19    2 endog       2
+318   21    3 endog       3
+318   17    4 endog       4
+318   11    5 endog       5
+319   24    0 endog       0
+319   18    1 endog       1
+319   10    2 endog       2
+319   14    3 endog       3
+319   13    4 endog       4
+319   12    5 endog       5
+322   28    0 endog       0
+322   21    1 endog       1
+322   25    2 endog       2
+322   32    3 endog       3
+322   34    4 endog       4
+322   NA    5 endog       5
+327   17    0 nonen       0
+327   18    1 nonen       0
+327   15    2 nonen       0
+327    8    3 nonen       0
+327   19    4 nonen       0
+327   17    5 nonen       0
+328   22    0 nonen       0
+328   24    1 nonen       0
+328   28    2 nonen       0
+328   26    3 nonen       0
+328   28    4 nonen       0
+328   29    5 nonen       0
+331   19    0 nonen       0
+331   21    1 nonen       0
+331   18    2 nonen       0
+331   16    3 nonen       0
+331   14    4 nonen       0
+331   10    5 nonen       0
+333   23    0 nonen       0
+333   20    1 nonen       0
+333   21    2 nonen       0
+333   20    3 nonen       0
+333   24    4 nonen       0
+333   14    5 nonen       0
+334   31    0 nonen       0
+334   25    1 nonen       0
+334   NA    2 nonen       0
+334    7    3 nonen       0
+334    8    4 nonen       0
+334   11    5 nonen       0
+335   21    0 nonen       0
+335   21    1 nonen       0
+335   18    2 nonen       0
+335   15    3 nonen       0
+335   12    4 nonen       0
+335   10    5 nonen       0
+337   27    0 nonen       0
+337   22    1 nonen       0
+337   23    2 nonen       0
+337   21    3 nonen       0
+337   12    4 nonen       0
+337   13    5 nonen       0
+338   22    0 nonen       0
+338   20    1 nonen       0
+338   22    2 nonen       0
+338   23    3 nonen       0
+338   19    4 nonen       0
+338   18    5 nonen       0
+339   27    0 endog       0
+339   NA    1 endog       1
+339   14    2 endog       2
+339   12    3 endog       3
+339   11    4 endog       4
+339   12    5 endog       5
+344   NA    0 endog       0
+344   21    1 endog       1
+344   12    2 endog       2
+344   13    3 endog       3
+344   13    4 endog       4
+344   18    5 endog       5
+345   29    0 nonen       0
+345   27    1 nonen       0
+345   27    2 nonen       0
+345   22    3 nonen       0
+345   22    4 nonen       0
+345   23    5 nonen       0
+346   25    0 endog       0
+346   24    1 endog       1
+346   19    2 endog       2
+346   23    3 endog       3
+346   14    4 endog       4
+346   21    5 endog       5
+347   18    0 endog       0
+347   15    1 endog       1
+347   14    2 endog       2
+347   10    3 endog       3
+347    8    4 endog       4
+347   NA    5 endog       5
+348   24    0 nonen       0
+348   21    1 nonen       0
+348   12    2 nonen       0
+348   13    3 nonen       0
+348   12    4 nonen       0
+348    5    5 nonen       0
+349   17    0 endog       0
+349   19    1 endog       1
+349   15    2 endog       2
+349   12    3 endog       3
+349    9    4 endog       4
+349   13    5 endog       5
+350   22    0 nonen       0
+350   25    1 nonen       0
+350   12    2 nonen       0
+350   16    3 nonen       0
+350   10    4 nonen       0
+350   16    5 nonen       0
+351   30    0 endog       0
+351   27    1 endog       1
+351   23    2 endog       2
+351   20    3 endog       3
+351   12    4 endog       4
+351   11    5 endog       5
+352   21    0 endog       0
+352   19    1 endog       1
+352   18    2 endog       2
+352   15    3 endog       3
+352   18    4 endog       4
+352   19    5 endog       5
+353   27    0 endog       0
+353   21    1 endog       1
+353   24    2 endog       2
+353   22    3 endog       3
+353   16    4 endog       4
+353   11    5 endog       5
+354   28    0 endog       0
+354   27    1 endog       1
+354   27    2 endog       2
+354   26    3 endog       3
+354   23    4 endog       4
+354   NA    5 endog       5
+355   22    0 endog       0
+355   26    1 endog       1
+355   20    2 endog       2
+355   13    3 endog       3
+355   10    4 endog       4
+355    7    5 endog       5
+357   27    0 endog       0
+357   22    1 endog       1
+357   24    2 endog       2
+357   25    3 endog       3
+357   19    4 endog       4
+357   19    5 endog       5
+360   21    0 endog       0
+360   28    1 endog       1
+360   27    2 endog       2
+360   29    3 endog       3
+360   28    4 endog       4
+360   33    5 endog       5
+361   30    0 endog       0
+361   22    1 endog       1
+361   11    2 endog       2
+361    8    3 endog       3
+361    7    4 endog       4
+361   19    5 endog       5
+