Hi,

I am computing the lambda statistics in R to check whether the distribution of my p-values follow a null distribution using: 

    pvalue <- runif(1000, min=0, max=1)
    chisq <- qchisq(1 - pvalue, 1)
    lambda <- median(chisq) / qchisq(0.5, 1)

Is there a procedure for constructing the confidence interval for this lambda genomic factor in R? 

Thanks

I encourage you to follow up on my points here with a statistician in your department. Always better to get more brains looking at the same thing. There are various ways to calculate confidence intervals for series of data.

You could start with a typical formula for the standard error of the mean (SE mean), which is the standard deviation divided by the square root of the sample n:

    σ / sqrt (n)


Lambda is a median value, so, we have to convert the SE mean to SE median:

    1.253 * (σ / sqrt (n))

(various articles on the WWW validate this)


As you are dealing with a Chi-squared quantile distribution, in order to calculate the confidence intervals, I believe you can do this with the following:

    pvalue <- runif(1000000, min=0, max=1)
    chisq <- qchisq(1 - pvalue, 1)
    lambda <- median(chisq) / qchisq(0.5, 1)
    
    lambda 
    [1] 0.9963626

    SE.median <- qchisq(0.975, 1) * (1.253 * ( sd(chisq) / sqrt( length(chisq) ) ) )
    
    SE.median
    [1] 0.008914307
     
    lambda + (SE.median / qchisq(0.5, 1))
    [1] 1.015957
     
    lambda - (SE.median / qchisq(0.5, 1))
    [1] 0.976768


I increased the number of random p-values to 1 million just to make the confidence intervals tighter. At 1000 (your original post), they were pretty wide.

Take this as a first step into this. I encourage you to look up all functions / formulae that I have used.

Kevin

