1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
|
#!/usr/bin/env python
#file test_alpha_diversity.py
from __future__ import division
from numpy import array, log, sqrt, exp
from math import e
from cogent.util.unit_test import TestCase, main
from cogent.maths.stats.alpha_diversity import expand_counts, counts, observed_species, singles, \
doubles, osd, margalef, menhinick, dominance, simpson, \
simpson_reciprocal, reciprocal_simpson,\
shannon, equitability, berger_parker_d, mcintosh_d, brillouin_d, \
strong, kempton_taylor_q, fisher_alpha, \
mcintosh_e, heip_e, simpson_e, robbins, robbins_confidence, \
chao1_uncorrected, chao1_bias_corrected, chao1, chao1_var, \
chao1_confidence, ACE, michaelis_menten_fit
__author__ = "Rob Knight"
__copyright__ = "Copyright 2007-2012, The Cogent Project"
__credits__ = ["Rob Knight","Justin Kuczynski"]
__license__ = "GPL"
__version__ = "1.5.3"
__maintainer__ = "Rob Knight"
__email__ = "rob@spot.colorado.edu"
__status__ = "Production"
class diversity_tests(TestCase):
"""Tests of top-level functions"""
def setUp(self):
"""Set up shared variables"""
self.TestData = array([0,1,1,4,2,5,2,4,1,2])
self.NoSingles = array([0,2,2,4,5,0,0,0,0,0])
self.NoDoubles = array([0,1,1,4,5,0,0,0,0,0])
def test_expand_counts(self):
"""expand_counts should return correct expanded array"""
c = array([2,0,1,2])
self.assertEqual(expand_counts(c), array([0,0,2,3,3]))
def test_counts(self):
"""counts should return correct array"""
c = array([5,0,1,1,5,5])
obs = counts(c)
exp = array([1,2,0,0,0,3])
self.assertEqual(obs, exp)
d = array([2,2,1,0])
obs = counts(d, obs)
exp = array([2,3,2,0,0,3])
self.assertEqual(obs, exp)
def test_singles(self):
"""singles should return correct # of singles"""
self.assertEqual(singles(self.TestData), 3)
self.assertEqual(singles(array([0,3,4])), 0)
self.assertEqual(singles(array([1])), 1)
def test_doubles(self):
"""doubles should return correct # of doubles"""
self.assertEqual(doubles(self.TestData), 3)
self.assertEqual(doubles(array([0,3,4])), 0)
self.assertEqual(doubles(array([2])), 1)
def test_osd(self):
"""osd should return correct # of observeds, singles, doubles"""
self.assertEqual(osd(self.TestData), (9,3,3))
def test_margalef(self):
"""margalef should match hand-calculated values"""
self.assertEqual(margalef(self.TestData), 8/log(22))
def test_menhinick(self):
"""menhinick should match hand-calculated values"""
self.assertEqual(menhinick(self.TestData), 9/sqrt(22))
def test_dominance(self):
"""dominance should match hand-calculated values"""
c = array([1,0,2,5,2])
self.assertFloatEqual(dominance(c), .34)
d = array([5])
self.assertEqual(dominance(d), 1)
def test_simpson(self):
"""simpson should match hand-calculated values"""
c = array([1,0,2,5,2])
self.assertFloatEqual(simpson(c), .66)
d = array([5])
self.assertFloatEqual(simpson(d), 0)
def test_reciprocal_simpson(self):
"""reciprocal_simpson should match hand-calculated results"""
c = array([1,0,2,5,2])
self.assertFloatEqual(reciprocal_simpson(c), 1/.66)
def test_simpson_reciprocal(self):
"""simpson_reciprocal should match 1/D results"""
c = array([1,0,2,5,2])
self.assertFloatEqual(simpson_reciprocal(c), 1./dominance(c))
def test_shannon(self):
"""shannon should match hand-calculated values"""
c = array([5])
self.assertFloatEqual(shannon(c), 0)
c = array([5,5])
self.assertFloatEqual(shannon(c), 1)
c = array([1,1,1,1,0])
self.assertEqual(shannon(c), 2)
def test_equitability(self):
"""equitability should match hand-calculated values"""
c = array([5])
self.assertFloatEqual(equitability(c), 0)
c = array([5,5])
self.assertFloatEqual(equitability(c), 1)
c = array([1,1,1,1,0])
self.assertEqual(equitability(c), 1)
def test_berger_parker_d(self):
"""berger-parker_d should match hand-calculated values"""
c = array([5])
self.assertFloatEqual(berger_parker_d(c), 1)
c = array([5,5])
self.assertFloatEqual(berger_parker_d(c), 0.5)
c = array([1,1,1,1,0])
self.assertEqual(berger_parker_d(c), 0.25)
def test_mcintosh_d(self):
"""mcintosh_d should match hand-calculated values"""
c = array([1,2,3])
self.assertFloatEqual(mcintosh_d(c), 0.636061424871458)
def test_brillouin_d(self):
"""brillouin_d should match hand-calculated values"""
c = array([1,2,3,1])
self.assertFloatEqual(brillouin_d(c), 0.86289353018248782)
def test_strong(self):
"""strong's dominance index should match hand-calculated values"""
c = array([1,2,3,1])
self.assertFloatEqual(strong(c), 0.214285714)
def test_kempton_taylor_q(self):
"""kempton_taylor_q should approximate Magurran 1998 calculation p143"""
c = array([2,3,3,3,3,3,4,4,4,6,6,7,7,9,9,11,14,15,15,20,29,33,34,
36,37,53,57,138,146,170])
self.assertFloatEqual(kempton_taylor_q(c), 14/log(34/4))
def test_fisher_alpha(self):
"""fisher alpha should match hand-calculated value."""
c = array([4,3,4,0,1,0,2])
obs = fisher_alpha(c)
self.assertFloatEqual(obs, 2.7823795367398798)
def test_mcintosh_e(self):
"""mcintosh e should match hand-calculated value."""
c = array([1,2,3,1])
num = sqrt(15)
den = sqrt(19)
exp = num/den
self.assertEqual(mcintosh_e(c), exp)
def test_heip_e(self):
"""heip e should match hand-calculated value"""
c = array([1,2,3,1])
h = shannon(c, base=e)
expected = exp(h-1)/3
self.assertEqual(heip_e(c), expected)
def test_simpson_e(self):
"""simpson e should match hand-calculated value"""
c = array([1,2,3,1])
s = simpson(c)
self.assertEqual((1/s)/4, simpson_e(c))
def test_robbins(self):
"""robbins metric should match hand-calculated value"""
c = array([1,2,3,0,1])
r = robbins(c)
self.assertEqual(r,2./7)
def test_robbins_confidence(self):
"""robbins CI should match hand-calculated value"""
c = array([1,2,3,0,1])
r = robbins_confidence(c, 0.05)
n = 7
s = 2
k = sqrt(8/0.05)
self.assertEqual(r, ((s-k)/(n+1), (s+k)/(n+1)))
def test_observed_species(self):
"""observed_species should return # observed species"""
c = array([4,3,4,0,1,0,2])
obs = observed_species(c)
exp = 5
self.assertEqual(obs, exp)
c = array([0,0,0])
obs = observed_species(c)
exp = 0
self.assertEqual(obs, exp)
self.assertEqual(observed_species(self.TestData), 9)
def test_chao1_bias_corrected(self):
"""chao1_bias_corrected should return same result as EstimateS"""
obs = chao1_bias_corrected(*osd(self.TestData))
self.assertEqual(obs, 9.75)
def test_chao1_uncorrected(self):
"""chao1_uncorrected should return same result as EstimateS"""
obs = chao1_uncorrected(*osd(self.TestData))
self.assertEqual(obs, 10.5)
def test_chao1(self):
"""chao1 should use right decision rules"""
self.assertEqual(chao1(self.TestData), 9.75)
self.assertEqual(chao1(self.TestData,bias_corrected=False),10.5)
self.assertEqual(chao1(self.NoSingles), 4)
self.assertEqual(chao1(self.NoSingles,bias_corrected=False),4)
self.assertEqual(chao1(self.NoDoubles), 5)
self.assertEqual(chao1(self.NoDoubles,bias_corrected=False),5)
def test_chao1_var(self):
"""chao1_var should match observed results from EstimateS"""
#NOTE: EstimateS reports sd, not var, and rounds to 2 dp
self.assertFloatEqual(chao1_var(self.TestData), 1.42**2, eps=0.01)
self.assertFloatEqual(chao1_var(self.TestData,bias_corrected=False),\
2.29**2, eps=0.01)
self.assertFloatEqualAbs(chao1_var(self.NoSingles), 0.39**2, eps=0.01)
self.assertFloatEqualAbs(chao1_var(self.NoSingles, \
bias_corrected=False), 0.39**2, eps=0.01)
self.assertFloatEqualAbs(chao1_var(self.NoDoubles), 2.17**2, eps=0.01)
self.assertFloatEqualAbs(chao1_var(self.NoDoubles, \
bias_corrected=False), 2.17**2, eps=0.01)
def test_chao1_confidence(self):
"""chao1_confidence should match observed results from EstimateS"""
#NOTE: EstimateS rounds to 2 dp
self.assertFloatEqual(chao1_confidence(self.TestData), (9.07,17.45), \
eps=0.01)
self.assertFloatEqual(chao1_confidence(self.TestData, \
bias_corrected=False), (9.17,21.89), eps=0.01)
self.assertFloatEqualAbs(chao1_confidence(self.NoSingles),\
(4, 4.95), eps=0.01)
self.assertFloatEqualAbs(chao1_confidence(self.NoSingles, \
bias_corrected=False), (4,4.95), eps=0.01)
self.assertFloatEqualAbs(chao1_confidence(self.NoDoubles), \
(4.08,17.27), eps=0.01)
self.assertFloatEqualAbs(chao1_confidence(self.NoDoubles, \
bias_corrected=False), (4.08,17.27), eps=0.01)
def test_ACE(self):
"""ACE should match values calculated by hand"""
self.assertFloatEqual(ACE(array([2,0])), 1.0, eps=0.001)
# next: just returns the number of species when all are abundant
self.assertFloatEqual(ACE(array([12,0,9])), 2.0, eps=0.001)
self.assertFloatEqual(ACE(array([12,2,8])), 3.0, eps=0.001)
self.assertFloatEqual(ACE(array([12,2,1])), 4.0, eps=0.001)
self.assertFloatEqual(ACE(array([12,1,2,1])), 7.0, eps=0.001)
self.assertFloatEqual(ACE(array([12,3,2,1])), 4.6, eps=0.001)
self.assertFloatEqual(ACE(array([12,3,6,1,10])), 5.62749672, eps=0.001)
def test_michaelis_menten_fit(self):
""" michaelis_menten_fit should match hand values in limiting cases"""
res = michaelis_menten_fit([22])
self.assertFloatEqual(res,1.0,eps=.01)
res = michaelis_menten_fit([42])
self.assertFloatEqual(res,1.0,eps=.01)
res = michaelis_menten_fit([34],num_repeats=3,params_guess=[13,13])
self.assertFloatEqual(res,1.0,eps=.01)
res = michaelis_menten_fit([70,70],num_repeats=5)
self.assertFloatEqual(res,2.0,eps=.01)
if __name__ == '__main__':
main()
|
