1
$\begingroup$

In order to solve heteroscedasticity in my data, I ran a regression with heteroscedasticity-consistent ("robust") standard errors. I would also like to report i) standardized betas together with confidence intervals, ii) partial eta squared values, iii) Tolerance and VIF in one table, similar to the first one in this image: regression table with more data, but additionally with confidence intervals for beta (rather than the unstandardized coefficient).

Here is what I did so far:


> load(file = "Amman.rda")
> model <- lm(progress ~ opi + competence + integration + indegree + voterank, data = Amman)
> summary(model)

Call:
lm(formula = progress ~ opi + competence + integration + indegree + 
    voterank, data = Amman)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.150864 -0.028803  0.003795  0.026640  0.124062 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.496290   0.078140   6.351 2.69e-06 ***
opi         -0.424556   0.153329  -2.769   0.0115 *  
competence   0.010045   0.004462   2.252   0.0352 *  
integration -0.238163   0.099404  -2.396   0.0260 *  
indegree     0.023413   0.013545   1.729   0.0986 .  
voterank    -0.002266   0.003058  -0.741   0.4669    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06614 on 21 degrees of freedom
  (10 observations deleted due to missingness)
Multiple R-squared:  0.5663,    Adjusted R-squared:  0.4631 
F-statistic: 5.485 on 5 and 21 DF,  p-value: 0.002205

> coeffHC4 <- coeftest(model, vcov = vcovHC(model, type = "HC4"))
> coeffHC4

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  0.4962902  0.0808139  6.1411 4.299e-06 ***
opi         -0.4245559  0.1508714 -2.8140  0.010397 *  
competence   0.0100454  0.0044635  2.2506  0.035257 *  
integration -0.2381628  0.0805296 -2.9575  0.007517 ** 
indegree     0.0234134  0.0098380  2.3799  0.026873 *  
voterank    -0.0022659  0.0023894 -0.9483  0.353756    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> conf_ints <- confint (coeffHC4, level = 0.95) |> 
     as.data.frame () |> 
     tibble::rownames_to_column("Variables") |>
     `colnames<-`(c("Variables", "95%CI_low", "95%CI_hi"))
> conf_ints

    Variables     95%CI_low     95%CI_hi
1 (Intercept)  0.3282284541  0.664351977
2         opi -0.7383101090 -0.110801602
3  competence  0.0007629805  0.019327796
4 integration -0.4056331442 -0.070692390
5    indegree  0.0029541503  0.043872647
6    voterank -0.0072348011  0.002703083

> VIF_tol = vif(model) |>
     as.data.frame () |>
     tibble::rownames_to_column("Variables") |>
     mutate (Tolerance = 1/`vif(model)`) |>
     `colnames<-`(c("Variables", "VIF", "Tolerance"))
> VIF_tol
           Variables      VIF Tolerance
1         scale(opi) 1.493807 0.6694305
2  scale(competence) 1.175722 0.8505409
3 scale(integration) 1.392198 0.7182888
4    scale(indegree) 1.050531 0.9518994
5    scale(voterank) 1.414064 0.7071815

> model_summary <- summary(model)
> output <- model_summary$coefficients |> 
     as.data.frame () |>
     tibble::rownames_to_column("Variables") |>
     left_join (conf_ints, by = "Variables") |> 
     left_join (VIF_tol, by = "Variables") |> 
     mutate (across(c(2:4, 6:9), .fns = function(x) {format(round(x, 5), nsmall = 5)})) |> 
     relocate (`95%CI_low`, .after = Estimate) |>
     relocate (`95%CI_hi`, .after = `95%CI_low`)
> output[ ,7] <- format.pval(output[ ,7], eps = .001, digits = 4)
> output

    Variables Estimate 95%CI_low 95%CI_hi Std. Error  t value Pr(>|t|)     VIF Tolerance
1 (Intercept)  0.49629   0.32823  0.66435    0.07814  6.35133  < 0.001      NA        NA
2         opi -0.42456  -0.73831 -0.11080    0.15333 -2.76892  0.01150 1.49381   0.66943
3  competence  0.01005   0.00076  0.01933    0.00446  2.25151  0.03519 1.17572   0.85054
4 integration -0.23816  -0.40563 -0.07069    0.09940 -2.39592  0.02597 1.39220   0.71829
5    indegree  0.02341   0.00295  0.04387    0.01354  1.72862  0.09855 1.05053   0.95190
6    voterank -0.00227  -0.00723  0.00270    0.00306 -0.74103  0.46688 1.41406   0.70718

# looks great, but the columns for SE, t and p values are pulled from the simple linear regression, not the one with HC4. I think there is a straightforward way to paste these in, will try it in another spare while. And then I will also try to add:

> etaSquared(model)
                eta.sq eta.sq.part
opi         0.15832698  0.26744774
competence  0.10468472  0.19445477
integration 0.11854396  0.21467219
indegree    0.06170723  0.12456734
voterank    0.01133971  0.02548222

# I will need to learn to add the eta.sq.part column to the output

> library(rempsyc)
> table <- nice_table(output)
> flextable::save_as_docx(table, path = "table.docx")

