Title:  Weighted Linear Fixed Effects Regression Models for Causal Inference 

Description:  Provides a computationally efficient way of fitting weighted linear fixed effects estimators for causal inference with various weighting schemes. Weighted linear fixed effects estimators can be used to estimate the average treatment effects under different identification strategies. This includes stratified randomized experiments, matching and stratification for observational studies, first differencing, and differenceindifferences. The package implements methods described in Imai and Kim (2017) "When should We Use Linear Fixed Effects Regression Models for Causal Inference with Longitudinal Data?", available at <https://imai.fas.harvard.edu/research/FEmatch.html>. 
Authors:  In Song Kim [aut, cre], Kosuke Imai [aut] 
Maintainer:  In Song Kim <[email protected]> 
License:  GPL(>= 2) 
Version:  1.9.1 
Built:  20240923 05:44:07 UTC 
Source:  https://github.com/insongkim/wfe 
pwfe
is used to fit weighted fixed effects model for causal
inference after transforming outcome variable based on estimated
propensity score. pwfe
also derives the regression weights for
different causal quantity of interest.
pwfe(formula, treat = "treat.name", outcome, data, pscore = NULL, unit.index, time.index = NULL, method = "unit", within.unit = TRUE, qoi = c("ate", "att"), estimator = NULL, C.it = NULL, White = TRUE, White.alpha = 0.05, hetero.se = TRUE, auto.se = TRUE, unbiased.se = FALSE, verbose = TRUE)
pwfe(formula, treat = "treat.name", outcome, data, pscore = NULL, unit.index, time.index = NULL, method = "unit", within.unit = TRUE, qoi = c("ate", "att"), estimator = NULL, C.it = NULL, White = TRUE, White.alpha = 0.05, hetero.se = TRUE, auto.se = TRUE, unbiased.se = FALSE, verbose = TRUE)
formula 
a symbolic description of the model for estimating propensity score. The formula should not include dummmies for fixed effects. The details of model specifications are given under ‘Details’. 
treat 
a character string indicating the name of treatment variable used in the models. The treatment should be binary indicator (integer with 0 for the control group and 1 for the treatment group). 
outcome 
a character string indicating the name of outcome variable. 
data 
data frame containing the variables in the model. 
pscore 
an optional character string indicating the name of estimated propensity score. Note that prespecified propensity score should be bounded away from zero and one. 
unit.index 
a character string indicating the name of unit variable used in the models. The index of unit should be factor. 
time.index 
a character string indicating the name of time variable used in the models. The index of time should be factor. 
method 
method for weighted fixed effects regression, either

