reg_to_mean.R

set.seed(42)

# sample size
N = 10e3

# generate two groups
group = rbinom( N , 1 , 0.5 )

# generate group-specific environments
shared_env = rnorm(N,group*-15,1)

# generate random individual level effects (no group differences)
g_sib1 = rnorm(N , 100 , 15 )
g_sib2 = rnorm(N , 100 , 15 )

# add the individual effect and the shraed environmental effect
y_sib1 = g_sib1 + shared_env
# do the same to generate siblings
y_sib2 = g_sib2 + shared_env

# find people in the two groups that have similar outcomes
midpoint_group0 = (y_sib1 > 85 & y_sib1 < 95 & group == 0)
midpoint_group1 = (y_sib1 > 85 & y_sib1 < 95 & group == 1)

# now compare to their sibs
group0_change = c( mean( y_sib1[midpoint_group0] ) , mean( y_sib2[midpoint_group0] ) )
group0_se = c( sd( y_sib1[midpoint_group0] )/sqrt(sum(midpoint_group0)) , sd( y_sib2[midpoint_group0] )/sqrt(sum(midpoint_group0)) )

group1_change = c( mean( y_sib1[midpoint_group1] ) , mean( y_sib2[midpoint_group1] ) )
group1_se = c( sd( y_sib1[midpoint_group1] )/sqrt(sum(midpoint_group1)) , sd( y_sib2[midpoint_group1] )/sqrt(sum(midpoint_group1)) )

# --- now visualize
clr = c("#d73027","#1a9850")

plot( 0 , 0 , type="n" , ylim=c(80,110) , xlim=c(0.5,2.5) , las=1 , xaxt="n" , xlab="" , ylab="IQ" , main="regression to the mean with no genetic differences")

arrows(1:2,group0_change-1.96*group0_se,1:2,group0_change+1.96*group0_se,len=0, lwd=3 , col=paste(clr[1],"75",sep='') )
points( 1:2 , group0_change , pch=19 , col=clr[1] )
lines( 1:2 , group0_change , col=clr[1] , lwd=2 )

arrows(1:2,group1_change-1.96*group1_se,1:2,group1_change+1.96*group1_se,len=0, lwd=3 , col=paste(clr[2],"75",sep='') )
points( 1:2 , group1_change , pch=19  , col=clr[2] )
lines( 1:2 , group1_change , col=clr[2] , lwd=2 )

axis(side=1,at=1:2,labels=c("Reference Sibling","Comparison Sibling"))
legend("topleft",legend=c("Group 1","Group 2") , pch=19 , col=clr , bty="n" )