As some readers like to see standardized betas as well, I've generated those, too:

> model_std <- lm(scale(progress) ~ scale(opi) + scale(competence) + scale(integration) + scale(indegree) + scale(voterank), data = Amman)
> coeffHC4_std <- coeftest(model_std, vcov = vcovHC(model_std, type = "HC4"))
> conf_ints <- confint (coeffHC4_std, level = 0.95) |> 
     as.data.frame () |> 
     tibble::rownames_to_column("Variables") |>
     `colnames<-`(c("Variables", "95%CI_low", "95%CI_hi"))
> model_std_summary <- summary(model_std)
output_std <- model_std_summary$coefficients |> 
     as.data.frame () |>
     tibble::rownames_to_column("Variables") |>
     left_join (conf_ints, by = "Variables") |> 
     mutate (across(c(2:4, 6:7), .fns = function(x) {format(round(x, 5), nsmall = 5)})) |> 
     relocate (`95%CI_low`, .after = Estimate) |>
     relocate (`95%CI_hi`, .after = `95%CI_low`)
> output_std[ ,7] <- format.pval(output_std[ ,7], eps = .001, digits = 4)
> output_std

           Variables Estimate 95%CI_low 95%CI_hi Std. Error  t value Pr(>|t|)
1        (Intercept)  0.17556  -0.14697  0.49809    0.15728  1.11623  0.27693
2         scale(opi) -0.47667  -0.82894 -0.12440    0.17215 -2.76892  0.01150
3  scale(competence)  0.38072   0.02892  0.73252    0.16909  2.25151  0.03519
4 scale(integration) -0.49318  -0.83997 -0.14639    0.20584 -2.39592  0.02597
5    scale(indegree)  0.26905   0.03395  0.50415    0.15564  1.72862  0.09855
6    scale(voterank) -0.13614  -0.43469  0.16241    0.18372 -0.74103  0.46688

And finally I'll need to knit the first three columns from the above into the previous table, together with the partial eta squared, and it's done.

I would like to generate one regression table with the HC4 estimator so that the columns for SE, t and p values are pulled from the regression with robust SEs (HC4) not the regular linear regression, include beta values + their confidence intervals, as well as VIF and Tolerance all in one table. Almost there thanks to the great help form @RGrowinG!

Here is the dataset if it helps: https://docs.google.com/spreadsheets/d/13U7xDK7otj-vBOcx6yTpsguAd-wlFNad/edit?usp=sharing&ouid=113244913961284427801&rtpof=true&sd=true

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

Getting from SPSS to R can by hard at first, when you are used to getting the coefficients you need to conduct regression diagnostics straight from the default output. But the more time you are spending using R, the more you recognize the freedom it is providing by just offering a wealth of possibilities and options.

Here is an approach to get the results you are looking for:

#------------------------
# Load required packages
#------------------------
library (dplyr)
library (car)


#------------------------
# Fit a linear regression model (using the iris dataset in this example)
# Sepal.Length is the dependent variable
# Sepal.Width, Petal.Length, Petal.Width are the independent variables (predictors)
#------------------------
fit <- iris |>
  lm (Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data=_) # mind using placeholder "_" for base pipe here


#------------------------
# Create an object containing the model summary information
#------------------------
mod_summary <- summary (fit)


#------------------------
# Create a data.frame containing 95%CIs for Estimate
#------------------------
conf_ints <- confint (fit, level = 0.95) |> 
  as.data.frame () |> 
  tibble::rownames_to_column("Variables") |> # make rownames a column called "Variables"
  `colnames<-`(c("Variables", "95%CI_low", "95%CI_hi")) # set column names

#------------------------
# Create a data.frame containing VIF & Tolerance
#------------------------
VIF_tol = vif(fit) |> # calculate VIF
  as.data.frame () |>
  tibble::rownames_to_column("Variables") |>
  mutate (Tolerance = 1/`vif(fit)`) |> # calculate Tolerance given by 1/VIF
  `colnames<-`(c("Variables", "VIF", "Tolerance"))


#------------------------
# Put it all together using left_join by "Variables", format & arrange columns
#------------------------
output <- mod_summary$coefficients |> # use the coefficients table from the model summary as starting point
  as.data.frame () |>
  tibble::rownames_to_column("Variables") |>
  left_join (conf_ints, by = "Variables") |> # left_join with conf_ints data.frame
  left_join (VIF_tol, by = "Variables") |> # left_join with VIF_tol data.frame
  mutate (across(c(2:4, 6:9), .fns = function(x) {format(round(x, 2), nsmall = 2)})) |> # set decimal numbers for all columns (except p-values) to 2 
  relocate (`95%CI_low`, .after = Estimate) |>
  relocate (`95%CI_hi`, .after = `95%CI_low`)

output[ ,7] <- format.pval(output[ ,7], eps = .001, digits = 3) # format small p-values < 0.001 nicely

output

     Variables Estimate 95%CI_low 95%CI_hi Std. Error t value Pr(>|t|)   VIF Tolerance
