2
$\begingroup$

I have a couple arrays of dN/dS scores, and I would like to calculate the confidence interval for each array of data. dN/dS scores are not normally distributed but are log-normally distrbuted, so I want some non-parametric approach for finding confidence intervals.

An answer in the CrossValidated forum suggested that in R you can use a Wilcoxon test to find the median-based confidence interval, a function that I do not think the python scipy.stats.wilcoxon() package does. Is there anyway of doing something similar in Python for the confidence intervals of log-normally distributed data?

For a example short array of dN/dS scores, here is the dN/dS scores of ATPases:

array([0.02288081, 0.44170839, 0.10549733, 0.17515196, 0.09449279,
       0.07110412, 0.00893079, 0.23485109, 0.14533192, 0.05449631,
       0.10281173, 0.02113355, 0.05157087, 0.01705113, 0.02651948,
       0.09352947, 0.05828018, 0.20528157, 0.02873843, 0.0141598 ,
       0.11262881, 0.06337332, 0.13689815, 0.10557703, 0.04452507,
       0.05046484, 0.40241456, 0.05939199, 0.24423569, 0.05042892,
       0.164257  , 0.03408215, 0.04737262, 0.01037463, 0.01246448,
       0.02170375, 0.0241773 , 0.05136936, 0.02393366, 0.18913979,
       0.15781334, 0.06448557, 0.04355384, 0.02821125, 0.08015629,
       0.10985432, 0.06074574, 0.15775976, 0.05678278, 0.03749782,
       0.05518756, 0.00770479, 0.21248167, 0.08005044, 0.16307954,
       0.05783565, 0.05907416, 0.07044622, 0.13227131, 0.01627556,
       0.10859962, 0.08149819, 0.05600647, 0.16098728, 0.15183062,
       0.05202344, 0.01769589, 0.00789287, 0.07777749, 0.1324942 ,
       0.12734709, 0.17146938, 0.03890857, 0.1296019 , 0.085146  ,
       0.05602965, 0.02708089, 0.11807038, 0.07848828, 0.03291032,
       0.36302533, 0.07800343, 0.06551307, 0.04676282, 0.04765273,
       0.08060882, 0.06339636, 0.03349833, 0.01224308, 0.1481316 ,
       0.31738452, 0.15690855, 0.0693822 , 0.020425  , 0.02909208,
       0.03499225, 0.03019904, 0.13722717, 0.14403507, 0.01257245,
       0.02223452, 0.07068784, 0.0544813 , 0.08738558, 0.02884046,
       0.10549474, 0.06695546, 0.01341142, 0.09440411, 0.11840834,
       0.08558889, 0.2688645 , 0.28313546, 0.15127967, 0.01463191,
       0.1728421 ])
$\endgroup$
1
  • $\begingroup$ I am not at sure about that, I would be amazed if SciPy couldn't do this. However in R is a base function. How many columns, comparisons? What are you measuring? $\endgroup$
    – M__
    Oct 13 '19 at 12:15
2
$\begingroup$

I guess you mean confidence interval of the mean (or median). I looked at the array you provided, even after log it does not look normal-distributed. Use the bootstrap for non-parametric values. I modified the code byJames Triveri, see example below, and you can use the mean or median. The textbook by Efron is a good reference for the method, which is also used in R package boot. There's also awebsite which should be ok.


import numpy as np

def bootstrap(data, n=1000, func=np.mean,p=0.95):

    sample_size = len(data)
    simulations = [func(np.random.choice(data, size=sample_size, replace=True)) for i in range(n)]
    simulations.sort()
    u_pval = (1+p)/2.
    l_pval = (1-u_pval)
    l_indx = int(np.floor(n*l_pval))
    u_indx = int(np.floor(n*u_pval))
    return(simulations[l_indx],simulations[u_indx])

x = np.array([0.02288081, 0.44170839, 0.10549733, 0.17515196, 0.09449279,
       0.07110412, 0.00893079, 0.23485109, 0.14533192, 0.05449631,
       0.10281173, 0.02113355, 0.05157087, 0.01705113, 0.02651948,
       0.09352947, 0.05828018, 0.20528157, 0.02873843, 0.0141598 ,
       0.11262881, 0.06337332, 0.13689815, 0.10557703, 0.04452507,
       0.05046484, 0.40241456, 0.05939199, 0.24423569, 0.05042892,
       0.164257  , 0.03408215, 0.04737262, 0.01037463, 0.01246448,
       0.02170375, 0.0241773 , 0.05136936, 0.02393366, 0.18913979,
       0.15781334, 0.06448557, 0.04355384, 0.02821125, 0.08015629,
       0.10985432, 0.06074574, 0.15775976, 0.05678278, 0.03749782,
       0.05518756, 0.00770479, 0.21248167, 0.08005044, 0.16307954,
       0.05783565, 0.05907416, 0.07044622, 0.13227131, 0.01627556,
       0.10859962, 0.08149819, 0.05600647, 0.16098728, 0.15183062,
       0.05202344, 0.01769589, 0.00789287, 0.07777749, 0.1324942 ,
       0.12734709, 0.17146938, 0.03890857, 0.1296019 , 0.085146  ,
       0.05602965, 0.02708089, 0.11807038, 0.07848828, 0.03291032,
       0.36302533, 0.07800343, 0.06551307, 0.04676282, 0.04765273,
       0.08060882, 0.06339636, 0.03349833, 0.01224308, 0.1481316 ,
       0.31738452, 0.15690855, 0.0693822 , 0.020425  , 0.02909208,
       0.03499225, 0.03019904, 0.13722717, 0.14403507, 0.01257245,
       0.02223452, 0.07068784, 0.0544813 , 0.08738558, 0.02884046,
       0.10549474, 0.06695546, 0.01341142, 0.09440411, 0.11840834,
       0.08558889, 0.2688645 , 0.28313546, 0.15127967, 0.01463191,
       0.1728421 ])

# for mean,95% confidence interval, do
bootstrap(x,1000,np.mean,0.95)
# for median, 95% confidence interval, do 
bootstrap(x,1000,np.median,0.95)

$\endgroup$
3
  • $\begingroup$ This is a good answer and it would certainly work. This is the 95% CI of the mean and it's not the Wilcoxan. If the OP is using W it suggests the data is separated into columns and they are trying to determine the difference because different groups. $\endgroup$
    – M__
    Oct 13 '19 at 16:14
  • 1
    $\begingroup$ Essentially the question isn't clear enough at present $\endgroup$
    – M__
    Oct 13 '19 at 16:15
  • 1
    $\begingroup$ I agree with you @MichaelG. You can replace np.mean with np.median, its the same concept. Yes, let's hope VeritatemAmo provides more info. $\endgroup$
    – StupidWolf
    Oct 13 '19 at 16:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.