R/get_sdy.R

In gap: Genetic Analysis Package

#' Get sd(y) from AF, n, b, se
#'
#' @param f Allele frequency.
#' @param n Sample size.
#' @param b effect size.
#' @param se standard error.
#' @param method method of averaging: "mean" or "median".
#' @param ... argument(s) passed to method
#'
#' @details
#' This function obtains standard error of a continuous outcome.
#'
#' @export
#' @return sd(y).
#' @examples
#' \dontrun{
#' set.seed(1)
#' X1 <- matrix(rbinom(1200,1,0.4),ncol=2)
#' X2 <- matrix(rbinom(1000,1,0.6),ncol=2)
#' colnames(X1) <- colnames(X2) <- c("f1","f2")
#' Y1 <- rnorm(600,apply(X1,1,sum),2)
#' Y2 <- rnorm(500,2*apply(X2,1,sum),5)
#' summary(lm1 <- lm(Y1~f1+f2,data=as.data.frame(X1)))
#' summary(lm2 <- lm(Y2~f1+f2,data=as.data.frame(X2)))
#' b1 <- coef(lm1)
#' b2 <- coef(lm2)
#' v1 <- vcov(lm1)
#' v2 <- vcov(lm2)
#' require(coloc)
#' ## Bayesian approach, esp. when only p values are available
#' abf <- coloc.abf(list(beta=b1, varbeta=diag(v1), N=nrow(X1), sdY=sd(Y1), type="quant"),
#'                  list(beta=b2, varbeta=diag(v2), N=nrow(X2), sdY=sd(Y2), type="quant"))
#' abf
#' # sdY
#' cat("sd(Y)=",sd(Y1),"==> Estimates:",sqrt(diag(var(X1)*b1[-1]^2+var(X1)*v1[-1,-1]*nrow(X1))),"\n")
#' for(k in 1:2)
#' {
#'   k1 <- k + 1
#'   cat("Based on b",k," sd(Y1) = ",sqrt(var(X1[,k])*(b1[k1]^2+nrow(X1)*v1[k1,k1])),"\n",sep="")
#' }
#' cat("sd(Y)=",sd(Y2),"==> Estimates:",sqrt(diag(var(X2)*b2[-1]^2+var(X2)*v2[-1,-1]*nrow(X2))),"\n")
#' for(k in 1:2)
#' {
#'   k1 <- k + 1
#'   cat("Based on b",k," sd(Y2) = ",sqrt(var(X2[,k])*(b2[k1]^2+nrow(X2)*v2[k1,k1])),"\n",sep="")
#' }
#' get_sdy(0.6396966,23991,0.04490488,0.009504684)
#' }

get_sdy <- function(f,n,b,se,method="mean",...)
{
  z <- sqrt(2*f*(1-f)*(b^2+n*se^2))
  ifelse(method=="mean",mean(z,...),median(z,...))
}