library(ggplot2)
mean_response <- expand.grid(treatment = c('placebo', 'active'),
                         gender = c('male', 'female'))
mean_response$mean_response <- ifelse(mean_response$treatment == 'active', 0.75, 0.25)
mean_response$mean_response <- mean_response$mean_response + ifelse(mean_response$gender == 'female', 0.25, 0)
ggplot(mean_response, aes(treatment, mean_response, fill = gender)) + 
  geom_bar(position = 'dodge', stat = 'identity') +
  ggtitle("Assumed mean response pattern")

set.seed(20130514)
nobs <- 120
sim_data <- cbind(mean_response, data.frame(id = 1:nobs))
sim_data$noise <- rnorm(n=nobs)
sim_data$Y <- with(sim_data, mean_response + noise)
t.test(Y~treatment, data = sim_data)

	Welch Two Sample t-test

data:  Y by treatment
t = -2, df = 116, p-value = 0.06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.70  0.02
sample estimates:
mean in group placebo  mean in group active 
                 0.48                  0.83 

summary(lm(Y~treatment + gender, data = sim_data))

Call:
lm(formula = Y ~ treatment + gender, data = sim_data)

Residuals:
   Min     1Q Median     3Q    Max 
-2.441 -0.678  0.159  0.611  2.858 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)        0.326      0.156    2.09    0.039 *
treatmentactive    0.340      0.181    1.88    0.062 .
genderfemale       0.317      0.181    1.76    0.082 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.99 on 117 degrees of freedom
Multiple R-squared:  0.0537,	Adjusted R-squared:  0.0375 
F-statistic: 3.32 on 2 and 117 DF,  p-value: 0.0396

ggplot(sim_data, aes(treatment, Y, fill = gender)) + 
  geom_boxplot(position = 'dodge') +
  ggtitle("Example simulated dataset")

set.seed(20130514)
nsim <- 10000
t_test_result <- lm_result <- matrix(NA, nsim, 2)
for(i in 1:nsim){
  sim_data$noise <- rnorm(n=nobs)
  sim_data$Y <- with(sim_data, mean_response + noise)
  tt1 <- t.test(Y~treatment, data = sim_data)
  lm1 <- summary(lm(Y~treatment + gender, data = sim_data))
  t_test_result[i, ] <- c(tt1$estimate %*% c(-1,1), tt1$p.value)
  lm_result[i, ] <- lm1$coef['treatmentactive', c('Estimate', 'Pr(>|t|)')]
}