Linear combinations of coefficients in R

Aim to reproduce the function of the “lincom” Stata command as in pg 96 of the notes.

Load the required libraries:

library(foreign)
library(multcomp)

## Loading required package: mvtnorm Loading required package: survival
## Loading required package: splines

Load the dataset and generate the required Age variable:

bp <- read.dta("~/Dropbox/Biostats/LMR/Module1/sbp_reg.dta")
bp$agehi <- bp$age >= 40
bp$agehi <- factor(bp$agehi)

Generate the model with interactions

mod <- lm(sbp ~ gender * agehi, data = bp)

## 
## Call:
## lm(formula = sbp ~ gender * agehi, data = bp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.686  -6.684  -0.383   5.423  12.518 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            131.11       3.36   39.07   <2e-16 ***
## genderMale               3.74       5.31    0.71    0.489    
## agehiTRUE               11.24       4.98    2.26    0.035 *  
## genderMale:agehiTRUE    -2.27       7.01   -0.32    0.749    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.22 on 20 degrees of freedom
## Multiple R-squared:  0.343,  Adjusted R-squared:  0.245 
## F-statistic: 3.49 on 3 and 20 DF,  p-value: 0.0349

Use “glht” to form the linear combination

names(coef(mod))  #Extract the coefficient names (b0,b1,b2,b3)

## [1] "(Intercept)"          "genderMale"           "agehiTRUE"           
## [4] "genderMale:agehiTRUE"

# linfct specifies the required combination In this case we want b1+b3=0
# (effect of gender in age>40)
summary(glht(mod, linfct = c("genderMale + genderMale:agehiTRUE = 0")))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = sbp ~ gender * agehi, data = bp)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error t value
## genderMale + genderMale:agehiTRUE == 0     1.47       4.58    0.32
##                                        Pr(>|t|)
## genderMale + genderMale:agehiTRUE == 0     0.75
## (Adjusted p values reported -- single-step method)

# Show the confidence interval
mod.lh <- glht(mod, linfct = c("genderMale + genderMale:agehiTRUE = 0"))
confint(mod.lh)

## 
##   Simultaneous Confidence Intervals
## 
## Fit: lm(formula = sbp ~ gender * agehi, data = bp)
## 
## Quantile = 2.086
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                                        Estimate lwr    upr   
## genderMale + genderMale:agehiTRUE == 0  1.469   -8.096 11.033