8.4 Weighted linear regression

8.4.1 Fitting the model

To carry out a linear regression that incorporates a survey design, use svyglm() with family=gaussian(). “Gaussian” means “normally distributed” so this is specifying a model with an outcome that, given the predictors, is normally distributed, which is equivalent to specifying a linear regression with normally distributed residuals. See, for example, Lumley and Scott (2017) and for further reading about fitting regression models to complex survey data.

After fitting the model, use summary(), confint(), and car::Anova(, type = 3, test.statistic = "F") to view the regression coefficients, their 95% confidence intervals, and p-values. Each of these functions has an option to change the DF to the design DF. Finally, table the regression results using tbl_regression(). There is an option to change the DF in tbl_regression() but it only applies to the overall p-values for categorical predictors, so the confidence intervals and remaining p-values in the table are based on default DF rather than the design DF. See the “Details” section in ?summary.svyglm for information regarding the default DF.

Example 8.1 (continued): Regress the outcome fasting glucose (LBDGLUSI) on the predictors waist circumference (BMXWAIST), smoking status (smoker) , age (RIDAGEYR), gender (RIAGENDR), race/ethnicity (RIDRETH3), and income (income), correctly incorporating the NHANES survey design to obtain estimates representative of the population. Use the NHANES 2017-2018 fasting subset survey design we have been working with in this chapter (design.FST.nomiss).

fit.ex8.1 <- svyglm(LBDGLUSI ~ BMXWAIST + smoker + RIDAGEYR +
                      RIAGENDR + race_eth + income,
                 family=gaussian(), design=design.FST.nomiss)

When calling summary(), confint(), and car::Anova() to summarize the results, use the df.resid, ddf, and error.df options, respectively, to specify the design DF for confidence intervals and hypothesis tests. To ensure you are using the same design as was used in the regression, extract the design from the svyglm object as fit$survey.design.

round(
  cbind(
    summary(fit.ex8.1, df.resid = degf(fit.ex8.1$survey.design))$coef,
    confint(fit.ex8.1, ddf      = degf(fit.ex8.1$survey.design))
  )
, 4)
##                            Estimate Std. Error t value Pr(>|t|)   2.5 %  97.5 %
## (Intercept)                  3.3890     0.3542  9.5670   0.0000  2.6340  4.1441
## BMXWAIST                     0.0221     0.0031  7.0544   0.0000  0.0154  0.0288
## smokerPast                   0.1914     0.0961  1.9916   0.0649 -0.0134  0.3963
## smokerCurrent                0.0055     0.1103  0.0499   0.9609 -0.2297  0.2407
## RIDAGEYR                     0.0210     0.0020 10.2867   0.0000  0.0166  0.0253
## RIAGENDRFemale              -0.3331     0.0684 -4.8694   0.0002 -0.4789 -0.1873
## race_ethNon-Hispanic White  -0.3448     0.1461 -2.3599   0.0322 -0.6561 -0.0334
## race_ethNon-Hispanic Black  -0.2577     0.1565 -1.6470   0.1203 -0.5912  0.0758
## race_ethNon-Hispanic Other  -0.1095     0.1381 -0.7928   0.4402 -0.4040  0.1849
## income$25,000 to < $55,000  -0.1437     0.1414 -1.0162   0.3256 -0.4452  0.1577
## income$55,000+              -0.1224     0.0895 -1.3674   0.1917 -0.3132  0.0684
car::Anova(fit.ex8.1, type = 3, test.statistic = "F", 
           error.df = degf(fit.ex8.1$survey.design))
## Analysis of Deviance Table (Type III tests)
## 
## Response: LBDGLUSI
##             Df      F      Pr(>F)    
## (Intercept)  1  91.53 0.000000089 ***
## BMXWAIST     1  49.77 0.000003909 ***
## smoker       2   2.09      0.1584    
## RIDAGEYR     1 105.82 0.000000034 ***
## RIAGENDR     1  23.71      0.0002 ***
## race_eth     3   3.00      0.0636 .  
## income       2   1.03      0.3793    
## Residuals   15                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

To table these results, use the tbl_regression() function in gtsummary after adding test.statistic = "F" in add_global_p(). The results are shown in Table 8.3.

