(from Richard McElreath, Statistical Rethinking, 2nd edition)

Assumption

• 20 people born each year
• Uniform happyness at birth, never changes
• At 18 years old, eligible to marry. Probability of mariage in each year proportional to happiness.
• Married people remain married until death.
• At age 65, move to south coast of Spain

Purpose of simulation: show that "controlling" for marriage status creates a correlation between age and happiness, even though we know there is no correlation.

Pseudocode:

for t in 1:N: for p in 1:20: u ~ uniform(0, 1) add person (yob=t, id=p, happyness=u, marriageYear=-1, active=1) to data end for (y, i, h, m, a) in data: if a == 1 and t-y >= 65: a <- 0 end if a == 1 and t-y >= 18 and m == -1: x ~ uniform(0, 1) if x < h m = t end end end end

options(repr.plot.width=5, repr.plot.height=5)
N = 1000
d <- c()
for (t in 1:N) {
    for (p in 1:20) {
        u <- runif(1)
        d <- rbind(d, c(t, p, u, -1, 1))
    }

    for (k in 1:nrow(d)) {
        if (d[k, 5] == 1 & t-d[k, 1] >= 65) {
            d[k, 5] <- 0
        }
        if (d[k, 5] ==1 & t-d[k, 1] >= 18 & d[k, 4] == -1) {
            x <- runif(1)
            if (x < d[k, 3]) {
                d[k, 4] <- t
            }
        }
    }
}
d <- data.frame(d)
colnames(d) = c('yearOfBirth', 'id', 'happyness', 'yearOfMarriage', 'active')

d

A data.frame: 20000 × 5

yearOfBirth	id	happyness	yearOfMarriage	active
<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	1	0.62603998	21	0
1	2	0.17610604	29	0
1	3	0.63755080	20	0
1	4	0.19059768	19	0
1	5	0.62739305	19	0
1	6	0.88193941	19	0
1	7	0.78229004	19	0
1	8	0.94091459	19	0
1	9	0.39370558	19	0
1	10	0.56506059	20	0
1	11	0.42473782	19	0
1	12	0.91578508	19	0
1	13	0.02277007	27	0
1	14	0.24671927	22	0
1	15	0.26506390	22	0
1	16	0.78611475	19	0
1	17	0.54984230	22	0
1	18	0.25997598	19	0
1	19	0.28386308	19	0
1	20	0.50323644	21	0
2	1	0.11223783	27	0
2	2	0.53179581	20	0
2	3	0.83149464	20	0
2	4	0.30098566	20	0
2	5	0.04900386	-1	0
2	6	0.22119243	28	0
2	7	0.10302163	36	0
2	8	0.37679712	20	0
2	9	0.29198968	22	0
2	10	0.07820483	28	0
⋮	⋮	⋮	⋮	⋮
999	11	0.49704276	-1	1
999	12	0.02731464	-1	1
999	13	0.59909935	-1	1
999	14	0.32947182	-1	1
999	15	0.93678503	-1	1
999	16	0.58662415	-1	1
999	17	0.30133592	-1	1
999	18	0.90417594	-1	1
999	19	0.47295345	-1	1
999	20	0.36885327	-1	1
1000	1	0.10986075	-1	1
1000	2	0.72820126	-1	1
1000	3	0.30867509	-1	1
1000	4	0.16785673	-1	1
1000	5	0.46957143	-1	1
1000	6	0.62324436	-1	1
1000	7	0.11293581	-1	1
1000	8	0.10722321	-1	1
1000	9	0.33942030	-1	1
1000	10	0.23653725	-1	1
1000	11	0.16829029	-1	1
1000	12	0.57695470	-1	1
1000	13	0.80887565	-1	1
1000	14	0.80032018	-1	1
1000	15	0.40082283	-1	1
1000	16	0.37687157	-1	1
1000	17	0.30728460	-1	1
1000	18	0.17809839	-1	1
1000	19	0.97013598	-1	1
1000	20	0.37275293	-1	1

d$ageOfMarriage <- d$yearOfMarriage - d$yearOfBirth
d$ageOfMarriage[d$yearOfMarriage == -1] <- -1

hist(d[d$yearOfMarriage > 0,]$ageOfMarriage, freq=F, nclass=30, 
     xlab='age of marriage', main='Distribution of age of marriage')

d contains all relevant information. Let us draw some random data points from this log.

drawSample <- function() {
    yob <- sample(N, 1) ### uniformly distributed random integer from 1..N
    p <- sample(20, 1) ### uniformly distributed random integer from 1..20
    s <- d[(d$yearOfBirth == yob) & (d$id == p),]
    age <- sample(0:65, 1) ### uniformly distributed random integer from 0..65
    married <- 0
    if ((yob+age > as.numeric(s$yearOfMarriage)) & (as.numeric(s$yearOfMarriage) > -1)) {
        married <- 1
    }
    return(c(age, s$happyness, c('single', 'married')[1+married]))
}

df <- data.frame(t(replicate(10000, drawSample())))
colnames(df) <- c('age', 'happyness', 'married')
df$age <- as.numeric(as.character(df$age))
df$happyness <- as.numeric(as.character(df$happyness))

plot(happyness~married, df)

plot(married ~ happyness, df)

plot(age~married, df)

plot(married~age, df)

plot(happyness~age, df, cex=.3)

summary(lm(happyness ~ age, df))

Call:
lm(formula = happyness ~ age, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.50269 -0.24986  0.00556  0.25416  0.49694 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.035e-01  5.691e-03  88.462   <2e-16 ***
age         -9.763e-06  1.512e-04  -0.065    0.949    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2881 on 9998 degrees of freedom
Multiple R-squared:  4.172e-07,	Adjusted R-squared:  -9.96e-05 
F-statistic: 0.004171 on 1 and 9998 DF,  p-value: 0.9485

summary(lm(happyness ~ age + married, df))

Call:
lm(formula = happyness ~ age + married, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.5738 -0.2495 -0.0028  0.2475  0.6427 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.6532495  0.0100565   64.96   <2e-16 ***
age           -0.0029129  0.0002199  -13.25   <2e-16 ***
marriedsingle -0.1572672  0.0087688  -17.93   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2836 on 9997 degrees of freedom
Multiple R-squared:  0.03117,	Adjusted R-squared:  0.03098 
F-statistic: 160.8 on 2 and 9997 DF,  p-value: < 2.2e-16

library(ggplot2)

ggplot(df, aes(x=age, y=happyness, col=married)) + 
    geom_point(size=.1) +
    geom_smooth(method='lm') +
    geom_smooth(col='black', method='lm')