anomalies

anomalies.R

n <- 100000
nsnp1 <- 10
nsnp2 <- 10
G1 <- matrix(rbinom(n*nsnp1, 2, 0.5), n, nsnp1)
G2 <- matrix(rbinom(n*nsnp2, 2, 0.5), n, nsnp2)

G2[,1] <- G1[,1]

eff1 <- rnorm(nsnp1)
eff2 <- rnorm(nsnp2)

prs1 <- scale(G1 %*% eff1)
prs2 <- scale(G2 %*% eff2)

x1 <- prs1 + rnorm(n)
x2 <- prs2 + rnorm(n)

bx1y <- 5
bx2y <- -5

y <- x1 * bx1y + x2 * bx2y + rnorm(n)



effx1h <- sex1 <- array(0, nsnp1)
effy1h <- sey1 <- array(0, nsnp1)
effx2h <- sex2 <- array(0, nsnp2)
effy2h <- sey2 <- array(0, nsnp2)

for(i in 1:nsnp1)
{
	m <- summary(lm(x1 ~ G1[,i]))$coefficients
	effx1h[i] <- m[2,1]
	sex1[i] <- m[2,2]
	m <- summary(lm(y ~ G1[,i]))$coefficients
	effy1h[i] <- m[2,1]
	sey1[i] <- m[2,2]
}

for(i in 1:nsnp2)
{
	m <- summary(lm(x2 ~ G2[,i]))$coefficients
	effx2h[i] <- m[2,1]
	sex2[i] <- m[2,2]
	m <- summary(lm(y ~ G2[,i]))$coefficients
	effy2h[i] <- m[2,1]
	sey2[i] <- m[2,2]
}



flipplot <- function(bexp, bout)
{
	index <- bout < 0
	bout[index] <- abs(bout[index])
	bexp[index] <- bexp[index] * -1
	plot(bout ~ bexp)
}

flipplot(effx1h, effy1h)

plot(effy1h ~ effx1h)
plot(effy2h ~ effx2h)

effx1h
effy1h
effx2h
effy2h