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")

def _compute_scaled_schoenfeld(self, X, T, E, weights, index=None):
        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")


        """

        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