\begin{Schunk}
\begin{Sinput}
 R <- 0.30 #"medium" in Cohen's opinion
 ## Start with a correlation matrix of predictors
 ## If you do this the standard R way
 IQ.cor <- matrix(c(1, R, R, 1), nrow = 2, ncol = 2,
                   dimnames = list(c("dadIQ", "momIQ"), 
                                c("dadIQ", "momIQ")))
 IQ.cor 
\end{Sinput}
\begin{Soutput}
      dadIQ momIQ
dadIQ   1.0   0.3
momIQ   0.3   1.0
\end{Soutput}
\begin{Sinput}
 ## My equivalent method in rockchalk
 IQ.cor <- lazyCor(X = R, d = 2)
 IQ.cor
\end{Sinput}
\begin{Soutput}
     [,1] [,2]
[1,]  1.0  0.3
[2,]  0.3  1.0
\end{Soutput}
\begin{Sinput}
 ## these are the assumed means of momIQ and dadIQ
 IQ.M <- c(dadIQ = 100, momIQ = 100)
 ## mvrnorm will take these as column names in the 
 ## output. That's why those are named
 
 ## these are the standard devations of momIQ and dadIQ
 IQ.SD <- c(15, 15)
 ## The diagonal matrix we need will be
 diag(IQ.SD) 
\end{Sinput}
\begin{Soutput}
     [,1] [,2]
[1,]   15    0
[2,]    0   15
\end{Soutput}
\begin{Sinput}
 ## Create the covariance matrix
 IQ.cov <- diag(IQ.SD) %*% IQ.cor %*% diag(IQ.SD)
 IQ.cov
\end{Sinput}
\begin{Soutput}
      [,1]  [,2]
[1,] 225.0  67.5
[2,]  67.5 225.0
\end{Soutput}
\begin{Sinput}
 ## Use R's builtin cov2cor to double-check thc
 ## correlations
 cov2cor(IQ.cov)
\end{Sinput}
\begin{Soutput}
     [,1] [,2]
[1,]  1.0  0.3
[2,]  0.3  1.0
\end{Soutput}
\begin{Sinput}
 N <- 100
 set.seed(123)
 dat <- mvrnorm(n = N, mu = IQ.M, Sigma = IQ.cov)
 dat <- as.data.frame(round(dat))
 head(dat)
\end{Sinput}
\begin{Soutput}
  dadIQ momIQ
1    95    91
2   118   119
3    86   117
4   117    94
5    96    88
6   112   118
\end{Soutput}
\begin{Sinput}
 ## Do observed means and correlations reflect 
 ## the population parameters? Rounding is not too harmful
 colMeans(dat)
\end{Sinput}
\begin{Soutput}
 dadIQ  momIQ 
100.23  99.86 
\end{Soutput}
\begin{Sinput}
 cor(dat)
\end{Sinput}
\begin{Soutput}
          dadIQ     momIQ
dadIQ 1.0000000 0.2170863
momIQ 0.2170863 1.0000000
\end{Soutput}
\begin{Sinput}
 ## Now we will once again draw random birth-orders
 dat$first <- rbinom(n = N, size = 1, prob = .4)
 ## Now that we have our multiple predictors, we 
 ## specify a model to generate outcomes (child's IQ) 
 ## with random sampling error.
 
 stde <- 9
 b <- c(-3, 5, .5, .5)
 ## parameters designed so child IQ is average of parents
 dat$IQnoe <- b[1] + b[2]*dat$first + b[3]*dat$dadIQ + b[4]*dat$momIQ 
 dat$IQ <- dat$IQnoe + rnorm(N, m = 0, sd = stde)
 dat$IQ <- round(dat$IQ)
 head(dat)
\end{Sinput}
\begin{Soutput}
  dadIQ momIQ first IQnoe  IQ
1    95    91     1  95.0  92
2   118   119     0 115.5 110
3    86   117     1 103.5 100
4   117    94     0 102.5 103
5    96    88     0  89.0 103
6   112   118     0 112.0 111
\end{Soutput}
\begin{Sinput}
 ## Do sample statistics match data generator parameters?
 mod0 <- lm(IQ ~ first + momIQ + dadIQ, data = dat)
 summary(mod0)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = IQ ~ first + momIQ + dadIQ, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.0270  -5.5010   0.1397   6.3073  19.2770 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.73128    8.59794   0.318   0.7514    
first        4.26077    1.86205   2.288   0.0243 *  
momIQ        0.46511    0.07340   6.336 7.60e-09 ***
dadIQ        0.49168    0.06181   7.955 3.53e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.074 on 96 degrees of freedom
Multiple R-squared:  0.5867,	Adjusted R-squared:  0.5738 
F-statistic: 45.42 on 3 and 96 DF,  p-value: < 2.2e-16
\end{Soutput}
\end{Schunk}
