Homework 2

2.1 Simulation Run

Parameter Definitions

## set seed to make simulations reproducible
set.seed(20210108)

## let's start with some parameter definitions
nsamp = 100
beta = 1 #effect size 
h2 = 0.1 #heritability
sig2X = h2
sig2epsi = (1 - sig2X) * beta^2
sigX = sqrt(sig2X) #genetic standard deviation
sigepsi = sqrt(sig2epsi)

#simulate X
X = rnorm(nsamp,mean=0, sd= sigX) #SNP
epsi = rnorm(nsamp,mean=0, sd=sigepsi) #environmental effects

#simulate Ynull
Ynull = rnorm(nsamp, mean=0, sd=beta) #phenotype

#simulate Yalt
Yalt = X * beta + epsi #phenotype

First, defining fastlm to repeat linear regression efficiently

fastlm = function(xx,yy)
{
  ## compute betahat (regression coef) and pvalue with Ftest
  ## for now it does not take covariates
  
  df1 = 2
  df0 = 1
  ind = !is.na(xx) & !is.na(yy)
  xx = xx[ind]
  yy = yy[ind]
  n = sum(ind)
  xbar = mean(xx)
  ybar = mean(yy)
  xx = xx - xbar
  yy = yy - ybar
  
  SXX = sum( xx^2 )
  SYY = sum( yy^2 )
  SXY = sum( xx * yy )
  
  betahat = SXY / SXX
  
  RSS1 = sum( ( yy - xx * betahat )^2 )
  RSS0 = SYY
  
  fstat = ( ( RSS0 - RSS1 ) / ( df1 - df0 ) )  / ( RSS1 / ( n - df1 ) )
  pval = 1 - pf(fstat, df1 = ( df1 - df0 ), df2 = ( n - df1 ))
  res = list(betahat = betahat, pval = pval)
  
  return(res)
}

nsim = 10000
## simulate normally distributed X and epsi
Xmat = matrix(rnorm(nsim * nsamp,mean=0, sd= sigX), nsamp, nsim)
epsimat = matrix(rnorm(nsim * nsamp,mean=0, sd=sigepsi), nsamp, nsim)

## generate Yalt (X has an effect on Yalt)
Ymat_alt = Xmat * beta + epsimat

## generate Ynull (X has no effect on Ynull)
Ymat_null = matrix(rnorm(nsim * nsamp, mean=0, sd=beta), nsamp, nsim)

## let's look at the dimensions of the simulated matrices
print(dim(Ymat_null))

[1]   100 10000

print(dim(Ymat_alt))

[1]   100 10000

The matrix dimensions for Ynull and Yalt is 100 x 10000

#Running simulation 10,000 times for Ynull
pvec_null = rep(NA,nsim)
bvec_null = rep(NA,nsim)

for(ss in 1:nsim)
{
  fit = fastlm(Xmat[,ss], Ymat_null[,ss]) #linear regression
  pvec_null[ss] = fit$pval  #find p-values and save them
  bvec_null[ss] = fit$betahat #find coefficients and save them
}

summary(pvec_null)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.244   0.491   0.494   0.740   1.000 

#This is for Yalt
pvec_alt = rep(NA,nsim)
bvec_alt = rep(NA,nsim)

for(ss in 1:nsim)
{
  fit = fastlm(Xmat[,ss], Ymat_alt[,ss]) #linear regression
  pvec_alt[ss] = fit$pval    #find p-values and save them
  bvec_alt[ss] = fit$betahat #find coefficients and save them
}

summary(pvec_alt)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0001  0.0013  0.0205  0.0098  0.9727 

2.2 P-value Distributions

Distribution of p-values under the Null

hist(pvec_null,xlab="p-value",main="Histogram of p-values under Null")

The histogram of the p-values is uniformly distributed. This is because if the null hypothesis is true, then all possible p-values are equally likely across the interval (0,1).

Distribution of p-values under the Alternative

hist(pvec_alt,xlab="p-value",main="Histogram of p-values under Alt")

Under the alternative hypothesis, we simulated the data so that there WOULD be association between X and Y. As such, we expect to see strong evidence against the null. This means many p-values will be small and close to 0, leading to a left-ward skew.

2.3 Mixed P-values

is_null = sample(c(rep(TRUE, nsim * 0.9), rep(FALSE, nsim * 0.1)))

p_mixed = ifelse(is_null, pvec_null, pvec_alt)

#histogram
hist(p_mixed,xlab="p-value",main="Histogram of p-values under mixture of null and alt")

In the histogram above, we see a mixture between the histogram of p-values under the null distribution and the histogram under the alternative distribution. We would expect to see a non-uniform distribution because some of the p-values were observed under the alternative hypothesis, where there is an association between X and Y. There is a greater frequency of significant p-values (less than 0.05) compared to all other p-values, as shown by the spike on the left.

