A specific data generating process

In this simulation, we will generate the observed data O in the following way.

  1. Draw \(W1\) from a Bernoulli(\(p = 0.20\))

  2. Given \(W1\), draw \(W2\) from a \(Bernoulli(p = expit(0.5 * W1))\)

  3. Given \(W1\), \(W2\) draw \(A\) from a \(Bernoulli(p = expit[W1 * W2])\)

  4. Given \(W1\), \(W2\), \(A\) draw \(C\) from a \(Bernoulli(p = expit[-A + 0.3 * W1 - W2])\)

  5. Given \(W1\), \(W2\), \(A\), \(C\), draw \(Y\) from a \(Normal(\mu = (4.4 * A + 0.7 * W1 - 2 *A * W2 - 2 * A * C), \sigma = 0.3^2)\) The random errors are independent.


1. Use Monte Carlo simulation to evaluate the true value of the following parameter of the observed data distribution :

\[ \phi = E_{0,W1,W2,C} [ \mathrm{E}(Y|A=1, W1, W2, C) - \mathrm{E}(Y|A=0, W1, W2, C)] \]

n <- 5000
R <- 10000 # 10K simulations
phi <- rep(NA, R)

for (r in 1:R) {
  W1 <- rbinom(n, 1, prob=0.20)
  W2 <- rbinom(n, 1, prob=plogis(0.5*W1))
  A <- rbinom(n, 1,  prob=plogis(W1*W2))
  C <- rbinom(n, 1, prob=plogis(-A+0.3*W1-W2))
  Y <- rnorm(n=n, mean = 4.4*A + 0.7*W1 - 2*A*W2 - 2*A*C, sd = 0.3)
  
  Y_1 <- 4.4 * 1 + 0.7 * W1 - 2 * 1 * W2 - 2 * 1 * C
  Y_0 <- 4.4 * 0 + 0.7 * W1 - 2 * 0 * W2 - 2 * 0 * C
  phi[r] <- mean(Y_1 - Y_0)
}

mean(phi)
## [1] 2.769436


Based on 10000 Monte Carlo simulations, the true value of the following parameter of the observed data distribution, was 2.769. \(\phi = E_{0,W1,W2,C} [ \mathrm{E}(Y|A=1, W1, W2, C) - \mathrm{E}(Y|A=0, W1, W2, C)] =\) 2.769.


2. Use Monte Carlo simulation to evaluate the true value of the average treatment effect:

\[ \theta^{*} = \mathrm{E}^*(Y_1)- \mathrm{E}^*(Y_0) \]

estimates <- rep(NA, R)

for (r in 1:R) {
  W1 <- rbinom(n, 1, prob=0.20)
  W2 <- rbinom(n, 1, prob=plogis(0.5*W1))
  A <- rbinom(n, 1,  prob=plogis(W1*W2))
  C <- rbinom(n, 1, prob=plogis(-A+0.3*W1-W2))
  Y <- rnorm(n=n, mean = 4.4*A + 0.7*W1 - 2*A*W2 - 2*A*C, sd = 0.3)
  
    
  Obs <- data.frame(W1, W2, A, C, Y)
    
  txt <- control <- Obs
  
  txt$A <- 1
  control$A <- 0
  
  reg1 <- glm(Y ~ A + W1 + W2 + C, data=Obs)
  
  # reg1
  pred1.txt <- predict(reg1, newdata = txt)
  pred1.ctrl <- predict(reg1, newdata = control)
    
  theta_star <- mean(pred1.txt-pred1.ctrl)
    
  estimates[r] <- theta_star
}

mean(estimates)
## [1] 2.794629


Based on 10000 Monte Carlo simulations, the true value of the average treatment effect, was 2.795. \(\theta^{*}(P^{*}) = 2.795\) is the difference in the counterfactual expected weight gain if all children were given RUTF and the counterfactual expected weight gain if all children were given the standard supplement.


End of Document