Anonymous user
P-value correction: Difference between revisions
added Ruby
(added Ruby) |
|||
Line 4,706:
Hommel</pre>
=={{header|Ruby}}==
{{trans|Perl}}
<lang ruby>def pmin(array)
x = 1
pmin_array = []
array.each_index do |i|
pmin_array[i] = [array[i], x].min
abort if pmin_array[i] > 1
end
pmin_array
end
def cummin(array)
cumulative_min = array[0]
arr_cummin = []
array.each do |p|
cumulative_min = [p, cumulative_min].min
arr_cummin.push(cumulative_min)
end
arr_cummin
end
def cummax(array)
cumulative_max = array[0]
arr_cummax = []
array.each do |p|
cumulative_max = [p, cumulative_max].max
arr_cummax.push(cumulative_max)
end
arr_cummax
end
# decreasing variable is optional
def order(array, decreasing = false)
if decreasing == false
array.sort.map { |n| array.index(n) }
else
array.sort.map { |n| array.index(n) }.reverse
end
end
def p_adjust(arr_pvalues, method = 'Holm')
lp = arr_pvalues.size
n = lp
if method.casecmp('hochberg').zero?
arr_o = order(arr_pvalues, true)
arr_cummin_input = []
(0..n).each do |index|
arr_cummin_input[index] = (index + 1) * arr_pvalues[arr_o[index].to_i]
end
arr_cummin = cummin(arr_cummin_input)
arr_pmin = pmin(arr_cummin)
arr_ro = order(arr_o)
return arr_pmin.values_at(*arr_ro)
elsif method.casecmp('bh').zero? || method.casecmp('benjamini-hochberg').zero?
arr_o = order(arr_pvalues, true)
arr_cummin_input = []
(0..(n - 1)).each do |i|
arr_cummin_input[i] = (n / (n - i).to_f) * arr_pvalues[arr_o[i]]
end
arr_ro = order(arr_o)
arr_cummin = cummin(arr_cummin_input)
arr_pmin = pmin(arr_cummin)
return arr_pmin.values_at(*arr_ro)
elsif method.casecmp('by').zero? || method.casecmp('benjamini-yekutieli').zero?
q = 0.0
arr_o = order(arr_pvalues, true)
arr_ro = order(arr_o)
(1..n).each do |index|
q += 1.0 / index
end
arr_cummin_input = []
(0..(n - 1)).each do |i|
arr_cummin_input[i] = q * (n / (n - i).to_f) * arr_pvalues[arr_o[i]]
end
arr_cummin = cummin(arr_cummin_input)
arr_pmin = pmin(arr_cummin)
return arr_pmin.values_at(*arr_ro)
elsif method.casecmp('bonferroni').zero?
arr_qvalues = []
(0..(n - 1)).each do |i|
q = arr_pvalues[i] * n
if (q >= 0) && (q < 1)
arr_qvalues[i] = q
elsif q >= 1
arr_qvalues[i] = 1.0
else
puts "Falied to get Bonferroni adjusted p for #{arr_pvalues[i]}"
end
end
return arr_qvalues
elsif method.casecmp('holm').zero?
o = order(arr_pvalues)
cummax_input = []
(0..(n - 1)).each do |index|
cummax_input[index] = (n - index) * arr_pvalues[o[index]]
end
ro = order(o)
arr_cummax = cummax(cummax_input)
arr_pmin = pmin(arr_cummax)
return arr_pmin.values_at(*ro)
elsif method.casecmp('hommel').zero?
o = order(arr_pvalues)
arr_p = arr_pvalues.values_at(*o)
ro = order(o)
q = []
pa = []
min = n * arr_p[0]
(0..(n - 1)).each do |index|
temp = n * arr_p[index] / (index + 1)
min = [min, temp].min
end
(0..(n - 1)).each do |index|
pa[index] = min
q[index] = min
end
j = n - 1
while j >= 2
ij = Array 0..(n - j)
i2_length = j - 1
i2 = []
(0..(i2_length - 1)).each do |i|
i2[i] = n - j + 2 + i - 1 # R's indices are 1-based, C's are 0-based
end
q1 = j * arr_p[i2[0]] / 2.0
(1..(i2_length - 1)).each do |i|
temp_q1 = j * arr_p[i2[i]] / (2 + i)
q1 = [temp_q1, q1].min
end
(0..(n - j)).each do |i|
tmp = j * arr_p[ij[i]]
q[ij[i]] = [tmp, q1].min
end
(0..(i2_length - 1)).each do |i|
q[i2[i]] = q[n - j]
end
(0..(n - 1)).each do |i|
pa[i] = q[i] if pa[i] < q[i]
end
j -= 1
end
return pa.values_at(*ro)
else
puts "#{method} isn't accepted."
abort
end
end
pvalues =
[4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02,
1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01,
4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01,
3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01,
4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02,
3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01,
9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02,
7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04,
1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04,
8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04,
1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03,
8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05,
2.900813e-02, 5.735490e-03]
correct_answers = {
'Benjamini-Hochberg' => [6.126681e-01, 8.521710e-01, 1.987205e-01,
1.891595e-01, 3.217789e-01, 9.301450e-01,
4.870370e-01, 9.301450e-01, 6.049731e-01,
6.826753e-01, 6.482629e-01, 7.253722e-01,
5.280973e-01, 8.769926e-01, 4.705703e-01,
9.241867e-01, 6.049731e-01, 7.856107e-01,
4.887526e-01, 1.136717e-01, 4.991891e-01,
8.769926e-01, 9.991834e-01, 3.217789e-01,
9.301450e-01, 2.304958e-01, 5.832475e-01,
3.899547e-02, 8.521710e-01, 1.476843e-01,
1.683638e-02, 2.562902e-03, 3.516084e-02,
6.250189e-02, 3.636589e-03, 2.562902e-03,
2.946883e-02, 6.166064e-03, 3.899547e-02,
2.688991e-03, 4.502862e-04, 1.252228e-05,
7.881555e-02, 3.142613e-02, 4.846527e-03,
2.562902e-03, 4.846527e-03, 1.101708e-03,
7.252032e-02, 2.205958e-02],
'Benjamini-Yekutieli' => [1.000000e+00, 1.000000e+00, 8.940844e-01,
8.510676e-01, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 5.114323e-01, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00,
1.754486e-01, 1.000000e+00, 6.644618e-01,
7.575031e-02, 1.153102e-02, 1.581959e-01,
2.812089e-01, 1.636176e-02, 1.153102e-02,
1.325863e-01, 2.774239e-02, 1.754486e-01,
1.209832e-02, 2.025930e-03, 5.634031e-05,
3.546073e-01, 1.413926e-01, 2.180552e-02,
1.153102e-02, 2.180552e-02, 4.956812e-03,
3.262838e-01, 9.925057e-02],
'Bonferroni' => [1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 7.019185e-01,
1.000000e+00, 1.000000e+00, 2.020365e-01, 1.516674e-02,
5.625735e-01, 1.000000e+00, 2.909271e-02, 1.537741e-02,
4.125636e-01, 6.782670e-02, 6.803480e-01, 1.882294e-02,
9.005725e-04, 1.252228e-05, 1.000000e+00, 4.713920e-01,
4.395577e-02, 1.088915e-02, 4.846527e-02, 3.305125e-03,
1.000000e+00, 2.867745e-01],
'Hochberg' => [9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 4.632662e-01,
9.991834e-01, 9.991834e-01, 1.575885e-01, 1.383967e-02,
3.938014e-01, 7.600230e-01, 2.501973e-02, 1.383967e-02,
3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01,
3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03,
8.992520e-01, 2.179486e-01],
'Holm' => [1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
1.000000e+00, 1.000000e+00, 1.000000e+00, 4.632662e-01,
1.000000e+00, 1.000000e+00, 1.575885e-01, 1.395341e-02,
3.938014e-01, 7.600230e-01, 2.501973e-02, 1.395341e-02,
3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01,
3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03,
8.992520e-01, 2.179486e-01],
'Hommel' => [9.991834e-01, 9.991834e-01, 9.991834e-01, 9.987624e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.595180e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
9.991834e-01, 9.991834e-01, 9.991834e-01, 4.351895e-01,
9.991834e-01, 9.766522e-01, 1.414256e-01, 1.304340e-02,
3.530937e-01, 6.887709e-01, 2.385602e-02, 1.322457e-02,
2.722920e-01, 5.426136e-02, 4.218158e-01, 1.581127e-02,
8.825610e-04, 1.252228e-05, 8.743649e-01, 3.016908e-01,
3.516461e-02, 9.582456e-03, 3.877222e-02, 3.172920e-03,
8.122276e-01, 1.950067e-01]
}
# correct_answers.each do |method, answers|
methods = ['Benjamini-Yekutieli', 'Benjamini-Hochberg', 'Hochberg',
'Bonferroni', 'Holm', 'Hommel']
methods.each do |method|
puts method
error = 0.0
arr_q = p_adjust(pvalues, method)
arr_q.each_index do |p|
error += (correct_answers[method][p] - arr_q[p])
end
puts "total error for #{method} = #{error}"
end
</lang>
{{out}}
<pre>Benjamini-Yekutieli
total error for Benjamini-Yekutieli = -1.7373780825929845e-07
Benjamini-Hochberg
total error for Benjamini-Hochberg = -1.4736877299143964e-08
Hochberg
total error for Hochberg = -1.2354999978398105e-07
Bonferroni
total error for Bonferroni = 4.49999999152751e-08
Holm
total error for Holm = -1.163499997815258e-07
Hommel
total error for Hommel = 1.1483094955369324e-07
</pre>
=={{header|SAS}}==
| |||