QTM 385 - Experimental Methods

Lecture 13 - Introduction to Quarto and Declare Design

Danilo Freire

Emory University

Hello, everyone! 😊

Brief recap 📚

Last class, we discussed…

  • The historical development of research ethics (Nuremberg Code, Tuskegee Study, Belmont Report)
  • Core ethical principles: respect for persons, beneficence, and justice
  • Informed consent requirements and procedures
  • Considerations for anonymity and confidentiality in research data
  • How Institutional Review Boards (IRBs) protect human subjects
  • Ethical controversies in recent studies like Facebook’s emotional contagion study
  • Adaptive research designs that maximise benefits and minimize harm
  • Ethical responsibilities toward research staff and communities
  • How ethical research practices improve data quality and research outcomes
  • The principle of equipoise in determining when experiments are justified
  • The need for special protections for vulnerable populations
  • How ethical practices help prevent attrition and strengthen research validity

Today’s plan 📅

Introduction to Quarto and DeclareDesign

Two useful tools for experimental research (and beyond)

  • Quarto is a great tool to write and publish reproducible documents
  • It fits well with R Markdown and Jupyter Notebooks, so there is very little learning curve for you!
  • Quarto also meets all the requirements for a good research workflow and literate programming
    • It is text-based, so you can use version control (e.g., Git)
    • It supports code execution and visualisation
    • It allows you to write in Markdown, LaTeX, and other formats

Introduction to Quarto and DeclareDesign

Two useful tools for experimental research (and beyond)

  • DeclareDesign is a package for R that helps you design and analyse experiments
  • We have already used it before, but today we will see more of its capabilities
  • As we will see, DeclareDesign is a powerful tool to plan and execute experiments
    • It helps you define your research question and hypotheses
    • It guides you in choosing the right design and sample size
    • It provides tools to analyse your data and draw conclusions
  • And it is simulation-based, so you can test your design before running the experiment
  • I hope you can use it in your pre-analysis plans! 😉

Quarto 📝

What is Quarto?

  • Quarto (.qmd) is a document format that combines the best features of R Markdown and Jupyter Notebooks
  • It is designed for reproducibility and literate programming, as mentioned before
  • Literate programming is a programming paradigm that combines code and documentation in a single document
    • This makes your code more readable and easier to maintain
    • It also helps you communicate your results more effectively
  • If you use R Markdown or Jupyter Notebooks, you already do literate programming!
  • Why is literate programming important for experimental research?
    • It helps you keep track of your code and results
    • It makes it easier to share your work with others
    • It allows you to reproduce your analysis and results

Why Quarto?

Can’t I just use R Markdown or Jupyter Notebooks instead?

  • Yes, you can!
  • But Quarto has some advantages over R Markdown and Jupyter Notebooks
  • For example, Quarto supports multiple programming languages (not just R and Python)
  • It also has built-in support for interactive visualisations and data tables
  • And it is designed to work well with version control systems like Git
  • You can also use it to publish your documents online (like this one!), create slideshows, websites, and books
  • Learn it once, use it everywhere! 🤓

How to get started with Quarto?

Click on the image to visit the Quarto installation page

Writing a Quarto document 📝

How to write a Quarto document?

  • Quarto documents are directly influenced by R Markdown
  • All documents have three main components:
    • YAML header with metadata
    • Markdown content (text)
    • Code chunks (R, Python, Julia, etc.)
  • YAML stands for “YAML Ain’t Markup Language”
    • It is a human-readable data serialisation format
    • It is often used for configuration files and metadata
    • In Quarto, the YAML header contains information about the document (title, author, date, etc.)
---
title: My Quarto Document
subtitle: A simple example
author: Danilo Freire
date: 2025-03-03
format: html
editor:
  render-on-save: true
---

Markdown

  • Headings: #, ##, ###, etc.
  • Emphasis: *italic* or _italics_, **boldface** or __boldface__, ~~strikethrough~~
  • Lists: *, 1., -
  • Links: [text](url)
  • Images: ![alt text](url)
  • Code: `code`
  • Tables:
