code: def _compute_scaled_schoenfeld(self, X, T, E, weights, index=None):
          r
  
          n_deaths = self.event_observed.sum()
          scaled_schoenfeld_resids = n_deaths * self._compute_schoenfeld(X, T, E, weights, index).dot(
              self.variance_matrix_
          )
  
          scaled_schoenfeld_resids.columns = self.hazards_.index
          return scaled_schoenfeld_resids
docs: r"""
          Let s_k be the kth schoenfeld residuals. Then E[s_k] = 0.
          For tests of proportionality, we want to test if \beta_i(t) is \beta_i (constant) or not.
  
          Let V_k be the contribution to the information matrix at time t_k. A main result from Grambsch and Therneau is that
  
          \beta(t) = E[s_k*V_k^{-1} + \hat{beta}]
  
          so define s_k^* = s_k*V_k^{-1} + \hat{beta} as the scaled schoenfeld residuals.
  
          We can approximate V_k with Hessian/d, so the inverse of Hessian/d is (d * variance_matrix_)
  
          Notes
          -------
          lifelines does not add the coefficients to the final results, but R does when you call residuals(c, "scaledsch")