p_unif <- runif(nsim) #simulating random uniform variables 
qqplot(-log10(p_unif), -log10(p_mixed), 
       main = "QQ Plot of P-values",
       xlab = "Expected -log10(P)",
       ylab = "Observed -log10(P)",
       pch = 1)

abline(0, 1, col = "red", lwd = 2, lty = 2) #identity line (y=x)

Deviations from the identity line shows a departure from the expected distribution. In other words, if the p-values are from the null hypothesis (uniform), they should lie along the identity line. The plot indicates a mix of null and non-null effects.The initial proximity to the identity line demonstrates some p-values align with the null distribution. Then, there is an upwards stretch that moves away from the red line, indicating that other p-values are significant under the alternative distribution.

2.4 Confusion Matrix & Reality Check

alpha = 0.05

#TN: null is true and test does not reject
TN = sum(is_null & (p_mixed >= alpha))

#FP: null is true but test rejected
FP = sum(is_null & (p_mixed < alpha)) 

#FN: null is false but test fails to reject
FN = sum(!is_null & (p_mixed >= alpha))

#TP: null is false and test fails to reject
TP = sum(!is_null & (p_mixed < alpha))

cat("true negatives:", TN, "\n")

true negatives: 8521 

cat("false positives:", FP, "\n")

false positives: 479 

cat("false negatives:", FN, "\n")

false negatives: 98 

cat("true positives:", TP, "\n")

true positives: 902 

This table cannot be generated in real data analysis. In a simulation, we create the data, therefore we know the “ground truth” (i.e. which events are TRULY significant/not significant). However, in real data analysis, we can only observe which tests are called significant or not based on p-values. We cannot directly determine which hypotheses are truly null or alternative because the true state (ground truth) of each hypothesis is unknown.

2.5 Bonferroni Correction

bonferroni <- 0.05 / nsim 
cat("Bonferroni-adjusted threshold:", bonferroni, "\n")

Bonferroni-adjusted threshold: 0.000005 

sig_bonferroni <- sum(p_mixed < bonferroni)
cat("Number of significant tests:", sig_bonferroni, "\n")

Number of significant tests: 83 

TN_bonferroni <- sum(is_null & (p_mixed >= bonferroni)) 
FP_bonferroni <- sum(is_null & (p_mixed < bonferroni))  
FN_bonferroni <- sum(!is_null & (p_mixed >= bonferroni))
TP_bonferroni <- sum(!is_null & (p_mixed < bonferroni)) 

cat("true negatives:", TN_bonferroni, "\n")

true negatives: 8999 

cat("false positives:", FP_bonferroni, "\n")

false positives: 1 

cat("false negatives:", FN_bonferroni, "\n")

false negatives: 918 

cat("true positives:", TP_bonferroni, "\n")

true positives: 82 

Without the correction, there is a higher rate of false positives in exchange for more true positives, making it less conservative. With the correction, the number of false positives drops dramatically from 479 to 1, indicating a much more strict threshold for calling something significant. As such, the correction minimizes the number of false positives but simultaneously identifies fewer true positives (902 vs 82). The trade-off with the correction is fewer false positives but also fewer true positives.

2.6 FDR Control with qvalue

qres_mix = qvalue::qvalue(p_mixed) #calculates qvalues for pvalues in mix

qres_mix$pi0

[1] 0.8601

When looking at the pi0 in the mixed simulations, we see quite a high proportion of true nulls. I simulated a mixture of null and alternative distributions, with 10% alternative hypotheses (prop_alt = 0.10). This means that 90% are null tests. 0.86 from the qvalues is close to 0.9, which confirms that the qvalue method is accurately estimating the proportion of true nulls.

#count how many tests are significant at an FDR threshold of 0.1
sig_tests_FDR <- sum(qres$qvalues <= 0.1)
cat("Number of significant tests at FDR threshold 0.1:", sig_tests_FDR, "\n")

Number of significant tests at FDR threshold 0.1: 867 

#count how TRUE positves
sig_tests_FDR_true <- sum((qres$qvalues <= 0.1) & !is_null)
cat("Number of true positives using FDR threshold", sig_tests_FDR_true, "\n")

Number of true positives using FDR threshold 762 

The number of true positives in the uncorrected is 902. With the correction, the number of true positives is 82. The FDR threshold has 867 positives, with 762 being true positives. The number of true positives using FDR is between that of the uncorrected and of the Bonferroni correction. FDR balances the uncorrected with the perhaps overly stringent threshold of the Bonferroni. FDR is useful for finding more true positives without excessively increasing false positives (such as in uncorrected).

The conceptual difference between FWER such as Bonferroni and FDR is that FWER strictly controls the probability of making any false positive errors while FDR controls the proportion of false positives but less strictly.