In bsaul/inferference: Methods for Causal Inference with Interference
JSS article
For a complete overview of the software, please refer to the Journal of Statistical Software paper: "A Recipe for inferference: Start with Causal Inference. Add Interference. Mix Well with R."
Using inferference
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 \citep{formula2010}. 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) = \mbox{logit}^{-1}(\mathbf{X}_{ij}\theta_x + b_i)$, $b_i$ is a group-level random effect, $f_b$ is a $N(0, \theta_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 \eqref{eq:propensitywhere $h_{ij}(b_i)$ is scaled by a constant $r$ term. If no random effect is included in theformula}, 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' [@lme4]. 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) $\times$ (# 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) $\times$ (# of allocations) $\times$ (# of parameters) array of weights computed using derivatives of the propensity function with respect to each parameter;scores: (# of groups) $\times$ (# 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.
Example: Vaccine Study
This section illustrates the use of inferference with an example drawn from vaccine research. The package includes a single dataset based on the same set of parameters used in the simulation study by @carolina2014. The vaccinesim dataset consists of 3000 units in 250 groups and contains two covariates (X1 = age in decades and X2 = distance to river), a vaccination indicator (A), a participation indicator (B), a binary outcome (Y) indicating cholera infection (1 yes, 0 no), and the unit's group.
library(knitr)
opts_chunk$set(prompt=TRUE)
options(continue ="+ ")
library(inferference)
head(vaccinesim)
Like the original study [@ali2005] that inspired the simulation, individuals were randomized to vaccine with a known probability of $2/3$, but subjects could opt to not participate in the trial. In essence, there are both experimental and observational aspects to the data. The interference function handles this design when logit_integrand's randomization argument is used and a participation variable is included in the formula.
example1 <- interference(
formula = Y | A | B ~ X1 + X2 + (1|group) | group,
allocations = c(.3, .45, .6),
data = vaccinesim,
randomization = 2/3,
method = 'simple')
The only arguments required for interference to run are formula, allocations, and data. When using the `robust' method (the default) to compute the variance, the internal workings call numDeriv::grad [@numDeriv2012] and stats::integrate frequently. The option method = 'simple' greatly speeds up the numDeriv::grad function. For more accurate derivatives, leave out this option. See ?numDeriv::grad for more options.
The print.interference function provides an overview of the causal effect estimates, estimated standard errors, and Wald-type confidence intervals. In the output, alpha1 and alpha2 refer to $\alpha$ and $\alpha'$, while trt1 and trt2 refer to $a_1$ and $a_2$, respectively.
print(example1)
The utility functions return selected effect estimates.
direct_effect(example1, .3)
ie(example1, .3)
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)
Example: Voting Experiment
The preceding example used the default logit_integrand function to define the group-level propensities. The following example demonstrates how to customize the propensity score function.
@nickerson2008 reported an experiment on voter behavior to examine peer-to-peer indirect effects on voting participation. The experiment randomized households with only two registered voters in Denver and Minneapolis to receive one of three assignments: voting encouragement, recycling encouragement, or nothing. Canvassers knocked on doors of households randomized to the voting or recycling groups a week before the 2002 primary. If a registered voter answered the door, the canvassers delivered a scripted message about voting (treatment) or recycling (control). The researchers used voter turnout records to determine if each member of the household then voted in the election. Nickerson was interested in the potential spillover effect of the voting encouragement to the untreated individual via the treated individual. From analysis of the observed data, he concluded there was a "secondary effect" where the household members not contacted by the canvassers voted more often in the treatment groups compared to the control groups.
The dataset voters contains information for 3861 households, 2549 in Denver and 1312 in Minneapolis, including covariates such as age, gender, previous voting history, and party affiliation. Our estimand of interest involves average voting outcomes when households receive voting encouragement compared to when household receive the recycling message, hence we exclude households not contacted by canvassers. We also exclude the single household where a canvasser appears to have contacted both registered voters.
voters <- within(voters, {
treated = (treatment == 1 & reached == 1) * 1
c_age = (age - mean(age))/10
})
reach_cnt <- tapply(voters$reached, voters$family, sum)
voters <- voters[!(voters$family %in% names(reach_cnt[reach_cnt > 1])), ]
voters <- voters[voters$hsecontact == 1, ]
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, $\mathbf{A}i = (A{i1}, A_{i2})$ is the treatment allocation for group $i$ and $\mathbf{X}i = (\mathbf{X}{i1}, \mathbf{X}{i2})$ is the matrix of individuals' covariate matrices for group $i$. Let $\mathbf{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, $\mathbf{B}i \in {(1, 0), (0, 1)}$ and $\Pr(B{i1} = 1| \mathbf{X}{i}) + \Pr(B{i2} = 1| \mathbf{X}{i}) = 1$. Let $h{ij} = \Pr(B_{ij} = 1| \mathbf{X}{i}; \theta) = \mbox{logit}^{-1}(\mathbf{X}{i}\theta_x)$. Let $G_i \in {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} | \mathbf{X}i; \theta) = 1 - \Pr(A{i2} | \mathbf{X}i; \theta)$, $\Pr(\mathbf{A}_i | \mathbf{X}_i; \theta)$ 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, \mathbf{X}i; \theta) = h{i1}^{A_{i1}}(1 - h_{i1})^{1 - A_{i1}}$. Thus, the group-level propensity can be expressed:
\begin{align}
\Pr(\mathbf{A}i | \mathbf{X}_i ; \theta) &= \Pr(\mathbf{A}_i | G_i = 1, \mathbf{X}_i ; \theta)\Pr(G_i = 1) + \Pr(\mathbf{A}_i | G_i = 0, \mathbf{X}_i ; \theta)\Pr(G_i = 0) \
&= 0.5 { \Pr(\mathbf{A}_i | G_i = 1, \mathbf{X}_i ; \theta) + \Pr(\mathbf{A}_i | G_i = 0, \mathbf{X}_i ; \theta) } \
&= \begin{cases} .5 & \text{if $A{i1} = 0, A_{i2} = 0$ and $G_i = 0$} \
.5 h_{i1}^{A_{i1}}(1 - h_{i1})^{1 - A_{i1}} & \text{if ($A_{i1} = 0, A_{i2} = 1$ or $A_{i1} = 1, A_{i2} = 0$) and $G_i = 1$} \
0 & \text{otherwise} \end{cases} \
\end{align}
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 $i$th groupA: the vector of treatment indicators for the $i$th groupparameters: 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')]
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 @nickerson2008.
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')]
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 $\alpha = 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)
Computational Issues with IPW Estimators
We show in this section how computation of the group-level weights may affect estimation as the number of individuals in groups grows. To illustrate, consider the IPW estimator of the overall effect, which weights individual outcomes in group $i$ with:
[
w_{1i} = \frac{\pi(\textbf{A}i ; \alpha)}{f{\mathbf{A}i | \mathbf{X}_i}(\mathbf{A}_i | \mathbf{X}_i; \hat{\boldsymbol{\theta}})} = \frac{\prod{j=1}^{n_i} \alpha^{A_{ij}} (1 - \alpha)^{1 - A_{ij}}}{ \int \prod_{j=1}^{n_i} h_{ij}^{A_{ij}} (1 - h_{ij})^{1 - A_{ij}} f_b (b_i ; \hat{\theta}_s) db_i } ]
or equivalently,
[
w_{2i} = \left{ \int \prod_{j=1}^{n_i} \left(\frac{h_{ij}}{\alpha}\right)^{A_{ij}} \left(\frac{1 - h_{ij}}{1 - \alpha}\right)^{1 - A_{ij}} f_b (b_i ; \hat{\theta}_s) db_i \right}^{-1} ]
or,
[
w_{3i} = \left{ \int \exp\left[ \sum_{j=1}^{n_i} \left{ A_{ij} \log \left( \frac{h_{ij}}{\alpha}\right) + (1 - A_{ij})\log \left(\frac{1 - h_{ij}}{1 - \alpha}\right) \right} \right] f_b (b_i ; \hat{\theta}_s) db_i \right}^{-1}.
]
While mathematically equivalent, these weights may be computationally dissimilar. In the case of $w_{1i}$, the product term within the integral entails multiplying probabilities and thus will tend to 0 as $n_i$ increases, causing the denominator of $w_{1i}$ to get arbitrarily large. In contrast, the product term in $w_{2i}$ entails multiplying values which may be less than or greater than 1 and thus tends to be less susceptible to numerical instability. Summing $\log(h_{ij}/\alpha)$ or $\log(1 - h_{ij})/(1 - \alpha))$ and then exponentiating the result may provide greater numerical stability. Internally, inferference uses $w_{3i}$.
When group sizes are small, the differences between these weights tend to be infinitesimal, but as group sizes grow the differences become important. To be specific, consider the code below comparing $w_{1i}$, $w_{2i}$, and $w_{3i}$ for increasing group sizes where $\alpha = 0.5$, all units are treated, there is no random effect, and $h_{ij}$ is fixed at $0.5$.
compare_weights <- function(n, alpha = .5, h = .5){
pi <- rep(alpha, n)
PrA <- rep(h, n)
c(w1 = prod(pi)/prod(PrA),
w2 = 1/prod(PrA/pi),
w3 = 1/exp(sum(log(PrA/pi))))
}
n <- c(50, 100, 500, 1074, 1075, 10000)
cbind(n, t(sapply(n, FUN = compare_weights)))
For group sizes up to $1074$ there is no difference, but when $n$ reaches $1075$, $w_{1i}$ returns NaN while $w_{2i}$ and $w_{3i}$ correctly return 1.
@carolina2014 used $w_{1i}$ to calculate weights, but 15 groups in their analysis had over 1000 subjects. These groups had missing values for weights for all the values of $\alpha$ considered and were excluded from computing the average IPW estimate. Rather than computing the average IPW across 700 groups, they inadvertently took the average across 685 groups. Correcting the estimates by using $w_{2i}$ or $w_{3i}$ did not alter the conclusions in this case, but analysts should be aware of this issue when dealing with large groups.
References
bsaul/inferference documentation built on April 21, 2021, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
JSS article
For a complete overview of the software, please refer to the Journal of Statistical Software paper: "A Recipe for inferference: Start with Causal Inference. Add Interference. Mix Well with R."
Using inferference
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.formulahas 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 \citep{formula2010}. Theexposure ~ covariatespiece is passed as a single formula to the chosenmodel_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.
- The following includes a random effect for the group:
propensity_integrand: a function, which may be created by the user, used to compute the IP weights. This defaults to the functionlogit_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) = \mbox{logit}^{-1}(\mathbf{X}_{ij}\theta_x + b_i)$, $b_i$ is a group-level random effect, $f_b$ is a $N(0, \theta_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 \eqref{eq:propensitywhere $h_{ij}(b_i)$ is scaled by a constant $r$ term. If no random effect is included in theformula},logit_integrand()ignores the random effect. IP weights are computed by numerically integratingpropensity_integrandover the random effect distribution using bystats::integrate()to which arguments may be passed via...(see below). The defaultlogit_integrand()also takes the following argument that can be passed via the...argument ininterference():randomization: a scalar. This is the $r$ in the formula just above. It defaults to 1 in the case that aparticipationvector 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 uselme4's'glmer'[@lme4].logit_integrandonly 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 ofinterference, use the'oracle'option. Seemodel_optionsfor details on how to pass the oracle parameters. Defaults to'glmer'}.model_options: a list of options passed to the function inmodel_method. Defaults tolist(family = binomial(link = 'logit')). Whenmodel_method = 'oracle', the list must have two elements,fixed.effectsandrandom.effects. If the model does not include random effects, setrandom.effects = NULL.causal_estimation_method: currently only supports and defaults to'ipw'.causal_estmation_options: a list with a single itemvariance_estimation, which is either'naive'or'robust'. Defaults to'robust'.conf.level: level for confidence intervals. Defaults to0.95.rescale.factor: a scalar multiplication factor by which to rescale outcomes and effects. Defaults to1.integrate_allocation: indicator of whether the integrand function uses the allocation parameter. Defaults to TRUE....: used to pass additional arguments to internal functions such asnumDeriv::grad()orstats::integrate(). Arguments can also be passed to thepropensity_integrandandloglihood_integrandfunctions.
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: theglmorglmerobject;summary: a list of objects summarizing the causal model such as the number of groups, number of allocations, and the formula used in theinterferencecall;weights: (# of groups) $\times$ (# 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) $\times$ (# of allocations) $\times$ (# of parameters) array of weights computed using derivatives of the propensity function with respect to each parameter;scores: (# of groups) $\times$ (# 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.
Example: Vaccine Study
This section illustrates the use of inferference with an example drawn from vaccine research. The package includes a single dataset based on the same set of parameters used in the simulation study by @carolina2014. The vaccinesim dataset consists of 3000 units in 250 groups and contains two covariates (X1 = age in decades and X2 = distance to river), a vaccination indicator (A), a participation indicator (B), a binary outcome (Y) indicating cholera infection (1 yes, 0 no), and the unit's group.
library(knitr) opts_chunk$set(prompt=TRUE) options(continue ="+ ")
library(inferference) head(vaccinesim)
Like the original study [@ali2005] that inspired the simulation, individuals were randomized to vaccine with a known probability of $2/3$, but subjects could opt to not participate in the trial. In essence, there are both experimental and observational aspects to the data. The interference function handles this design when logit_integrand's randomization argument is used and a participation variable is included in the formula.
example1 <- interference( formula = Y | A | B ~ X1 + X2 + (1|group) | group, allocations = c(.3, .45, .6), data = vaccinesim, randomization = 2/3, method = 'simple')
The only arguments required for interference to run are formula, allocations, and data. When using the `robust' method (the default) to compute the variance, the internal workings call numDeriv::grad [@numDeriv2012] and stats::integrate frequently. The option method = 'simple' greatly speeds up the numDeriv::grad function. For more accurate derivatives, leave out this option. See ?numDeriv::grad for more options.
The print.interference function provides an overview of the causal effect estimates, estimated standard errors, and Wald-type confidence intervals. In the output, alpha1 and alpha2 refer to $\alpha$ and $\alpha'$, while trt1 and trt2 refer to $a_1$ and $a_2$, respectively.
print(example1)
The utility functions return selected effect estimates.
direct_effect(example1, .3) ie(example1, .3)
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)
Example: Voting Experiment
The preceding example used the default logit_integrand function to define the group-level propensities. The following example demonstrates how to customize the propensity score function.
@nickerson2008 reported an experiment on voter behavior to examine peer-to-peer indirect effects on voting participation. The experiment randomized households with only two registered voters in Denver and Minneapolis to receive one of three assignments: voting encouragement, recycling encouragement, or nothing. Canvassers knocked on doors of households randomized to the voting or recycling groups a week before the 2002 primary. If a registered voter answered the door, the canvassers delivered a scripted message about voting (treatment) or recycling (control). The researchers used voter turnout records to determine if each member of the household then voted in the election. Nickerson was interested in the potential spillover effect of the voting encouragement to the untreated individual via the treated individual. From analysis of the observed data, he concluded there was a "secondary effect" where the household members not contacted by the canvassers voted more often in the treatment groups compared to the control groups.
The dataset voters contains information for 3861 households, 2549 in Denver and 1312 in Minneapolis, including covariates such as age, gender, previous voting history, and party affiliation. Our estimand of interest involves average voting outcomes when households receive voting encouragement compared to when household receive the recycling message, hence we exclude households not contacted by canvassers. We also exclude the single household where a canvasser appears to have contacted both registered voters.
voters <- within(voters, { treated = (treatment == 1 & reached == 1) * 1 c_age = (age - mean(age))/10 }) reach_cnt <- tapply(voters$reached, voters$family, sum) voters <- voters[!(voters$family %in% names(reach_cnt[reach_cnt > 1])), ] voters <- voters[voters$hsecontact == 1, ]
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, $\mathbf{A}i = (A{i1}, A_{i2})$ is the treatment allocation for group $i$ and $\mathbf{X}i = (\mathbf{X}{i1}, \mathbf{X}{i2})$ is the matrix of individuals' covariate matrices for group $i$. Let $\mathbf{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, $\mathbf{B}i \in {(1, 0), (0, 1)}$ and $\Pr(B{i1} = 1| \mathbf{X}{i}) + \Pr(B{i2} = 1| \mathbf{X}{i}) = 1$. Let $h{ij} = \Pr(B_{ij} = 1| \mathbf{X}{i}; \theta) = \mbox{logit}^{-1}(\mathbf{X}{i}\theta_x)$. Let $G_i \in {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} | \mathbf{X}i; \theta) = 1 - \Pr(A{i2} | \mathbf{X}i; \theta)$, $\Pr(\mathbf{A}_i | \mathbf{X}_i; \theta)$ 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, \mathbf{X}i; \theta) = h{i1}^{A_{i1}}(1 - h_{i1})^{1 - A_{i1}}$. Thus, the group-level propensity can be expressed:
\begin{align} \Pr(\mathbf{A}i | \mathbf{X}_i ; \theta) &= \Pr(\mathbf{A}_i | G_i = 1, \mathbf{X}_i ; \theta)\Pr(G_i = 1) + \Pr(\mathbf{A}_i | G_i = 0, \mathbf{X}_i ; \theta)\Pr(G_i = 0) \ &= 0.5 { \Pr(\mathbf{A}_i | G_i = 1, \mathbf{X}_i ; \theta) + \Pr(\mathbf{A}_i | G_i = 0, \mathbf{X}_i ; \theta) } \ &= \begin{cases} .5 & \text{if $A{i1} = 0, A_{i2} = 0$ and $G_i = 0$} \ .5 h_{i1}^{A_{i1}}(1 - h_{i1})^{1 - A_{i1}} & \text{if ($A_{i1} = 0, A_{i2} = 1$ or $A_{i1} = 1, A_{i2} = 0$) and $G_i = 1$} \ 0 & \text{otherwise} \end{cases} \ \end{align}
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 theintegratefunction 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 $i$th groupA: the vector of treatment indicators for the $i$th groupparameters: vector of estimated parameters from themodel_method*allocation: the allocation level for which the propensity is currently being estimated...: other arguments can be passed via the ellipsis ininterferenceNow 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')]
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 @nickerson2008.
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')]
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 $\alpha = 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)
Computational Issues with IPW Estimators
We show in this section how computation of the group-level weights may affect estimation as the number of individuals in groups grows. To illustrate, consider the IPW estimator of the overall effect, which weights individual outcomes in group $i$ with:
[ w_{1i} = \frac{\pi(\textbf{A}i ; \alpha)}{f{\mathbf{A}i | \mathbf{X}_i}(\mathbf{A}_i | \mathbf{X}_i; \hat{\boldsymbol{\theta}})} = \frac{\prod{j=1}^{n_i} \alpha^{A_{ij}} (1 - \alpha)^{1 - A_{ij}}}{ \int \prod_{j=1}^{n_i} h_{ij}^{A_{ij}} (1 - h_{ij})^{1 - A_{ij}} f_b (b_i ; \hat{\theta}_s) db_i } ]
or equivalently, [ w_{2i} = \left{ \int \prod_{j=1}^{n_i} \left(\frac{h_{ij}}{\alpha}\right)^{A_{ij}} \left(\frac{1 - h_{ij}}{1 - \alpha}\right)^{1 - A_{ij}} f_b (b_i ; \hat{\theta}_s) db_i \right}^{-1} ]
or, [ w_{3i} = \left{ \int \exp\left[ \sum_{j=1}^{n_i} \left{ A_{ij} \log \left( \frac{h_{ij}}{\alpha}\right) + (1 - A_{ij})\log \left(\frac{1 - h_{ij}}{1 - \alpha}\right) \right} \right] f_b (b_i ; \hat{\theta}_s) db_i \right}^{-1}. ]
While mathematically equivalent, these weights may be computationally dissimilar. In the case of $w_{1i}$, the product term within the integral entails multiplying probabilities and thus will tend to 0 as $n_i$ increases, causing the denominator of $w_{1i}$ to get arbitrarily large. In contrast, the product term in $w_{2i}$ entails multiplying values which may be less than or greater than 1 and thus tends to be less susceptible to numerical instability. Summing $\log(h_{ij}/\alpha)$ or $\log(1 - h_{ij})/(1 - \alpha))$ and then exponentiating the result may provide greater numerical stability. Internally, inferference uses $w_{3i}$.
When group sizes are small, the differences between these weights tend to be infinitesimal, but as group sizes grow the differences become important. To be specific, consider the code below comparing $w_{1i}$, $w_{2i}$, and $w_{3i}$ for increasing group sizes where $\alpha = 0.5$, all units are treated, there is no random effect, and $h_{ij}$ is fixed at $0.5$.
compare_weights <- function(n, alpha = .5, h = .5){ pi <- rep(alpha, n) PrA <- rep(h, n) c(w1 = prod(pi)/prod(PrA), w2 = 1/prod(PrA/pi), w3 = 1/exp(sum(log(PrA/pi)))) } n <- c(50, 100, 500, 1074, 1075, 10000) cbind(n, t(sapply(n, FUN = compare_weights)))
For group sizes up to $1074$ there is no difference, but when $n$ reaches $1075$, $w_{1i}$ returns NaN while $w_{2i}$ and $w_{3i}$ correctly return 1.
@carolina2014 used $w_{1i}$ to calculate weights, but 15 groups in their analysis had over 1000 subjects. These groups had missing values for weights for all the values of $\alpha$ considered and were excluded from computing the average IPW estimate. Rather than computing the average IPW across 700 groups, they inadvertently took the average across 685 groups. Correcting the estimates by using $w_{2i}$ or $w_{3i}$ did not alter the conclusions in this case, but analysts should be aware of this issue when dealing with large groups.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.