[This article was first published on Econometrics by Simulation, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I have received a number of requests for demonstration code on how to perform a power analysis using simulation in R.  I have already demonstrated howto do this in Stata but lacked the easy to use Stata command “simulate” that I preferred.  However, in a recent post I have written up a command very similar to simulate in R called SimpleSim.

This command takes is capable of taking a vector of parameters and feeding them into a function which returns a vector of results.  The results are turned into a table which can be used to show the rates of type 1 and type 2 error given a particular sample size and underlying effect.

Example:
Let’s imagine that we are interested in testing the effect that bednet usage has on the rate of infection of malaria.  We are concerned that the decision to purchase and use bednets might endogenous as a result of those more likely to come down with malaria being more prone to seeing bednets as a valuable investment.

So, in response we are going to design an randomized control trial (RCT) in which we either supply fully paid bednet coupons or half discounted coupons each to different thirds of the population leaving the remaining third as a control.  We are concerned that if supply the bednets at no cost then the recipients will not value them so that is why we give coupons that only reduce the cost of the bednets by 1/2.

After distributing the coupons, we would like to use the random distribution of the coupons as an instrument for use of a bednet.  Ultimately we would like to do two levels of analysis.

1. Using instrumental variables (IV) to see how effective bednets are at preventing Malaria in an active population. This is actually an interesting question because often there is imperfect compliance because few individuals spend all of the hours at dusk and dawn under a bednet.  In this case I should use an estimator that takes the form of the first stage being a binary regressor followed by a second stage binary regressor.  However, since linear models seem to be pretty good at approximating average partial effects I am just going to use a IV estimator.
2. Using ordinary least squares (OLS) we would like to see how effective each of the coupon programs are at increasing use of bednets within the homes.

require(“AER”) # We will use the ivreg function later

# We define a function that both simulates a population and estimates
# results and returns and returns the results.
MalariaIV <- function(
nsim = 100,
npop = 1000,    # Define the sampling population
treat0 = 1/3,   # Define the proportion which is control
treat1 = 1/3,   # Define the proportion which is treated with free bednets
treat2 = 1/3,   # Define the proportion which is treated with 50% cost
t0comp = .05,   # Bednet useage among control group
t1comp = .85,   # Bednet useage if recieving free net
t2comp = .5 ,   # Bednet useage if recieving 50% cost
malariaRT = .85,# Rate of getting malaria without bednet
netRT = .45,    # Rate of getting malaria with bednet
Tdetect = 1,    # Likihood of detecting malaria if present
Fdetect = 0,    # Likihood of detecting malaria if not present
alpha = .05     # The alpha level that a p-value must be below
) {
# Define how the function works
# First, define a vector to store results
pvalues <- NULL

#
for (i in 1:nsim) { # Repeat the simulation a number of times
# Generate the population by technology useage
simdata <- data.frame(
control=c(rep(1,npop*treat0),rep(0,npop*treat1), rep(0,npop*treat2)),
treat1 =c(rep(0,npop*treat0),rep(1,npop*treat1), rep(0,npop*treat2)),
treat2 =c(rep(0,npop*treat0),rep(0,npop*treat1), rep(1,npop*treat2)))
# Calculate the actual population generated (should be 999 in this case)
npeople <- nrow(simdata)
# Generate the rate of bednet usage.
simdata$bednet <- simdata$control*rbinom(npeople,1,t0comp) + # Bednet usage control
simdata$treat1 *rbinom(npeople,1,t1comp) + # Free simdata$treat2 *rbinom(npeople,1,t2comp)   # 50% cost
# Now let's generate the rate of malaria
simdata$malaria <- (simdata$bednet==0)*rbinom(npeople,1,malariaRT)+
(simdata$bednet==1)*rbinom(npeople,1,netRT) # Finally generate the rate of malaria detection as a function # of true parasite levels. simdata$mdetect <-
(simdata$malaria==0)*rbinom(npeople,1,Fdetect)+ # False detection (simdata$malaria==1)*rbinom(npeople,1,Tdetect)  # True detection

# Time to do our simple estimation of the effect of treatment on bednet use
lmcoef <- summary(lm(bednet~treat1+treat2,data=simdata))$coefficients # Now let's try the 2SLS to estimate the effect of bednets on contraction # of malaria. ivregest <- ivreg(mdetect~bednet | treat1+treat2, data=simdata) ivregcoef <- summary(ivregest)$coefficients

# Save the rejection of the null rates
pvalues <- rbind(pvalues,c(treat1=lmcoef[2,4]# Thus we want to ensure our study has at least 5000 participants.

# I think there might be some additional considerations for the
# number of participants required to ensure that there is no
# false rejection of the null.  However, I am no expert on Power
# Analysis so this is what I know.Formatted by Pretty R at inside-R.org

var vglnk = {key: '949efb41171ac6ec1bf7f206d57e90b8'};
(function(d, t) {
var s = d.createElement(t);
s.type = 'text/javascript';
s.async = true;
// s.defer = true;
//          s.src = '//cdn.viglink.com/api/vglnk.js';
var r = d.getElementsByTagName(t)[0];
r.parentNode.insertBefore(s, r);
}(document, 'script'));

Related
ShareTweet