simulation tools

simtools.r

fastAssoc <- function(y, x)
{
	index <- is.finite(y) & is.finite(x)
	n <- sum(index)
	y <- y[index]
	x <- x[index]
	vx <- var(x)
	vy <- var(y)
	bhat <- cov(y, x) / vx
	ahat <- mean(y) - bhat * mean(x)
	# fitted <- ahat + x * bhat
	# residuals <- y - fitted
	# SSR <- sum((residuals - mean(residuals))^2)
	# SSF <- sum((fitted - mean(fitted))^2)

	rsq <- (bhat * vx)^2 / (vx * vy)
	fval <- rsq * (n-2) / (1-rsq)
	tval <- sqrt(fval)
	se <- abs(bhat / tval)

	# Fval <- (SSF) / (SSR/(n-2))
	# pval <- pf(Fval, 1, n-2, lowe=F)
	p <- pf(fval, 1, n-2, lowe=F)
	return(list(
		ahat=ahat, bhat=bhat, se=se, fval=fval, pval=p, n=n, rsq=bhat^2 * vx / vy
	))
}

getFittedVals <- function(y, x)
{
	n <- length(x)
	bhat <- cov(y, x) / var(x)
	ahat <- mean(y) - bhat * mean(x)
	fitted <- ahat + x * bhat
	return(fitted)
}

bhat = cov(x,y) / var(x)
rsq = cov(x,y)^2 / (var(x) * var(y))


runGwas <- function(x, y)
{
	require(dplyr)
	y <- as.matrix(y)
	x <- as.matrix(x)
	ntrait <- ncol(y)
	nsnp <- ncol(x)
	l <- list()
	for(i in 1:ntrait)
	{
		message(i)
		res <- array(0, c(nsnp, 7))
		res[,1] <- 1:nsnp
		for(j in 1:nsnp)
		{
			res[j,2:7] <- unlist(fastAssoc(y[,i], x[,j]))
		}
		l[[i]] <- data.frame(res)
		names(l[[i]]) <- c("snp", "a", "b", "se", "fval", "pval", "n")
		l[[i]]$trait <- paste("trait", i)

	}
	return(bind_rows(l))
}


makeGeno <- function(nid, nsnp, maf=0.5)
{
	matrix(rbinom(nid * nsnp, 2, maf), nid, nsnp)
}


makePhen <- function(effs, indep, vy=1, vx=rep(1, length(effs)))
{
	if(is.null(dim(indep))) indep <- cbind(indep)
	stopifnot(ncol(indep) == length(effs))
	stopifnot(length(vx) == length(effs))
	cors <- effs * vx / sqrt(vx) / sqrt(vy)
	stopifnot(sum(cors^2) <= 1)
	cors <- c(cors, sqrt(1-sum(cors^2)))
	indep <- t(t(scale(cbind(indep, rnorm(nrow(indep))))) * cors * c(vx, 1))
	y <- drop(scale(rowSums(indep)) * sqrt(vy))
	return(y)
}

chooseEffects <- function(nsnp, totvar, sqrt=TRUE)
{
	eff <- rnorm(nsnp)
	aeff <- abs(eff)
	sc <- sum(aeff) / totvar
	out <- eff / sc
	if(sqrt)
	{
		out <- sqrt(abs(out)) * sign(out)
	}
	return(out)
}

spreadGwas <- function(gwas, col)
{
	require(tidyr)
	l <- list()
	stopifnot(all(col %in% names(gwas)))
	for(i in 1:length(col))
	{
		x <- select(gwas, snp, trait, get(col[i]))
		l[[i]] <- spread_(x, "trait", col[i])
	}

	if(length(col)==1) 
	{
		return(l[[i]])
	} else {
		names(l) <- col
		return(l)
	}
}





# g <- makeGeno(1000, 20)
# u <- makePhen(chooseEffects(20, 0.3), g)

# x <- makePhen(sqrt(0.6), u)
# y <- makePhen(sqrt(0.2), u)

# gwas <- runGwas(g, cbind(x,y))
# out <- spreadGwas(gwas, c("b"))

# plot(out[,3] ~ out[,2])




# g <- makeGeno(10000, 20)
# x <- makePhen(chooseEffects(20, 0.3), g)
# u <- rnorm(10000)
# y <- makePhen(sqrt(c(0.4, 0.02)), cbind(x,u))

# gwas <- runGwas(g, cbind(x,y))
# out <- spreadGwas(gwas, c("b"))

# plot(out[,3] ~ out[,2])

# gwas