Documentation

library(bios6312)

The following function are available:

• get_robust_variance
• get_studentized_residuals
• qqnorm_studentized_residuals
• studentized_residuals_vs_fitted_plot
• studentized_residuals_vs_predictors_plot
• testparm
• lincom

To get help on any function, you can access the documentation with ?example_function.

We will use the built-in mtcars dataset to demonstrate each function.

head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

my_model <- lm(mpg~wt, data=mtcars)

summary(my_model)

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

get_robust_variance(my_model)

## Robust Variance 
##             (Intercept)         wt
## (Intercept)    4.817955 -1.3927648
## wt            -1.392765  0.4283488
## 
## 
## Robust Standard Error 
## (Intercept)          wt 
##   2.1949842   0.6544836

get_studentized_residuals(my_model)

##           Mazda RX4       Mazda RX4 Wag          Datsun 710      Hornet 4 Drive 
##         -0.77883244         -0.31251235         -0.71741861          0.43990454 
##   Hornet Sportabout             Valiant          Duster 360           Merc 240D 
##         -0.06792331         -0.23530952         -1.32710250          1.41185562 
##            Merc 230            Merc 280           Merc 280C          Merc 450SE 
##          0.79688521          0.10176282         -0.37335834          0.29772794 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood Lincoln Continental 
##         -0.01711622         -0.64204006          0.42995227          0.78252353 
##   Chrysler Imperial            Fiat 128         Honda Civic      Toyota Corolla 
##          2.20946174          2.37349819          0.62044614          2.25373119 
##       Toyota Corona    Dodge Challenger         AMC Javelin          Camaro Z28 
##         -0.89418062         -1.00955292         -1.26474650         -1.18201197 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2        Lotus Europa 
##          0.84140193          0.12446808          0.05262766          0.42952735 
##      Ford Pantera L        Ferrari Dino       Maserati Bora          Volvo 142E 
##         -1.54054834         -0.94625660         -1.08923170         -0.34956654

qqnorm_studentized_residuals(my_model)

studentized_residuals_vs_fitted_plot(my_model)

studentized_residuals_vs_predictors_plot(my_model)

For the next functions, we will need to use a model with interactions.

my_interaction_model <- lm(mpg ~ wt*cyl, data=mtcars)

summary(my_interaction_model)

## 
## Call:
## lm(formula = mpg ~ wt * cyl, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2288 -1.3495 -0.5042  1.4647  5.2344 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  54.3068     6.1275   8.863 1.29e-09 ***
## wt           -8.6556     2.3201  -3.731 0.000861 ***
## cyl          -3.8032     1.0050  -3.784 0.000747 ***
## wt:cyl        0.8084     0.3273   2.470 0.019882 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.368 on 28 degrees of freedom
## Multiple R-squared:  0.8606, Adjusted R-squared:  0.8457 
## F-statistic: 57.62 on 3 and 28 DF,  p-value: 4.231e-12

# Replicates the `testparm` function in Stata
# Test if beta2 + beta3 == 0
testparm(par=c(3,4),
         model=my_interaction_model,
         type="F")

##            F            P 
## 1.182702e+01 1.891263e-04

# Replicates the `lincom` function in Stata
# Estimates beta2 + 3.2beta3
lincom(par=c(3,4),
       mult=c(1,3.2),
       model=my_interaction_model)

##          EST        CI.LO        CI.HI            P 
## -1.216355561 -1.967601441 -0.465109682  0.002530604