robust() computes robust standard error for regression models. This method calls one of the vcov*()-functions from the sandwich-package for robust covariance matrix estimators. Results are returned as tidy data frame.

svy() is intended to compute standard errors for survey designs (complex samples) fitted with regular lm or glm functions, as alternative to the survey-package. It simulates sampling weights by adjusting the residual degrees of freedom based on the precision weights used to fit x, and then calls robust() with the adjusted model.

robust(x, vcov.fun = "vcovHC", vcov.type = c("HC3", "const", "HC",
  "HC0", "HC1", "HC2", "HC4", "HC4m", "HC5"), vcov.args = NULL,
  conf.int = FALSE, exponentiate = FALSE)

svy(x, vcov.fun = "vcovHC", vcov.type = c("HC1", "const", "HC", "HC0",
  "HC3", "HC2", "HC4", "HC4m", "HC5"), vcov.args = NULL,
  conf.int = FALSE, exponentiate = FALSE)

Arguments

x

A fitted model of any class that is supported by the vcov*()-functions from the sandwich package. For svy(), x must be lm object, fitted with weights.

vcov.fun

String, indicating the name of the vcov*()-function from the sandwich-package, e.g. vcov.fun = "vcovCL".

vcov.type

Character vector, specifying the estimation type for the robust covariance matrix estimation (see vcovHC for details).

vcov.args

List of named vectors, used as additional arguments that are passed down to vcov.fun.

conf.int

Logical, TRUE if confidence intervals based on robust standard errors should be included.

exponentiate

Logical, whether to exponentiate the coefficient estimates and confidence intervals (typical for logistic regression).

Value

A summary of the model, including estimates, robust standard error, p-value and - optionally - the confidence intervals.

Note

svy() simply calls robust(), but first adjusts the residual degrees of freedom based on the model weights. Hence, for svy(), x should be fitted with weights. This simulates sampling weights like in survey designs, though lm and glm implement precision weights. The results from svy() are usually more accurate than simple weighted standard errors for complex samples. However, results from the survey package are still more exactly, especially regarding the estimates.

vcov.type for svy() defaults to "HC1", because standard errors with this estimation type come closest to the standard errors from the survey-package.

Currently, svy() only works for objects of class lm.

Examples

data(efc) fit <- lm(barthtot ~ c160age + c12hour + c161sex + c172code, data = efc) summary(fit)
#> #> Call: #> lm(formula = barthtot ~ c160age + c12hour + c161sex + c172code, #> data = efc) #> #> Residuals: #> Min 1Q Median 3Q Max #> -74.639 -15.246 4.251 19.009 73.327 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 90.06448 6.17237 14.592 <2e-16 *** #> c160age -0.22156 0.07111 -3.116 0.0019 ** #> c12hour -0.27810 0.01865 -14.915 <2e-16 *** #> c161sex -0.26178 2.08649 -0.125 0.9002 #> c172code -0.76215 1.41971 -0.537 0.5915 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Residual standard error: 25.35 on 816 degrees of freedom #> (87 observations deleted due to missingness) #> Multiple R-squared: 0.2696, Adjusted R-squared: 0.266 #> F-statistic: 75.28 on 4 and 816 DF, p-value: < 2.2e-16 #>
robust(fit)
#> term estimate std.error statistic p.value #> 1 (Intercept) 90.0644792 5.58174936 16.1355291 4.832660e-51 #> 2 c160age -0.2215581 0.07106417 -3.1177189 1.886499e-03 #> 3 c12hour -0.2781004 0.02119537 -13.1208077 8.012747e-36 #> 4 c161sex -0.2617783 1.91864853 -0.1364389 8.915080e-01 #> 5 c172code -0.7621525 1.39456606 -0.5465159 5.848608e-01
confint(fit)
#> 2.5 % 97.5 % #> (Intercept) 77.9488902 102.18006829 #> c160age -0.3611297 -0.08198647 #> c12hour -0.3146997 -0.24150107 #> c161sex -4.3573007 3.83374416 #> c172code -3.5488594 2.02455439
robust(fit, conf.int = TRUE)
#> term estimate std.error conf.low conf.high statistic #> 1 (Intercept) 90.0644792 5.58174936 79.1082006 101.02075786 16.1355291 #> 2 c160age -0.2215581 0.07106417 -0.3610482 -0.08206799 -3.1177189 #> 3 c12hour -0.2781004 0.02119537 -0.3197042 -0.23649650 -13.1208077 #> 4 c161sex -0.2617783 1.91864853 -4.0278463 3.50428977 -0.1364389 #> 5 c172code -0.7621525 1.39456606 -3.4995120 1.97520692 -0.5465159 #> p.value #> 1 4.832660e-51 #> 2 1.886499e-03 #> 3 8.012747e-36 #> 4 8.915080e-01 #> 5 5.848608e-01
robust(fit, vcov.type = "HC1", conf.int = TRUE) # "HC1" should be Stata default
#> term estimate std.error conf.low conf.high statistic #> 1 (Intercept) 90.0644792 5.55391758 79.1628309 100.96612755 16.2163874 #> 2 c160age -0.2215581 0.07066514 -0.3602650 -0.08285124 -3.1353241 #> 3 c12hour -0.2781004 0.02103734 -0.3193941 -0.23680669 -13.2193675 #> 4 c161sex -0.2617783 1.90939743 -4.0096876 3.48613101 -0.1370999 #> 5 c172code -0.7621525 1.38747118 -3.4855856 1.96128056 -0.5493105 #> p.value #> 1 1.788307e-51 #> 2 1.778136e-03 #> 3 2.730414e-36 #> 4 8.909856e-01 #> 5 5.829427e-01
library(sjmisc) # dichtomozize service usage by "service usage yes/no" efc$services <- sjmisc::dicho(efc$tot_sc_e, dich.by = 0) fit <- glm(services ~ neg_c_7 + c161sex + e42dep, data = efc, family = binomial(link = "logit")) robust(fit)
#> term estimate std.error statistic p.value #> 1 (Intercept) -0.5198021 0.38048112 -1.366171 0.172232067 #> 2 neg_c_7 0.0419026 0.02035790 2.058296 0.039853669 #> 3 c161sex -0.2189336 0.16263979 -1.346126 0.178606278 #> 4 e42dep 0.2134784 0.07970607 2.678320 0.007536269
robust(fit, conf.int = TRUE, exponentiate = TRUE)
#> term estimate std.error conf.low conf.high statistic p.value #> 1 (Intercept) 0.5946382 0.38048112 0.2818017 1.254764 -1.366171 0.172232067 #> 2 neg_c_7 1.0427929 0.02035790 1.0019492 1.085302 2.058296 0.039853669 #> 3 c161sex 0.8033751 0.16263979 0.5838345 1.105470 -1.346126 0.178606278 #> 4 e42dep 1.2379767 0.07970607 1.0587020 1.447609 2.678320 0.007536269