sashagusev

## An implementation of Gaussian Process regression in R with examples of fitting and plotting with multiple kernels.
## Input: Does not require any input
## Output: Generates multiple SVG plots

## This is work based on Chapter 2 of Rasmussen and Williams, available at [ http://www.gaussianprocess.org/ ].
## It was extremely helpful to follow along with the implementation at [ http://www.jameskeirstead.ca/blog/gaussian-process-regression-with-r/ with some change in notation ]

# this library is necessary for sampling from the Multivariate Normal
require(MASS)

sashagusev

require(MASS)

## This is work based on Chapter 2 of Rasmussen and Williams, availabler at [ http://www.gaussianprocess.org/ ].

### GP solution and plotting:

gp_solve_basis = function( x.train , y.train , x.pred , h.train, h.pred , kernel , sigma2e = 0 ) {
	# Args:
	#	x.train:	matrix of training data
	#	y.train:	matrix of training labels

sashagusev

# Num SNPS
M = 50
# Num POP
N_POP = 10000
# Num GWAS
N_GWAS = 5000
# Num Trios
N_trio = 2000

# GWAS heritability

sashagusev

library("NMF")
library("ggplot2")
library("reshape2")

N_cells = 500
N_genes = 100

N_cells_program = 20
N_genes_program = 20

sashagusev

library("survival")
library("cmprsk")
library("reshape2")
library("ggplot2")
library("RColorBrewer")

# Standard function for simulating survival with a Weibull baseline hazard
# https://stats.stackexchange.com/questions/135124/how-to-create-a-toy-survival-time-to-event-data-with-right-censoring

# baseline hazard: Weibull

sashagusev

library("survival")
library("cmprsk")
library("RColorBrewer")

set.seed(66)

# Standard function for simulating survival with a Weibull baseline hazard
# https://stats.stackexchange.com/questions/135124/how-to-create-a-toy-survival-time-to-event-data-with-right-censoring

# baseline hazard: Weibull

sashagusev

library("RColorBrewer")
clr = brewer.pal(3,"Set1")

# unrelated samples
N = 50000
# heritability parameter
h2g = 0.4
# shared environment parameter
h2c = 0.4
# environment parameter

sashagusev

library("RColorBrewer")

# function for generating twis
make_all_twins = function(N,rge,rgc,h2g,h2c,h2e,h2gxe=0,h2gxc=0,pop_strat=0,env_strat=0){
  # --- generate MZs
  # generate MZ twins
  gv_MZ1 = scale(rnorm(N))
  gv_MZ2 = scale(gv_MZ1)
  # generate rGE
  ev_MZ1 = scale(rge * gv_MZ1 + sqrt(1 - rge^2) * scale(rnorm(N)))

sashagusev

N = 50000

get_correlated_trait = function(input,r) {
  r * scale(input) + sqrt(1-r^2) * rnorm(length(input),0,1)
}
get_noisy_trait = function(input,noise) {
  sqrt(1-noise) * scale(input) + sqrt(noise) * rnorm(length(input),0,1)
}

# good parameters from a 2-dof grid search

sashagusev

get_correlated_trait = function(input,r) {
  r * (input) + sqrt(1-r^2) * rnorm(length(input),0,1)
}
get_noisy_trait = function(input,noise) {
  sqrt(1-noise) * (input) + sqrt(noise) * rnorm(length(input),0,1)
}

simulate_non_genetic = function(N,AM,ENV,err) {
  set.seed(42)