HomogeneousModel is not producing what is expected?

[damingli09]

We used the lines of code below to generate a connectivity matrix based on Oh method. However, it produces a matrix that is wrong: zero correlation with the new one based on voxel model, which does not make sense.

cache = VoxelModelCache() voxel_data = cache.get_experiment_data(cre=False, injection_structure_ids=iso_list, projection_structure_ids=iso_list, injection_hemisphere_id = 2, projection_hemisphere_id = 2)

iso_list is the list of ids of the 43 isocortex ROIs, right hemisphere only (we tried both hemispheres as well)

regional_data = RegionalData.from_voxel_data(voxel_data) reg = HomogeneousModel() reg.fit(regional_data.injections, regional_data.projections) homogeneous_connectivity = np.array(reg.weights)

[kamdh]

I'd like to echo Daming's concerns here. Is anybody working on this?

[kamdh]

@damingli91 may want to try this code from an unresolved pull request: https://github.com/AllenInstitute/mouse_connectivity_models/pull/25/commits/89092e6dd9e992cb1db7d625852c3b8b57318834

[DrJigsaw]

The issue seems to be the order of the regions (sources and targets). They are sometimes in "ontology order" and sometimes in ascending numerical order. I was careful to separate the two in the code below and I found a correlation coefficient of ~0.6 between the regionalized voxel model and the homogeneous model in isocortex.

from mcmodels.core import VoxelModelCache, RegionalData
from mcmodels.models.voxel import RegionalizedModel
from mcmodels.models import HomogeneousModel
import scipy.stats as st
import numpy as np
import pandas as pd

# set up cache and get list of isocortex structures
cache = VoxelModelCache()
structure_tree = cache.get_structure_tree()
summary_structures = structure_tree.get_structures_by_set_id([687527945])
iso = structure_tree.get_structures_by_acronym(['Isocortex'])[0]['id']
iso_desc = structure_tree.descendant_ids([iso])[0]
summary_structure_ids = [structure['id'] for structure in summary_structures]
iso_list = [item for item in iso_desc if item in summary_structure_ids]
str_iso = [str(strid) for strid in iso_list]

#%% regionalized voxel model method 1
normalized_connection_density = cache.get_normalized_connection_density()
ncd_ipsi = normalized_connection_density['ipsi']
# ontology order
ncd_method1_sorted = ncd_ipsi.loc[iso_list, str_iso]
# ascending order
ncd_method1 = ncd_method1_sorted.sort_index()
ncd_method1 = ncd_method1[[str(item) for item in sorted(iso_list)]]

#%% regionalized voxel model method 2
voxel_array, source_mask, target_mask = cache.get_voxel_connectivity_array()
source_key = source_mask.get_key(structure_ids=iso_list)
target_key = target_mask.get_key(structure_ids=iso_list)
regionalized_model = RegionalizedModel.from_voxel_array(
        voxel_array, source_key, target_key)
ncd_method2 = pd.DataFrame(regionalized_model.normalized_connection_density,
                           columns = regionalized_model.source_regions, 
                           index = regionalized_model.target_regions)
sorter = dict(zip(iso_list, np.arange(len(iso_list))))
# make a copy and sort by ontology order
ncd_method2_sorted = ncd_method2.copy()
ncd_method2_sorted['sort'] = [sorter[structure] for structure in ncd_method2_sorted.index]
ncd_method2_sorted.sort_values(by='sort', inplace = True)
ncd_method2_sorted.drop(columns = {'sort'}, inplace = True)
ncd_method2_sorted = ncd_method2_sorted[iso_list]

#%% homogeneous model  - regions in ascending numerical order
voxel_data = cache.get_experiment_data(cre=False, 
                                       injection_structure_ids=iso_list, 
                                       projection_structure_ids=iso_list,
                                       injection_hemisphere_id = 2, 
                                       projection_hemisphere_id = 2)
regional_data = RegionalData.from_voxel_data(voxel_data)
reg = HomogeneousModel()
reg.fit(regional_data.injections, regional_data.projections)

#%% compare methods 1 and 2 with each other and with homogeneous model
# all in ascending order

spearmanr_2methods, pval_2methods = st.spearmanr(ncd_method1.values.flatten(),
                                                 ncd_method2.values.flatten())
print('correlation methods 1 and 2', spearmanr_2methods)
spearmanr1, pval1 = st.spearmanr(ncd_method1.values.flatten(), 
                                reg.weights.flatten())
print('correlation method 1 and homogeneous', spearmanr1)
spearmanr2, pval2 = st.spearmanr(ncd_method2.values.flatten(),
                                reg.weights.flatten())
print('correlation method 2 and homogeneous', spearmanr2)

[smih]

Based on the answer in the previous comment, the issue looks solved.