Gottesman_h2

require(MASS)

# fixed a2 and e2 for the entire script:
a <- .87 # additive genetic variance
e <- .13 # environmental variance

# make ZM covariance, i.e. the cov is a, the var = 1
sigma_mz <- matrix(c(1,a,a,1),2,2) 

# sample size (is big becuase rare traits)
n <- 100000

# make the outcome liabilities (which are correlated a between MZ)
liabilities <- mvrnorm(n,c(0,0),sigma_mz)

# make a rare (1% prevalence) binary trait/outcome 
binary <- liabilities > qnorm(0.99)

# concordance in 1st degree relatives:
tab <- table(binary[,1],binary[,2])/n # ± 50% concordance
tab[2,2]/(tab[2,2]+tab[2,1]) # ± 50% concordance

# repeat for 1st degree relatives

# correalted 0.5*a
sigma_sib_dz <- matrix(c(1,0.5*a,0.5*a,1),2,2)

# again continous trait risk
liabilities <- mvrnorm(n,c(0,0),sigma_sib_dz )

# binary trait (rare, still 1% obv)
binary <- liabilities > qnorm(0.99)


tab <- table(binary[,1],binary[,2])/n # ± 10% concordance
tab[2,2]/(tab[2,2]+tab[2,1]) # ± 10% concordance


# 2nd degree relatives:

# correlated 0.25*a
sigma_hz_cous <- matrix(c(1,0.25*a,0.25*a,1),2,2)

# again with the liability...
liabilities <- mvrnorm(n,c(0,0),sigma_hz_cous )

# binary trait (rare, still 1% obv)
binary <- liabilities > qnorm(0.99)


tab <- table(binary[,1],binary[,2])/n # ± 3/4% concordance
tab[2,2]/(tab[2,2]+tab[2,1]) # ± 3.4% concordance