effect and se from z

eff_from_z.R

get_beta_se_from_p_z_n_vary <- function(z, n, vy, maf)
{
    # qval <- qchisq(pnorm(abs(z), low=FALSE)*2, 1, low=F) / (qchisq(pnorm(abs(z), low=FALSE)*2, n-2, low=F)/(n-2))
    # r <- sqrt(sum(qval / (n - qval)))
    r <- sqrt(abs(z)^2 / (abs(z)^2 + n-2))
    b <- sign(z) * sqrt(r^2 * vy / (2*maf*(1-maf)))
    se <- b / z
    return(list(b=b, se=se))
}

nsim <- 500

emp_b <- array(0, nsim)
emp_se <- array(0, nsim)
est_b <- array(0, nsim)
est_se <- array(0, nsim)

for(i in 1:nsim)
{
	message(i)
	n <- round(runif(1, min=1000, max=10000))
	# n <- 1000
	p <- runif(1, min=0.1, max=0.5)
	# p <- 0.1
	a <- rbinom(n, 2, p)
	eff <- 1
	b <- eff * a + rnorm(n)*runif(1, min=1, max=10)
	b1 <- scale(b)
	mod <- coefficients(summary(lm(b ~ a)))
	mod1 <- coefficients(summary(lm(b1 ~ a)))
	emp_b[i] <- mod1[2,1]
	emp_se[i] <- mod1[2,2]
	jitter <- runif(1, min=0.5, max=1.5)

	z <- mod[2,3]

	temp <- get_beta_se_from_p_z_n_vary(z, n*1.5, 1, mean(a)/2 * jitter)
	est_b[i] <- temp$b
	est_se[i] <- temp$se
}


summary(lm(est_b ~ emp_b))
summary(lm(est_se ~ emp_se))

plot(est_b ~ emp_b)
plot(est_se ~ emp_se)


p <- runif(1000)
z <- qnorm(p)
t <- qt(p, 5000)

n <- round(runif(1, min=1000, max=10000))
p <- runif(1, min=0.1, max=0.5)
a <- rbinom(n, 2, p)
eff <- 1
b <- eff * a + rnorm(n)
mod <- coefficients(summary(lm(b ~ a)))
z <- mod[2,3]

f <- anova(lm(b ~ a))$F[1]
f_p <- qf(2*pt(z, n-1, lower=FALSE), 1, n-1, lower=FALSE)
f_p <- qchisq(pnorm(abs(z), low=FALSE)*2, 1, low=F) / (qchisq(pt(abs(z), n-1, low=FALSE)*2, n-2, low=F)/(n-2))

f
f_p

nsim <- 1000
f <- rep(0, nsim)
t <- rep(0, nsim)
z <- rep(0, nsim)

for(i in 1:nsim)
{
	message(i)
	p <- runif(1)
	n <- runif(1, min=1000, max=10000)
	f[i] <- qf(p, 1, n-1, low=FALSE)
	t[i] <- qt(p, n-1, low=FALSE)
	z[i] <- qnorm(p, low=FALSE)
}
dat <- data.frame(f,t,z,sqrt(f))

pairs(dat)