within.unit 
a logical value indicating whether propensity score
is estimated within unit. The default is 
qoi 
one of 
estimator 
an optional character string 
C.it 
an optional nonnegative numeric vector specifying relative weights for each unit of analysis. 
White 
a logical value indicating whether White misspecification
statistics should be calculated. The default is 
White.alpha 
level of functional specification test. See White
(1980) and Imai . The default is 
hetero.se 
a logical value indicating whether heteroskedasticity
across units is allowed in calculating standard errors. The default
is 
auto.se 
a logical value indicating whether arbitrary
autocorrelation is allowed in calculating standard errors. The
default is 
unbiased.se 
logical. If 
verbose 
logical. If 
To fit the weighted unit (time) fixed effects model with propensity
score weighting, use the syntax for the formula, ~ x1 + x2
,
where x1
and x2
are unit (time) varying
covariates.
One can provide his/her own estimated pscore
which can be used
to transform the outcome varialbe. If so, one does not need to specify
formula
.
If pscore
is not provided, bayesglm
will be used to
estimate propensity scores. If within.unit = TRUE
, propensity
score will be separately estimated within time (unit) when
method
is unit
(time
). Otherwise, propensity
score will be estimated on entire data at once.
The estimated propensity scores will be used to transform the
outcome
variable as described in Imai and Kim (2018).
pwfe
calculates weights based on different underlying causal
quantity of interest: Average Treatment Effect (qoi = "ate"
) or
Average Treatment Effect for the Treated (qoi = "att"
).
One can further set estimating methods: FirstDifference
(estimator ="fd"
) or Differenceindifferences (estimator
= "did"
).
To specify different exante weights for each unit of analysis, use
nonnegative weights C.it
. For instance, using the survey
weights for C.it
enables the estimation fo the average
treatement effect for the target population.
pwfe
returns an object of class "pwfe", a list that contains the
components listed below.
The function summary
(i.e., summary.pwfe
) can be used to
obtain a table of the results.
coefficients 
a named vector of coefficients 
residuals 
the residuals, that is respons minus fitted values 
df 
the degree of freedom 
W 
weight matrix calculated from the model. Row and column indices can be found from unit.name, time.name. 
call 
the matched call 
causal 
causal quantity of interest 
estimator 
the estimating method 
unit.name 
a vector containing unique unit names 
unit.index 
a vector containing unique unit index number 
time.name 
a vector containing unique time names 
time.index 
a vector containing unique time index number 
method 
call of the method used 
vcov 
the variance covariance matrix 
White.alpha 
the alpha level for White specification test 
White.pvalue 
the pvalue for White specification test 
White.stat 
the White statistics 
x 
the design matrix 
y 
the response vector 
mf 
the model frame 
In Song Kim, Massachusetts Institute of Technology, [email protected] and Kosuke Imai, Princeton University, [email protected]
Imai, Kosuke and In Song Kim. (2018) “When Should We Use Unit Fixed Effects Regression Models for Causal Inference with Longitudinal Data?" American Journal of Political Science, Forthcoming.
Stock, James and Mark Watson. (2008) “HeteroskedasticityRobust Standard Errors for Fixed Effect Panel Data Regression” Econometrica, 76, 1.
White, Halbert. (1980) 'Using Least Squares to Approximate Unknown Regression Functions.” International Economic Review, 21, 1, 149–170.
wfe
for fitting weighted fixed effect models.
### NOTE: this example illustrates the use of wfe function with randomly ### generated panel data with arbitrary number of units and time. ## generate panel data with number of units = N, number of time = Time ## Not run: N < 10 # number of distinct units Time < 15 # number of distinct time ## generate treatment variable treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) ## make sure at least one observation is treated for each unit while ((sum(apply(treat, 2, mean) == 0) > 0)  (sum(apply(treat, 2, mean) == 1) > 0)  (sum(apply(treat, 1, mean) == 0) > 0)  (sum(apply(treat, 1, mean) == 1) > 0)) { treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) } treat.vec < c(treat) ## unit fixed effects alphai < rnorm(N, mean = apply(treat, 2, mean)) ## geneate two random covariates x1 < matrix(rnorm(N*Time, 0.5,1), ncol=N) x2 < matrix(rbeta(N*Time, 5,1), ncol=N) pscore < matrix(runif(N*Time, 0,1), ncol=N) x1.vec < c(x1) x2.vec < c(x2) pscore < c(pscore) ## generate outcome variable y < matrix(NA, ncol = N, nrow = Time) for (i in 1:N) { y[, i] < alphai[i] + treat[, i] + x1[,i] + x2[,i] + rnorm(Time) } y.vec < c(y) ## generate unit and time index unit.index < rep(1:N, each = Time) time.index < rep(1:Time, N) Data.str < as.data.frame(cbind(y.vec, treat.vec, unit.index, x1.vec, x2.vec)) colnames(Data.str) < c("y", "tr", "strata.id", "x1", "x2") Data.obs < as.data.frame(cbind(y.vec, treat.vec, unit.index, time.index, x1.vec, x2.vec, pscore)) colnames(Data.obs) < c("y", "tr", "unit", "time", "x1", "x2", "pscore") ############################################################ # Example 1: Stratified Randomized Experiments ############################################################ ## run the weighted fixed effect regression with strata fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. ### Average Treatment Effect ps.ate < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.str, unit.index = "strata.id", method = "unit", within.unit = TRUE, qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.ate) ### Average Treatment Effect for the Treated ps.att < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.str, unit.index = "strata.id", method = "unit", within.unit = TRUE, qoi = "att", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.att) ############################################################ # Example 2: Observational Studies with Unit Fixedeffects ############################################################ ## run the weighted fixed effect regression with unit fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. ### Average Treatment Effect ps.obs < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.obs, unit.index = "unit", time.index = "time", method = "unit", within.unit = TRUE, qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.obs) ## extracting weigths summary(ps.obs)$Weights ### Average Treatment Effect with Firstdifference ps.fd < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.obs, unit.index = "unit", time.index = "time", method = "unit", within.unit = TRUE, qoi = "ate", estimator = "fd", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.fd) ############################################################ # Example 3: Estimation with prespecified propensity score ############################################################ ### Average Treatment Effect with Prespecified Propensity Scores mod.ps < pwfe(treat = "tr", outcome = "y", data = Data.obs, pscore = "pscore", unit.index = "unit", time.index = "time", method = "unit", within.unit = TRUE, qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(mod.ps) ## End(Not run)
### NOTE: this example illustrates the use of wfe function with randomly ### generated panel data with arbitrary number of units and time. ## generate panel data with number of units = N, number of time = Time ## Not run: N < 10 # number of distinct units Time < 15 # number of distinct time ## generate treatment variable treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) ## make sure at least one observation is treated for each unit while ((sum(apply(treat, 2, mean) == 0) > 0)  (sum(apply(treat, 2, mean) == 1) > 0)  (sum(apply(treat, 1, mean) == 0) > 0)  (sum(apply(treat, 1, mean) == 1) > 0)) { treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) } treat.vec < c(treat) ## unit fixed effects alphai < rnorm(N, mean = apply(treat, 2, mean)) ## geneate two random covariates x1 < matrix(rnorm(N*Time, 0.5,1), ncol=N) x2 < matrix(rbeta(N*Time, 5,1), ncol=N) pscore < matrix(runif(N*Time, 0,1), ncol=N) x1.vec < c(x1) x2.vec < c(x2) pscore < c(pscore) ## generate outcome variable y < matrix(NA, ncol = N, nrow = Time) for (i in 1:N) { y[, i] < alphai[i] + treat[, i] + x1[,i] + x2[,i] + rnorm(Time) } y.vec < c(y) ## generate unit and time index unit.index < rep(1:N, each = Time) time.index < rep(1:Time, N) Data.str < as.data.frame(cbind(y.vec, treat.vec, unit.index, x1.vec, x2.vec)) colnames(Data.str) < c("y", "tr", "strata.id", "x1", "x2") Data.obs < as.data.frame(cbind(y.vec, treat.vec, unit.index, time.index, x1.vec, x2.vec, pscore)) colnames(Data.obs) < c("y", "tr", "unit", "time", "x1", "x2", "pscore") ############################################################ # Example 1: Stratified Randomized Experiments ############################################################ ## run the weighted fixed effect regression with strata fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. ### Average Treatment Effect ps.ate < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.str, unit.index = "strata.id", method = "unit", within.unit = TRUE, qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.ate) ### Average Treatment Effect for the Treated ps.att < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.str, unit.index = "strata.id", method = "unit", within.unit = TRUE, qoi = "att", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.att) ############################################################ # Example 2: Observational Studies with Unit Fixedeffects ############################################################ ## run the weighted fixed effect regression with unit fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. ### Average Treatment Effect ps.obs < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.obs, unit.index = "unit", time.index = "time", method = "unit", within.unit = TRUE, qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.obs) ## extracting weigths summary(ps.obs)$Weights ### Average Treatment Effect with Firstdifference ps.fd < pwfe(~ x1+x2, treat = "tr", outcome = "y", data = Data.obs, unit.index = "unit", time.index = "time", method = "unit", within.unit = TRUE, qoi = "ate", estimator = "fd", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(ps.fd) ############################################################ # Example 3: Estimation with prespecified propensity score ############################################################ ### Average Treatment Effect with Prespecified Propensity Scores mod.ps < pwfe(treat = "tr", outcome = "y", data = Data.obs, pscore = "pscore", unit.index = "unit", time.index = "time", method = "unit", within.unit = TRUE, qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(mod.ps) ## End(Not run)
wfe
is used to fit weighted fixed effects model for causal
inference. wfe
also derives the regression weights for
different causal quantity of interest.
wfe(formula, data, treat = "treat.name", unit.index, time.index = NULL, method = "unit", dyad1.index = NULL, dyad2.index = NULL, qoi = "ate", estimator = NULL, C.it = NULL, hetero.se = TRUE, auto.se = TRUE, dyad.se = FALSE, White = TRUE, White.alpha = 0.05, verbose = TRUE, unbiased.se = FALSE, unweighted = FALSE, store.wdm = FALSE, maxdev.did = NULL, tol = sqrt(.Machine$double.eps))
wfe(formula, data, treat = "treat.name", unit.index, time.index = NULL, method = "unit", dyad1.index = NULL, dyad2.index = NULL, qoi = "ate", estimator = NULL, C.it = NULL, hetero.se = TRUE, auto.se = TRUE, dyad.se = FALSE, White = TRUE, White.alpha = 0.05, verbose = TRUE, unbiased.se = FALSE, unweighted = FALSE, store.wdm = FALSE, maxdev.did = NULL, tol = sqrt(.Machine$double.eps))
formula 
a symbolic description of the model to be fitted. The formula should not include dummmies for fixed effects. The details of model specifications are given under ‘Details’. 
data 
data frame containing the variables in the model. 
treat 
a character string indicating the name of treatment variable used in the models. The treatment should be binary indicator (integer with 0 for the control group and 1 for the treatment group). 
unit.index 
a character string indicating the name of unit variable used in the models. The index of unit should be factor. 
time.index 
a character string indicating the name of time variable used in the models. The index of time should be factor. 
method 
method for weighted fixed effects regression, either

