Consider the following multiple linear regression model \begin{equation} y=\beta_0 + \beta_1 x_1 + ... + \beta_k x_k +u \end{equation} If the variance of the error is not constant for any value of the regressors, we have heteroskedasticity and Gauss-Markov assumption of homoskedasticity is violeted, i.e., $Var(u_i|x_{i1},...,x_{ik})= \sigma^2_i$.
A valid estimator of $Var(\hat{\beta}_j)$, under Gauss-Markov with homoskedasticity not necessarily holding, is \begin{equation} \widehat{Var}(\hat{\beta}_j) = \frac{ \sum_{i=1}^{n} \hat{r}^2_{ij} \hat{u}^2_i}{RSS^2_j} = \frac{ \sum_{i=1}^{n} \hat{r}^2_{ij} \hat{u}^2_i}{[TSS_j (1-R^2_j)]^2} \label{varbeta3} \end{equation} where $\hat{r}_{ij}$ denotes the $i$ residual from regressing $x_j$ on all the other independent variables, and $RSS_j$ is the residual sum of squares from this regression. Also, $TSS_j$ is the total sample variation in $x_j$ and $R^2_j$ is the R-squared from regressing $x_j$ on all the other independent variables(including the intercept).
I tried to compute the heteroskedasticity robust st.error in R studio, by appling the following:
\begin{equation} se(\hat{\beta}_j)= \frac{ \sqrt{\frac{\sum_{i=1}^{n} \hat{u}_i^2 \hat{r}^2_{ij} }{n-k-1}} }{ n^{-1} \sum_{i=1}^{n} (x_{ij} - \bar{x}_j)^2 (1-R^2_j)} \end{equation}
However, I got a wrong heteroskedasticity robust st.error. What am I doing wrong?