| Header 1 | Header 2 |
|----------|----------|
| Cell 1   | Cell 2   |

Heading

Italics, boldface, and code.

List:

  • Item 1
  • Item 2
    • Subitem 1

This is a link.

Header 1 Header 2
Cell 1 Cell 2

Markdown

  • Blockquotes: > quote
  • Horizontal rules: ---
  • Footnotes: [^1]
  • Superscript: 2^10^
  • Subscript: H~2~O
  • Math: $\LaTeX$
  • For equations, use $ or $$
  • Check the LaTeX Wiki for math symbols
  • More at Markdown Guide

This is a blockquote.

This is a paragraph1.

Equation: \(\mu = \frac{1}{n} \sum_{i=1}^{n} x_i\)

210 is 1024, and H2O is water.


This is a horizontal rule.

Quarto engines

knitr

---
title: "ggplot2 demo"
  format: 
    html:
      code-fold: true
---

## Meet Quarto

Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see <https://quarto.org>.

```{r}
#| label: plot-penguins
#| echo: false
#| message: false
#| warning: false

library(tidyverse)
library(palmerpenguins)

ggplot(penguins, 
       aes(x = flipper_length_mm, y = bill_length_mm)) +
geom_point(aes(color = species, shape = species)) +
scale_color_manual(values = c("darkorange","purple","cyan4")) +
  labs(
    title = "Flipper and bill length",
    subtitle = "Dimensions for penguins at Palmer Station LTER",
    x = "Flipper length (mm)", y = "Bill length (mm)",
    color = "Penguin species", shape = "Penguin species"
  ) +
theme_minimal()

Jupyter

---
title: "Palmer Penguins Demo"
  format: 
    html:
      code-fold: true
jupyter: python3
---


## Meet Quarto

Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see <https://quarto.org>.

```{python}
#| echo: false
#| message: false

import pandas as pd
import seaborn as sns 
from palmerpenguins import load_penguins
sns.set_style('whitegrid')

penguins = load_penguins()

g = sns.lmplot(x="flipper_length_mm",
               y="body_mass_g",
               hue="species",
               height=7,
               data=penguins,
               palette=['#FF8C00','#159090','#A034F0']);
g.set_xlabels('Flipper Length');
g.set_ylabels('Body Mass');

PDFs

  • Most articles and reports are in PDF format
  • You will need a LaTeX distribution to create PDFs
  • You can install it with quarto:
    • quarto install tinytex
  • If you use R, you can use tinytex with:
    • install.packages("tinytex")
    • tinytex::install_tinytex()
  • This is strongly recommended as a LaTeX distribution is huge and takes a long time to install
  • More information here

