docs: def fit_iSet(Y, U_R=None, S_R=None, covs=None, Xr=None, n_perms=0, Ie=None,
               strat=False, verbose=True):
      """
      Args:
          Y:          [N, P] phenotype matrix
          S_R:        N vector of eigenvalues of R
          U_R:        [N, N] eigenvector matrix of R
          covs:       [N, K] matrix for K covariates
          Xr:         [N, S] genotype data of the set component
          n_perms:    number of permutations to consider
          Ie:         N boolean context indicator
          strat:      if True, the implementation with stratified designs is considered
      """
      factr=1e7 # remove?
      if strat:
          assert Ie is not None, 'Ie must be specified for stratification analyses'
          assert Y.shape[1]==1, 'Y must be Nx1 for stratification analysis'
      else:
          assert covs==None, 'Covariates are not supported for analysis of fully observed phenotypes'
  
      if verbose:     print('fittng iSet')
      if strat:
          mtSetGxE = ISet_Strat(Y, Ie, Xr, covs=covs)
          RV = {}
          RV['null'] = mtSetGxE.fitNull()
          RV['rank2'] = mtSetGxE.fitFullRank()
          RV['rank1'] = mtSetGxE.fitLowRank()
          RV['block'] = mtSetGxE.fitBlock()
          RV['var'] = mtSetGxE.getVC()
      else:
          mtSetGxE = ISet_Full(Y=Y, S_R=S_R, U_R=U_R, Xr=Xr, factr=factr)
          RV = {}
          RV['null'] = mtSetGxE.fitNull()
          RV['rank2'] = mtSetGxE.fitFullRank()
          RV['rank1'] = mtSetGxE.fitLowRank()
          LLR = RV['rank1']['NLLAlt'] - RV['rank2']['NLLAlt']
          if LLR<-1e-6:
              RV['rank2'] = mtSetGxE.fitFullRank(init_method='lr')
          try:
              RV['block'] = mtSetGxE.fitBlock()
          except:
              try:
                  RV['block'] = mtSetGxE.fitBlock(init_method='null')
              except:
                  RV['block'] = mtSetGxE.fitBlock(init_method='null_no_opt')
          RV['var'] = mtSetGxE.getVC()
  
      if n_perms>0:
          RVperm = {}
          nulls = ['null', 'block', 'rank1']
          tests = ['mtSet', 'iSet', 'iSet-het']
          for test in tests:
              RVperm[test+' LLR0'] = sp.zeros(n_perms)
          for seed_i in range(n_perms):
              if verbose:     print('permutation %d / %d' % (seed_i, n_perms))
              for it, test in enumerate(tests):
                  if test=='mtSet':
                      idxs = sp.random.permutation(Xr.shape[0])
                      _Xr = Xr[idxs, :]
                      df0 = fit_iSet(Y, U_R=U_R, S_R=S_R, covs=covs, Xr=_Xr, n_perms=0, Ie=Ie, strat=strat, verbose=False)
                  else:
                      Y0 = mtSetGxE._sim_from(set_covar=nulls[it])
                      Y0 -= Y0.mean(0)
                      df0 = fit_iSet(Y0, U_R=U_R, S_R=S_R, covs=covs, Xr=Xr, n_perms=0, Ie=Ie, strat=strat, verbose=False)
                  RVperm[test+' LLR0'][seed_i]  = df0[test+' LLR'][0]
  
      # output
      LLR_mtSet = RV['null']['NLL']-RV['rank2']['NLL']
      LLR_iSet = RV['block']['NLL']-RV['rank2']['NLL']
      LLR_iSet_het = RV['rank1']['NLL']-RV['rank2']['NLL']
  
      if strat:   var_keys = ['var_r_full', 'var_c', 'var_n']
      else:       var_keys = ['var_r_full', 'var_g', 'var_n']
      varT = sp.sum([RV['var'][key] for key in var_keys])
      var_pers = RV['var']['var_r_block'] / varT
      var_resc = (RV['var']['var_r_rank1'] - RV['var']['var_r_block']) / varT
      var_het = (RV['var']['var_r_full'] - RV['var']['var_r_rank1']) / varT
      conv = RV['null']['conv']
      conv*= RV['block']['conv']
      conv*= RV['rank1']['conv']
      conv*= RV['rank2']['conv']
  
      M = sp.array([LLR_mtSet, LLR_iSet, LLR_iSet_het, var_pers, var_resc, var_het, conv]).T
      columns = ['mtSet LLR', 'iSet LLR', 'iSet-het LLR',
                  'Persistent Var', 'Rescaling-GxC Var', 'Heterogeneity-GxC var', 'Converged']
      df = pd.DataFrame(M, columns=columns)
  
      if n_perms>0:
          return df, pd.DataFrame(RVperm)
      return df