code: def linear_regression(
      X,
      y,
      add_intercept=True,
      weights=None,
      coef_only=False,
      alpha=0.05,
      as_dataframe=True,
      remove_na=False,
      relimp=False,
  ):
      """(Multiple) Linear regression.
  
      Parameters
      ----------
      X : array_like
          Predictor(s), of shape *(n_samples, n_features)* or *(n_samples)*.
      y : array_like
          Dependent variable, of shape *(n_samples)*.
      add_intercept : bool
          If False, assume that the data are already centered. If True, add a
          constant term to the model. In this case, the first value in the
          output dict is the intercept of the model.
  
          .. note:: It is generally recommended to include a constant term
              (intercept) to the model to limit the bias and force the residual
              mean to equal zero. The intercept coefficient and p-values
              are however rarely meaningful.
      weights : array_like
          An optional vector of sample weights to be used in the fitting
          process, of shape *(n_samples)*. Missing or negative weights are not
          allowed. If not null, a weighted least squares is calculated.
  
          .. versionadded:: 0.3.5
      coef_only : bool
          If True, return only the regression coefficients.
      alpha : float
          Alpha value used for the confidence intervals.
          :math:`\\text{CI} = [\\alpha / 2 ; 1 - \\alpha / 2]`
      as_dataframe : bool
          If True, returns a pandas DataFrame. If False, returns a dictionnary.
      remove_na : bool
          If True, apply a listwise deletion of missing values (i.e. the entire
          row is removed). Default is False, which will raise an error if missing
          values are present in either the predictor(s) or dependent
          variable.
      relimp : bool
          If True, returns the relative importance (= contribution) of
          predictors. This is irrelevant when the predictors are uncorrelated:
          the total :math:`R^2` of the model is simply the sum of each univariate
          regression :math:`R^2`-values. However, this does not apply when
          predictors are correlated. Instead, the total :math:`R^2` of the model
          is partitioned by averaging over all combinations of predictors,
          as done in the `relaimpo
          <https://cran.r-project.org/web/packages/relaimpo/relaimpo.pdf>`_
          R package (``calc.relimp(type="lmg")``).
  
          .. warning:: The computation time roughly doubles for each
              additional predictor and therefore this can be extremely slow for
              models with more than 12-15 predictors.
  
          .. versionadded:: 0.3.0
  
      Returns
      -------
      stats : :py:class:`pandas.DataFrame` or dict
          Linear regression summary:
  
          * ``'names'``: name of variable(s) in the model (e.g. x1, x2...)
          * ``'coef'``: regression coefficients
          * ``'se'``: standard errors
          * ``'T'``: T-values
          * ``'pval'``: p-values
          * ``'r2'``: coefficient of determination (:math:`R^2`)
          * ``'adj_r2'``: adjusted :math:`R^2`
          * ``'CI[2.5%]'``: lower confidence intervals
          * ``'CI[97.5%]'``: upper confidence intervals
          * ``'relimp'``: relative contribution of each predictor to the final\
                          :math:`R^2` (only if ``relimp=True``).
          * ``'relimp_perc'``: percent relative contribution
  
          In addition, the output dataframe comes with hidden attributes such as
          the residuals, and degrees of freedom of the model and residuals, which
          can be accessed as follow, respectively:
  
          >>> lm = pg.linear_regression() # doctest: +SKIP
          >>> lm.residuals_, lm.df_model_, lm.df_resid_ # doctest: +SKIP
  
          Note that to follow scikit-learn convention, these hidden atributes end
          with an "_". When ``as_dataframe=False`` however, these attributes
          are no longer hidden and can be accessed as any other keys in the
          output dictionary.
  
          >>> lm = pg.linear_regression() # doctest: +SKIP
          >>> lm['residuals'], lm['df_model'], lm['df_resid'] # doctest: +SKIP
  
          When ``as_dataframe=False`` the dictionary also contains the
          processed ``X`` and ``y`` arrays (i.e, with NaNs removed if
          ``remove_na=True``) and the model's predicted values ``pred``.
  
          >>> lm['X'], lm['y'], lm['pred'] # doctest: +SKIP
  
          For a weighted least squares fit, the weighted ``Xw`` and ``yw``
          arrays are included in the dictionary.
  
          >>> lm['Xw'], lm['yw'] # doctest: +SKIP
  
      See also
      --------
      logistic_regression, mediation_analysis, corr
  
      Notes
      -----
      The :math:`\\beta` coefficients are estimated using an ordinary least
      squares (OLS) regression, as implemented in the
      :py:func:`scipy.linalg.lstsq` function. The OLS method minimizes
      the sum of squared residuals, and leads to a closed-form expression for
      the estimated :math:`\\beta`:
  
      .. math:: \\hat{\\beta} = (X^TX)^{-1} X^Ty
  
      It is generally recommended to include a constant term (intercept) to the
      model to limit the bias and force the residual mean to equal zero.
      Note that intercept coefficient and p-values are however rarely meaningful.
  
      The standard error of the estimates is a measure of the accuracy of the
      prediction defined as:
  
      .. math:: \\sigma = \\sqrt{\\text{MSE} \\cdot (X^TX)^{-1}}
  
      where :math:`\\text{MSE}` is the mean squared error,
  
      .. math::
  
          \\text{MSE} = \\frac{SS_{\\text{resid}}}{n - p - 1}
           = \\frac{\\sum{(\\text{true} - \\text{pred})^2}}{n - p - 1}
  
      :math:`p` is the total number of predictor variables in the model
      (excluding the intercept) and :math:`n` is the sample size.
  
      Using the :math:`\\beta` coefficients and the standard errors,
      the T-values can be obtained:
  
      .. math:: T = \\frac{\\beta}{\\sigma}
  
      and the p-values approximated using a T-distribution with
      :math:`n - p - 1` degrees of freedom.
  
      The coefficient of determination (:math:`R^2`) is defined as:
  
      .. math:: R^2 = 1 - (\\frac{SS_{\\text{resid}}}{SS_{\\text{total}}})
  
      The adjusted :math:`R^2` is defined as:
  
      .. math:: \\overline{R}^2 = 1 - (1 - R^2) \\frac{n - 1}{n - p - 1}
  
      The relative importance (``relimp``) column is a partitioning of the
      total :math:`R^2` of the model into individual :math:`R^2` contribution.
      This is calculated by taking the average over average contributions in
      models of different sizes. For more details, please refer to
      `Groemping et al. 2006 <http://dx.doi.org/10.18637/jss.v017.i01>`_
      and the R package `relaimpo
      <https://cran.r-project.org/web/packages/relaimpo/relaimpo.pdf>`_.
  
      Note that Pingouin will automatically remove any duplicate columns
      from :math:`X`, as well as any column with only one unique value
      (constant), excluding the intercept.
  
      Results have been compared against sklearn, R, statsmodels and JASP.
  
      Examples
      --------
      1. Simple linear regression using columns of a pandas dataframe
  
      In this first example, we'll use the tips dataset to see how well we
      can predict the waiter's tip (in dollars) based on the total bill (also
      in dollars).
  
      >>> import numpy as np
      >>> import pingouin as pg
      >>> df = pg.read_dataset('tips')
      >>> # Let's predict the tip ($) based on the total bill (also in $)
      >>> lm = pg.linear_regression(df['total_bill'], df['tip'])
      >>> lm.round(2)
              names  coef    se      T  pval    r2  adj_r2  CI[2.5%]  CI[97.5%]
      0   Intercept  0.92  0.16   5.76   0.0  0.46    0.45      0.61       1.23
      1  total_bill  0.11  0.01  14.26   0.0  0.46    0.45      0.09       0.12
  
      It comes as no surprise that total bill is indeed a significant predictor
      of the waiter's tip (T=14.26, p<0.05). The :math:`R^2` of the model is 0.46
      and the adjusted :math:`R^2` is 0.45, which means that our model roughly
      explains ~45% of the total variance in the tip amount.
  
      2. Multiple linear regression
  
      We can also have more than one predictor and run a multiple linear
      regression. Below, we add the party size as a second predictor of tip.
  
      >>> # We'll add a second predictor: the party size
      >>> lm = pg.linear_regression(df[['total_bill', 'size']], df['tip'])
      >>> lm.round(2)
              names  coef    se      T  pval    r2  adj_r2  CI[2.5%]  CI[97.5%]
      0   Intercept  0.67  0.19   3.46  0.00  0.47    0.46      0.29       1.05
      1  total_bill  0.09  0.01  10.17  0.00  0.47    0.46      0.07       0.11
      2        size  0.19  0.09   2.26  0.02  0.47    0.46      0.02       0.36
  
      The party size is also a significant predictor of tip (T=2.26, p=0.02).
      Note that adding this new predictor however only improved the :math:`R^2`
      of our model by ~1%.
  
      This function also works with numpy arrays:
  
      >>> X = df[['total_bill', 'size']].to_numpy()
      >>> y = df['tip'].to_numpy()
      >>> pg.linear_regression(X, y).round(2)
             names  coef    se      T  pval    r2  adj_r2  CI[2.5%]  CI[97.5%]
      0  Intercept  0.67  0.19   3.46  0.00  0.47    0.46      0.29       1.05
      1         x1  0.09  0.01  10.17  0.00  0.47    0.46      0.07       0.11
      2         x2  0.19  0.09   2.26  0.02  0.47    0.46      0.02       0.36
  
      3. Get the residuals
  
      >>> # For clarity, only display the first 9 values
      >>> np.round(lm.residuals_, 2)[:9]
      array([-1.62, -0.55,  0.31,  0.06, -0.11,  0.93,  0.13, -0.81, -0.49])
  
      Using pandas, we can show a summary of the distribution of the residuals:
  
      >>> import pandas as pd
      >>> pd.Series(lm.residuals_).describe().round(2)
      count    244.00
      mean      -0.00
      std        1.01
      min       -2.93
      25%       -0.55
      50%       -0.09
      75%        0.51
      max        4.04
      dtype: float64
  
      5. No intercept and return only the regression coefficients
  
      Sometimes it may be useful to remove the constant term from the regression,
      or to only return the regression coefficients without calculating the
      standard errors or p-values. This latter can potentially save you a lot of
      time if you need to calculate hundreds of regression and only care about
      the coefficients!
  
      >>> pg.linear_regression(X, y, add_intercept=False, coef_only=True)
      array([0.1007119 , 0.36209717])
  
      6. Return a dictionnary instead of a dataframe
  
      >>> lm_dict = pg.linear_regression(X, y, as_dataframe=False)
      >>> lm_dict.keys()
      dict_keys(['names', 'coef', 'se', 'T', 'pval', 'r2', 'adj_r2', 'CI[2.5%]',
                 'CI[97.5%]', 'df_model', 'df_resid', 'residuals', 'X', 'y',
                 'pred'])
  
      7. Remove missing values
  
      >>> X[4, 1] = np.nan
      >>> y[7] = np.nan
      >>> pg.linear_regression(X, y, remove_na=True, coef_only=True)
      array([0.65749955, 0.09262059, 0.19927529])
  
      8. Get the relative importance of predictors
  
      >>> lm = pg.linear_regression(X, y, remove_na=True, relimp=True)
      >>> lm[['names', 'relimp', 'relimp_perc']]
             names    relimp  relimp_perc
      0  Intercept       NaN          NaN
      1         x1  0.342503    73.045583
      2         x2  0.126386    26.954417
  
      The ``relimp`` column is a partitioning of the total :math:`R^2` of the
      model into individual contribution. Therefore, it sums to the :math:`R^2`
      of the full model. The ``relimp_perc`` is normalized to sum to 100%. See
      `Groemping 2006 <https://www.jstatsoft.org/article/view/v017i01>`_
      for more details.
  
      >>> lm[['relimp', 'relimp_perc']].sum()
      relimp           0.468889
      relimp_perc    100.000000
      dtype: float64
  
      9. Weighted linear regression
  
      >>> X = [1, 2, 3, 4, 5, 6]
      >>> y = [10, 22, 11, 13, 13, 16]
      >>> w = [1, 0.1, 1, 1, 0.5, 1]  # Array of weights. Must be >= 0.
      >>> lm = pg.linear_regression(X, y, weights=w)
      >>> lm.round(2)
             names  coef    se     T  pval    r2  adj_r2  CI[2.5%]  CI[97.5%]
      0  Intercept  9.00  2.03  4.42  0.01  0.51    0.39      3.35      14.64
      1         x1  1.04  0.50  2.06  0.11  0.51    0.39     -0.36       2.44
      """
      # Extract names if X is a Dataframe or Series
      if isinstance(X, pd.DataFrame):
          names = X.keys().tolist()
      elif isinstance(X, pd.Series):
          names = [X.name]
      else:
          names = []
  
      # Convert input to numpy array
      X = np.asarray(X)
      y = np.asarray(y)
      assert y.ndim == 1, "y must be one-dimensional."
      assert 0 < alpha < 1
  
      if X.ndim == 1:
          # Convert to (n_samples, n_features) shape
          X = X[..., np.newaxis]
  
      # Check for NaN / Inf
      if remove_na:
          X, y = rm_na(X, y[..., np.newaxis], paired=True, axis="rows")
          y = np.squeeze(y)
      y_gd = np.isfinite(y).all()
      X_gd = np.isfinite(X).all()
      assert y_gd, (
          "Target (y) contains NaN or Inf. Please remove them "
          "manually or use remove_na=True."
      )
      assert X_gd, (
          "Predictors (X) contain NaN or Inf. Please remove them "
          "manually or use remove_na=True."
      )
  
      # Check that X and y have same length
      assert (
          y.shape[0] == X.shape[0]
      ), "X and y must have same number of samples"
  
      if not names:
          names = ["x" + str(i + 1) for i in range(X.shape[1])]
  
      if add_intercept:
          # Add intercept
          X = np.column_stack((np.ones(X.shape[0]), X))
          names.insert(0, "Intercept")
  
      # FINAL CHECKS BEFORE RUNNING LEAST SQUARES REGRESSION
      # 1. Let's remove column(s) with only zero, otherwise the regression fails
      n_nonzero = np.count_nonzero(X, axis=0)
      idx_zero = np.flatnonzero(
          n_nonzero == 0
      )  # Find columns that are only 0
      if len(idx_zero):
          X = np.delete(X, idx_zero, 1)
          names = np.delete(names, idx_zero)
  
      # 2. We also want to make sure that there is no more than one constant
      # column (= intercept), otherwise the regression fails
      # This is equivalent, but much faster, to pd.DataFrame(X).nunique()
      idx_unique = np.where(np.all(X == X[0, :], axis=0))[0]
      if len(idx_unique) > 1:
          # We remove all but the first "Intercept" column.
          X = np.delete(X, idx_unique[1:], 1)
          names = np.delete(names, idx_unique[1:])
      # Is there a constant in our predictor matrix? Useful for dof and R^2.
      constant = 1 if len(idx_unique) > 0 else 0
  
      # 3. Finally, we want to remove duplicate columns
      if X.shape[1] > 1:
          idx_duplicate = []
          for pair in itertools.combinations(range(X.shape[1]), 2):
              if np.array_equal(X[:, pair[0]], X[:, pair[1]]):
                  idx_duplicate.append(pair[1])
          if len(idx_duplicate):
              X = np.delete(X, idx_duplicate, 1)
              names = np.delete(names, idx_duplicate)
  
      # 4. Check that we have enough samples / features
      n, p = X.shape[0], X.shape[1]
      assert n >= 3, "At least three valid samples are required in X."
      assert p >= 1, "X must have at least one valid column."
  
      # 5. Handle weights
      if weights is not None:
          if relimp:
              raise ValueError(
                  "relimp = True is not supported when using "
                  "weights."
              )
          w = np.asarray(weights)
          assert w.ndim == 1, "weights must be a 1D array."
          assert w.size == n, "weights must be of shape n_samples."
          assert not np.isnan(
              w
          ).any(), "Missing weights are not accepted."
          assert not (w < 0).any(), "Negative weights are not accepted."
          # Do not count weights == 0 in dof
          # This gives similar results as R lm() but different from statsmodels
          n = np.count_nonzero(w)
          # Rescale (whitening)
          wts = np.diag(np.sqrt(w))
          Xw = wts @ X
          yw = wts @ y
      else:
          # Set all weights to one, [1, 1, 1, ...]
          w = np.ones(n)
          Xw = X
          yw = y
  
      # FIT (WEIGHTED) LEAST SQUARES REGRESSION
      coef, ss_res, rank, _ = lstsq(Xw, yw, cond=None)
      ss_res = ss_res[0] if ss_res.shape == (1,) else ss_res
      if coef_only:
          return coef
      calc_ss_res = False
      if rank < Xw.shape[1]:
          # in this case, ss_res is of shape (0,), i.e., an empty array
          warnings.warn(
              "Design matrix supplied with `X` parameter is rank "
              f"deficient (rank {rank} with {Xw.shape[1]} columns). "
              "That means that one or more of the columns in `X` "
              "are a linear combination of one of more of the "
              "other columns."
          )
          calc_ss_res = True
  
      # Degrees of freedom
      df_model = rank - constant
      df_resid = n - rank
  
      # Calculate predicted values and (weighted) residuals
      pred = Xw @ coef
      resid = yw - pred
      if calc_ss_res:
          # In case we did not get ss_res from lstsq due to rank deficiency
          ss_res = (resid**2).sum()
  
      # Calculate total (weighted) sums of squares and R^2
      ss_tot = yw @ yw
      ss_wtot = np.sum(w * (y - np.average(y, weights=w)) ** 2)
      if constant:
          r2 = 1 - ss_res / ss_wtot
      else:
          r2 = 1 - ss_res / ss_tot
      adj_r2 = 1 - (1 - r2) * (n - constant) / df_resid
  
      # Compute mean squared error, variance and SE
      mse = ss_res / df_resid
      beta_var = mse * (np.linalg.pinv(Xw.T @ Xw).diagonal())
      beta_se = np.sqrt(beta_var)
  
      # Compute T and p-values
      T = coef / beta_se
      pval = 2 * t.sf(np.fabs(T), df_resid)
  
      # Compute confidence intervals
      crit = t.ppf(1 - alpha / 2, df_resid)
      marg_error = crit * beta_se
      ll = coef - marg_error
      ul = coef + marg_error
  
      # Rename CI
      ll_name = "CI[%.1f%%]" % (100 * alpha / 2)
      ul_name = "CI[%.1f%%]" % (100 * (1 - alpha / 2))
  
      # Create dict
      stats = {
          "names": names,
          "coef": coef,
          "se": beta_se,
          "T": T,
          "pval": pval,
          "r2": r2,
          "adj_r2": adj_r2,
          ll_name: ll,
          ul_name: ul,
      }
  
      # Relative importance
      if relimp:
          data = pd.concat(
              [
                  pd.DataFrame(y, columns=["y"]),
                  pd.DataFrame(X, columns=names),
              ],
              sort=False,
              axis=1,
          )
          if "Intercept" in names:
              # Intercept is the first column
              reli = _relimp(
                  data.drop(columns=["Intercept"]).cov(
                      numeric_only=True
                  )
              )
              reli["names"] = ["Intercept"] + reli["names"]
              reli["relimp"] = np.insert(reli["relimp"], 0, np.nan)
              reli["relimp_perc"] = np.insert(
                  reli["relimp_perc"], 0, np.nan
              )
          else:
              reli = _relimp(data.cov(numeric_only=True))
          stats.update(reli)
  
      if as_dataframe:
          stats = _postprocess_dataframe(pd.DataFrame(stats))
          stats.df_model_ = df_model
          stats.df_resid_ = df_resid
          stats.residuals_ = 0  # Trick to avoid Pandas warning
          stats.residuals_ = resid  # Residuals is a hidden attribute
      else:
          stats["df_model"] = df_model
          stats["df_resid"] = df_resid
          stats["residuals"] = resid
          stats["X"] = X
          stats["y"] = y
          stats["pred"] = pred
          if weights is not None:
              stats["yw"] = yw
              stats["Xw"] = Xw
      return stats
  
file_path: pingouin/pingouin/regression.py
imports: import itertools
  import warnings
  import numpy
  import pandas
  import pandas_flavor
  from scipy.stats import t, norm
  from scipy.linalg import pinvh, lstsq
  from pingouin.config import options
  from pingouin.utils import remove_na
  from pingouin.utils import _flatten_list
  from pingouin.utils import _postprocess_dataframe
label_desc: Use image node_id 1 for calling a global function with example usage: linear_regression(X, y, add_intercept, weights, coef_only, alpha, as_dataframe, remove_na, relimp) and returns: stats, coef