Warning: There is no option in tbl_regression() to change the DF to make the confidence intervals and p-values for the individual regression coefficients be based on the design DF. There is an option in add_global_p() to change the DF for the Type III tests (commented out in the code below), however using that option will lead to p-values that do not match between the overall (Type III) test for a binary predictor and the test for the single term in the model for that binary predictor. The code below uses R’s default DF throughout. If, instead, you want to use the design DF throughout, uncomment the error.df row below, assign the table to an object, and export (see Sections 3.3.2 and 3.3.3). Then, in a word processing program, edit the confidence intervals and p-values for each regression coefficient, replacing them with those you computed above using the design DF.

library(gtsummary)
fit.ex8.1 %>%
  tbl_regression(intercept = T,
                 estimate_fun = function(x) style_sigfig(x, digits = 3),
                 pvalue_fun   = function(x) style_pvalue(x, digits = 3),
                 label  = list(BMXWAIST ~ "Waist Circumference (cm)",
                               smoker   ~ "Smoking Status",
                               RIDAGEYR ~ "Age (years)",
                               RIAGENDR ~ "Gender",
                               race_eth ~ "Race/ethnicity",
                               income   ~ "Annual Income")) %>%
  add_global_p(keep = T, test.statistic = "F"
  # Uncomment to use design df for Type III test p-values
  # (but then they will not match the regression term p-values for
  #  binary predictors)
  #             , error.df = degf(fit.ex8.1$survey.design)
                                 ) %>%
  modify_caption("Weighted linear regression results for fasting glucose (mmol/L)")
Table 8.3: Weighted linear regression results for fasting glucose (mmol/L)
Characteristic Beta 95% CI1 p-value
(Intercept) 3.39 2.48, 4.30 <0.001
Waist Circumference (cm) 0.022 0.014, 0.030 <0.001
Smoking Status

0.219
    Never
    Past 0.191 -0.056, 0.438 0.103
    Current 0.006 -0.278, 0.289 0.962
Age (years) 0.021 0.016, 0.026 <0.001
Gender

0.005
    Male
    Female -0.333 -0.509, -0.157 0.005
Race/ethnicity

0.134
    Hispanic
    Non-Hispanic White -0.345 -0.720, 0.031 0.065
    Non-Hispanic Black -0.258 -0.660, 0.145 0.160
    Non-Hispanic Other -0.110 -0.465, 0.246 0.464
Annual Income

0.421
    < $25,000
    $25,000 to < $55,000 -0.144 -0.507, 0.220 0.356
    $55,000+ -0.122 -0.353, 0.108 0.230
1 CI = Confidence Interval

Conclusion: After adjusting for age, gender, race/ethnicity, income, and each other, waist circumference is significantly associated with fasting glucose (B = 0.022; 95% CI = 0.014, 0.030; p < .001) but smoking status is not (p = .219).

8.4.2 Prediction

As was the case with glm(), predict() applied to a svyglm() object does not return a CI. Instead, use svycontrast_design_df() (in Functions_rmph.R, which you loaded at the beginning of this chapter), a modified version of svycontrast() with a DF option, which is similar to gmodels::estimable(). See Chapter 6.8 for how to specify the second argument to get the prediction of interest.

Example 8.1 (continued): What is the estimated population mean fasting glucose among individuals with a waist circumference of 120 cm, who have never smoked, are age 35 years, female, and Hispanic, and have an annual income from $25,000 to < $55,000?

# Always include the intercept for prediction.
# Specify a 1 for the intercept, a # for each continuous predictor
# and a 1 for each non-reference level of a categorical variable.
# If a predictor is at its reference level, it should be left out.
svycontrast_design_df(fit.ex8.1,
                      c("(Intercept)"                = 1,
                        "BMXWAIST"                   = 120,
                        "RIDAGEYR"                   = 35,
                        "RIAGENDRFemale"             = 1,
                        "income$25,000 to < $55,000" = 1))
##     est lower upper
## 1 6.302 6.018 6.585

Conclusion: The estimated population mean fasting glucose among such individuals is 6.30 mmol/L (95% CI = 6.02, 6.59).

8.4.3 Interactions

The syntax for including an interaction in svyglm() is exactly as in lm(): add X:Z to include an interaction between X and Z.

Example 8.2: Does the association between the outcome fasting glucose (LBDGLUSI) and waist circumference (BMXWAIST) depend on gender (RIAGENDR)? Correctly incorporate the NHANES survey design to obtain an estimate of the waist circumference slope for each gender that is generalizable to the population.

fit.ex8.2 <- svyglm(LBDGLUSI ~ cBMXWAIST + RIAGENDR + cBMXWAIST:RIAGENDR,
                 family=gaussian(), design=design.FST.nomiss)

