code: def group_reproducibility_metrics(
      mask_images,
      contrast_images,
      variance_images,
      thresholds,
      ngroups,
      method,
      cluster_threshold=10,
      number_of_samples=10,
      sigma=6.0,
      do_clusters=True,
      do_voxels=True,
      do_peaks=True,
      swap=False,
  ):
      """
      Main function to perform reproducibility analysis, including nifti1 io
  
      Parameters
      ----------
      threshold: list or 1-d array,
                 the thresholds to be tested
  
      Returns
      -------
      cluster_rep_results: dictionary,
                           results of cluster-level reproducibility analysis
      voxel_rep_results: dictionary,
                         results of voxel-level reproducibility analysis
      peak_rep_results: dictionary,
                        results of peak-level reproducibility analysis
      """
      from nibabel import load
  
      from ..mask import intersect_masks
  
      if len(variance_images) == 0 and method != "crfx":
          raise ValueError("Variance images are necessary")
  
      nsubj = len(contrast_images)
  
      # compute the group mask
      affine = get_affine(load(mask_images[0]))
      mask = intersect_masks(mask_images, threshold=0) > 0
      domain = grid_domain_from_binary_array(mask, affine)
  
      # read the data
      group_con = []
      group_var = []
      for s in range(nsubj):
          group_con.append(load(contrast_images[s]).get_fdata()[mask])
          if len(variance_images) > 0:
              group_var.append(
                  load(variance_images[s]).get_fdata()[mask]
              )
  
      group_con = np.squeeze(np.array(group_con)).T
      group_con[np.isnan(group_con)] = 0
      if len(variance_images) > 0:
          group_var = np.squeeze(np.array(group_var)).T
          group_var[np.isnan(group_var)] = 0
          group_var = np.maximum(group_var, 1.0e-15)
  
      # perform the analysis
      voxel_rep_results = {}
      cluster_rep_results = {}
      peak_rep_results = {}
  
      for ng in ngroups:
          if do_voxels:
              voxel_rep_results.update({ng: {}})
          if do_clusters:
              cluster_rep_results.update({ng: {}})
          if do_peaks:
              peak_rep_results.update({ng: {}})
          for th in thresholds:
              kappa = []
              cls = []
              pk = []
              kwargs = {"threshold": th, "csize": cluster_threshold}
  
              for i in range(number_of_samples):
                  if do_voxels:
                      kappa.append(
                          voxel_reproducibility(
                              group_con,
                              group_var,
                              domain,
                              ng,
                              method,
                              swap,
                              **kwargs,
                          )
                      )
                  if do_clusters:
                      cls.append(
                          cluster_reproducibility(
                              group_con,
                              group_var,
                              domain,
                              ng,
                              sigma,
                              method,
                              swap,
                              **kwargs,
                          )
                      )
                  if do_peaks:
                      pk.append(
                          peak_reproducibility(
                              group_con,
                              group_var,
                              domain,
                              ng,
                              sigma,
                              method,
                              swap,
                              **kwargs,
                          )
                      )
  
              if do_voxels:
                  voxel_rep_results[ng].update({th: np.array(kappa)})
              if do_clusters:
                  cluster_rep_results[ng].update({th: np.array(cls)})
              if do_peaks:
                  peak_rep_results[ng].update({th: np.array(cls)})
  
      return voxel_rep_results, cluster_rep_results, peak_rep_results
  
imports: import numpy
  from nipy.io.nibcompat import get_affine
  from nipy.labs.spatial_models.discrete_domain import grid_domain_from_binary_array
label_desc: Use image node_id 18 for calling a global function with example usage: group_reproducibility_metrics(mask_images, contrast_images, variance_images, thresholds, ngroups, method, cluster_threshold, number_of_samples, sigma, do_clusters, do_voxels, do_peaks, swap) and returns: voxel_rep_results,  cluster_rep_results,  peak_rep_results