Simulation

A good way to test your intuitions about how the stats works is to run simulations.

# The basics

You'll need to generate random numbers: rnorm, rbin, rpois, will all come in handy.

# Some examples

## Generating the null distribution of Pearson's r for various sample sizes

Simulate a series of correlational studies each of which has n participants where there is no correlation between x and y, i.e. gives an estimate of the null sampling distribution of r.

``````studies = 1000
n = 50
r1 = c()

for (i in 1:studies) {
x = rnorm(n,0,1)
y = rnorm(n,0,1)
r1[i] = cor(x,y)
}

n = 200
r2 = c()

for (i in 1:studies) {
x = rnorm(n,0,1)
y = rnorm(n,0,1)
r2[i] = cor(x,y)
}

par(mfrow = c(1,2))
hist(r1, main="n = 50", col = "lightblue",xlim=c(-.7,.7), xlab = "r")
hist(r2, main="n = 200", col = "lightblue",xlim=c(-.7,.7), xlab = "r")
par(mfrow = c(1,1))```
```

Picture below: ## Generating data with different means and variances

``````# covariate values
x <- runif(100, -5, 5)

# Example 1
# model for variance
sds <- 0.5 + abs(x) * 2

# model for responses
y <- 2 + 3 * x + rnorm(100, sd = sds)

fm <- lm(y ~ x)
plot(fm, 1)
plot(fm, 3)

# Example 2
sds <- 1 + (abs(min(x)) + x) * 2

# model for responses
y <- 2 + 3 * x + rnorm(100, sd = sds)

fm <- lm(y ~ x)
plot(fm, 1)
plot(fm, 3)```
```

Thank you to Dimitris Rizopoulos

## Simulating mixture distributions

Here we simulate an experiment with two (between-subject) conditions. In both cases, there's a mixture distribution (intuitively there are two groups: very fast responses and slightly slower, perhaps caused by different systems), however in the second (represented by x2), the upper distribution is shifted upwards.

``````x1a  = rnorm(2000,100,5)
x1b  = rnorm(5000,110,5)
x1  = c(x1a,x1b)
x2a = rnorm(2000,100,5)
x2b = rnorm(5000,115,5)
x2  = c(x2a,x2b)

par(mfrow =c(1,2))
plot(density(x1),main="", xlab = "RT")
lines(density(x2), lty=2)
plot(density(x1b),main="", xlab = "RT")
lines(density(x2b), lty=2)
par(mfrow = c(1,1))

t.test(x1,x2)
t.test(x1b,x2b)```
```

The t-test on the mixture gives a significant result, but note the means are skewed because of the presence of the lower distro. (Mixture difference: 107 vs 111; pure effect difference: 110 vs. 115)

The density plot below shows both the mixture distributions and the pure effect bit of the distribution. ## Positive correlation isn't transitive

Let U, V, and W be independent random variables. Define X = U+V, Y = V+W, and Z = W-U. Then the correlation between X and Y is positive, Y and Z is positive, but the correlation between X and Z is negative. (See Langford, E.; Schwertman, N. & Owens, M. (2001) [Is the Property of Being Positively Correlated Transitive? The American Statistician, 55, 322-325.])

``````U = rnorm(100,0,1)
V = rnorm(100,0,1)
W = rnorm(100,0,1)

X = U+V
Y = V+W
Z = W-U

cor.test(X,Y)
cor.test(Y,Z)
cor.test(X,Z)

par(mfrow = c(2,2))

plot(Y~X,xlim=c(-5,5),ylim=c(-5,5))
abline(lm(Y~X),col="red")
plot(Y~Z,xlim=c(-5,5),ylim=c(-5,5))
abline(lm(Y~Z),col="red")
plot(Z~X,xlim=c(-5,5),ylim=c(-5,5))
abline(lm(Z~X),col="red")

par(mfrow = c(1,1))```
``` ## Median splits are a bad idea…

``````> x = rnorm(100,100,15)
> y = 2*x + 5 + rnorm(100,0,35)
>
> summary(lm(y ~ x))

Estimate Std. Error t value Pr(>|t|)
(Intercept)   22.442     21.503    1.04      0.3
x              1.885      0.212    8.88  3.3e-14 ***

Multiple R-squared: 0.446,      Adjusted R-squared: 0.44
F-statistic: 78.8 on 1 and 98 DF,  p-value: 3.32e-14

> summary(lm(y ~ cut(x,2)))

Estimate Std. Error t value Pr(>|t|)
(Intercept)           178.45       6.58   27.12  < 2e-16 ***
cut(x, 2)(94.2,135]    49.26       8.10    6.08  2.3e-08 ***

Multiple R-squared: 0.274,      Adjusted R-squared: 0.267
F-statistic:   37 on 1 and 98 DF,  p-value: 2.30e-08```
```

(Note the R-squares and the t's and the p's.)

## Sums of uniformly distributed random variables

See here for the background.

``````x1 = runif(1000,1,6)
x2 = runif(1000,1,6)
x3 = runif(1000,1,6)
x4 = runif(1000,1,6)
x5 = runif(1000,1,6)
x6 = runif(1000,1,6)
x7 = runif(1000,1,6)
x8 = runif(1000,1,6)
x9 = runif(1000,1,6)

par(mfrow=c(3,3))
hist(x1, main="",xlab="x")
hist(x1+x2, main="",xlab="x")
hist(x1+x2+x3, main="",xlab="x")
hist(x1+x2+x3+x4, main="",xlab="x")
hist(x1+x2+x3+x4+x5, main="",xlab="x")
hist(x1+x2+x3+x4+x5+x6, main="",xlab="x")
hist(x1+x2+x3+x4+x5+x6+x7, main="",xlab="x")
hist(x1+x2+x3+x4+x5+x6+x7+x8, main="",xlab="x")
hist(x1+x2+x3+x4+x5+x6+x7+x8+x9, main="",xlab="x")
par(mfrow=c(1,1))```
``` page revision: 26, last edited: 13 Jun 2008 10:52