Adding references

  • You can add academic references to your documents using bibtex
  • BibTeX is a reference management software developed for LaTeX
  • A bibtex file is a plain text file with .bib extension
  • You can add references to your YAML header
  • Use @key to cite a reference, [@key] to cite in parentheses
  • The key is the first line of the reference in the .bib file (e.g., @article{nash1950equilibrium,)
  • More information here
  • The easiest way to manage references is with a reference manager like Zotero or download a BibTex reference from Google Scholar
  • Just search for a paper and click on the cite button
  • Then click on BibTeX and copy the reference to your .bib file

Example

@article{nash1950equilibrium,
  title={Equilibrium points in n-person games},
  author={Nash Jr, John F},
  journal={Proceedings of the national academy of sciences},
  volume={36},
  number={1},
  pages={48--49},
  year={1950},
  publisher={National Acad Sciences}
}
  • Save it on your references.bib file (or any other name, as long as you reference it in your YAML header and the file ends with .bib)
  • Quarto will automatically format the references
  • More information about BibTeX here

Example

Then you add the references to your YAML header:

---
title: "Matplotlib Demo"
author: "Norah Jones"
date: "2024-09-30"
date-format: "MMMM D, YYYY"
format: 
  pdf:
    fig-width: 3
    fig-height: 2
    margins: 2in
bibliography: references.bib
---

DeclareDesign 📊

What is DeclareDesign?

A package for designing and analysing data

  • DeclareDesign is an R package that helps you design and analyse experiments, quasi-experiments, and other observational studies
  • Its main goal is to provide tools for declaration, diagnosis, and redesign in code, that is, to help you plan and execute your research before data analysis
  • DeclareDesign can be thought of in 2 parts
    • The data fabrication step (this is the bulk of the code)
    • The evaluation step, where users can change initial assumptions and see how well their design actually worked (this is the computationally intensive part of the code)
  • The great thing about DeclareDesign is that, as far as I understand, you can simulate any kind of design you can think of (e.g., experiments, regression discontinuity, difference-in-differences, even qualitative research!)
  • However, most studies are not that complex, so you will rarely need to use all of DeclareDesign’s capabilities
  • The authors have written a book that explains the package in detail
  • And they have also created a companion package, DesignLibrary, with pre-built designs that you can use right away! 🤓

Six components of a DeclareDesign study

  • Population: Set of units about which inferences are sought and their characteristics
    • Where should the theory apply? Who are the units of interest?
  • Potential outcomes: Outcomes each unit might exhibit depending on how causal process changes the world
    • Rooted in theory; non-compliance, spillovers and attrition affect potential outcomes
  • Sampling strategy: Strategy used to select units to include in study
    • How are we selecting units to analyse?
  • Assignment: Manner in which units are assigned to reveal one potential outcome or another
    • Randomisation strategy used
  • Estimand: Quantities that we want to learn about in the world, in terms of potential outcomes
    • What are we trying to estimate? ATE? Difference in means? Other quantities?
  • Estimator: Procedure for generating estimates of the quantities we want to learn about
    • This is your statistical model (e.g., OLS, IV, etc.)

DeclareDesign in action

  • First, you need to install DeclareDesign with install.packages("DeclareDesign") and DesignLibrary with install.packages("DesignLibrary")
  • Under the hood, DeclareDesign uses fabricatr to generate data, randomizr to randomise units, and estimatr to estimate causal effects
  • All functions start with declare_, and you can chain them together with +
  • Once you have your design, you can run it with diagnose_design and draw_data
  • Then you can also use ggplot2 to visualise your results or dplyr to summarise them
  • The authors have written a tutorial that explains how to use DeclareDesign
  • And they have also created a gallery with examples of different designs
Design component Function Description
Model declare_model() background variables and potential outcomes
Inquiry declare_inquiry() research questions
Data strategy declare_sampling() sampling procedures
declare_assignment() assignment procedures
declare_measurement() measurement procedures
Answer strategy declare_estimator() estimation procedures
declare_test() testing procedures

Chapter 18 of the book explains how to use DeclareDesign for experimental causal inference. Lots of examples! Some we are going to see today 🤓

Example

Two-arm experiment

library(DeclareDesign)
library(tidyverse)
set.seed(385)

sample_size <- 50

two_arm <-
  declare_model(
    N = 10000,
    X = rep(c(0, 1), each = N / 2),
    U = rnorm(N, sd = 0.25),
    potential_outcomes(Y ~ 0.2 * Z + X + U)
  ) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_sampling(S = complete_rs(N = N, n = sample_size)) +
  declare_assignment(Z = complete_ra(N = N, m = sample_size/2)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z, .method = lm_robust, inquiry = "ATE")

diagnosis_two_arm <- diagnose_design(two_arm)
diagnosis_two_arm

Research design diagnosis based on 500 simulations. Diagnosis completed in 2 secs. Diagnosand estimates with bootstrapped standard errors in parentheses (100 replicates).

  Design Inquiry Estimator Outcome Term N Sims Mean Estimand Mean Estimate
 two_arm     ATE estimator       Y    Z    500          0.20          0.21
                                                      (0.00)        (0.01)
   Bias SD Estimate   RMSE  Power Coverage
   0.01        0.16   0.16   0.24     0.96
 (0.01)      (0.00) (0.00) (0.02)   (0.01)

Power is a bit too low…

Let’s redesign the experiment

set.seed(385)

two_arm_power <-
  two_arm |>
  redesign(sample_size = c(50, 100, 250, 500)) |>
  diagnose_designs() |>
  tidy() |>
  filter(diagnosand == "power")

two_arm_power |>
  select(design, sample_size, diagnosand, estimate)
    design sample_size diagnosand estimate
1 design_1          50      power    0.236
2 design_2         100      power    0.398
3 design_3         250      power    0.788
4 design_4         500      power    0.986
ggplot(two_arm_power, aes(sample_size, estimate)) + 
  geom_point() +
  geom_smooth(method = "loess", color = "blue") +
  geom_hline(yintercept = 0.8, color = "red") +
  labs(x = "Sample Size", y = "Statistical Power") +
  theme_minimal()

Two-arm experiments with covariates

set.seed(385)

model <-
  declare_model(
    N = 10000,
    X = rep(c(0, 1), each = N / 2),
    U = rnorm(N, sd = 0.25),
    potential_outcomes(Y ~ 0.2 * Z + X + U)
  ) 
inquiry <-
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0))
sampling <-
  declare_sampling(S = complete_rs(N = N, n = sample_size))