1  (Intercept)     1.86      1.36     2.35       0.25    7.40   <0.001    NA        NA
2  Sepal.Width     0.65      0.52     0.78       0.07    9.77   <0.001  1.27      0.79
3 Petal.Length     0.71      0.60     0.82       0.06   12.50   <0.001 15.10      0.07
4  Petal.Width    -0.56     -0.81    -0.30       0.13   -4.36   <0.001 14.23      0.07

Regarding your edited question, please see if that helps. The results for this example are nearly identical but may vary in case of your data.

#------------------------
# Load required packages
#------------------------
library (dplyr)
library (car)
library (lmtest)
library (sandwich)


#------------------------
# Fit a linear regression model (using the iris dataset in this example)
# Sepal.Length is the dependent variable
# Sepal.Width, Petal.Length, Petal.Width are the independent variables (predictors)
#------------------------
fit <- iris |>
  lm (Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data= _)


#------------------------
# Create an object containing the model summary information with robust estimators
#------------------------
#mod_summary <- summary (fit)


#------------------------
# Create an coefficient tables (t-Test) using HC3 and HC4 types
#------------------------
coeffs1 <- coeftest(fit, vcov = vcovHC(fit, type = "HC3")) # HC3
coeffs2 <- coeftest(fit, vcov = vcovHC(fit, type = "HC4")) # HC4

list_coeffs <- list (coeffs1, coeffs2) # make a list containing the coefficient tables for lapply processing


#------------------------
# Create a list of data.frames containing 95%CIs & VIF (plus Tolerance)
#------------------------
list_results <- lapply (list_coeffs, function (x) {

 
## 95%CIs   
conf_ints <- confint (x, level = 0.95) |> 
    as.data.frame () |> 
    tibble::rownames_to_column("Variables") |> # make rownames a column called "Variables"
    `colnames<-`(c("Variables", "95%CI_low", "95%CI_hi")) # set column names


## VIF & Tolerance
VIF_tol = vif(fit) |> # calculate VIF
    as.data.frame () |>
    tibble::rownames_to_column("Variables") |>
    mutate (Tolerance = 1/`vif(fit)`) |> # calculate Tolerance given by 1/VIF
    `colnames<-`(c("Variables", "VIF", "Tolerance"))
  
  
## Put it all together using left_join by "Variables" & format table nicely
  output <- x[,] |> # use the robust coefficients table 
    as.data.frame () |>
    tibble::rownames_to_column("Variables") |>
    left_join (conf_ints, by = "Variables") |> # left_join with conf_ints data.frame
    left_join (VIF_tol, by = "Variables") |> # left_join with VIF_tol data.frame
    mutate (across(c(2:4, 6:9), .fns = function(x) {format(round(x, 2), nsmall = 2)})) |> # set decimal numbers for all columns (except p-values) to 2 
    relocate (`95%CI_low`, .after = Estimate) |>
    relocate (`95%CI_hi`, .after = `95%CI_low`)
  
  output[ ,7] <- format.pval(output[ ,7], eps = .001, digits = 3) # format small p-values < 0.001 nicely
  
  output
    
})

list_results 

[[1]]
     Variables Estimate 95%CI_low 95%CI_hi Std. Error t value Pr(>|t|)   VIF Tolerance
1  (Intercept)     1.86      1.40     2.32       0.23    7.99   <0.001    NA        NA
2  Sepal.Width     0.65      0.53     0.77       0.06   10.48   <0.001  1.27      0.79
3 Petal.Length     0.71      0.58     0.84       0.06   11.04   <0.001 15.10      0.07
4  Petal.Width    -0.56     -0.85    -0.26       0.15   -3.74   <0.001 14.23      0.07

[[2]]
     Variables Estimate 95%CI_low 95%CI_hi Std. Error t value Pr(>|t|)   VIF Tolerance
1  (Intercept)     1.86      1.40     2.32       0.23    7.99   <0.001    NA        NA
2  Sepal.Width     0.65      0.53     0.77       0.06   10.50   <0.001  1.27      0.79
3 Petal.Length     0.71      0.58     0.84       0.06   11.03   <0.001 15.10      0.07
4  Petal.Width    -0.56     -0.85    -0.26       0.15   -3.74   <0.001 14.23      0.07

From the answer in the following question, it seems to me that VIF does not need to be adapted. https://stackoverflow.com/q/59041445/11809780

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks so much, @GRowinG! These seem to work for the regular regression. What to do if I need one with "robust" (heteroscedasticity-consistent) standard errors? I've edited the original post to show the code used so far. $\endgroup$
    – mbp
    Commented Apr 14 at 15:02
  • $\begingroup$ Thanks so much again, @RGrowinG! This is great. I've just updated my question above. The only issue so far is that the columns for SE, t and p values are pulled from the simple linear regression, not the one with HC4 (which I eventually decided on). I'll also try to add in partial eta squared values as well as standardized betas + CIs. Generated them, will just try to add them into one table at one go rather than manually in post-editing. $\endgroup$
    – mbp
    Commented Apr 18 at 17:27
  • $\begingroup$ BTW, @RGrowinG, if you leave your full name, I'd like to include you in the acknowledgments section of our paper. $\endgroup$
    – mbp
    Commented Jun 4 at 8:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.