| Title: | Companion Datasets and Functions for Research Design in the Social Sciences |
|---|---|
| Description: | Helper functions to accompany the Blair, Coppock, and Humphreys (2023) "Research Design in the Social Sciences: Declaration, Diagnosis, and Redesign" <https://book.declaredesign.org>. 'rdss' includes datasets, helper functions, and plotting components to enable use and replication of the book. |
| Authors: | Graeme Blair [aut] (ORCID: <https://orcid.org/0000-0001-9164-2102>), Alexander Coppock [aut, cre] (ORCID: <https://orcid.org/0000-0002-5733-2386>), Macartan Humphreys [aut] (ORCID: <https://orcid.org/0000-0001-7029-2326>) |
| Maintainer: | Alexander Coppock <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.16 |
| Built: | 2026-10-27 14:53:11 UTC |
| Source: | https://github.com/declaredesign/rdss |
Add parentheses around standard error estimates
add_parens(x, digits = 3)add_parens(x, digits = 3)
x |
A numeric vector. |
digits |
The number of decimal places. Defaults to 3. |
A character vector, e.g. "(0.120)".
std.error <- c(0.12, 0.001, 1.2) add_parens(std.error)std.error <- c(0.12, 0.001, 1.2) add_parens(std.error)
An inquiry handler. For each covariate, regresses the true unit-level effect tau on the covariate cut into intervals, and reports which covariate gives the highest R-squared. Pair it with causal_forest_handler(), whose var_imp column is the estimate of the same quantity.
best_predictor(data, covariate_names, cuts = 20)best_predictor(data, covariate_names, cuts = 20)
data |
A data frame with a column |
covariate_names |
A character vector of covariates to assess. |
cuts |
Either a numeric vector of two or more unique cut points or a single number (greater than or equal to 2) giving the number of intervals into which each covariate is to be cut. Defaults to 20. |
See https://book.declaredesign.org/complex-designs.html#discovery-using-causal-forests
A data frame with one row: inquiry ("best_predictor") and estimand, the position in covariate_names of the best predictor.
set.seed(343) dat <- data.frame(A = rnorm(500), B = rnorm(500)) dat$tau <- 1 + dat$A # The effect varies with A, so A (position 1) is the best predictor best_predictor(dat, covariate_names = c("A", "B"), cuts = 4)set.seed(343) dat <- data.frame(A = rnorm(500), B = rnorm(500)) dat$tau <- 1 + dat$A # The effect varies with A, so A (position 1) is the best predictor best_predictor(dat, covariate_names = c("A", "B"), cuts = 4)
One row per survey respondent. Z is the randomly assigned message (a factor with four levels) and blm_support the outcome; the remaining columns are respondent characteristics.
bonilla_tillerybonilla_tillery
A data frame with 849 rows and 10 columns: female, lgbtq, age, religiosity, income, college, blm_familiarity, linked_fate, Z, and blm_support.
Bonilla, Tabitha, and Alvin B. Tillery. 2020. "Which Identity Frames Boost Support for and Mobilization in the #BlackLivesMatter Movement? An Experimental Test." American Political Science Review 114(4): 947-962. doi:10.1017/S0003055420000544. Data from Dataverse, doi:10.7910/DVN/IUZDQI.
estimatr::lm_robust(blm_support ~ Z, data = bonilla_tillery)estimatr::lm_robust(blm_support ~ Z, data = bonilla_tillery)
Splits the data at random into a training and a test set, fits grf::causal_forest() on the training set, and predicts each unit's treatment effect. The data must contain the outcome as Y and the treatment as Z.
causal_forest_handler(data, covariate_names, share_train = 0.5, ...)causal_forest_handler(data, covariate_names, share_train = 0.5, ...)
data |
A data frame with columns |
covariate_names |
A character vector naming the covariate columns the forest may split on. |
share_train |
The share of units assigned to the training set. Defaults to 0.5. |
... |
Further arguments passed to |
See https://book.declaredesign.org/complex-designs.html#discovery-using-causal-forests
data with four columns added: pred (the predicted treatment effect), var_imp (the position in covariate_names of the most important covariate), and the logical indicators train and test.
library(DeclareDesign) library(ggplot2) dat <- fabricate( N = 1000, A = rnorm(N), B = rnorm(N), Z = complete_rs(N), Y = A*Z + rnorm(N)) # note: remove num.threads = 1 to use more processors estimates <- causal_forest_handler(data = dat, covariate_names = c("A", "B"), num.threads = 1) ggplot(data = estimates, aes(A, pred)) + geom_point()library(DeclareDesign) library(ggplot2) dat <- fabricate( N = 1000, A = rnorm(N), B = rnorm(N), Z = complete_rs(N), Y = A*Z + rnorm(N)) # note: remove num.threads = 1 to use more processors estimates <- causal_forest_handler(data = dat, covariate_names = c("A", "B"), num.threads = 1) ggplot(data = estimates, aes(A, pred)) + geom_point()
From "Estimating the Impact of the Hajj: Religion and Tolerance in Islam's Global Gathering," The Quarterly Journal of Economics 124(3): 1133-1170. One row per applicant to the Hajj visa lottery. success is 1 for lottery winners; the views_ columns are attitudes toward each group, and views is their sum.
clingingsmith_etalclingingsmith_etal
A data frame with 958 rows and 8 columns: success, views_saudi, views_indonesian, views_turkish, views_african, views_chinese, views_european, and views.
Clingingsmith, David, Asim Ijaz Khwaja, and Michael Kremer. 2009. "Estimating the Impact of the Hajj: Religion and Tolerance in Islam's Global Gathering." The Quarterly Journal of Economics 124(3): 1133-1170. doi:10.1162/qjec.2009.124.3.1133.
estimatr::difference_in_means(views ~ success, data = clingingsmith_etal)estimatr::difference_in_means(views ~ success, data = clingingsmith_etal)
For each attribute, assigns every row one of that attribute's levels by complete random assignment, so the levels appear in as equal numbers as possible.
conjoint_assignment(data, levels_list)conjoint_assignment(data, levels_list)
data |
A data frame with one row per profile. |
levels_list |
A named list with one element per attribute, each a vector of that attribute's levels, e.g. |
See https://book.declaredesign.org/experimental-descriptive.html#conjoint-experiments
data with one column per attribute, named as in levels_list.
set.seed(343) levels_list <- list(gender = c("woman", "man"), age = c("young", "old")) # 25 subjects, each choosing between two profiles in each of four tasks profiles <- expand.grid(profile = 1:2, task = 1:4, subject = 1:25) head(conjoint_assignment(profiles, levels_list))set.seed(343) levels_list <- list(gender = c("woman", "man"), age = c("young", "old")) # 25 subjects, each choosing between two profiles in each of four tasks profiles <- expand.grid(profile = 1:2, task = 1:4, subject = 1:25) head(conjoint_assignment(profiles, levels_list))
Computes the true average marginal component effect (AMCE) of every attribute level relative to that attribute's first level. Each AMCE is the average change in the probability that the second profile in a task is chosen when the attribute moves from its first level to the given level, with the other attributes drawn at random by conjoint_assignment().
conjoint_inquiries(data, levels_list, utility_fn)conjoint_inquiries(data, levels_list, utility_fn)
data |
A data frame with one row per profile and columns |
levels_list |
A named list of attribute levels, as in |
utility_fn |
A function that adds the utility column |
See https://book.declaredesign.org/experimental-descriptive.html#conjoint-experiments
A data frame with one row per non-reference level: attribute, reference, level, inquiry, and estimand.
set.seed(343) levels_list <- list(gender = c("woman", "man"), age = c("young", "old")) # 25 subjects, each choosing between two profiles in each of four tasks profiles <- expand.grid(profile = 1:2, task = 1:4, subject = 1:25) utility_fn <- function(data) { data$U <- 0.5 * (data$gender == "woman") - 0.2 * (data$age == "old") + rnorm(nrow(data)) data } conjoint_inquiries(profiles, levels_list, utility_fn)set.seed(343) levels_list <- list(gender = c("woman", "man"), age = c("young", "old")) # 25 subjects, each choosing between two profiles in each of four tasks profiles <- expand.grid(profile = 1:2, task = 1:4, subject = 1:25) utility_fn <- function(data) { data$U <- 0.5 * (data$gender == "woman") - 0.2 * (data$age == "old") + rnorm(nrow(data)) data } conjoint_inquiries(profiles, levels_list, utility_fn)
Applies utility_fn and, within each subject and task, compares the two profiles and sets choice to 1 for the one with the higher utility U (the second profile on a tie).
conjoint_measurement(data, utility_fn)conjoint_measurement(data, utility_fn)
data |
A data frame with one row per profile and columns |
utility_fn |
A function that takes |
See https://book.declaredesign.org/experimental-descriptive.html#conjoint-experiments
data with U and choice added.
set.seed(343) levels_list <- list(gender = c("woman", "man"), age = c("young", "old")) # 25 subjects, each choosing between two profiles in each of four tasks profiles <- expand.grid(profile = 1:2, task = 1:4, subject = 1:25) profiles <- conjoint_assignment(profiles, levels_list) utility_fn <- function(data) { data$U <- 0.5 * (data$gender == "woman") - 0.2 * (data$age == "old") + rnorm(nrow(data)) data } head(conjoint_measurement(profiles, utility_fn))set.seed(343) levels_list <- list(gender = c("woman", "man"), age = c("young", "old")) # 25 subjects, each choosing between two profiles in each of four tasks profiles <- expand.grid(profile = 1:2, task = 1:4, subject = 1:25) profiles <- conjoint_assignment(profiles, levels_list) utility_fn <- function(data) { data$U <- 0.5 * (data$gender == "woman") - 0.2 * (data$age == "old") + rnorm(nrow(data)) data } head(conjoint_measurement(profiles, utility_fn))
Based on Karthik Ram's wesanderson package (https://github.com/karthik/wesanderson)
dd_palette(name, n)dd_palette(name, n)
name |
The palette name, as a string. One of the names listed below. |
n |
The number of colors to return, from the start of the palette. Defaults to all of them. |
Palettes:
three_color_palette: light blue, orange, pink
grey_palette: light blue, orange, pink, light gray
quilt_palette: light gray, pink, purple, light blue, orange
two_color_palette: dark blue, pink
quilt_three_color_palette: light gray, translucent light blue, light blue
two_color_gray: dark blue, light gray
Single colors: dd_dark_blue ("#3564ED"), dd_light_blue ("#72B4F3"), dd_orange ("#F38672"), dd_purple ("#7E43B6"), dd_gray (gray(0.2)), dd_pink ("#C6227F"), dd_light_gray (gray(0.8)), and the translucent dd_dark_blue_alpha and dd_light_blue_alpha.
A character vector of hex color codes.
dd_palette("three_color_palette") dd_palette("quilt_palette", n = 3)dd_palette("three_color_palette") dd_palette("quilt_palette", n = 3)
Runs DIDmultiplegt::did_multiplegt(), the de Chaisemartin and D'Haultfoeuille estimator, and returns its estimates as a data frame.
did_multiplegt_tidy(data, ...)did_multiplegt_tidy(data, ...)
data |
A data frame, passed to |
... |
Further arguments to |
Use mode = "dyn". It estimates the effect l periods after a group first switches treatment, for l from 1 to effects; effects = 1 is the effect in the period of the switch. This mode runs did_multiplegt_dyn() from the 'DIDmultiplegtDYN' package, which also needs the 'polars' package. polars is not on CRAN; install it with
install.packages("polars", repos = "https://rpolars.r-universe.dev")
and attach it with library(polars).
mode = "old" is still passed through, but in 'DIDmultiplegt' 2.1.0 it returns NaN: it takes first differences with stats::lag(), which does not shift the data, so no group is ever seen to switch. The function warns when that happens.
See https://book.declaredesign.org/observational-causal.html#difference-in-differences
For mode = "dyn", a data frame with one row per effect and columns term ("Effect_1", "Effect_2", ...), estimate, std.error, conf.low, conf.high, n (the number of group-period cells used), and n_switchers (the number of switching groups). For mode = "old", a data frame with one column, estimate.
library(polars) set.seed(343) # 20 units in 10 periods; units switch into treatment at different times panel <- data.frame(unit = rep(1:20, times = 10), period = rep(1:10, each = 20)) panel$D <- as.numeric(panel$period > (panel$unit %% 5) + 4) panel$Y <- 0.5 * panel$D + rnorm(nrow(panel)) did_multiplegt_tidy( panel, mode = "dyn", outcome = "Y", group = "unit", time = "period", treatment = "D", effects = 2 )library(polars) set.seed(343) # 20 units in 10 periods; units switch into treatment at different times panel <- data.frame(unit = rep(1:20, times = 10), period = rep(1:10, each = 20)) panel$D <- as.numeric(panel$period > (panel$unit %% 5) + 4) panel$Y <- 0.5 * panel$D + rnorm(nrow(panel)) did_multiplegt_tidy( panel, mode = "dyn", outcome = "Y", group = "unit", time = "period", treatment = "D", effects = 2 )
Runs the Aronow and Samii estimator from the 'interference' package and returns Horvitz-Thompson and Hajek estimates of the total, direct, and indirect effects. Exposure is defined by one-hop neighbours in adj_matrix. The data must contain the treatment as Z and the outcome as Y.
estimator_AS_tidy( data, permutatation_matrix = NULL, adj_matrix, obs_prob_exposure = NULL )estimator_AS_tidy( data, permutatation_matrix = NULL, adj_matrix, obs_prob_exposure = NULL )
data |
A data frame with columns |
permutatation_matrix |
(The spelling is the published one and is kept so existing code runs.) A matrix of random treatment assignments with one
row per permutation of the treatment vector and one column per unit, as returned
by |
adj_matrix |
An adjacency matrix defining the network structure. This can be created, for example, as follows:
adjacency <- fairfax |>
as("Spatial") |>
spdep::poly2nb(queen = TRUE) |>
spdep::nb2mat(style = "B", zero.policy = TRUE)
|
obs_prob_exposure |
Optional exposure probabilities. They depend only on the permutation matrix and the network, so in a simulation they are the same on every draw; computing them once and passing them here avoids recomputing them each time. For example: prob_exposure <- interference::make_exposure_prob( potential_tr_vector = permutatation_matrix, adj_matrix = adjacency, exposure_map_fn = interference::make_exposure_map_AS, exposure_map_fn_add_args = list(hop = 1) ) When |
The function requires the 'interference' package, which is not available on CRAN.
To use this function, install it with remotes::install_github('szonszein/interference')
Without it the function returns nothing and explains why, so a design that includes it still declares and diagnoses.
See https://book.declaredesign.org/experimental-causal.html#experiments-over-networks
A data frame with six rows (three inquiries by two estimators) and columns term, inquiry, estimator, estimate, std.error, p.value, conf.low, and conf.high. Each effect is a contrast with the no exposure condition. The Hajek rows use the Hajek variance estimator. The Horvitz-Thompson rows use the conservative estimator of Aronow and Samii (2017), the larger of the Horvitz-Thompson variance and the constant-effects variance, because the Horvitz-Thompson variance of a difference can be negative. The p-values and 95 percent confidence intervals use a normal approximation. When an exposure condition is rare the variance estimates are noisy and the intervals can undercover: in 300 simulations of the book's Fairfax design with 10 percent of units treated, the total effect's intervals covered 82 percent of the time (Horvitz-Thompson) and 90 percent (Hajek); the direct and indirect effects' covered 94 to 98 percent.
set.seed(343) # A ring of 20 units, each adjacent to the units on either side of it N <- 20 adjacency <- matrix(0, N, N) adjacency[cbind(1:N, c(2:N, 1))] <- 1 adjacency <- adjacency + t(adjacency) Z <- randomizr::complete_ra(N) # 100 possible assignments, one per row permutations <- t(replicate(100, randomizr::complete_ra(N))) dat <- data.frame(Z = Z, Y = Z + rnorm(N)) estimator_AS_tidy(dat, permutatation_matrix = permutations, adj_matrix = adjacency)set.seed(343) # A ring of 20 units, each adjacent to the units on either side of it N <- 20 adjacency <- matrix(0, N, N) adjacency[cbind(1:N, c(2:N, 1))] <- 1 adjacency <- adjacency + t(adjacency) Z <- randomizr::complete_ra(N) # 100 possible assignments, one per row permutations <- t(replicate(100, randomizr::complete_ra(N))) dat <- data.frame(Z = Z, Y = Z + rnorm(N)) estimator_AS_tidy(dat, permutatation_matrix = permutations, adj_matrix = adjacency)
Boundaries of the 238 voting precincts in Fairfax County, Virginia, with each precinct's ID, name, district, and polling place. One row per precinct. The book uses it to build the adjacency matrix for the experiment over networks; see estimator_AS_tidy().
fairfaxfairfax
An sf object with 238 rows and 10 columns:
Fairfax County Government, voting precincts layer (REGMGR.VOTING_PRECINCTS).
names(fairfax)names(fairfax)
One row per voter. treat is the randomly assigned treatment, marked_register_2014 the turnout outcome, and ward and street (a hashed identifier) the levels at which voters are clustered.
foos_etalfoos_etal
A data frame with 8,375 rows and 5 columns: marked_register_2014, ward, street, treat, and weights.
Foos, Florian, Peter John, Christian Müller, and Kevin Cunningham. 2021. "Social Mobilization in Partisan Spaces." The Journal of Politics 83(3): 1190-1197. doi:10.1086/710970. Data from Dataverse, doi:10.7910/DVN/NDPXND.
estimatr::lm_robust(marked_register_2014 ~ treat, clusters = street, data = foos_etal)estimatr::lm_robust(marked_register_2014 ~ treat, clusters = street, data = foos_etal)
Rounds to digits decimal places and keeps trailing zeros, so a column of numbers lines up in a table.
format_num(x, digits = 3)format_num(x, digits = 3)
x |
A numeric vector. |
digits |
The number of decimal places. Defaults to 3. |
A character vector.
std.error <- c(0.12, 0.001, 1.2) format_num(std.error)std.error <- c(0.12, 0.001, 1.2) format_num(std.error)
Converts the exposure map from interference::make_exposure_map_AS(), a 0/1 matrix with one row per unit and one column per exposure condition, into a vector naming each unit's condition (e.g. "dir_ind1" or "no").
get_exposure_AS(obs_exposure)get_exposure_AS(obs_exposure)
obs_exposure |
An exposure map, as returned by |
See https://book.declaredesign.org/experimental-causal.html#experiments-over-networks
A character vector with one element per unit.
set.seed(343) # A ring of 20 units, each adjacent to the units on either side of it N <- 20 adjacency <- matrix(0, N, N) adjacency[cbind(1:N, c(2:N, 1))] <- 1 adjacency <- adjacency + t(adjacency) Z <- randomizr::complete_ra(N) exposure <- interference::make_exposure_map_AS( adj_matrix = adjacency, tr_vector = Z, hop = 1 ) get_exposure_AS(exposure)set.seed(343) # A ring of 20 units, each adjacent to the units on either side of it N <- 20 adjacency <- matrix(0, N, N) adjacency[cbind(1:N, c(2:N, 1))] <- 1 adjacency <- adjacency + t(adjacency) Z <- randomizr::complete_ra(N) exposure <- interference::make_exposure_map_AS( adj_matrix = adjacency, tr_vector = Z, hop = 1 ) get_exposure_AS(exposure)
Downloads one of the book's design declarations or diagnosis objects from its Dataverse archive and returns it as an R object. Requires an internet connection.
get_rdss_file(name, verbose = TRUE)get_rdss_file(name, verbose = TRUE)
name |
The name of the file, as a string, e.g. |
verbose |
If |
See https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/HYVPO5 for further details and the code used to create these files.
The available names include:
Design declaration objects:
declaration_9.5
declaration_2.1
declaration_2.2
declaration_4.1
declaration_5.1
declaration_7.1
declaration_9.1
declaration_9.2
declaration_9.3
declaration_9.4
declaration_9.6
declaration_9.7
declaration_10.1
declaration_10.2
declaration_10.3
declaration_10.4
declaration_10a
declaration_11.1
declaration_11.2
declaration_11.3
declaration_11.4
declaration_11.5
declaration_12.1a
declaration_12.1b
declaration_12.1c
declaration_12.1d
declaration_13.1
declaration_13.2
declaration_15.1
declaration_15.2
declaration_15.3a
declaration_15.3b
declaration_15.3c
declaration_15.4
declaration_15.5
declaration_15.6
declaration_16.1a
declaration_16.1b
declaration_16.2
declaration_16.3
declaration_16.4
declaration_16.5
declaration_16.6
declaration_17.1
declaration_17.2
declaration_17.3
declaration_17.4
declaration_17.5
declaration_17.6_a
declaration_17.6_b
declaration_18.1
declaration_18.2
declaration_18.3
declaration_18.4
declaration_18.5
declaration_18.6
declaration_18.7
declaration_18.8
declaration_18.9a
declaration_18.9b
declaration_18.9c
declaration_18.10
declaration_18.11
declaration_18.12
declaration_18.13
declaration_19.1
declaration_19.2
declaration_19.3
declaration_19.4
declaration_23.1a
declaration_23.1b
declaration_23.1c
declaration_23.1d
Diagnosis objects:
diagnosis_2.1
diagnosis_4.1
diagnosis_9.1
diagnosis_9.2
diagnosis_9.3
diagnosis_9.4
diagnosis_9.5
diagnosis_9.6
diagnosis_9.7
simulation_10.1
diagnosis_10.1
diagnosis_10.2
diagnosis_10.3
diagnosis_10.4
diagnosis_10.5
diagnosis_10a
diagnosis_11.1
diagnosis_11.2
diagnosis_11.3
diagnosis_11.4
diagnosis_11.5
diagnosis_12.1
diagnosis_12.2
diagnosis_13.1
diagnosis_15.1
diagnosis_15.2
diagnosis_15.3
diagnosis_15.4
diagnosis_15.5
diagnosis_16.1
diagnosis_16.2
diagnosis_16.3
diagnosis_16.4
diagnosis_16.5
diagnosis_17.1
diagnosis_17.2
diagnosis_17.3
diagnosis_17.4
diagnosis_17.5
diagnosis_18.1
diagnosis_18.10_encouragment
diagnosis_18.10_placebo
diagnosis_18.11
diagnosis_18.12
diagnosis_18.13
diagnosis_18.2
diagnosis_18.3
diagnosis_18.4
diagnosis_18.5
diagnosis_18.6
diagnosis_18.7
diagnosis_18.8
diagnosis_18.9
diagnosis_19.1
diagnosis_19.2
diagnosis_19.3
diagnosis_19.4
diagnosis_19a
diagnosis_21a
diagnosis_21b
diagnosis_23.1
diagnosis_23a
The requested object: a design declaration, or a diagnosis or simulation data frame.
## Not run: # Requires internet access if(curl::has_internet()) { diagnosis_2.1 <- get_rdss_file("diagnosis_2.1") diagnosis_2.1 } ## End(Not run)## Not run: # Requires internet access if(curl::has_internet()) { diagnosis_2.1 <- get_rdss_file("diagnosis_2.1") diagnosis_2.1 } ## End(Not run)
Add alpha transparency to a color defined in hexadecimal
hex_add_alpha(col, alpha)hex_add_alpha(col, alpha)
col |
A color as a six-digit hex code, e.g. |
alpha |
Opacity, from 0 (transparent) to 1 (opaque). |
The color as an eight-digit hex code.
hex_add_alpha("#72B4F3", 0.5)hex_add_alpha("#72B4F3", 0.5)
Party registration, age, census tract, and 2012 turnout for 1,000 registered voters sampled at random from the Los Angeles County voter file. One row per voter.
la_voter_filela_voter_file
A data frame with 1000 rows and 4 variables:
political party registration
age of voter in years
US Census tract number
1 if the voter voted in the 2012 general election, 0 otherwise
California Secretary of State.
table(la_voter_file$party)table(la_voter_file$party)
Lags x within each group, ordering by order_by rather than by row order. Write the variables as bare column names inside fabricatr::fabricate() or dplyr::mutate(); see the README for an example.
lag_by_group(x, groups, n = 1, order_by, default = NA)lag_by_group(x, groups, n = 1, order_by, default = NA)
x |
The variable to lag. |
groups |
The grouping variable, e.g. a unit identifier. |
n |
A positive integer: the number of periods to lag by. Defaults to 1. |
order_by |
Ordering variable within group (e.g., time) |
default |
The value for rows with no earlier period in their group. Defaults to |
See https://book.declaredesign.org/observational-causal.html#difference-in-differences
A vector the length of x, in the original row order.
# Three units observed in four periods; unit u is treated after period u panel <- expand.grid(unit = 1:3, period = 1:4) panel$D <- as.numeric(panel$period > panel$unit) dplyr::mutate(panel, D_lag = lag_by_group(D, groups = unit, order_by = period))# Three units observed in four periods; unit u is treated after period u panel <- expand.grid(unit = 1:3, period = 1:4) panel$D <- as.numeric(panel$period > panel$unit) dplyr::mutate(panel, D_lag = lag_by_group(D, groups = unit, order_by = period))
Used in the book's student exercises. These data were resampled with replacement from the LAPOP data, to 10,000 rows, for a subset of variables. These data cannot be used for scientific inferences, and are only useful for teaching purposes. ID numbers were scrambled so that individuals and municipalities cannot easily be identified.
lapop_brazillapop_brazil
A data frame with 10,000 rows and 10 columns: ID, municipality, and eight survey items (trust_police, govt_responsive, ideology, govt_pride, self_efficacy_political, trust_military, trust_supreme_court, and support_political_system).
Download the original data from https://www.vanderbilt.edu/cgd/data-access/
See https://www.vanderbilt.edu/lapop/core-surveys.php for the survey questionnaire.
LAPOP AmericasBarometer 2018/19 survey of Brazil (fieldwork January to March 2019), Vanderbilt University, https://www.vanderbilt.edu/cgd/data-access/.
summary(lapop_brazil$trust_police)summary(lapop_brazil$trust_police)
Format confidence intervals for nice printing
make_interval_entry(conf.low, conf.high, digits = 2)make_interval_entry(conf.low, conf.high, digits = 2)
conf.low |
A numeric vector of lower bounds. |
conf.high |
A numeric vector of upper bounds. |
digits |
The number of decimal places. Defaults to 2. |
A character vector, e.g. "[-0.17, 0.31]".
conf.low <- c(-0.1652, 0.00304, -6.352) conf.high <- c(0.3052, 0.00696, -1.648) make_interval_entry(conf.low, conf.high)conf.low <- c(-0.1652, 0.00304, -6.352) conf.high <- c(0.3052, 0.00696, -1.648) make_interval_entry(conf.low, conf.high)
Format estimates and standard errors for nice printing
make_se_entry(estimate, std.error, digits = 2)make_se_entry(estimate, std.error, digits = 2)
estimate |
A numeric vector of estimates. |
std.error |
A numeric vector of standard errors, the same length as |
digits |
The number of decimal places. Defaults to 2. |
A character vector, e.g. "0.07 (0.12)".
estimate <- c(0.07, 0.005, -4) std.error <- c(0.12, 0.001, 1.2) make_se_entry(estimate, std.error)estimate <- c(0.07, 0.005, -4) std.error <- c(0.12, 0.001, 1.2) make_se_entry(estimate, std.error)
Predicts each row's outcome from a fitted model, then averages the predictions within each group using the post-stratification weights. In the book, the rows are demographic cells and the groups are states.
post_stratification_helper(model_fit, data, group, weights)post_stratification_helper(model_fit, data, group, weights)
model_fit |
A fitted model that |
data |
A data frame of post-stratification cells to predict for. |
group |
The bare column name of the grouping variable. |
weights |
The bare column name of the post-stratification weights, e.g. each cell's population share. |
See https://book.declaredesign.org/observational-descriptive.html#multi-level-regression-and-poststratification
A data frame with one row per group and columns for the group and estimate.
set.seed(343) # Post-stratification cells: high school graduates and non-graduates in ten # states, weighted by each state's share of graduates cells <- expand.grid(state = state.name[1:10], HS = 0:1) prob_HS <- rep(seq(0.5, 0.9, length.out = 10), 2) cells$PS_weight <- ifelse(cells$HS == 1, prob_HS, 1 - prob_HS) survey <- cells[sample(nrow(cells), 1000, replace = TRUE), ] survey$Y <- rbinom(1000, 1, prob = plogis(0.5 * survey$HS)) model_fit <- lme4::glmer(Y ~ HS + (1 | state), data = survey, family = binomial) post_stratification_helper(model_fit, data = cells, group = state, weights = PS_weight)set.seed(343) # Post-stratification cells: high school graduates and non-graduates in ten # states, weighted by each state's share of graduates cells <- expand.grid(state = state.name[1:10], HS = 0:1) prob_HS <- rep(seq(0.5, 0.9, length.out = 10), 2) cells$PS_weight <- ifelse(cells$HS == 1, prob_HS, 1 - prob_HS) survey <- cells[sample(nrow(cells), 1000, replace = TRUE), ] survey$Y <- rbinom(1000, 1, prob = plogis(0.5 * survey$HS)) model_fit <- lme4::glmer(Y ~ HS + (1 | state), data = survey, family = binomial) post_stratification_helper(model_fit, data = cells, group = state, weights = PS_weight)
Answers a causal query about one case using a causal model from 'CausalQueries', once for each process tracing strategy. A strategy names the nodes whose values are observed, so the estimate shows how the answer changes as more of the case's data is examined.
process_tracing_estimator(causal_model, query, data, strategies)process_tracing_estimator(causal_model, query, data, strategies)
causal_model |
A model from |
query |
The causal query, in 'CausalQueries' syntax, e.g. |
data |
A one-row data frame giving the case's value on each observed node. Nodes that are |
strategies |
A character vector of strategies, each listing the nodes observed, separated by hyphens, e.g. |
See https://book.declaredesign.org/observational-causal.html#process-tracing
A data frame with one row per strategy: term (the strategy), XY (the case's values of X and Y, e.g. "X1Y1"), and estimate (the posterior mean of the query).
# Simple example showing ambiguity in attribution process_tracing_estimator( causal_model = CausalQueries::make_model("X -> Y"), query = "Y[X=1] > Y[X=0]", data = data.frame(X=1, Y = 1), strategies = "X-Y") # Example where M=1 acts as a hoop test process_tracing_estimator( causal_model = CausalQueries::make_model("X -> M -> Y") |> CausalQueries::set_restrictions("Y[M=1] < Y[M=0]") |> CausalQueries::set_restrictions("M[X=1] < M[X=0]"), query = "Y[X=1] > Y[X=0]", data = data.frame(X=1, Y = 1, M = 0), strategies = c("Y", "X-Y", "X-M-Y"))# Simple example showing ambiguity in attribution process_tracing_estimator( causal_model = CausalQueries::make_model("X -> Y"), query = "Y[X=1] > Y[X=0]", data = data.frame(X=1, Y = 1), strategies = "X-Y") # Example where M=1 acts as a hoop test process_tracing_estimator( causal_model = CausalQueries::make_model("X -> M -> Y") |> CausalQueries::set_restrictions("Y[M=1] < Y[M=0]") |> CausalQueries::set_restrictions("M[X=1] < M[X=0]"), query = "Y[X=1] > Y[X=0]", data = data.frame(X=1, Y = 1, M = 0), strategies = c("Y", "X-Y", "X-M-Y"))
declare_estimator
Lets rdrobust::rdrobust(), which takes vectors, be called with a data frame and column names. Pass it as .method in declare_estimator(); tidy.rdrobust() then tidies the fit.
rdrobust_helper(data, y, x, subset = NULL, ...)rdrobust_helper(data, y, x, subset = NULL, ...)
data |
A data frame. |
y |
The bare column name of the outcome. |
x |
The bare column name of the running variable. |
subset |
An optional logical expression in terms of the columns of |
... |
Further arguments to |
See https://book.declaredesign.org/observational-causal.html#regression-discontinuity-designs
rdrobust model fit object
set.seed(343) rd <- data.frame(X = runif(1000, min = -1, max = 1)) rd$Y <- 0.5 * (rd$X > 0) + rd$X + rnorm(1000, sd = 0.5) fit <- rdrobust_helper(rd, y = Y, x = X, c = 0) tidy(fit)set.seed(343) rd <- data.frame(X = runif(1000, min = -1, max = 1)) rd$Y <- 0.5 * (rd$X > 0) + rd$X + rnorm(1000, sd = 0.5) fit <- rdrobust_helper(rd, y = Y, x = X, c = 0) tidy(fit)
Companion datasets and functions for the book Research Design in the Social Sciences: Declaration, Diagnosis, and Redesign by Graeme Blair, Alexander Coppock, and Macartan Humphreys (Princeton University Press, 2023), free to read at https://book.declaredesign.org.
Most functions fill one step of a design declaration in the book, and each help page links to the section that uses it. The Get started article on the package website, https://declaredesign.org/r/rdss/articles/rdss.html, lists them chapter by chapter.
Maintainer: Alexander Coppock [email protected] (ORCID)
Authors:
Alexander Coppock [email protected] (ORCID)
Graeme Blair [email protected] (ORCID)
Macartan Humphreys [email protected] (ORCID)
Useful links:
Report bugs at https://github.com/DeclareDesign/rdss/issues
See https://book.declaredesign.org/complex-designs.html#meta-analysis
rma_helper(data, yi, sei, method = "REML", ...)rma_helper(data, yi, sei, method = "REML", ...)
data |
A data frame with one row per study. |
yi |
The bare column name of the study estimates. |
sei |
The bare column name of their standard errors. |
method |
character string to specify whether a fixed- or a random/mixed-effects model should be fitted. A fixed-effects model (with or without moderators) is fitted when using method = "FE". Random/mixed-effects models are fitted by setting method equal to one of the following: "DL", "HE", "SJ", "ML", "REML", "EB", "HS", "HSk", or "GENQ". Default is "REML". |
... |
Further arguments to |
Fits metafor::rma() to study-level estimates and standard errors. If the fit fails, it returns the error as an object rma_mu_tau() recognizes, so one failed simulation does not stop a diagnosis.
An rma.uni fit; pass it to rma_mu_tau() to tidy it.
set.seed(343) studies <- data.frame( est = rnorm(10, mean = 0.2, sd = 0.1), se = runif(10, min = 0.05, max = 0.15) ) fit <- rma_helper(studies, yi = est, sei = se) fitset.seed(343) studies <- data.frame( est = rnorm(10, mean = 0.2, sd = 0.1), se = runif(10, min = 0.05, max = 0.15) ) fit <- rma_helper(studies, yi = est, sei = se) fit
A tidier for fits from rma_helper(): mu is the average effect across studies and tau_sq the variance of the true effects across studies.
rma_mu_tau(fit)rma_mu_tau(fit)
fit |
A fit from |
See https://book.declaredesign.org/complex-designs.html#meta-analysis
A data frame with two rows, mu and tau_sq, with estimates, standard errors, and confidence intervals (for tau_sq, intervals only when method = "REML"). If the fit failed, one row with estimate = NA and error = TRUE.
set.seed(343) studies <- data.frame( est = rnorm(10, mean = 0.2, sd = 0.1), se = runif(10, min = 0.05, max = 0.15) ) fit <- rma_helper(studies, yi = est, sei = se) rma_mu_tau(fit)set.seed(343) studies <- data.frame( est = rnorm(10, mean = 0.2, sd = 0.1), se = runif(10, min = 0.05, max = 0.15) ) fit <- rma_helper(studies, yi = est, sei = se) rma_mu_tau(fit)
ggplot2::theme_minimal() with light gray major gridlines only, no axis ticks, and no legend. The book labels series directly, so it hides the legend; add theme(legend.position = "right") after it to get one back.
theme_dd()theme_dd()
A ggplot2 theme.
library(ggplot2) ggplot(mtcars, aes(wt, mpg)) + geom_point() + theme_dd()library(ggplot2) ggplot(mtcars, aes(wt, mpg)) + geom_point() + theme_dd()
Deprecated: use broom.mixed::tidy(), which this function now calls with the same arguments.
tidy_stan(x, conf.int = FALSE, conf.level = 0.95, exponentiate = FALSE, ...)tidy_stan(x, conf.int = FALSE, conf.level = 0.95, exponentiate = FALSE, ...)
x |
A fit from |
conf.int |
Logical indicating whether or not to include a confidence interval in the tidied output. Defaults to FALSE. |
conf.level |
The confidence level to use for the confidence interval if conf.int = TRUE. Must be strictly greater than 0 and less than 1. Defaults to 0.95, which corresponds to a 95 percent confidence interval. |
exponentiate |
Logical indicating whether or not to exponentiate the coefficient estimates. Defaults to FALSE. |
... |
Further arguments to |
See https://book.declaredesign.org/choosing-an-answer-strategy.html#bayesian-formalizations
A data frame with one row per coefficient.
fit <- rstanarm::stan_glm(mpg ~ wt, data = mtcars, chains = 1, iter = 1000, refresh = 0, seed = 343) # Prefer broom.mixed::tidy(fit, conf.int = TRUE) tidy_stan(fit, conf.int = TRUE)fit <- rstanarm::stan_glm(mpg ~ wt, data = mtcars, chains = 1, iter = 1000, refresh = 0, seed = 343) # Prefer broom.mixed::tidy(fit, conf.int = TRUE) tidy_stan(fit, conf.int = TRUE)
A generics::tidy() method for fits from cjoint::amce(): one row per attribute level, with normal-approximation confidence intervals.
## S3 method for class 'amce' tidy(x, alpha = 0.05, ...)## S3 method for class 'amce' tidy(x, alpha = 0.05, ...)
x |
A fit from |
alpha |
The significance level for the confidence intervals. Defaults to 0.05, which gives 95 percent intervals. |
... |
Not used. |
See https://book.declaredesign.org/experimental-descriptive.html#conjoint-experiments
A data frame with columns attribute, level, estimate, std.error, statistic, p.value, conf.low, and conf.high.
library(cjoint) data(immigrationconjoint) data(immigrationdesign) # Run AMCE estimator using all attributes in the design results <- amce(Chosen_Immigrant ~ Gender + Education + `Language Skills` + `Country of Origin` + Job + `Job Experience` + `Job Plans` + `Reason for Application` + `Prior Entry`, data = immigrationconjoint, cluster = TRUE, respondent.id = "CaseID", design = immigrationdesign) tidy(results)library(cjoint) data(immigrationconjoint) data(immigrationdesign) # Run AMCE estimator using all attributes in the design results <- amce(Chosen_Immigrant ~ Gender + Education + `Language Skills` + `Country of Origin` + Job + `Job Experience` + `Job Plans` + `Reason for Application` + `Prior Entry`, data = immigrationconjoint, cluster = TRUE, respondent.id = "CaseID", design = immigrationdesign) tidy(results)
A generics::tidy() method for fits from rdrobust::rdrobust(), so they can be used in declare_estimator(). Use rdrobust_helper() to fit the model.
## S3 method for class 'rdrobust' tidy(x, ...)## S3 method for class 'rdrobust' tidy(x, ...)
x |
A fit from |
... |
Not used. |
See https://book.declaredesign.org/observational-causal.html#regression-discontinuity-designs
A data frame with one row per estimate rdrobust reports (conventional, bias-corrected, and robust) and columns term, estimate, std.error, statistic, p.value, conf.low, conf.high, and cutoff.
set.seed(343) rd <- data.frame(X = runif(1000, min = -1, max = 1)) rd$Y <- 0.5 * (rd$X > 0) + rd$X + rnorm(1000, sd = 0.5) fit <- rdrobust_helper(rd, y = Y, x = X, c = 0) tidy(fit)set.seed(343) rd <- data.frame(X = runif(1000, min = -1, max = 1)) rd$Y <- 0.5 * (rd$X > 0) + rd$X + rnorm(1000, sd = 0.5) fit <- rdrobust_helper(rd, y = Y, x = X, c = 0) tidy(fit)