assignment <-
  declare_assignment(Z = complete_ra(N = N, m = sample_size/2))
measurement <-
  declare_measurement(Y = reveal_outcomes(Y ~ Z))
answer_strategy <-
  declare_estimator(Y ~ Z, .method = lm_robust, inquiry = "ATE", label = "OLS")

two_arms_a <-
  model +
  inquiry +
  sampling +
  assignment +
  measurement +
  answer_strategy

diagnosis_two_arms_a <- diagnose_design(two_arms_a)
diagnosis_two_arms_a

Research design diagnosis based on 500 simulations. Diagnosis completed in 2 secs. Diagnosand estimates with bootstrapped standard errors in parentheses (100 replicates).

     Design Inquiry Estimator Outcome Term N Sims Mean Estimand Mean Estimate
 two_arms_a     ATE       OLS       Y    Z    500          0.20          0.21
                                                         (0.00)        (0.01)
   Bias SD Estimate   RMSE  Power Coverage
   0.01        0.16   0.16   0.24     0.96
 (0.01)      (0.00) (0.00) (0.02)   (0.01)

Two-arm experiments with covariates

set.seed(385)

two_arm_cov <-
  model +
  inquiry +
  sampling +
  assignment +
  measurement +
  declare_estimator(Y ~ Z, covariates = ~ X, .method = lm_lin, inquiry = "ATE", label = "OLS with covariates")

diagnosis_two_arm_cov <- diagnose_design(two_arm_cov)
diagnosis_two_arm_cov

Research design diagnosis based on 500 simulations. Diagnosis completed in 2 secs. Diagnosand estimates with bootstrapped standard errors in parentheses (100 replicates).

      Design Inquiry           Estimator Outcome Term N Sims Mean Estimand
 two_arm_cov     ATE OLS with covariates       Y    Z    500          0.20
                                                                    (0.00)
 Mean Estimate   Bias SD Estimate   RMSE  Power Coverage
          0.20  -0.00        0.07   0.07   0.81     0.95
        (0.00) (0.00)      (0.00) (0.00) (0.02)   (0.01)

Two-arm experiments without covariates

simulations_ta <- get_simulations(diagnosis_two_arm)

# first create summary for vertical lines
summary_ta <- 
  simulations_ta |> 
  group_by(estimator) |> 
  summarize(estimand = mean(estimand))

# then plot simulations
ggplot(simulations_ta) +
  geom_histogram(aes(estimate),
                 bins = 40, fill = "#72B4F3") +
  geom_vline(data = summary_ta,
             aes(xintercept = estimand),
             lty = "dashed", color = "#C6227F") +
  annotate("text", y = 300, x = 0.18, label = "Estimand",
           color = "#C6227F", hjust = 1) +
  facet_wrap(~ estimator) +
  labs(x = "Estimate", y = "Count of simulations") +
  theme_minimal()

Two-arm experiments with covariates

simulations_tac <- get_simulations(diagnosis_two_arm_cov)

summary_tac <- 
  simulations_tac |> 
  group_by(estimator) |> 
  summarize(estimand = mean(estimand))

