sashagusev

set.seed(42)

### --- REML
# Original code:
# https://github.com/sashagusev/SKK-REML-sim/blob/master/func_reml.R

library('msm')

# Utility function for calculating log of determinant
`logdet` <-

sashagusev

# Simulations and visualization of a heritable liability threshold model across generations and in families
# ---
# Liability functions (from Baselmans et al.)
# See: https://github.com/BartBaselmans/CHARRGe

risk_from_liability = function( prev , h2 ) {
  # Threshold
  Tx = -qnorm(prev,0,1)
  z = dnorm(Tx)
  i = z/prev

sashagusev

set.seed(1)

# Library for GLMM computations
library(PQLseq2)
# Library for mixed model / Wald test
library("GMMAT")

# convenience function to suppress prints
quiet <- function(x) { 
  sink(tempfile()) 

sashagusev

set.seed(42)

# Library for GLMM computations
library(PQLseq2)
# Library for mixed model / Wald test
library("GMMAT")

# convenience function to suppress prints
quiet <- function(x) { 
  sink(tempfile()) 

sashagusev

set.seed(42)

# Library for GLMM computations
library(PQLseq2)

quiet <- function(x) { 
  sink(tempfile()) 
  on.exit(sink()) 
  invisible(force(x)) 
} 

sashagusev

# ---
# Code to estimate and simulate selection on healthy individuals (controls) for a threshold trait
# Updated with alternative simulation setup (producing the same result) and averaging to get smoother curves
# ---

library("VGAM")

# --- Estimate using the Breeder's Equation for threshold traits (Walsh & Lynch 2018; Dempster & Lerner 1950)
# We will select on being "healthy" where 1-prev of the population is healthy
ThresholdSelection <- function( FractionHealthy, Heritability ) {

sashagusev

library("RColorBrewer")
library("umap")
library("Rtsne")

# Number of markers
M = 2e3
# Number of people
N = 1e3

# frequencies in the ancestral population

sashagusev

par(mfrow=c(1,5))
par(oma=c(0.5,0.5,0.5,0.5))
par(mar=c(0.5,0.5,1,0.5))

# fst parameters
# frequency in one population
p_b = 0.2
# scan of frequencies in the second population
p_a = seq(0.01,0.99,by=0.01)
# estimate using Hudson's Fst

sashagusev

par(mfrow=c(1,2))
seeds = 10
x = seq(1,10,by=0.001)
all_u_coefs = seq(-5,5,by=0.1)

# lin reg
all_bs = rep(NA,length(all_u_coefs))
for ( i in 1:length(all_u_coefs) ) {
  u_coef = all_u_coefs[i]
  cur_r = rep(NA,seeds)

sashagusev

# simple factor analysis example
# ---

# lavaan for factor analysis
library('lavaan')

# --- parameters
# number of people
N = 1e3
# number of factors