Skip to content

Commit

Permalink
Cerberus validation of observations; sd->std
Browse files Browse the repository at this point in the history
  • Loading branch information
@rgerkin
rgerkin committed Aug 9, 2018
1 parent 52528e5 commit 3223517
Show file tree
Hide file tree
Showing 13 changed files with 310 additions and 336 deletions.
48 changes: 24 additions & 24 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
| [![Travis](https://travis-ci.org/scidash/neuronunit.svg?branch=master)](https://travis-ci.org/scidash/neuronunit) | [![Travis](https://travis-ci.org/scidash/neuronunit.svg?branch=dev)](https://travis-ci.org/scidash/neuronunit) |
| [![RTFD](https://readthedocs.org/projects/neuronunit/badge/?version=master)](http://neuronunit.readthedocs.io/en/latest/?badge=master) | [![RTFD](https://readthedocs.org/projects/neuronunit/badge/?version=dev)](http://neuronunit.readthedocs.io/en/latest/?badge=dev) |
| [![Coveralls](https://coveralls.io/repos/github/scidash/neuronunit/badge.svg?branch=master)](https://coveralls.io/github/scidash/neuronunit?branch=master) | [![Coveralls](https://coveralls.io/repos/github/scidash/neuronunit/badge.svg?branch=dev)](https://coveralls.io/github/scidash/neuronunit?branch=dev) |
| [![Requirements](https://requires.io/github/scidash/neuronunit/requirements.svg?branch=master)](https://requires.io/github/scidash/neuronunit/requirements/?branch=master) | [![Requirements](https://requires.io/github/scidash/neuronunit/requirements.svg?branch=dev)](https://requires.io/github/scidash/neuronunit/requirements/?branch=dev) |
| [![Requirements](https://requires.io/github/scidash/neuronunit/requirements.svg?branch=master)](https://requires.io/github/scidash/neuronunit/requirements/?branch=master) | [![Requirements](https://requires.io/github/scidash/neuronunit/requirements.svg?branch=dev)](https://requires.io/github/scidash/neuronunit/requirements/?branch=dev) |
| [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/scidash/neuronunit/master) |


Expand Down Expand Up @@ -44,13 +44,13 @@ model = ChannelModel(channel_file_path, channel_index=0, name=channel_model_name
# Get the experiment data from ChannelWorm and instantiate the test
doi = '10.1083/jcb.200203055'
fig = '2B'
sample_data = GraphData.objects.get(graph__experiment__reference__doi=doi,
sample_data = GraphData.objects.get(graph__experiment__reference__doi=doi,
graph__figure_ref_address=fig)
voltage, current_per_farad = sample_data.asunitedarray()
patch_capacitance = pq.Quantity(1e-13,'F') # Assume recorded patch had this capacitance;
patch_capacitance = pq.Quantity(1e-13,'F') # Assume recorded patch had this capacitance;
# an arbitrary scaling factor.
current = current_per_farad * patch_capacitance
observation = {'v':voltage,
observation = {'v':voltage,
'i':current}
test = IVCurvePeakTest(observation)

Expand All @@ -63,7 +63,7 @@ score.plot(rd['v'],rd['i_pred'],same_fig=True,color='r',label='Predicted (model)
![png](https://raw.githubusercontent.com/scidash/assets/master/figures/SCU_IVCurve_Model_6_0.png)

```
score.summarize()
score.summarize()
""" OUTPUT:
Model EGL-19.channel (ChannelModel) achieved score Fail on test 'IV Curve Test (IVCurvePeakTest)'. ===
"""
Expand All @@ -90,27 +90,27 @@ neurolex_id = 'nifext_128' # Cerebellar Granule Cell
# Specify reference data for a test of resting potential for a granule cell.
reference_data = neuroelectro.NeuroElectroSummary(
neuron = {'nlex_id':neurolex_id}, # Neuron type.
ephysprop = {'name':'Resting Membrane Potential'}) # Electrophysiological property name.
# Get and verify summary data for the combination above from neuroelectro.org.
ephysprop = {'name':'Resting Membrane Potential'}) # Electrophysiological property name.
# Get and verify summary data for the combination above from neuroelectro.org.
reference_data.get_values()
vm_test = tests.RestingPotentialTest(
observation = {'mean':reference_data.mean,
'std':reference_data.std},
'sd':reference_data.std},
name = 'Resting Potential')
# Specify reference data for a test of action potential width.
reference_data = neuroelectro.NeuroElectroSummary(
neuron = {'nlex_id':neurolex_id}, # Neuron type.
ephysprop = {'name':'Spike Half-Width'}) # Electrophysiological property name.
# Get and verify summary data for the combination above from neuroelectro.org.
ephysprop = {'name':'Spike Half-Width'}) # Electrophysiological property name.
# Get and verify summary data for the combination above from neuroelectro.org.
reference_data.get_values()
spikewidth_test = tests.InjectedCurrentAPWidthTest(
observation = {'mean':reference_data.mean,
'std':reference_data.std},
'sd':reference_data.std},
name = 'Spike Width',
params={'injected_square_current':{'amplitude':5.3*pq.pA,
'delay':50.0*pq.ms,
'duration':500.0*pq.ms}})
'duration':500.0*pq.ms}})
# 5.3 pA of injected current in a 500 ms square pulse.
# Create a test suite from these two tests.
Expand All @@ -122,7 +122,7 @@ for model_name in model_names # Iterate through a list of models downloaded from
model = nc_models.OSBModel(*model_info)
models.append(model) # Add to the list of models to be tested.
score_matrix = suite.judge(models,stop_on_error=True)
score_matrix = suite.judge(models,stop_on_error=True)
score_matrix.view()
```
### Score Matrix for Test Suite 'Neuron Tests'
Expand Down Expand Up @@ -153,7 +153,7 @@ NeuronUnit is based on [SciUnit](http://github.com/scidash/sciunit), a disciplin
2. Check that the model has the capabilities required to take the test.
1. Make the model take the test.
2. Generate a score from that test run.
1. Bind the score to the specific model/test combination and any related data from test execution.
1. Bind the score to the specific model/test combination and any related data from test execution.
1. Visualize the score (i.e. print or display the result of the test).
Here, we will break down how this is accomplished in NeuronUnit. Although NeuronUnit contains several model and test classes that make it easy to work with standards in neuron modeling and electrophysiology data reporting, here we will use toy model and test classes constructed on-the-fly so the process of model and test construction is fully transparent.
Expand All @@ -180,7 +180,7 @@ Let's see what the `neuronunit.capabilities.ProducesMembranePotential` capabilit
```python
class ProducesMembranePotential(Capability):
"""Indicates that the model produces a somatic membrane potential."""

def get_membrane_potential(self):
"""Must return a neo.core.AnalogSignal."""
raise NotImplementedError()
Expand All @@ -206,12 +206,12 @@ Now we can then construct a simple test to use on this model or any other test t
```python
class ToyAveragePotentialTest(sciunit.Test):
"""Tests the average membrane potential of a neuron."""

def __init__(self,
observation={'mean':None,'std':None},
observation={'mean':None,'sd':None},
name="Average potential test"):
"""Takes the mean and standard deviation of reference membrane potentials."""

sciunit.Test.__init__(self,observation,name) # Call the base constructor.
self.required_capabilities += (neuronunit.capabilities.ProducesMembranePotential,)
# This test will require a model to express the above capabilities
Expand All @@ -220,18 +220,18 @@ class ToyAveragePotentialTest(sciunit.Test):
score_type = sciunit.scores.ZScore # The test will return this kind of score.

def validate_observation(self, observation):
"""An optional method that makes sure an observation to be used as
"""An optional method that makes sure an observation to be used as
reference data has the right form"""
try:
assert type(observation['mean']) is quantities.Quantity # From the 'quantities' package
assert type(observation['std']) is quantities.Quantity
assert type(observation['sd']) is quantities.Quantity
except Exception as e:
raise sciunit.ObservationError(("Observation must be of the form "
"{'mean':float*mV,'std':float*mV}"))
"{'mean':float*mV,'sd':float*mV}"))

def generate_prediction(self, model):
"""Implementation of sciunit.Test.generate_prediction."""
vm = model.get_median_vm() # If the model has the capability 'ProducesMembranePotential',
vm = model.get_median_vm() # If the model has the capability 'ProducesMembranePotential',
# then it implements this method
prediction = {'mean':vm}
return prediction
Expand All @@ -248,8 +248,8 @@ The test constructor takes an observation to parameterize the test, e.g.:

```python
from quantities import mV
my_observation = {'mean':-60.0*mV,
'std':3.5*mV}
my_observation = {'mean':-60.0*mV,
'sd':3.5*mV}
my_average_potential_test = ToyAveragePotentialTest(my_observation, name='my_average_potential_test')
```

Expand Down
62 changes: 31 additions & 31 deletions neuronunit/neuroelectro.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -317,7 +317,7 @@ def get_observation(self, params=None, show=False):
class NeuroElectroPooledSummary(NeuroElectroDataMap):
"""Class for getting summary values from neuroelectro.org.
Values are computed by pooling each report's mean and SD across reports.
Values are computed by pooling each report's mean and std across reports.
"""

def get_values(self, params=None, quiet=False):
Expand Down Expand Up @@ -348,7 +348,7 @@ def get_values(self, params=None, quiet=False):
self.neuron_name = data[0]['ncm']['n']['name']
self.ephysprop_name = data[0]['ecm']['e']['name']

# Pool each paper by weighing each mean by the paper's N and SD
# Pool each paper by weighing each mean by the paper's N and std
stats = self.get_pooled_stats(data, quiet)

self.mean = stats['mean']
Expand Down Expand Up @@ -383,7 +383,7 @@ def get_pooled_stats(self, data, quiet=True):
lines = []
means = []
sems = []
sds = []
stds = []
ns = []
sources = []

Expand All @@ -393,54 +393,54 @@ def get_pooled_stats(self, data, quiet=True):
# Collect raw values for each paper from NeuroElectro
for item in data:

err_is_sem = item['error_type'] == "sem" # SEM or SD
err_is_sem = item['error_type'] == "sem" # SEM or std
err = (item['err_norm'] if item['err_norm'] is not None
else item['err'])

sem = err if err_is_sem else None
sd = err if not err_is_sem else None
std = err if not err_is_sem else None
mean = (item['val_norm'] if item['val_norm'] is not None
else item['val'])
n = item['n']
source = item['source']

means.append(mean)
sems.append(sem)
sds.append(sd)
stds.append(std)
ns.append(n)
sources.append(source)

if not quiet:
print({'mean': mean, 'std': sd, 'sem': sem, 'n': n})
print({'mean': mean, 'std': std, 'sem': sem, 'n': n})

# Fill in missing values
self.fill_missing_ns(ns)
self.fill_missing_sems_sds(sems, sds, ns)
self.fill_missing_sems_stds(sems, stds, ns)

if not quiet:
print("---------------------------------------------------")
print("Filled in Values (computed or median where missing)")

for i, _ in enumerate(means):
line = {'mean': means[i], 'sd': sds[i], 'sem': sems[i],
line = {'mean': means[i], 'std': stds[i], 'sem': sems[i],
'n': ns[i], "source": sources[i]}
lines.append(line)
print(line)

# Compute the weighted grand_mean
# grand_mean = SUM( N[i]*Mean[i] ) / SUM(N[i])
# grand_sd = SQRT( SUM( (N[i]-1)*SD[i]^2 ) / SUM(N[i]-1) )
# grand_std = SQRT( SUM( (N[i]-1)*std[i]^2 ) / SUM(N[i]-1) )

ns_np = np.array(ns)
means_np = np.array(means)
sds_np = np.array(sds)
stds_np = np.array(stds)
n_sum = ns_np.sum()

grand_mean = np.sum(ns_np * means_np) / n_sum
grand_sd = np.sqrt(np.sum((ns_np-1)*(sds_np**2))/np.sum(ns_np-1))
grand_sem = grand_sd / np.sqrt(n_sum)
grand_std = np.sqrt(np.sum((ns_np-1)*(stds_np**2))/np.sum(ns_np-1))
grand_sem = grand_std / np.sqrt(n_sum)

return {'mean': grand_mean, 'sem': grand_sem, 'std': grand_sd,
return {'mean': grand_mean, 'sem': grand_sem, 'std': grand_std,
'n': n_sum, 'items': lines}

def fill_missing_ns(self, ns):
Expand All @@ -456,33 +456,33 @@ def fill_missing_ns(self, ns):
if ns[i] is None:
ns[i] = n_median

def fill_missing_sems_sds(self, sems, sds, ns):
"""Fill in computable sems/sds."""
def fill_missing_sems_stds(self, sems, stds, ns):
"""Fill in computable sems/stds."""
for i, _ in enumerate(sems):

# Check if sem or sd is computable
if sems[i] is None and sds[i] is not None:
sems[i] = sds[i] / np.sqrt(ns[i])
# Check if sem or std is computable
if sems[i] is None and stds[i] is not None:
sems[i] = stds[i] / np.sqrt(ns[i])

if sds[i] is None and sems[i] is not None:
sds[i] = sems[i] * np.sqrt(ns[i])
if stds[i] is None and sems[i] is not None:
stds[i] = sems[i] * np.sqrt(ns[i])

# Fill in the remaining missing using median sd
none_free_sds = np.array(sds)[sds != np.array(None)]
# Fill in the remaining missing using median std
none_free_stds = np.array(stds)[stds != np.array(None)]

if none_free_sds:
sd_median = np.median(none_free_sds)
if none_free_stds:
std_median = np.median(none_free_stds)

else: # If no SDs or SEMs reported at all, raise error
else: # If no stds or SEMs reported at all, raise error

# Perhaps the median SD of all cells for this property could
# Perhaps the median std of all cells for this property could
# be used. However, NE API nes interface does not support summary
# prop values without specifying the neuron id
msg = 'No StDevs or SEMs reported for "%s" property "%s"'
msg = msg % (self.neuron_name, self.ephysprop_name)
raise NotImplementedError(msg)

for i, _ in enumerate(sds):
if sds[i] is None:
sds[i] = sd_median
sems[i] = sd_median / np.sqrt(ns[i])
for i, _ in enumerate(stds):
if stds[i] is None:
stds[i] = std_median
sems[i] = std_median / np.sqrt(ns[i])
Loading

0 comments on commit 3223517

Please sign in to comment.