# then plot simulations
ggplot(simulations_tac) +
  geom_histogram(aes(estimate),
                 bins = 40, fill = "#72B4F3") +
  geom_vline(data = summary_tac,
             aes(xintercept = estimand),
             lty = "dashed", color = "#C6227F") +
  annotate("text", y = 300, x = 0.18, label = "Estimand",
           color = "#C6227F", hjust = 1) +
  facet_wrap(~ estimator) +
  labs(x = "Estimate", y = "Count of simulations") +
  theme_minimal()

Block randomisation two-arm experiment

set.seed(385)

block_randomisation <-
  model +
  inquiry +
  declare_sampling(S = strata_rs(strata = X, n = 50)) +
  declare_assignment(Z = block_ra(blocks = X)) +
  measurement +
  declare_estimator(Y ~ Z, covariates = ~ X, .method = lm_lin, inquiry = "ATE", label = "Lin")

diagnosis_block_randomisation <- diagnose_design(block_randomisation)
diagnosis_block_randomisation

Research design diagnosis based on 500 simulations. Diagnosis completed in 4 secs. Diagnosand estimates with bootstrapped standard errors in parentheses (100 replicates).

              Design Inquiry Estimator Outcome Term N Sims Mean Estimand
 block_randomisation     ATE       Lin       Y    Z    500          0.20
                                                                  (0.00)
 Mean Estimate   Bias SD Estimate   RMSE  Power Coverage
          0.20   0.00        0.05   0.05   0.97     0.94
        (0.00) (0.00)      (0.00) (0.00) (0.01)   (0.01)
simulations_br <- get_simulations(diagnosis_block_randomisation)

summary_br <- 
  simulations_br |> 
  group_by(estimator) |> 
  summarize(estimand = mean(estimand))

# then plot simulations
ggplot(simulations_br) +
  geom_histogram(aes(estimate),
                 bins = 40, fill = "#72B4F3") +
  geom_vline(data = summary_br,
             aes(xintercept = estimand),
             lty = "dashed", color = "#C6227F") +
  annotate("text", y = 300, x = 0.18, label = "Estimand",
           color = "#C6227F", hjust = 1) +
  facet_wrap(~ estimator) +
  labs(x = "Estimate", y = "Count of simulations") +
  theme_minimal()

Blocked and clustered randomisation

set.seed(385)

ICC <- 0.9

block_cluster <-
  declare_model(
    cluster =
      add_level(
        N = 10,
        cluster_size = rep(seq(10, 50, 10), 2),
        cluster_shock =
          scale(cluster_size + rnorm(N, sd = 5)) * sqrt(ICC),
        cluster_tau = rnorm(N, sd = sqrt(ICC))
      ),
    individual =
      add_level(
        N = cluster_size,
        individual_shock = rnorm(N, sd = sqrt(1 - ICC)),
        individual_tau = rnorm(N, sd = sqrt(1 - ICC)),
        Y_Z_0 = cluster_shock + individual_shock,
        Y_Z_1 = Y_Z_0 + cluster_tau + individual_tau
      )
  ) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_assignment(Z = block_and_cluster_ra(clusters = cluster, blocks = cluster_size)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z,
                    clusters = cluster,
                    inquiry = "ATE")

diagnosis_block_cluster <- diagnose_design(block_cluster)
diagnosis_block_cluster

Research design diagnosis based on 500 simulations. Diagnosis completed in 2 secs. Diagnosand estimates with bootstrapped standard errors in parentheses (100 replicates).

        Design Inquiry Estimator Outcome Term N Sims Mean Estimand
 block_cluster     ATE estimator       Y    Z    500          0.01
                                                            (0.01)
 Mean Estimate   Bias SD Estimate   RMSE  Power Coverage
         -0.02  -0.03        0.52   0.41   0.00     1.00
        (0.02) (0.02)      (0.01) (0.01) (0.00)   (0.00)

Redesigning over values of ICC

set.seed(385)

block_cluster_icc <-
  block_cluster |>
  redesign(ICC = seq(0.1, 0.9, by = 0.4)) |> 
  diagnose_designs() |>
  tidy()

