\begin{Schunk}
\begin{Sinput}
 ## Group mean differences
 
 ## There is evidence that birth order (social more than biological)
 ## affects IQ (doi:10.1126/science.1141493).  Suppose that we only
 ## record whether a child is raised as first-born (1 = first-born, 0 =
 ## not), and that the true mean difference between those two
 ## populations is 5 IQ points.
 103 - 98
\end{Sinput}
\begin{Soutput}
[1] 5
\end{Soutput}
\begin{Sinput}
 ## Specifically, let's suppose that first-borns represent 40% of
 ## children and have a mean IQ = 103, whereas children raised with
 ## older siblings are 60% and have a mean IQ = 98, so the overall
 ## average IQ is 100.
 .4*103 + .6*98
\end{Sinput}
\begin{Soutput}
[1] 100
\end{Soutput}
\begin{Sinput}
 ## We can sample separately from two groups with these
 ## characteristics.
 set.seed(123)
 dat <- data.frame(first = rbinom(n = 100, size = 1, prob = .4))
 ## We have just assigned rows of data into groups labeled 1 and 0
 head(dat)
\end{Sinput}
\begin{Soutput}
  first
1     0
2     1
3     0
4     1
5     1
6     0
\end{Soutput}
\begin{Sinput}
 ## Now we sample IQ scores, taking group membership into account.
 ## Here we use vectorized inputs to the data generator
 dat$IQ <- rnorm(NROW(dat), m = 98 + 5 * dat$first, sd = 15)
 ## round to nearest whole number, since that is how IQ scores are
 ## reported
 dat$IQ <- round(dat$IQ)
 head(dat, 20)
\end{Sinput}
\begin{Soutput}
   first  IQ
1      0 102
2      1 103
3      0  97
4      1 124
5      1 100
6      0 121
7      0  75
8      1 112
9      0 100
10     0 101
11     1 109
12     0  90
13     1  98
14     0  83
15     0  82
16     1 108
17     0 105
18     0  99
19     0 112
20     1 134
\end{Soutput}
\begin{Sinput}
 ## Do the sample statistics match the data generator (population) parameters?
 ## Several convenient ways to retrieve the answers
 m1 <- lm(IQ ~ first, data = dat)
 summary(m1)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = IQ ~ first, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-33.425  -8.969  -0.817   8.727  33.183 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   97.817      1.884  51.920   <2e-16 ***
first          3.608      2.979   1.211    0.229    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.59 on 98 degrees of freedom
Multiple R-squared:  0.01475,	Adjusted R-squared:  0.004698 
F-statistic: 1.467 on 1 and 98 DF,  p-value: 0.2287
\end{Soutput}
\begin{Sinput}
 ## Or
 
 t.test(IQ ~ first, data = dat)
\end{Sinput}
\begin{Soutput}
	Welch Two Sample t-test

data:  IQ by first
t = -1.2193, df = 85.601, p-value = 0.2261
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.491625  2.274959
sample estimates:
mean in group 0 mean in group 1 
       97.81667       101.42500 
\end{Soutput}
\begin{Sinput}
 ## or
 aggregate(IQ ~ first, data = dat,
            FUN = function(x) c(M = mean(x), SD = sd(x)))
\end{Sinput}
\begin{Soutput}
  first      IQ.M     IQ.SD
1     0  97.81667  14.78030
2     1 101.42500  14.30597
\end{Soutput}
\begin{Sinput}
 ## We expect sampling variability, the observed group
 ## difference will not be exactly the same as the true
 ## difference in the data generator
 diff(aggregate(IQ ~ first, data = dat, mean)$IQ)
\end{Sinput}
\begin{Soutput}
[1] 3.608333
\end{Soutput}
\begin{Sinput}
 ## Let's use Monte Carlo methods to find out how much the
 ## mean-difference varies
 
 ## define a function for one replication
 getData <- function(rep, N, M1, M2, SD1, SD2) {
      dat <- data.frame(first = rbinom(n = N, size = 1, prob = .4))
      dat$IQ <- rnorm(N, m = M2 + dat$first * (M1 - M2), 
                         s = SD2 + dat$first * (SD1 - SD2))
      dat$IQ <- round(dat$IQ)
      dat
  }
 getDiff <- function(dframe){
      diff(aggregate(IQ ~ first, data = dframe, mean)$IQ)
  }
 ## try it on one replication first
 dat1 <- getData(1, N = 100, M1 = 103, M2 = 98, SD1 = 15, SD2 = 15)
 getDiff(dat1)
\end{Sinput}
\begin{Soutput}
[1] 7.899117
\end{Soutput}
\begin{Sinput}
 ##
 ## Combine into 1 step if we don't want to save the data
 getDiff(getData(1, N = 100, M1 = 103, M2 = 98, SD1 = 15, SD2 = 15))
\end{Sinput}
\begin{Soutput}
[1] 5.374122
\end{Soutput}
\begin{Sinput}
 ## now do it 2000 times
 
 ## Make a wrapper function
 oneSim <- function(rep, N = 100, M1 = 103, M2 = 98, SD1 = 15, SD2 = 15){
     getDiff(getData(1, N = N, M1 = M1, M2 = M2, SD1 = SD1, SD2 = SD2))
  }
 set.seed(123)
 myMeanDiffs <- vapply(1:2000, oneSim, N = 100,
                        M1 = 103, M2 = 98, SD1 = 15, SD2 = 15,
                        numeric(1))
 ## check results
 mean(myMeanDiffs)
\end{Sinput}
\begin{Soutput}
[1] 4.920373
\end{Soutput}
\begin{Sinput}
 sd(myMeanDiffs)
\end{Sinput}
\begin{Soutput}
[1] 3.073736
\end{Soutput}
\begin{Sinput}
 
 
 ## Try it yourself!
 ## Generate a dichotomous grouping variable using rbinom
 ## Use a linear equation to generate dependent variable scores
 ## Don't forget the error term!
\end{Sinput}
\end{Schunk}
