From fab62fa1e8a8825688ab9a644bf36f2898761557 Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 12:11:00 -0500 Subject: [PATCH 01/12] permuted Mantel tests will report all permuted scores --- pyssage/mantel.py | 74 ++++++++++++++++++++++++++++---------------- tests/test_mantel.py | 15 ++++----- 2 files changed, 55 insertions(+), 34 deletions(-) diff --git a/pyssage/mantel.py b/pyssage/mantel.py index b2102ff..bfc12d1 100644 --- a/pyssage/mantel.py +++ b/pyssage/mantel.py @@ -21,7 +21,7 @@ def check_tail(tail: str) -> None: def mantel(input_matrix1: numpy.ndarray, input_matrix2: numpy.ndarray, partial, permutations: int = 0, - tail: str = "both") -> Tuple[float, float, list, float, float, float, float]: + tail: str = "both") -> Tuple[float, float, list, float, float, float, list, float]: check_tail(tail) n = check_for_square_matrix(input_matrix1) if n != check_for_square_matrix(input_matrix2): @@ -81,6 +81,7 @@ def mantel(input_matrix1: numpy.ndarray, input_matrix2: numpy.ndarray, partial, p_value = 1 - p_value # perform permutation tests + permuted_r_list = [r] if permutations > 0: cumulative_left = 0 cumulative_right = 0 @@ -94,31 +95,49 @@ def mantel(input_matrix1: numpy.ndarray, input_matrix2: numpy.ndarray, partial, if len(partial) > 0: # if partial, calculate residuals for permuted matrix matrix1 = residuals_from_matrix_regression(matrix1, partial) - # If it is a two-tailed test, we need to calculate r, otherwise for one-tailed tests we can stick with - # Z which is faster - if tail == "both": - numerator = square_matrix_covariance(matrix1, matrix2) - if len(partial) > 0: - denominator = sqrt(square_matrix_covariance(matrix1, matrix1) * sq_cov2) - else: # for non-partial tests can save computation as denominator is fixed - denominator = sqxy - permuted_r = numerator / denominator - if permuted_r < r: - cumulative_left += 1 - elif permuted_r > r: - cumulative_right += 1 - else: - cumulative_equal += 1 - if abs(permuted_r) >= abs(r): - cumulative_total += 1 + + # if we want to save the permuted values, we cannot use the quick way for one-tailed tests + numerator = square_matrix_covariance(matrix1, matrix2) + if len(partial) > 0: + denominator = sqrt(square_matrix_covariance(matrix1, matrix1) * sq_cov2) + else: # for non-partial tests can save computation as denominator is fixed + denominator = sqxy + permuted_r = numerator / denominator + if permuted_r < r: + cumulative_left += 1 + elif permuted_r > r: + cumulative_right += 1 else: - permuted_z = numpy.sum(numpy.multiply(matrix1, matrix2)) - if permuted_z < observed_z: - cumulative_left += 1 - elif permuted_z > observed_z: - cumulative_right += 1 - else: - cumulative_equal += 1 + cumulative_equal += 1 + if abs(permuted_r) >= abs(r): + cumulative_total += 1 + permuted_r_list.append(permuted_r) + + # # If it is a two-tailed test, we need to calculate r, otherwise for one-tailed tests we can stick with + # # Z which is faster + # if tail == "both": + # numerator = square_matrix_covariance(matrix1, matrix2) + # if len(partial) > 0: + # denominator = sqrt(square_matrix_covariance(matrix1, matrix1) * sq_cov2) + # else: # for non-partial tests can save computation as denominator is fixed + # denominator = sqxy + # permuted_r = numerator / denominator + # if permuted_r < r: + # cumulative_left += 1 + # elif permuted_r > r: + # cumulative_right += 1 + # else: + # cumulative_equal += 1 + # if abs(permuted_r) >= abs(r): + # cumulative_total += 1 + # else: + # permuted_z = numpy.sum(numpy.multiply(matrix1, matrix2)) + # if permuted_z < observed_z: + # cumulative_left += 1 + # elif permuted_z > observed_z: + # cumulative_right += 1 + # else: + # cumulative_equal += 1 permuted_right_p = (cumulative_equal + cumulative_right) / permutations permuted_left_p = (cumulative_equal + cumulative_left) / permutations if tail == "both": @@ -141,8 +160,9 @@ def mantel(input_matrix1: numpy.ndarray, input_matrix2: numpy.ndarray, partial, permuted_left_p, permuted_right_p, permuted_two_p = 1, 1, 1 mantel_output = namedtuple("mantel_output", ["r", "p_value", "output_text", "permuted_left_p", "permuted_right_p", - "permuted_two_p", "z_score"]) - return mantel_output(r, p_value, output_text, permuted_left_p, permuted_right_p, permuted_two_p, z_score) + "permuted_two_p", "permuted_r_list", "z_score"]) + return mantel_output(r, p_value, output_text, permuted_left_p, permuted_right_p, permuted_two_p, permuted_r_list, + z_score) def mantel_moments(x: numpy.ndarray, y: numpy.ndarray) -> Tuple[float, float]: diff --git a/tests/test_mantel.py b/tests/test_mantel.py index 26b89a4..77495fa 100644 --- a/tests/test_mantel.py +++ b/tests/test_mantel.py @@ -24,7 +24,7 @@ def test_mantel(): coords = create_test_coords() distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) - r, p_value, output_text, _, _, _, _ = pyssage.mantel.mantel(distances, angles, []) + r, p_value, output_text, _, _, _, _, _ = pyssage.mantel.mantel(distances, angles, []) for line in output_text: print(line) assert round(r, 5) == 0.29830 @@ -40,7 +40,7 @@ def test_mantel_with_permutation(): coords = create_test_coords() distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) - r, p_value, output_text, _, _, _, _ = pyssage.mantel.mantel(distances, angles, [], permutations=100) + r, p_value, output_text, _, _, _, _, _ = pyssage.mantel.mantel(distances, angles, [], permutations=100) for line in output_text: print(line) @@ -69,7 +69,7 @@ def test_mantel_with_partial_single(): angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) - r, p_value, output_text, _, _, _, _ = pyssage.mantel.mantel(data_distances, distances, [angles]) + r, p_value, output_text, _, _, _, _, _ = pyssage.mantel.mantel(data_distances, distances, [angles]) for line in output_text: print(line) assert round(r, 5) == 0.26699 @@ -87,7 +87,7 @@ def test_mantel_with_partial_single_permuted(): angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) - r, p_value, output_text, _, _, _, _ = pyssage.mantel.mantel(data_distances, distances, [angles], permutations=100) + r, p_value, output_text, _, _, _, _, _ = pyssage.mantel.mantel(data_distances, distances, [angles], permutations=100) for line in output_text: print(line) @@ -116,7 +116,8 @@ def test_mantel_with_partial_multi(): data_distances1 = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) data_distances2 = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_manhattan) - r, p_value, output_text, _, _, _, _ = pyssage.mantel.mantel(data_distances1, data_distances2, [distances, angles]) + r, p_value, output_text, _, _, _, _, _ = pyssage.mantel.mantel(data_distances1, data_distances2, + [distances, angles]) for line in output_text: print(line) assert round(r, 5) == 0.89973 @@ -135,8 +136,8 @@ def test_mantel_with_partial_multi_permuted(): data_distances1 = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) data_distances2 = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_manhattan) - r, p_value, output_text, _, _, _, _ = pyssage.mantel.mantel(data_distances1, data_distances2, [distances, angles], - permutations=100) + r, p_value, output_text, _, _, _, _, _ = pyssage.mantel.mantel(data_distances1, data_distances2, + [distances, angles], permutations=100) for line in output_text: print(line) From adc29f4bbb85a5c0aabae7a5573fce769d56ffdf Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 12:11:44 -0500 Subject: [PATCH 02/12] tweak to create_output_table parameter type --- pyssage/anisotropy.py | 10 +++++----- pyssage/utils.py | 3 +-- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/pyssage/anisotropy.py b/pyssage/anisotropy.py index 8620a2e..356adea 100644 --- a/pyssage/anisotropy.py +++ b/pyssage/anisotropy.py @@ -34,7 +34,7 @@ def bearing_analysis(data: numpy.ndarray, distances: numpy.ndarray, angles: nump for a in range(nbearings): test_angle = a * angle_width b_matrix = distances * numpy.square(numpy.cos(angles - test_angle)) - r, p_value, _, _, _, rand_p, _ = pyssage.mantel.mantel(data, b_matrix, [], npermutations) + r, p_value, _, _, _, rand_p,_, _ = pyssage.mantel.mantel(data, b_matrix, [], npermutations) if npermutations > 0: output.append([a*180/nbearings, r, p_value, rand_p]) else: @@ -47,11 +47,11 @@ def bearing_analysis(data: numpy.ndarray, distances: numpy.ndarray, angles: nump output_text.append("Tested {} vectors".format(nbearings)) output_text.append("") if npermutations > 0: - col_headers = ("Bearing", "Correlation", "Prob", "RandProb") - col_formats = ("f", "f", "f", "f") + col_headers = ["Bearing", "Correlation", "Prob", "RandProb"] + col_formats = ["f", "f", "f", "f"] else: - col_headers = ("Bearing", "Correlation", "Prob") - col_formats = ("f", "f", "f") + col_headers = ["Bearing", "Correlation", "Prob"] + col_formats = ["f", "f", "f"] create_output_table(output_text, output, col_headers, col_formats) bearing_output = namedtuple("bearing_output", ["output_values", "output_text"]) return bearing_output(output, output_text) diff --git a/pyssage/utils.py b/pyssage/utils.py index 6685c89..ad7960d 100644 --- a/pyssage/utils.py +++ b/pyssage/utils.py @@ -56,7 +56,6 @@ def deflatten_without_diagonal(vector: numpy.ndarray, n: int) -> numpy.ndarray: return numpy.array(output) - # def euc_dist(x1: float, x2: float, y1: float = 0, y2: float = 0, z1: float = 0, z2: float = 0) -> float: # return sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2) @@ -100,7 +99,7 @@ def check_for_square_matrix(test_matrix: numpy.ndarray) -> int: return len(test_matrix) -def create_output_table(output_text: list, table_data: list, col_headers: tuple, col_formats: tuple, +def create_output_table(output_text: list, table_data: list, col_headers: list, col_formats: list, out_dec: int = OUT_DEC, sbc: int = 3) -> None: """ Create a well-formatted set of strings representing an output table, including headers and computationally From caedc4f62ccb40fa60f66b203a83bb1e8af9f4dc Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 12:12:01 -0500 Subject: [PATCH 03/12] permutation tests added to correlograms bearing and windrose still need to be tested/tweaked --- pyssage/correlogram.py | 124 +++++++++++++++++++++++++++----------- tests/test_correlogram.py | 43 +++++++++++++ 2 files changed, 132 insertions(+), 35 deletions(-) diff --git a/pyssage/correlogram.py b/pyssage/correlogram.py index 92a7414..dee8a13 100644 --- a/pyssage/correlogram.py +++ b/pyssage/correlogram.py @@ -6,6 +6,8 @@ from pyssage.utils import create_output_table, check_for_square_matrix import pyssage.mantel +__all__ = ["bearing_correlogram", "correlogram", "windrose_correlogram"] + def check_variance_assumption(x: Optional[str]) -> None: """ @@ -21,21 +23,35 @@ def check_variance_assumption(x: Optional[str]) -> None: def morans_i(y: numpy.ndarray, weights: Connections, alt_weights: Optional[numpy.ndarray] = None, - variance: Optional[str] = "random"): + variance: Optional[str] = "random", permutations: int = 0): check_variance_assumption(variance) n = len(y) - # mean_y = numpy.average(y) mean_y = numpy.mean(y) dev_y = y - mean_y # deviations from mean w = weights.as_binary() if alt_weights is not None: # multiply to create non-binary weights, if necessary w = w * alt_weights - sumyij = numpy.sum(numpy.outer(dev_y, dev_y) * w, dtype=numpy.float64) sumy2 = numpy.sum(numpy.square(dev_y), dtype=numpy.float64) # sum of squared deviations from mean sumw = numpy.sum(w, dtype=numpy.float64) # sum of weight matrix sumw2 = sumw**2 - moran = n * sumyij / (sumw * sumy2) expected = -1 / (n - 1) + sumyij = numpy.sum(numpy.outer(dev_y, dev_y) * w, dtype=numpy.float64) + moran = n * sumyij / (sumw * sumy2) + + # permutations + permuted_i_list = [moran] + perm_p = 1 + if permutations > 0: + rand_y = numpy.copy(dev_y) + for k in range(permutations-1): + numpy.random.shuffle(rand_y) + perm_sumyij = numpy.sum(numpy.outer(rand_y, rand_y) * w, dtype=numpy.float64) + perm_moran = n * perm_sumyij / (sumw * sumy2) + permuted_i_list.append(perm_moran) + if abs(perm_moran-expected) >= abs(moran-expected): + perm_p += 1 + perm_p /= permutations + if variance is None: sd, z, p = None, None, None else: @@ -51,24 +67,38 @@ def morans_i(y: numpy.ndarray, weights: Connections, alt_weights: Optional[numpy z = abs(moran - expected) / sd p = scipy.stats.norm.sf(z)*2 # two-tailed test - return weights.min_scale, weights.max_scale, weights.n_pairs(), expected, moran, sd, z, p + return weights.min_scale, weights.max_scale, weights.n_pairs(), expected, moran, sd, z, p, perm_p, permuted_i_list def gearys_c(y: numpy.ndarray, weights: Connections, alt_weights: Optional[numpy.ndarray] = None, - variance: Optional[str] = "random"): + variance: Optional[str] = "random", permutations: int = 0): check_variance_assumption(variance) n = len(y) - # mean_y = numpy.average(y) mean_y = numpy.mean(y) dev_y = y - mean_y # deviations from mean w = weights.as_binary() if alt_weights is not None: # multiply to create non-binary weights, if necessary w *= alt_weights - sumdif2 = numpy.sum(numpy.square(w * (dev_y[:, numpy.newaxis] - dev_y)), dtype=numpy.float64) sumy2 = numpy.sum(numpy.square(dev_y), dtype=numpy.float64) # sum of squared deviations from mean sumw = numpy.sum(w, dtype=numpy.float64) # sum of weight matrix sumw2 = sumw**2 + sumdif2 = numpy.sum(numpy.square(w * (dev_y[:, numpy.newaxis] - dev_y)), dtype=numpy.float64) geary = (n - 1) * sumdif2 / (2 * sumw * sumy2) + + # permutations + permuted_c_list = [geary] + perm_p = 1 + if permutations > 0: + rand_y = numpy.copy(dev_y) + for k in range(permutations-1): + numpy.random.shuffle(rand_y) + perm_sumdif2 = numpy.sum(numpy.square(w * (rand_y[:, numpy.newaxis] - rand_y)), dtype=numpy.float64) + perm_geary = (n - 1) * perm_sumdif2 / (2 * sumw * sumy2) + permuted_c_list.append(perm_geary) + if abs(perm_geary - 1) >= abs(geary - 1): + perm_p += 1 + perm_p /= permutations + if variance is None: sd, z, p = None, None, None else: @@ -86,11 +116,11 @@ def gearys_c(y: numpy.ndarray, weights: Connections, alt_weights: Optional[numpy z = abs(geary - 1) / sd p = scipy.stats.norm.sf(z)*2 # two-tailed test - return weights.min_scale, weights.max_scale, weights.n_pairs(), 1, geary, sd, z, p + return weights.min_scale, weights.max_scale, weights.n_pairs(), 1, geary, sd, z, p, perm_p, permuted_c_list def mantel_correl(y: numpy.ndarray, weights: Connections, alt_weights: Optional[numpy.ndarray] = None, - variance: Optional[str] = "random"): + variance=None, permutations: int = 0): """ in order to get the bearing version to work right, we have to use normal binary weights, then reverse the sign of the resulting Mantel correlation. if we use reverse binary weighting we end up multiplying the 'out of @@ -99,12 +129,12 @@ def mantel_correl(y: numpy.ndarray, weights: Connections, alt_weights: Optional[ w = weights.as_binary() if alt_weights is not None: # multiply to create non-binary weights, if necessary w *= alt_weights - r, p_value, _, _, _, _, z = pyssage.mantel.mantel(y, w, []) - return weights.min_scale, weights.max_scale, weights.n_pairs(), 0, -r, -z, p_value + r, p_value, _, _, _, permuted_two_p, permuted_rs, z = pyssage.mantel.mantel(y, w, [], permutations=permutations) + return weights.min_scale, weights.max_scale, weights.n_pairs(), 0, -r, -z, p_value, permuted_two_p, permuted_rs def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: morans_i, - variance: Optional[str] = "random"): + variance: Optional[str] = "random", permutations: int = 0): if metric == morans_i: metric_title = "Moran's I" exp_format = "f" @@ -119,7 +149,11 @@ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: moran exp_format = "" output = [] for dc in dist_class_connections: - output.append(metric(data, dc, variance=variance)) + *tmp_out, permuted_values = metric(data, dc, variance=variance, permutations=permutations) + if permutations > 0: + output.append(tmp_out) + else: + output.append(tmp_out[:len(tmp_out)-1]) # create basic output text output_text = list() @@ -130,18 +164,22 @@ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: moran output_text.append("Distribution assumption = {}".format(variance)) output_text.append("") if metric == mantel_correl: - col_headers = ("Min dist", "Max dist", "# pairs", "Expected", metric_title, "Z", "Prob") - col_formats = ("f", "f", "d", exp_format, "f", "f", "f") + col_headers = ["Min dist", "Max dist", "# pairs", "Expected", metric_title, "Z", "Prob"] + col_formats = ["f", "f", "d", exp_format, "f", "f", "f"] else: - col_headers = ("Min dist", "Max dist", "# pairs", "Expected", metric_title, "SD", "Z", "Prob") - col_formats = ("f", "f", "d", exp_format, "f", "f", "f", "f") + col_headers = ["Min dist", "Max dist", "# pairs", "Expected", metric_title, "SD", "Z", "Prob"] + col_formats = ["f", "f", "d", exp_format, "f", "f", "f", "f"] + if permutations > 0: + col_headers.append("PermProb") + col_formats.append("f") + create_output_table(output_text, output, col_headers, col_formats) return output, output_text def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angles: numpy.ndarray, n_bearings: int = 18, - metric=morans_i, variance: Optional[str] = "random"): + metric=morans_i, variance: Optional[str] = "random", permutations: int = 0): if metric == morans_i: metric_title = "Moran's I" exp_format = "f" @@ -166,9 +204,13 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle output = [] for i, b in enumerate(bearing_weights): for dc in dist_class_connections: - tmp_out = list(metric(data, dc, alt_weights=b, variance=variance)) + *tmp_out, permuted_values = metric(data, dc, alt_weights=b, variance=variance) + tmp_out = list(tmp_out) tmp_out.insert(2, degrees(bearings[i])) - output.append(tmp_out) + if permutations > 0: + output.append(tmp_out) + else: + output.append(tmp_out[:len(tmp_out)-1]) # create basic output text output_text = list() @@ -179,11 +221,14 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle output_text.append("Distribution assumption = {}".format(variance)) output_text.append("") if metric == mantel_correl: - col_headers = ("Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "Z", "Prob") - col_formats = ("f", "f", "f", "d", exp_format, "f", "f", "f") + col_headers = ["Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "Z", "Prob"] + col_formats = ["f", "f", "f", "d", exp_format, "f", "f", "f"] else: - col_headers = ("Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "SD", "Z", "Prob") - col_formats = ("f", "f", "f", "d", exp_format, "f", "f", "f", "f") + col_headers = ["Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "SD", "Z", "Prob"] + col_formats = ["f", "f", "f", "d", exp_format, "f", "f", "f", "f"] + if permutations > 0: + col_headers.append("PermProb") + col_formats.append("f") create_output_table(output_text, output, col_headers, col_formats) return output, output_text @@ -211,7 +256,8 @@ def create_windrose_connections(distances: numpy.ndarray, angles: numpy.ndarray, def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: numpy.ndarray, radius_c: float, radius_d: float, radius_e: float, segment_param: int = 4, - min_pairs: int = 21, metric=morans_i, variance: Optional[str] = "random"): + min_pairs: int = 21, metric=morans_i, variance: Optional[str] = "random", + permutations: int = 0): if metric == morans_i: metric_title = "Moran's I" exp_format = "f" @@ -243,12 +289,17 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: segment_param, radius_c, radius_d, radius_e) np = connection.n_pairs() if np >= min_pairs: - tmp_out = list(metric(data, connection, variance=variance)) + *tmp_out, permuted_values = metric(data, connection, variance=variance) + tmp_out = list(tmp_out) # add sector angles to output tmp_out.insert(2, degrees(min_ang)) tmp_out.insert(3, degrees(max_ang)) - output.append(tmp_out) - all_output.append(tmp_out) + if permutations > 0: + output.append(tmp_out) + all_output.append(tmp_out) + else: + output.append(tmp_out[:len(tmp_out) - 1]) + all_output.append(tmp_out[:len(tmp_out) - 1]) else: # using -1 for the probability as an indicator that nothing was calculated if metric == mantel_correl: @@ -273,13 +324,16 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: output_text.append("") if metric == mantel_correl: - col_headers = ("Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, - "Z", "Prob") - col_formats = ("f", "f", "f", "f", "d", exp_format, "f", "f", "f") + col_headers = ["Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, + "Z", "Prob"] + col_formats = ["f", "f", "f", "f", "d", exp_format, "f", "f", "f"] else: - col_headers = ("Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, "SD", - "Z", "Prob") - col_formats = ("f", "f", "f", "f", "d", exp_format, "f", "f", "f", "f") + col_headers = ["Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, "SD", + "Z", "Prob"] + col_formats = ["f", "f", "f", "f", "d", exp_format, "f", "f", "f", "f"] + if permutations > 0: + col_headers.append("PermProb") + col_formats.append("f") create_output_table(output_text, output, col_headers, col_formats) return output, output_text, all_output diff --git a/tests/test_correlogram.py b/tests/test_correlogram.py index 118c055..d889c45 100644 --- a/tests/test_correlogram.py +++ b/tests/test_correlogram.py @@ -75,6 +75,20 @@ def test_morans_i(): assert round(output[i][j], 5) == ans +def test_morans_i_perm(): + data, _ = create_test_scattered() + coords = create_test_coords() + distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) + dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) + dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) + output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.morans_i, + permutations=100) + pyssage.graph.draw_correlogram(numpy.array(output), "Moran's I", "Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) + for line in output_text: + print(line) + + def test_gearys_c(): # answer calculated from PASSaGE 2 and exported to 5 decimals # tests both random and normal variance assumption @@ -133,6 +147,20 @@ def test_gearys_c(): assert round(output[i][j], 5) == ans +def test_gearys_c_perm(): + data, _ = create_test_scattered() + coords = create_test_coords() + distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) + dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) + dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) + output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.gearys_c, + permutations=100) + pyssage.graph.draw_correlogram(numpy.array(output), "Geary's c", "Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) + for line in output_text: + print(line) + + def test_mantel_correl(): # The Passage 2 Mantel correlogram code appears to be buggy # Some manual tests do seem to indicate this is working correctly @@ -150,6 +178,21 @@ def test_mantel_correl(): print(line) +def test_mantel_correl_perm(): + data, _ = create_test_scattered() + coords = create_test_coords() + distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) + dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) + dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) + data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) + output, output_text = pyssage.correlogram.correlogram(data_distances, dc_con, pyssage.correlogram.mantel_correl, + variance=None, permutations=100) + pyssage.graph.draw_correlogram(numpy.array(output), "Mantel r", "Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) + for line in output_text: + print(line) + + def test_bearing_correlogram(): # answer calculated from PASSaGE 2 and exported to 5 decimals, but only matches perfectly to 2 decimals answer = [(0, 2.44706, 0, 4187, -0.00282, 0.48074, 0.01802, 26.83519, 0), From bcecaa3857e3a678b8eb51e94774d93047b9a384 Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 12:21:57 -0500 Subject: [PATCH 04/12] permutations for windrose and bearing correlograms --- pyssage/correlogram.py | 6 ++- tests/test_correlogram.py | 90 ++++++++++++++++++++++++++++++++++----- 2 files changed, 83 insertions(+), 13 deletions(-) diff --git a/pyssage/correlogram.py b/pyssage/correlogram.py index dee8a13..56358bb 100644 --- a/pyssage/correlogram.py +++ b/pyssage/correlogram.py @@ -204,7 +204,7 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle output = [] for i, b in enumerate(bearing_weights): for dc in dist_class_connections: - *tmp_out, permuted_values = metric(data, dc, alt_weights=b, variance=variance) + *tmp_out, permuted_values = metric(data, dc, alt_weights=b, variance=variance, permutations=permutations) tmp_out = list(tmp_out) tmp_out.insert(2, degrees(bearings[i])) if permutations > 0: @@ -289,7 +289,7 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: segment_param, radius_c, radius_d, radius_e) np = connection.n_pairs() if np >= min_pairs: - *tmp_out, permuted_values = metric(data, connection, variance=variance) + *tmp_out, permuted_values = metric(data, connection, variance=variance, permutations=permutations) tmp_out = list(tmp_out) # add sector angles to output tmp_out.insert(2, degrees(min_ang)) @@ -308,6 +308,8 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: else: tmp_out = [connection.min_scale, connection.max_scale, degrees(min_ang), degrees(max_ang), np, 0, 0, 0, 0, -1] + if permutations > 0: + tmp_out.append(-1) all_output.append(tmp_out) # create basic output text diff --git a/tests/test_correlogram.py b/tests/test_correlogram.py index d889c45..9d782dc 100644 --- a/tests/test_correlogram.py +++ b/tests/test_correlogram.py @@ -484,6 +484,59 @@ def test_bearing_correlogram(): assert round(output[i][j], 2) == round(ans, 2) +def test_bearing_correlogram_perm(): + data, _ = create_test_scattered() + coords = create_test_coords() + distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) + angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) + dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) + dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) + output, output_text = pyssage.correlogram.bearing_correlogram(data[:, 22], dc_con, angles, permutations=100) + # pyssage.graph.draw_bearing_correlogram_old(numpy.array(output), "Moran's I Bearing Correlogram") + pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Moran's I Bearing Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) + for line in output_text: + print(line) + + +def test_bearing_correlogram_mantel(): + # The Passage 2 Mantel correlogram code appears to be buggy, which effects the Mantel bearing correlogram as well + # There was also potentially a logic error in using "reverse binary weights" with non-binary weighting which + # also would have led to problematic results, even if the first error hadn't overridden everything + data, _ = create_test_scattered() + coords = create_test_coords() + distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) + angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) + data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) + dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) + dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) + + output, output_text = pyssage.correlogram.bearing_correlogram(data_distances, dc_con, angles, + metric=pyssage.correlogram.mantel_correl) + pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Mantel Bearing Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) + for line in output_text: + print(line) + + +def test_bearing_correlogram_mantel_perm(): + data, _ = create_test_scattered() + coords = create_test_coords() + distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) + angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) + data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) + dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) + dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) + + output, output_text = pyssage.correlogram.bearing_correlogram(data_distances, dc_con, angles, + metric=pyssage.correlogram.mantel_correl, + permutations=100) + pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Mantel Bearing Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) + for line in output_text: + print(line) + + def test_windrose_correlogram(): data, _ = create_test_scattered() coords = create_test_coords() @@ -508,7 +561,23 @@ def test_windrose_correlogram(): # pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Moran's I Windrose Correlogram") -def test_mantel_windrose_correlogram(): +def test_windrose_correlogram_perm(): + data, _ = create_test_scattered() + coords = create_test_coords() + distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) + angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) + output, output_text, all_output = pyssage.correlogram.windrose_correlogram(data[:, 22], distances, angles, + radius_c=3, radius_d=0, radius_e=0, + permutations=100) + pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Moran's I Windrose Correlogram Pair Counts", + show_counts=True, figoutput=pyssage.graph.FigOutput(figshow=True)) + pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Moran's I Windrose Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) + for line in output_text: + print(line) + + +def test_windrose_correlogram_mantel(): data, _ = create_test_scattered() coords = create_test_coords() distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) @@ -526,21 +595,20 @@ def test_mantel_windrose_correlogram(): print(line) -def test_bearing_correlogram_mantel(): - # The Passage 2 Mantel correlogram code appears to be buggy, which effects the Mantel bearing correlogram as well - # There was also potentially a logic error in using "reverse binary weights" with non-binary weighting which - # also would have led to problematic results, even if the first error hadn't overridden everything +def test_windrose_correlogram_mantel_perm(): data, _ = create_test_scattered() coords = create_test_coords() distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) - dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) - dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.bearing_correlogram(data_distances, dc_con, angles, - metric=pyssage.correlogram.mantel_correl) - pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Mantel Bearing Correlogram", - figoutput=pyssage.graph.FigOutput(figshow=True)) + output, output_text, all_output = pyssage.correlogram.windrose_correlogram(data_distances, distances, angles, + radius_c=3, radius_d=0, radius_e=0, + metric=pyssage.correlogram.mantel_correl, + permutations=100) + pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Mantel Windrose Correlogram Pair Counts", + show_counts=True, figoutput=pyssage.graph.FigOutput(figshow=True)) + pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Mantel Windrose Correlogram", + figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: print(line) From d2b14d1bd46b82503914804c623d9465d235e837 Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 13:01:09 -0500 Subject: [PATCH 05/12] fixing documentation --- pyssage/utils.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pyssage/utils.py b/pyssage/utils.py index ad7960d..fd90228 100644 --- a/pyssage/utils.py +++ b/pyssage/utils.py @@ -56,10 +56,6 @@ def deflatten_without_diagonal(vector: numpy.ndarray, n: int) -> numpy.ndarray: return numpy.array(output) -# def euc_dist(x1: float, x2: float, y1: float = 0, y2: float = 0, z1: float = 0, z2: float = 0) -> float: -# return sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2) - - def euclidean_angle(x1: float, y1: float, x2: float, y2: float, do360: bool = False) -> float: """ This will calculate the angle between the two points @@ -110,8 +106,8 @@ def create_output_table(output_text: list, table_data: list, col_headers: list, already have been instantiated :param table_data: a list of lists containing the data to appear in the table; each sublist represents a row of the table and must contain the same number of columns - :param col_headers: a tuple containing strings representing headers for each column in the table - :param col_formats: a tuple containing basic string formatting codes, generally expected to be "f" of "d" + :param col_headers: a list containing strings representing headers for each column in the table + :param col_formats: a list containing basic string formatting codes, generally expected to be "f" of "d" :param out_dec: the number of decimal places to output floating point numbers in the table :param sbc: the number of spaces to use between each column in the table """ From 3c793e144e904e1900678f5766a0bf63491cb803 Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 13:01:58 -0500 Subject: [PATCH 06/12] output tweak --- pyssage/mantel.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyssage/mantel.py b/pyssage/mantel.py index bfc12d1..d60a2c1 100644 --- a/pyssage/mantel.py +++ b/pyssage/mantel.py @@ -65,6 +65,8 @@ def mantel(input_matrix1: numpy.ndarray, input_matrix2: numpy.ndarray, partial, # OutputAddLine(outstr); # end; output_text.append("Matrices are {0} x {0}".format(n)) + if len(partial) > 0: + output_text.append("{} matrices held constant".format(len(partial))) output_text.append("") output_text.append("Observed Z = " + format(observed_z, OUT_FRMT)) output_text.append("Correlation = " + format(r, OUT_FRMT)) @@ -138,6 +140,7 @@ def mantel(input_matrix1: numpy.ndarray, input_matrix2: numpy.ndarray, partial, # cumulative_right += 1 # else: # cumulative_equal += 1 + permuted_right_p = (cumulative_equal + cumulative_right) / permutations permuted_left_p = (cumulative_equal + cumulative_left) / permutations if tail == "both": From 5430b649329522102851a77e5b326c018c436afd Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 13:07:18 -0500 Subject: [PATCH 07/12] added permuted correlogram values to function output --- pyssage/correlogram.py | 19 +++++++---- tests/test_correlogram.py | 68 ++++++++++++++++++++------------------- 2 files changed, 48 insertions(+), 39 deletions(-) diff --git a/pyssage/correlogram.py b/pyssage/correlogram.py index 56358bb..6b0a72b 100644 --- a/pyssage/correlogram.py +++ b/pyssage/correlogram.py @@ -134,7 +134,7 @@ def mantel_correl(y: numpy.ndarray, weights: Connections, alt_weights: Optional[ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: morans_i, - variance: Optional[str] = "random", permutations: int = 0): + variance: Optional[str] = "random", permutations: int = 0) -> Tuple[list, list, list]: if metric == morans_i: metric_title = "Moran's I" exp_format = "f" @@ -148,10 +148,12 @@ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: moran metric_title = "" exp_format = "" output = [] + all_permuted_values = [] for dc in dist_class_connections: *tmp_out, permuted_values = metric(data, dc, variance=variance, permutations=permutations) if permutations > 0: output.append(tmp_out) + all_permuted_values.append(permuted_values) else: output.append(tmp_out[:len(tmp_out)-1]) @@ -175,11 +177,12 @@ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: moran create_output_table(output_text, output, col_headers, col_formats) - return output, output_text + return output, output_text, all_permuted_values def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angles: numpy.ndarray, n_bearings: int = 18, - metric=morans_i, variance: Optional[str] = "random", permutations: int = 0): + metric=morans_i, variance: Optional[str] = "random", + permutations: int = 0) -> Tuple[list, list, list]: if metric == morans_i: metric_title = "Moran's I" exp_format = "f" @@ -202,6 +205,7 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle bearing_weights.append(numpy.square(numpy.cos(angles - a))) output = [] + all_permuted_values = [] for i, b in enumerate(bearing_weights): for dc in dist_class_connections: *tmp_out, permuted_values = metric(data, dc, alt_weights=b, variance=variance, permutations=permutations) @@ -209,6 +213,7 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle tmp_out.insert(2, degrees(bearings[i])) if permutations > 0: output.append(tmp_out) + all_permuted_values.append(permuted_values) else: output.append(tmp_out[:len(tmp_out)-1]) @@ -231,7 +236,7 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle col_formats.append("f") create_output_table(output_text, output, col_headers, col_formats) - return output, output_text + return output, output_text, all_permuted_values def windrose_sectors_per_annulus(segment_param: int, annulus: int) -> int: @@ -257,7 +262,7 @@ def create_windrose_connections(distances: numpy.ndarray, angles: numpy.ndarray, def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: numpy.ndarray, radius_c: float, radius_d: float, radius_e: float, segment_param: int = 4, min_pairs: int = 21, metric=morans_i, variance: Optional[str] = "random", - permutations: int = 0): + permutations: int = 0) -> Tuple[list, list, list, list]: if metric == morans_i: metric_title = "Moran's I" exp_format = "f" @@ -283,6 +288,7 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: all_output = [] # all_output is needed for graphing the output *if* we want to include those sectors with too few pairs, but # still more than zero + all_permuted_values = [] for annulus in range(n_annuli): for sector in range(windrose_sectors_per_annulus(segment_param, annulus)): connection, min_ang, max_ang = create_windrose_connections(distances, angles, annulus, sector, @@ -297,6 +303,7 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: if permutations > 0: output.append(tmp_out) all_output.append(tmp_out) + all_permuted_values.append(permuted_values) else: output.append(tmp_out[:len(tmp_out) - 1]) all_output.append(tmp_out[:len(tmp_out) - 1]) @@ -338,4 +345,4 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: col_formats.append("f") create_output_table(output_text, output, col_headers, col_formats) - return output, output_text, all_output + return output, output_text, all_output, all_permuted_values diff --git a/tests/test_correlogram.py b/tests/test_correlogram.py index 9d782dc..7bc5620 100644 --- a/tests/test_correlogram.py +++ b/tests/test_correlogram.py @@ -40,7 +40,7 @@ def test_morans_i(): distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.morans_i) + output, output_text, _ = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.morans_i) pyssage.graph.draw_correlogram(numpy.array(output), "Moran's I", "Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -65,8 +65,8 @@ def test_morans_i(): (14.96656, 16.90406, 4189, -0.00282, -0.1481, 0.01496, 9.71083, 0), (16.90406, 20.06888, 4189, -0.00282, 0.00519, 0.01468, 0.54596, 0.58509)] - output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.morans_i, - variance="normal") + output, output_text, _ = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.morans_i, + variance="normal") for line in output_text: print(line) @@ -81,8 +81,8 @@ def test_morans_i_perm(): distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.morans_i, - permutations=100) + output, output_text, _ = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.morans_i, + permutations=100) pyssage.graph.draw_correlogram(numpy.array(output), "Moran's I", "Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -112,7 +112,7 @@ def test_gearys_c(): distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.gearys_c) + output, output_text, _ = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.gearys_c) pyssage.graph.draw_correlogram(numpy.array(output), "Geary's c", "Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -137,8 +137,8 @@ def test_gearys_c(): (14.96656, 16.90406, 4189, 1, 1.00765, 0.03815, 0.20057, 0.84104), (16.90406, 20.06888, 4189, 1, 0.68403, 0.0538, 5.87292, 0)] - output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.gearys_c, - variance="normal") + output, output_text, _ = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.gearys_c, + variance="normal") for line in output_text: print(line) @@ -153,8 +153,8 @@ def test_gearys_c_perm(): distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.gearys_c, - permutations=100) + output, output_text, _ = pyssage.correlogram.correlogram(data[:, 0], dc_con, pyssage.correlogram.gearys_c, + permutations=100) pyssage.graph.draw_correlogram(numpy.array(output), "Geary's c", "Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -170,8 +170,8 @@ def test_mantel_correl(): dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) - output, output_text = pyssage.correlogram.correlogram(data_distances, dc_con, pyssage.correlogram.mantel_correl, - variance=None) + output, output_text, _ = pyssage.correlogram.correlogram(data_distances, dc_con, pyssage.correlogram.mantel_correl, + variance=None) pyssage.graph.draw_correlogram(numpy.array(output), "Mantel r", "Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -185,8 +185,8 @@ def test_mantel_correl_perm(): dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) - output, output_text = pyssage.correlogram.correlogram(data_distances, dc_con, pyssage.correlogram.mantel_correl, - variance=None, permutations=100) + output, output_text, _ = pyssage.correlogram.correlogram(data_distances, dc_con, pyssage.correlogram.mantel_correl, + variance=None, permutations=100) pyssage.graph.draw_correlogram(numpy.array(output), "Mantel r", "Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -472,7 +472,7 @@ def test_bearing_correlogram(): angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.bearing_correlogram(data[:, 22], dc_con, angles) + output, output_text, _ = pyssage.correlogram.bearing_correlogram(data[:, 22], dc_con, angles) # pyssage.graph.draw_bearing_correlogram_old(numpy.array(output), "Moran's I Bearing Correlogram") pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Moran's I Bearing Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) @@ -491,7 +491,7 @@ def test_bearing_correlogram_perm(): angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.bearing_correlogram(data[:, 22], dc_con, angles, permutations=100) + output, output_text, _ = pyssage.correlogram.bearing_correlogram(data[:, 22], dc_con, angles, permutations=100) # pyssage.graph.draw_bearing_correlogram_old(numpy.array(output), "Moran's I Bearing Correlogram") pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Moran's I Bearing Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) @@ -511,8 +511,8 @@ def test_bearing_correlogram_mantel(): dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.bearing_correlogram(data_distances, dc_con, angles, - metric=pyssage.correlogram.mantel_correl) + output, output_text, _ = pyssage.correlogram.bearing_correlogram(data_distances, dc_con, angles, + metric=pyssage.correlogram.mantel_correl) pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Mantel Bearing Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -528,9 +528,9 @@ def test_bearing_correlogram_mantel_perm(): dist_classes = pyssage.distances.create_distance_classes(distances, "determine pair count", 15) dc_con = pyssage.connections.distance_classes_to_connections(dist_classes, distances) - output, output_text = pyssage.correlogram.bearing_correlogram(data_distances, dc_con, angles, - metric=pyssage.correlogram.mantel_correl, - permutations=100) + output, output_text, _ = pyssage.correlogram.bearing_correlogram(data_distances, dc_con, angles, + metric=pyssage.correlogram.mantel_correl, + permutations=100) pyssage.graph.draw_bearing_correlogram(numpy.array(output), "Mantel Bearing Correlogram", figoutput=pyssage.graph.FigOutput(figshow=True)) for line in output_text: @@ -542,8 +542,8 @@ def test_windrose_correlogram(): coords = create_test_coords() distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) - output, output_text, all_output = pyssage.correlogram.windrose_correlogram(data[:, 22], distances, angles, - radius_c=3, radius_d=0, radius_e=0) + output, output_text, all_output, _ = pyssage.correlogram.windrose_correlogram(data[:, 22], distances, angles, + radius_c=3, radius_d=0, radius_e=0) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Moran's I Windrose Correlogram Pair Counts", show_counts=True, figoutput=pyssage.graph.FigOutput(figshow=True)) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Moran's I Windrose Correlogram", @@ -566,9 +566,9 @@ def test_windrose_correlogram_perm(): coords = create_test_coords() distances = pyssage.distances.euclidean_distance_matrix(coords[:, 0], coords[:, 1]) angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) - output, output_text, all_output = pyssage.correlogram.windrose_correlogram(data[:, 22], distances, angles, - radius_c=3, radius_d=0, radius_e=0, - permutations=100) + output, output_text, all_output, _ = pyssage.correlogram.windrose_correlogram(data[:, 22], distances, angles, + radius_c=3, radius_d=0, radius_e=0, + permutations=100) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Moran's I Windrose Correlogram Pair Counts", show_counts=True, figoutput=pyssage.graph.FigOutput(figshow=True)) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Moran's I Windrose Correlogram", @@ -584,9 +584,10 @@ def test_windrose_correlogram_mantel(): angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) - output, output_text, all_output = pyssage.correlogram.windrose_correlogram(data_distances, distances, angles, - radius_c=3, radius_d=0, radius_e=0, - metric=pyssage.correlogram.mantel_correl) + (output, output_text, + all_output, _) = pyssage.correlogram.windrose_correlogram(data_distances, distances, angles, radius_c=3, + radius_d=0, radius_e=0, + metric=pyssage.correlogram.mantel_correl) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Mantel Windrose Correlogram Pair Counts", show_counts=True, figoutput=pyssage.graph.FigOutput(figshow=True)) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Mantel Windrose Correlogram", @@ -602,10 +603,11 @@ def test_windrose_correlogram_mantel_perm(): angles = pyssage.distances.euclidean_angle_matrix(coords[:, 0], coords[:, 1]) data_distances = pyssage.distances.data_distance_matrix(data, pyssage.distances.data_distance_euclidean) - output, output_text, all_output = pyssage.correlogram.windrose_correlogram(data_distances, distances, angles, - radius_c=3, radius_d=0, radius_e=0, - metric=pyssage.correlogram.mantel_correl, - permutations=100) + (output, output_text, + all_output, _) = pyssage.correlogram.windrose_correlogram(data_distances, distances, angles, radius_c=3, + radius_d=0, radius_e=0, + metric=pyssage.correlogram.mantel_correl, + permutations=100) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Mantel Windrose Correlogram Pair Counts", show_counts=True, figoutput=pyssage.graph.FigOutput(figshow=True)) pyssage.graph.draw_windrose_correlogram(numpy.array(all_output), title="Mantel Windrose Correlogram", From bc72bdc18c04c5e793c508ac351e7f589c485b5b Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 13:14:50 -0500 Subject: [PATCH 08/12] code tweak/simplification --- pyssage/correlogram.py | 36 ++++++++++++++++-------------------- 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/pyssage/correlogram.py b/pyssage/correlogram.py index 6b0a72b..d8bae9e 100644 --- a/pyssage/correlogram.py +++ b/pyssage/correlogram.py @@ -165,12 +165,11 @@ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: moran if variance is not None: output_text.append("Distribution assumption = {}".format(variance)) output_text.append("") - if metric == mantel_correl: - col_headers = ["Min dist", "Max dist", "# pairs", "Expected", metric_title, "Z", "Prob"] - col_formats = ["f", "f", "d", exp_format, "f", "f", "f"] - else: - col_headers = ["Min dist", "Max dist", "# pairs", "Expected", metric_title, "SD", "Z", "Prob"] - col_formats = ["f", "f", "d", exp_format, "f", "f", "f", "f"] + col_headers = ["Min dist", "Max dist", "# pairs", "Expected", metric_title, "Z", "Prob"] + col_formats = ["f", "f", "d", exp_format, "f", "f", "f"] + if metric != mantel_correl: + col_headers.insert(5, "SD") + col_formats.insert(5, "f") if permutations > 0: col_headers.append("PermProb") col_formats.append("f") @@ -225,12 +224,11 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle if variance is not None: output_text.append("Distribution assumption = {}".format(variance)) output_text.append("") - if metric == mantel_correl: - col_headers = ["Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "Z", "Prob"] - col_formats = ["f", "f", "f", "d", exp_format, "f", "f", "f"] - else: - col_headers = ["Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "SD", "Z", "Prob"] - col_formats = ["f", "f", "f", "d", exp_format, "f", "f", "f", "f"] + col_headers = ["Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "Z", "Prob"] + col_formats = ["f", "f", "f", "d", exp_format, "f", "f", "f"] + if metric != mantel_correl: + col_headers.insert(6, "SD") + col_formats.insert(6, "f") if permutations > 0: col_headers.append("PermProb") col_formats.append("f") @@ -332,14 +330,12 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: output_text.append("Distribution assumption = {}".format(variance)) output_text.append("") - if metric == mantel_correl: - col_headers = ["Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, - "Z", "Prob"] - col_formats = ["f", "f", "f", "f", "d", exp_format, "f", "f", "f"] - else: - col_headers = ["Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, "SD", - "Z", "Prob"] - col_formats = ["f", "f", "f", "f", "d", exp_format, "f", "f", "f", "f"] + col_headers = ["Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, + "Z", "Prob"] + col_formats = ["f", "f", "f", "f", "d", exp_format, "f", "f", "f"] + if metric != mantel_correl: + col_headers.insert(7, "SD") + col_formats.insert(7, "f") if permutations > 0: col_headers.append("PermProb") col_formats.append("f") From 816a45dc855c7760972de19034713fca8cae314a Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 13:16:55 -0500 Subject: [PATCH 09/12] added permutation info to textual output of correlograms --- pyssage/correlogram.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pyssage/correlogram.py b/pyssage/correlogram.py index d8bae9e..e8b689d 100644 --- a/pyssage/correlogram.py +++ b/pyssage/correlogram.py @@ -164,6 +164,8 @@ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: moran output_text.append("# of data points = {}".format(len(data))) if variance is not None: output_text.append("Distribution assumption = {}".format(variance)) + if permutations > 0: + output_text.append("Permutation probability calculated from {} permutations".format(permutations)) output_text.append("") col_headers = ["Min dist", "Max dist", "# pairs", "Expected", metric_title, "Z", "Prob"] col_formats = ["f", "f", "d", exp_format, "f", "f", "f"] @@ -223,6 +225,8 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle output_text.append("# of data points = {}".format(len(data))) if variance is not None: output_text.append("Distribution assumption = {}".format(variance)) + if permutations > 0: + output_text.append("Permutation probability calculated from {} permutations".format(permutations)) output_text.append("") col_headers = ["Min dist", "Max dist", "Bearing", "# pairs", "Expected", metric_title, "Z", "Prob"] col_formats = ["f", "f", "f", "d", exp_format, "f", "f", "f"] @@ -328,6 +332,8 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: output_text.append("") if variance is not None: output_text.append("Distribution assumption = {}".format(variance)) + if permutations > 0: + output_text.append("Permutation probability calculated from {} permutations".format(permutations)) output_text.append("") col_headers = ["Min dist", "Max dist", "Min angle", "Max angle", "# pairs", "Expected", metric_title, From 6bc778d3ea372fdd7fc75be8557456435b2e3f91 Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 14:25:03 -0500 Subject: [PATCH 10/12] Update connections.py added missing function to __all__ --- pyssage/connections.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyssage/connections.py b/pyssage/connections.py index 74236e1..7f159de 100644 --- a/pyssage/connections.py +++ b/pyssage/connections.py @@ -5,8 +5,9 @@ from pyssage.classes import Point, Triangle, VoronoiEdge, VoronoiTessellation, VoronoiPolygon, Connections from pyssage.utils import euclidean_angle, check_for_square_matrix -__all__ = ["delaunay_tessellation", "relative_neighborhood_network", "gabriel_network", - "minimum_spanning_tree", "connect_distance_range", "least_diagonal_network", "nearest_neighbor_connections"] +__all__ = ["connect_distance_range", "delaunay_tessellation", "distance_classes_to_connections", + "gabriel_network", "least_diagonal_network", "minimum_spanning_tree", "nearest_neighbor_connections", + "relative_neighborhood_network"] def create_point_list(x: numpy.ndarray, y: numpy.ndarray) -> list: From 7b80e0477f4415f7e44929393025156314f4f2a2 Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Thu, 18 Feb 2021 14:25:25 -0500 Subject: [PATCH 11/12] converted correlogram outputs to namedtuples --- pyssage/correlogram.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/pyssage/correlogram.py b/pyssage/correlogram.py index e8b689d..201f7f0 100644 --- a/pyssage/correlogram.py +++ b/pyssage/correlogram.py @@ -1,5 +1,6 @@ from math import sqrt, pi, degrees from typing import Optional, Tuple +from collections import namedtuple import numpy import scipy.stats from pyssage.classes import Connections @@ -178,7 +179,8 @@ def correlogram(data: numpy.ndarray, dist_class_connections: list, metric: moran create_output_table(output_text, output, col_headers, col_formats) - return output, output_text, all_permuted_values + correlogram_output = namedtuple("correlogram_output", ["output_values", "output_text", "permuted_values"]) + return correlogram_output(output, output_text, all_permuted_values) def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angles: numpy.ndarray, n_bearings: int = 18, @@ -238,7 +240,8 @@ def bearing_correlogram(data: numpy.ndarray, dist_class_connections: list, angle col_formats.append("f") create_output_table(output_text, output, col_headers, col_formats) - return output, output_text, all_permuted_values + correlogram_output = namedtuple("correlogram_output", ["output_values", "output_text", "permuted_values"]) + return correlogram_output(output, output_text, all_permuted_values) def windrose_sectors_per_annulus(segment_param: int, annulus: int) -> int: @@ -347,4 +350,6 @@ def windrose_correlogram(data: numpy.ndarray, distances: numpy.ndarray, angles: col_formats.append("f") create_output_table(output_text, output, col_headers, col_formats) - return output, output_text, all_output, all_permuted_values + windrose_correlogram_output = namedtuple("windrose_correlogram_output", ["output_values", "output_text", + "all_output", "permuted_values"]) + return windrose_correlogram_output(output, output_text, all_output, all_permuted_values) From 1d2af09357fe3f690748cc898c7b08474249b99d Mon Sep 17 00:00:00 2001 From: msrosenberg Date: Fri, 19 Feb 2021 09:32:51 -0500 Subject: [PATCH 12/12] Update setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 7466d90..dcf5dc2 100644 --- a/setup.py +++ b/setup.py @@ -7,7 +7,7 @@ setuptools.setup( name="pyssage", - version="0.1.1", + version="0.1.2", author="Michael S. Rosenberg", author_email="msrosenberg@vcu.edu", description="python version of some of the analyses from PASSaGE 2",