block_cluster_icc |>
  select(design, ICC, diagnosand, estimate) |>
  filter(diagnosand %in% c("mean_estimate", "sd_estimate"))
    design ICC    diagnosand    estimate
1 design_1 0.1 mean_estimate  0.00151772
2 design_1 0.1   sd_estimate  0.22702412
3 design_2 0.5 mean_estimate -0.01250764
4 design_2 0.5   sd_estimate  0.41689454
5 design_3 0.9 mean_estimate  0.04037843
6 design_3 0.9   sd_estimate  0.53993123

DesignLibrary 📚

DesignLibrary

Making your life easier!

library(DesignLibrary)

# two_arm_designer(
  # N = 100,
  # assignment_prob = 0.5,
  # control_mean = 0,
  # control_sd = 1,
  # ate = 1,
  # treatment_mean = control_mean + ate,
  # treatment_sd = control_sd,
  # rho = 1,
  # args_to_fix = NULL
# )

two_arm_design <- two_arm_designer(ate = 2, N = 1000)

two_arm_design

Research design declaration summary

Step 1 (population): declare_population(N = N, u_0 = rnorm(N), u_1 = rnorm(n = N, mean = rho * u_0, sd = sqrt(1 - rho^2))) 

Step 2 (potential outcomes): declare_potential_outcomes(Y ~ (1 - Z) * (u_0 * control_sd + control_mean) + Z * (u_1 * treatment_sd + treatment_mean)) 

Step 3 (inquiry): declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) -------------------

Step 4 (assignment): declare_assignment(Z = complete_ra(N, prob = assignment_prob)) 

Step 5 (reveal): declare_reveal() ----------------------------------------------

Step 6 (estimator): declare_estimator(Y ~ Z, inquiry = estimand) ---------------

Run of the design:

 inquiry estimand estimator term estimate std.error statistic   p.value
     ATE        2 estimator    Z     1.94    0.0654      29.7 3.38e-139
 conf.low conf.high  df outcome
     1.81      2.07 998       Y

DesignLibrary

Clustered designs

  • You can create complex designs easily too!
# A design with varying block size and varying cluster size
cluster_sampling_design <- cluster_sampling_designer(
  N_blocks = 2,
  N_clusters_in_block = 100,
  N_i_in_cluster = 50,
  n_clusters_in_block = 10,
  n_i_in_cluster = 10,
  icc = 0.1)

cluster_sampling_design

Research design declaration summary

Step 1 (population): declare_population(data = fixed_pop) ----------------------

Step 2 (inquiry): declare_inquiry(mean(Y), label = "Ybar") ---------------------

Step 3 (sampling): declare_sampling(S1 = strata_and_cluster_rs(strata = block, clusters = cluster, n = n_clusters_in_block), filter = S1 == 1) 

Step 4 (sampling): declare_sampling(S2 = strata_rs(strata = cluster, n = n_i_in_cluster), filter = S2 == 1) 

Step 5 (estimator): declare_estimator(Y ~ 1, model = lm_robust, clusters = cluster, inquiry = estimand, label = "Clustered Standard Errors") 

Run of the design:
 inquiry estimand                 estimator        term estimate std.error
    Ybar     3.98 Clustered Standard Errors (Intercept)     3.82     0.168
 statistic  p.value conf.low conf.high df outcome
      22.8 3.01e-15     3.47      4.18 19       Y

Summary 📚

Today we learned…

  • The basics of Quarto and DeclareDesign
  • How to write a Quarto document that includes text, code, and visualisations
  • The importance of literate programming for reproducible research
  • How to use DeclareDesign to plan and execute experiments
  • The six components of a DeclareDesign study
  • How to use DesignLibrary to create common designs

Further readings

Quarto cheat sheet

Source: https://raw.githubusercontent.com/rstudio/cheatsheets/master/quarto.pdf

DeclareDesign cheat sheet

Source: https://raw.githubusercontent.com/rstudio/cheatsheets/master/declaredesign.pdf

And that’s all for today! 🎉

Have a great week! 🎉