User’s guide

The list below details the arguments for interference, the primary function in inferference. Special attention should be given to the propensity_integrand and formula arguments.

• formula: formula used to define the causal model. formula has a minimum of 4 parts, separated by | and ~ in a specific structure: outcome | exposure ~ covariates | group. The order matters, and the pipes (|) split the data frame into corresponding pieces . The exposure ~ covariates piece is passed as a single formula to the chosen model_method (defined below) used to estimate or fix propensity parameters.
  • The following includes a random effect for the group: outcome | exposure ~ covariates + (1|group) | group. In this instance, the group variable appears twice.
  • If the study design includes a “participation” variable (as in both examples below), this is easily added to the formula: outcome | exposure | participation ~ covariates | group.
• propensity_integrand: a function, which may be created by the user, used to compute the IP weights. This defaults to the function logit_integrand(), which calculates the product of inverse logits for individuals in a group: $\prod_{j = 1}^{n_i} \{r h_{ij}(b_i)\}^{A_{ij}}\{1 - r h_{ij}(b_i)\}^{1 - A_{ij}} f_b(b_i; \theta_s)$, where h_ij(b_i) = logit⁻¹(X_ijθ_x + b_i), b_i is a group-level random effect, f_b is a N(0, θ_s) density, and r is a known constant. In an observational study typically r = 1. The examples below include individual randomized experiments in which case r denotes the randomization probability among trial participants. logit_integrand() is the integrand of , logit_integrand() ignores the random effect. IP weights are computed by numerically integrating propensity_integrand over the random effect distribution using by stats::integrate() to which arguments may be passed via ... (see below). The default logit_integrand() also takes the following argument that can be passed via the ... argument in interference():
  • randomization: a scalar. This is the r in the formula just above. It defaults to 1 in the case that a participation vector is not included. The vaccine study example below demonstrates use of this argument.
• loglihood_integrand: a function, which may be created by the user, that defines the log likelihood of the propensity score model. This should generally be the same function as `propensity_integrand}, which is the default.
• allocations: a vector of values in [0, 1]. Increasing the number of elements of the allocation vector increases computation time; however, a larger number of allocations will make plotted effect estimates smoother. A minimum of two allocations is required.
• data: the analysis data frame. This must include all the variables defined in the `formula}.
• model_method: the method used to estimate or set the propensity model parameters. Must be one of 'glm', 'glmer', or 'oracle'. For a fixed effects only model use 'glm', or to include random effects use lme4’s 'glmer' (Bates et al. 2014). logit_integrand only supports a single random effect for the grouping variable, corresponding to b_i. When the propensity parameters are known (as in simulations) or if estimating parameters for the propensity model outside of interference, use the 'oracle' option. See model_options for details on how to pass the oracle parameters. Defaults to 'glmer'}.
• model_options: a list of options passed to the function in model_method. Defaults to list(family = binomial(link = 'logit')). When model_method = 'oracle', the list must have two elements, fixed.effects and random.effects. If the model does not include random effects, set random.effects = NULL.
• causal_estimation_method: currently only supports and defaults to 'ipw'.
• causal_estmation_options: a list with a single item variance_estimation, which is either 'naive' or 'robust'. Defaults to 'robust'.
• conf.level: level for confidence intervals. Defaults to 0.95.
• rescale.factor: a scalar multiplication factor by which to rescale outcomes and effects. Defaults to 1.
• integrate_allocation: indicator of whether the integrand function uses the allocation parameter. Defaults to TRUE.
• ...: used to pass additional arguments to internal functions such as numDeriv::grad() or stats::integrate(). Arguments can also be passed to the propensity_integrand and loglihood_integrand functions.

The interference object

An interference() call results in an S3 object of class interference which contains:

• estimates: a data frame of causal effect estimates;
• models$propensity_model: the glm or glmer object;
• summary: a list of objects summarizing the causal model such as the number of groups, number of allocations, and the formula used in the interference call;
• weights: (# of groups) × (# of allocations) matrix of group-level weights: $$ w_{i, k} = \frac{\pi_i(\mathbf{A}_i; \alpha_k)}{f_{\mathbf{A}_i|\mathbf{X}_i}(\mathbf{A}_i|\mathbf{X}_i; \hat{\boldsymbol{\theta}})}. $$

If variance_estimation = 'robust', then the object also includes:

• weightd: (# of groups) × (# of allocations) × (# of parameters) array of weights computed using derivatives of the propensity function with respect to each parameter;
• scores: (# of groups) × (# of parameters) matrix of derivatives of the log likelihood.

Utility functions

The package includes tools to extract effect estimates of interest from the S3 object. The functions direct_effect, indirect_effect (or ie), total_effect (or te), and overall_effect (or oe) select appropriate records from the estimates data frame in the interference object.

Plotting effect estimates

Plots of effect estimates over a range of α levels may be helpful in summarizing results. Perez-Heydrich et al. (2014) present several such graphical displays. Here we demonstrate how to generate similar plots of effect estimates using inferference.

First, we estimate the effects over a dense sequence of allocations so that lines will be smooth.

> example2 <- interference( formula = Y | A | B ~ X1 + X2 + (1|group) | group, 
+     allocations = seq(.2, .8, by = .1), 
+     data = vaccinesim, randomization = 2/3, method = 'simple')

In the plots, we present the direct and indirect effect estimates over this range of allocations. For direct effects, a simple scatterplot showing the point-wise confidence intervals suffices. One approach with indirect effects fixes α and plots estimates over a range of α′, whereas a contour plot displays all pairwise (α, α′) comparisons over a range of allocation strategies.

> deff <- direct_effect(example2)
> x <- deff$alpha1
> y <- as.numeric(deff$estimate)
> u <- as.numeric(deff$conf.high)
> l <- as.numeric(deff$conf.low)
> plot(c(min(x), max(x)), c(-.15, .25), type = 'n', bty = 'l',
+      xlab = expression(alpha), ylab = '' )
> title(ylab = expression(widehat(DE) * "(" * alpha * ")"),
+       line = 2)
> polygon(c(x, rev(x)), c(u, rev(l)), col = 'skyblue', border = NA)
> lines(x, y, cex = 2)

> ieff.4 <- ie(example2, allocation1 = .4)
> x <- ieff.4$alpha2
> y <- as.numeric(ieff.4$estimate)
> u <- as.numeric(ieff.4$conf.high)
> l <- as.numeric(ieff.4$conf.low)
> plot(c(min(x), max(x)),c(-.15, .25), type = 'n', bty = 'l',
+      xlab = expression(alpha * "'"), ylab = '')
> title(ylab = expression(widehat(IE) * "(" * 0.4 * "," * alpha * "'" * ")"),
+       line = 2)
> polygon(c(x, rev(x)), c(u, rev(l)), col = 'skyblue', border = NA)
> lines(x, y, cex = 2)

> ieff <- subset(example2[["estimates"]], effect == 'indirect')
> x <- sort(unique(ieff$alpha1))
> y <- sort(unique(ieff$alpha2))
> z <- xtabs(estimate ~ alpha1 + alpha2, data= ieff)
> contour(x, y, z, xlab = expression(alpha), 
+         ylab = expression(alpha * "'"), bty = 'l')

Diagnostics

IPW estimators are known to be unstable if the weights range greatly. The package includes a basic utility to check the performance of the group-level weights, w_i, k, for multiple allocations. The function diagnose_weights plots histograms of weights for chosen allocation levels. If the allocations argument is left NULL, the function plots histograms for five allocation levels used the interference call. The plot below shows the resulting histogram for a single allocation. The analyst should examine groups with extreme weights, which may unduly influence population-level estimates.

> diagnose_weights(example2, allocations = .5, breaks = 30)

Household-level propensity

Unlike the vaccine study example, in this data set randomization occurred at the group level but individual level treatment was not randomized. With only two subjects, A_i = (A_i1, A_i2) is the treatment allocation for group i and X_i = (X_i1, X_i2) is the matrix of individuals’ covariate matrices for group i. Let B_i = (B_i1, B_i2) be indicators of being reached by a canvasser in group i. Since we only consider households where someone answered the door, B_i ∈ {(1, 0), (0, 1)} and Pr (B_i1 = 1|X_i) + Pr (B_i2 = 1|X_i) = 1. Let h_ij = Pr (B_ij = 1|X_i; θ) = logit⁻¹(X_iθ_x). Let G_i ∈ {0, 1} be the indicator that group i is randomized to treatment (1) or control (0). By design, Pr (G_i = 1) = 0.5. Since Pr (A_i1|X_i; θ) = 1 − Pr (A_i2|X_i; θ), Pr (A_i|X_i; θ) can arbitrarily be defined in terms of either household member. By convention we use the first subject (subject one) of each group found in the dataset. Among treated groups, the probability of subject one being treated is the probability that a canvasser reached subject one. That is, Pr (A_i1|G_i = 1, X_i; θ) = h_i1^A_i1(1 − h_i1)^{1 − A_i1}. Thus, the group-level propensity can be expressed:

Thus, h_i1 is sufficient to determine the group-level propensity. If we know whether or not the first subject was reached by a canvasser, then we know if the second was. Therefore, we can estimate parameters for h_i1 with a dataset that includes only subject one from each group. To do this, we must estimate the parameters outside of inferference and use model_method = 'oracle'. We include centered age (in decades) in the propensity model for demonstration purposes.

> voters1 <- do.call(rbind, by(voters, voters[, 'family'], function(x) x[1, ]))
> coef.voters <- coef(glm(reached ~ c_age, data = voters1, 
+                      family = binomial(link = 'logit')))

Coding the propensity function

Custom propensity_integrand and loglihood_integrand functions must have at least one argument:

• b: the first argument is the variable for which the integrate function computes the integral. As in this example, the function can be written so that the integral evaluates to 1 and has no effect.

For example, the following function will fix the group-level propensity to 0.5 for all groups when variance_estimation = "naive":

> fixed_propensity <- function(b){
+   return(0.5 * dnorm(b))
+ }

For more realistic models, additional arguments may be passed to the custom function:

• X: the covariate matrix (determined by the `formula}) for the ith group
• A: the vector of treatment indicators for the ith group
• parameters: vector of estimated parameters from the model_method * allocation: the allocation level for which the propensity is currently being estimated
• ...: other arguments can be passed via the ellipsis in interference Now we have the pieces to write the propensity function for the voting example.

> household_propensity <- function(b, X, A, 
+                                  parameters, 
+                                  group.randomization = .5){
+   if(!is.matrix(X)){
+     X <- as.matrix(X)
+   }   
+   if(sum(A) == 0){ 
+     pr <-  group.randomization
+   } else { 
+     X.1 <- X[1 ,]; A.1 <- A[1] 
+     h   <- plogis(X.1 %*% parameters)
+     pr  <-  group.randomization * dbinom(A.1, 1, h)
+   }   
+   out <- pr * dnorm(b) 
+   out
+ }

Evidence of a peer influence effect

The influence of the door opener on the non-door opener’s voting behavior corresponds to an indirect effect. Though the Bernoulli-type parametrization of the estimands allows us to look at a range of allocations, $\widehat{IE}(0.5, 0)$ makes the sensible comparison between a world where individuals receive a voting message with probability 0.5 to a world where individuals have zero probability of receiving the voting message.

> example3 <- interference(
+   formula = voted02p | treated | reached ~ c_age | family,
+   propensity_integrand = 'household_propensity',
+   data = voters,
+   model_method = 'oracle', 
+   model_options = list(fixed.effects = coef.voters, random.effects = NULL),
+   allocations   = c(0, .5),
+   integrate_allocation = FALSE,
+   causal_estimation_options = list(variance_estimation = 'robust'),
+   conf.level = .9)

> ie(example3, .5, 0)[ , c('estimate', 'conf.low', 'conf.high')]

##     estimate    conf.low  conf.high
## 1 0.03151501 0.002290586 0.06073944

The point estimate suggests an individual receiving the voting encouragement increases the voting likelihood of the other household member by 3.2%. The 90% confidence interval excludes zero, indicating a significant indirect effect corroborating the analysis in Nickerson (2008).

For comparison, suppose that a flip of a fair coin determined which registered voter opened the door. We exclude age as a covariate and instead set h_i1 = 0.5. Here we assume to know the propensity score, so we use variance_estimation = 'naive'.

> example4 <- interference(
+   formula = voted02p | treated | reached ~ 1 | family,
+   propensity_integrand = 'household_propensity',
+   data = voters,
+   model_method = 'oracle',
+   model_options = list(fixed.effects = 0, random.effects = NULL),
+   allocations   = c(0, .5),
+   integrate_allocation = FALSE,
+   causal_estimation_options = list(variance_estimation = 'naive'),
+   conf.level = .9)

> ie(example4, .5, 0)[ , c('estimate', 'conf.low', 'conf.high')]

##     estimate    conf.low  conf.high
## 1 0.03144654 0.002209019 0.06068406

Examining the group-level weights may help diagnose coding errors in the propensity score function. In the case of a fixed probability as in `example4}, the propensity weights are easily computed by hand. For example, for α = 0.5,

$$ w_i = \frac{\pi(\mathbf{A}_i; 0.5)}{f_{\mathbf{A}_i | \mathbf{X}_i}(\mathbf{A}_i | \mathbf{X}_i; \theta = 0)} = \begin{cases} \frac{0.5^2}{.5} = .5 & \text{ if } G_i = 0 \\ \frac{0.5^2}{.5^2} = 1 & \text{ if } G_i = 1\\ \end{cases}, $$

which we can confirm the software computed.

> G <- tapply(voters[1:12, 'treated'], voters[1:12, 'family'], sum)
> W <- head(example4[["weights"]])[, 2]
> cbind(G, W)

##    G   W
## 2  1 1.0
## 4  0 0.5
## 5  0 0.5
## 6  0 0.5
## 9  1 1.0
## 10 0 0.5