round(
  cbind(
    summary(fit.ex8.2,
              df.resid=degf(fit.ex8.2$survey.design))$coef,
    confint(fit.ex8.2,
              ddf=degf(fit.ex8.2$survey.design))
  )
, 4)
##                          Estimate Std. Error t value Pr(>|t|)   2.5 %  97.5 %
## (Intercept)                6.2555     0.0758  82.513   0.0000  6.0939  6.4170
## cBMXWAIST                  0.0317     0.0043   7.400   0.0000  0.0226  0.0408
## RIAGENDRFemale            -0.3139     0.0600  -5.232   0.0001 -0.4418 -0.1861
## cBMXWAIST:RIAGENDRFemale  -0.0084     0.0042  -1.993   0.0648 -0.0174  0.0006

The p-value for cBMXWAIST:RIAGENDRFemale is larger than 0.05 so the interaction is not statistically significant.

To estimate the waist circumference slope at each level of gender, use svycontrast_design_df() (in Functions_rmph.R, which you loaded at the beginning of this chapter). See Chapter 5.9.7 for how to specify the second argument to get the estimates of interest.

# Estimate slope at each level of other variable, along with a CI
# using the design DF (the default)
rbind(
  "cBMXWAIST @ Male"   = svycontrast_design_df(fit.ex8.2,
                              c("cBMXWAIST"                = 1)),
  "cBMXWAIST @ Female" = svycontrast_design_df(fit.ex8.2,
                              c("cBMXWAIST"                = 1,
                                "cBMXWAIST:RIAGENDRFemale" = 1))
)
##                        est   lower   upper
## cBMXWAIST @ Male   0.03170 0.02257 0.04084
## cBMXWAIST @ Female 0.02332 0.01597 0.03066

Conclusion: After incorporating the NHANES complex survey design, the association between waist circumference and fasting glucose was not significantly different between males and females (p = 0.065; slope for males = 0.032; slope for females = 0.023).

8.4.4 Visualize the weighted unadjusted relationships

Use svyplot() to visualize the weighted unadjusted relationships between predictors and the outcome using bubble plots, with points representing more individuals in the population plotted with larger circles. This function works for both continuous and categorical predictors (Figures 8.4 and 8.5, respectively). For a continuous predictor, use svysmooth() to superimpose an estimate of the mean outcome vs. the predictor that does not assume linearity.

NOTE: svyplot() requires a formula of the form Y ~ X. This is different from the plot() function, which we used in previous chapters, which accepts a formula or the individual variables in the order X, Y.

# Plot for continuous predictor
svyplot(LBDGLUSI ~ BMXWAIST, design.FST.nomiss,
        ylab = "Fasting Glucose (mmol/L)",
        xlab = "Waist Circumference (cm)")
# Fit unadjusted model
fit.wc <- svyglm(LBDGLUSI ~ BMXWAIST,
                 family=gaussian(), design=design.FST.nomiss)
# Add line
abline(fit.wc, col = "red", lwd = 2)
# Add smoother
lines(svysmooth(fit.wc$formula,
                fit.wc$survey.design,
                df = 5),
      col = "blue", lwd = 2, lty = 2)
legend("topright", c("Regression line", "Smoother"),
       title = "Weighted estimate",
       col = c("red", "blue"), lwd = c(2,2), lty = c(1,2),
       bty = "n", seg.len = 5)
Weighted scatterplot of outcome vs. a continuous predictor

Figure 8.4: Weighted scatterplot of outcome vs. a continuous predictor

# Plot for categorical predictor
svyplot(LBDGLUSI ~ smoker, design.FST.nomiss,
        ylab = "Fasting Glucose (mmol/L)",
        xlab = "Smoking Status",
        xaxt = "n")
axis(1, at = 1:3, labels = levels(nhanes.mod$smoker))
# Fit unadjusted model
fit.smoker <- svyglm(LBDGLUSI ~ smoker,
                 family=gaussian(), design=design.FST.nomiss)
# Compute means
MEANS <- as.numeric(
  predict(fit.smoker,
          data.frame(smoker=levels(nhanes.mod$smoker)))
)
# Add points and lines
points(1:length(MEANS), MEANS, col = "red", pch = 20, cex = 2)
lines( 1:length(MEANS), MEANS, col = "red")
Weighted plot of outcome vs. a categorical predictor

Figure 8.5: Weighted plot of outcome vs. a categorical predictor

References

———. 2017. “Fitting Regression Models to Survey Data.” Statistical Science 32 (2): 265–78. https://doi.org/10.1214/16-STS605.