dyad1.index 
a character string indicating the variable name of first unit
of a given dyad. The default is 
dyad2.index 
a character string indicating the variable name of second unit
of a given dyad. The default is 
qoi 
one of 
estimator 
an optional character string indicating the
estimating method. One of 
C.it 
an optional nonnegative numeric vector specifying relative weights for each unit of analysis. If not specified, the weights will be calculated based on the estimator and quantity of interest. 
hetero.se 
a logical value indicating whether heteroskedasticity
across units is allowed in calculating standard errors. The default
is 
auto.se 
a logical value indicating whether arbitrary
autocorrelation is allowed in calculating standard errors. The
default is 
dyad.se 
a logical value indicating whether correlations across dyads exist. The
default is 
White 
a logical value indicating whether White misspecification
statistics should be calculated. The default is 
White.alpha 
level of functional specification test. See White
(1980) and Imai and Kim (2018). The default is 
verbose 
logical. If 
unbiased.se 
logical. If 
unweighted 
logical. If 
store.wdm 
logical. If 
maxdev.did 
an optional positive numeric value specifying the
maximum deviation in pretreatment outcome when 
tol 
a relative tolerance to detect zero singular values for generalized inverse. The default is sqrt(.Machine$double.eps) 
To fit the weighted unit (time) fixed effects model, use the syntax
for the formula, y ~ x1 + x2
, where y
is a dependent
variable and x1
and x2
are unit (time) varying
covariates.
wfe
calculates weights based on different underlying causal
quantity of interest: Average Treatment Effect (qoi = "ate"
) or
Average Treatment Effect for the Treated (qoi = "att"
).
One can further set estimating methods: FirstDifference
(estimator ="fd"
) or Differenceindifferences (estimator
= "did"
). For the twoway fixed effects model, set estimator
= "did"
To specify different exante weights for each unit of analysis, use
nonnegative weights C.it
. For instance, using the survey
weights for C.it
enables the estimation fo the average
treatement effect for the target population.
An object of class "wfe" contains vectors of unique unit(time) names and unique unit(time) indices.
wfe
returns an object of class "wfe", a list that contains the
components listed below.
The function summary
(i.e., summary.wfe
) can be used to
obtain a table of the results.
coefficients 
a named vector of coefficients 
residuals 
the residuals, that is respons minus fitted values 
df 
the degree of freedom 
W 
a dataframe containing unit and time indices along with the
weights used for the observation. If method= 
Num.nonzero 
Number of observations with nonzero weights 
call 
the matched call 
causal 
causal quantity of interest 
estimator 
the estimating method 
units 
a dataframe containing unit names used for 
times 
a dataframe containing time names used for 
method 
call of the method used 
vcov 
the variance covariance matrix 
White.alpha 
the alpha level for White specification test 
White.pvalue 
the pvalue for White specification test 
White.stat 
the White statistics 
X 
the design matrix 
Y 
the response vector 
X.wdm 
the demeaned design matrix 
Y.wdm 
the demeaned response vector 
mf 
the model frame where the last column is the weights used for the analysis 
In Song Kim, Massachusetts Institute of Technology, [email protected] and Kosuke Imai, Princeton University, [email protected]
Imai, Kosuke and In Song Kim. (2018) “When Should We Use Unit Fixed Effects Regression Models for Causal Inference with Longitudinal Data?" American Journal of Political Science, Forthcoming.
Aronow, Peter M., Cyrus Samii, and Valentina A. Assenova (2015) “Cluster–robust Variance Estimation for Dyadic Data." Political Analysis 23, no. 4, 564–577.
Stock, James and Mark Watson. (2008) “HeteroskedasticityRobust Standard Errors for Fixed Effect Panel Data Regression” Econometrica, 76, 1.
White, Halbert. (1980) “Using Least Squares to Approximate Unknown Regression Functions.” International Economic Review, 21, 1, 149–170.
pwfe
for fitting weighted fixed effects models with propensity
score weighting
### NOTE: this example illustrates the use of wfe function with randomly ### generated panel data with arbitrary number of units and time. ## generate panel data with number of units = N, number of time = Time N < 10 # number of distinct units Time < 15 # number of distinct time ## treatment effect beta < 1 ## generate treatment variable treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) ## make sure at least one observation is treated for each unit while ((sum(apply(treat, 2, mean) == 0) > 0)  (sum(apply(treat, 2, mean) == 1) > 0)  (sum(apply(treat, 1, mean) == 0) > 0)  (sum(apply(treat, 1, mean) == 1) > 0)) { treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) } treat.vec < c(treat) ## unit fixed effects alphai < rnorm(N, mean = apply(treat, 2, mean)) ## geneate two random covariates x1 < matrix(rnorm(N*Time, 0.5,1), ncol=N) x2 < matrix(rbeta(N*Time, 5,1), ncol=N) x1.vec < c(x1) x2.vec < c(x2) ## generate outcome variable y < matrix(NA, ncol = N, nrow = Time) for (i in 1:N) { y[, i] < alphai[i] + treat[, i] + x1[,i] + x2[,i] + rnorm(Time) } y.vec < c(y) ## generate unit and time index unit.index < rep(1:N, each = Time) time.index < rep(1:Time, N) Data.str < as.data.frame(cbind(y.vec, treat.vec, unit.index, x1.vec, x2.vec)) colnames(Data.str) < c("y", "tr", "strata.id", "x1", "x2") Data.obs < as.data.frame(cbind(y.vec, treat.vec, unit.index, time.index, x1.vec, x2.vec)) colnames(Data.obs) < c("y", "tr", "unit", "time", "x1", "x2") ############################################################ # Example 1: Stratified Randomized Experiments ############################################################ ## run the weighted fixed effect regression with strata fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. ### Average Treatment Effect mod.ate < wfe(y~ tr+x1+x2, data = Data.str, treat = "tr", unit.index = "strata.id", method = "unit", qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(mod.ate) ### Average Treatment Effect for the Treated mod.att < wfe(y~ tr+x1+x2, data = Data.str, treat = "tr", unit.index = "strata.id", method = "unit", qoi = "att", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(mod.att) ############################################################ # Example 2: Observational Studies with Unit Fixedeffects ############################################################ ## run the weighted fixed effect regression with unit fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. mod.obs < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05) ## summarize the results summary(mod.obs) ## extracting weigths summary(mod.obs)$W ## Not run: ################################################################### # Example 3: Observational Studies with differencesindifferences ################################################################### ## run differenceindifferences estimator. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. mod.did < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", estimator ="did", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05, verbose = TRUE) ## summarize the results summary(mod.did) ## extracting weigths summary(mod.did)$W ######################################################################### # Example 4: DID with Matching on Pretreatment Outcomes ######################################################################### ## implements matching on pretreatment outcomes where the maximum ## deviation is specified as 0.5 mod.Mdid < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", estimator ="Mdid", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05, maxdev.did = 0.5, verbose = TRUE) ## summarize the results summary(mod.Mdid) ## Note: setting the maximum deviation to infinity (or any value ## bigger than the maximum pairwise difference in the outcome) will ## return the same result as Example 3. dev < 1000+max(Data.obs$y)min(Data.obs$y) mod.did2 < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", estimator ="Mdid", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05, maxdev.did = dev, verbose = TRUE) ## summarize the results summary(mod.did2) mod.did2$coef[1] == mod.did$coef[1] ## End(Not run)
### NOTE: this example illustrates the use of wfe function with randomly ### generated panel data with arbitrary number of units and time. ## generate panel data with number of units = N, number of time = Time N < 10 # number of distinct units Time < 15 # number of distinct time ## treatment effect beta < 1 ## generate treatment variable treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) ## make sure at least one observation is treated for each unit while ((sum(apply(treat, 2, mean) == 0) > 0)  (sum(apply(treat, 2, mean) == 1) > 0)  (sum(apply(treat, 1, mean) == 0) > 0)  (sum(apply(treat, 1, mean) == 1) > 0)) { treat < matrix(rbinom(N*Time, size = 1, 0.25), ncol = N) } treat.vec < c(treat) ## unit fixed effects alphai < rnorm(N, mean = apply(treat, 2, mean)) ## geneate two random covariates x1 < matrix(rnorm(N*Time, 0.5,1), ncol=N) x2 < matrix(rbeta(N*Time, 5,1), ncol=N) x1.vec < c(x1) x2.vec < c(x2) ## generate outcome variable y < matrix(NA, ncol = N, nrow = Time) for (i in 1:N) { y[, i] < alphai[i] + treat[, i] + x1[,i] + x2[,i] + rnorm(Time) } y.vec < c(y) ## generate unit and time index unit.index < rep(1:N, each = Time) time.index < rep(1:Time, N) Data.str < as.data.frame(cbind(y.vec, treat.vec, unit.index, x1.vec, x2.vec)) colnames(Data.str) < c("y", "tr", "strata.id", "x1", "x2") Data.obs < as.data.frame(cbind(y.vec, treat.vec, unit.index, time.index, x1.vec, x2.vec)) colnames(Data.obs) < c("y", "tr", "unit", "time", "x1", "x2") ############################################################ # Example 1: Stratified Randomized Experiments ############################################################ ## run the weighted fixed effect regression with strata fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. ### Average Treatment Effect mod.ate < wfe(y~ tr+x1+x2, data = Data.str, treat = "tr", unit.index = "strata.id", method = "unit", qoi = "ate", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(mod.ate) ### Average Treatment Effect for the Treated mod.att < wfe(y~ tr+x1+x2, data = Data.str, treat = "tr", unit.index = "strata.id", method = "unit", qoi = "att", hetero.se=TRUE, auto.se=TRUE) ## summarize the results summary(mod.att) ############################################################ # Example 2: Observational Studies with Unit Fixedeffects ############################################################ ## run the weighted fixed effect regression with unit fixed effect. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. mod.obs < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05) ## summarize the results summary(mod.obs) ## extracting weigths summary(mod.obs)$W ## Not run: ################################################################### # Example 3: Observational Studies with differencesindifferences ################################################################### ## run differenceindifferences estimator. ## Note: the quantity of interest is Average Treatment Effect ("ate") ## and the standard errors allow heteroskedasticity and arbitrary ## autocorrelation. mod.did < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", estimator ="did", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05, verbose = TRUE) ## summarize the results summary(mod.did) ## extracting weigths summary(mod.did)$W ######################################################################### # Example 4: DID with Matching on Pretreatment Outcomes ######################################################################### ## implements matching on pretreatment outcomes where the maximum ## deviation is specified as 0.5 mod.Mdid < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", estimator ="Mdid", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05, maxdev.did = 0.5, verbose = TRUE) ## summarize the results summary(mod.Mdid) ## Note: setting the maximum deviation to infinity (or any value ## bigger than the maximum pairwise difference in the outcome) will ## return the same result as Example 3. dev < 1000+max(Data.obs$y)min(Data.obs$y) mod.did2 < wfe(y~ tr+x1+x2, data = Data.obs, treat = "tr", unit.index = "unit", time.index = "time", method = "unit", qoi = "ate", estimator ="Mdid", hetero.se=TRUE, auto.se=TRUE, White = TRUE, White.alpha = 0.05, maxdev.did = dev, verbose = TRUE) ## summarize the results summary(mod.did2) mod.did2$coef[1] == mod.did$coef[1] ## End(Not run)