In order to solve heteroscedasticity in my data, I ran a regression with heteroscedasticity-consistent ("robust") standard errors. I used SPSS, but to my current knowledge that is only able to output the following values (columns): heteroscedasticity-consistent SEs I would also like to report i) standardized betas together with confidence intervals, ii) Tolerance andpartial eta squared values, iii) Tolerance and VIF. Is there a way to get these values with R? The R outputs for parameter estimates with robust SEs e.g. using the sandwich package do not seem to readily provide them. I would like to be able to generate a 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).
> library(dplyr)
> libraryload(ggplot2)
>file library(car)
>= data("Amman""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
# The above yield the standard regression. My data needs "robust" SEs. Trying two different packages:
# sandwich
> library(lmtest)
> library(sandwich)
> coeftest(model, vcov = vcovHC(model, type = "HC3"))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4962902 0.0857433 5.7881 9.575e-06 ***
opi -0.4245559 0.1660452 -2.5569 0.01837 *
competence 0.0100454 0.0049919 2.0123 0.05719 .
integration -0.2381628 0.0847770 -2.8093 0.01051 *
indegree 0.0234134 0.0102126 2.2926 0.03230 *
voterank -0.0022659 0.0026184 -0.8654 0.39662 > coeffHC4
--<-
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 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
#> sandwichconf_ints affords<- bothconfint HC3(coeffHC4, andlevel HC4,= but0.95) does|> not
as.data.frame (automatically?) provide|> confidence
intervals tibble::rownames_to_column("Variables") |>
# estimatr `colnames<-`(c("Variables", "95%CI_low", "95%CI_hi"))
> conf_ints
> library Variables 95%CI_low 95%CI_hi
1 (estimatrIntercept) 0.3282284541 0.664351977
>2 lmrobust <- lm_robust(progress ~ opi +-0.7383101090 -0.110801602
3 competence + 0.0007629805 0.019327796
4 integration +-0.4056331442 indegree-0.070692390
5 + voterank, data =indegree Amman, se_type0.0029541503 = "HC3")0.043872647
>6 summary(lmrobust) voterank -0.0072348011 0.002703083
Standard> errorVIF_tol type= vif(model) |>
as.data.frame () |>
tibble::rownames_to_column("Variables") |>
HC3 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
Coefficients:
> model_summary <- summary(model)
> output <- model_summary$coefficients |>
Estimate Stdas. Error t value data.frame Pr(>|t|) CI|>
Lower CI Upper DF
tibble::rownames_to_column(Intercept"Variables") 0.496290|>
0.085743 5.7881left_join 9.575e-06(conf_ints, by 0.3179773= "Variables") 0.674603|> 21
opi -0.424556 left_join 0.166045(VIF_tol, -2.5569by 1.837e-02= -0.7698657"Variables") -0.079246|> 21
competence 0.010045 0.004992 mutate (across(c(2.0123:4, 5.719e-026:9), -0.0003358 fns 0.020427= 21
integrationfunction(x) -0.238163{format(round(x, 5), nsmall 0.084777= -2.80935)})) 1.051e-02|> -0.4144662
-0.061859 21
indegree relocate (`95%CI_low`, 0.023413after = Estimate) 0.010213|>
2.2926 3.230e-02 0.0021751relocate (`95%CI_hi`, 0.044652after 21
voterank= `95%CI_low`)
> output[ ,7] <-0.002266 format.pval(output[ ,7], 0.002618eps -0.8654= 3.966e-01001, -0.0077111digits = 0.0031794)
> 21output
Multiple R 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 -squared:0.42456 -0.566373831 -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 Adjustedcolumns Rfor 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 <-squared 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.463195) |>
F as.data.frame () |>
tibble::rownames_to_column("Variables") |>
`colnames<-statistic`(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), .651fns on= function(x) {format(round(x, 5), andnsmall 21= DF5)})) |>
relocate (`95%CI_low`, .after = pEstimate) |>
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.000739817556 -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.
The estimatr package provides confidence intervals, but doesn't seem to offer the HC4 option. II would like to generate regressionsone regression table with HC3 andthe HC4 estimator so that the columns for SE, including betat and p values are pulled from the regression with robust SEs (if they differ fromHC4) not the regular linear regression), include beta values + their confidence intervals for beta, as well as VIF and Tolerance valuesall in one table. Almost there thanks to the great help form @RGrowinG!
