library(mgcv)
## Warning: package 'mgcv' was built under R version 4.1.3
## Loading required package: nlme
## This is mgcv 1.8-40. For overview type 'help("mgcv-package")'.
library(ggplot2)
## Warning in register(): Can't find generic `scale_type` in package ggplot2 to
## register S3 method.
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.1.6     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.1.1     v forcats 0.5.1
## v purrr   0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::collapse() masks nlme::collapse()
## x dplyr::filter()   masks stats::filter()
## x dplyr::lag()      masks stats::lag()
library(dplyr)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some

Splines Using the BRFSS

For this short report, the 2020 Behavioral Risk Factor Surveillance Survey (BRFSS) is used to compare a generalized linear model (GLM) to generalized additive model (GAM) of access to healthcare by race/ethnicity, educational attainment, insurance, income, and age.

brfss20 = readRDS(url("https://github.com/coreysparks/DEM7283/blob/master/data/brfss20sm.rds?raw=true"))

names(brfss20) = tolower(gsub(pattern = "_",replacement =  "",x =  names(brfss20)))

Recoding the Variables

A binary outcome variable “hlthcare” was defined by dichotomizing “hlthpln1,” a variable from the 2020 BRFSS survey data which measures whether a survey respondent has health care coverage, including health insurance, prepaid plans or government plans.

Does access to healthcare coverage vary on the basis of race/ethnicity or education?

Two categorical predictor variables were recoded and reveled for this analysis: “race_eth” and “educ.”

#Access to healthcare
brfss20$hlthcare = Recode(brfss20$hlthpln1, recodes = "1=1; 2=0; else=NA")

#Race/ethnicity
brfss20$black = Recode(brfss20$racegr3,
                      recodes="2=1; 9=NA; else=0")

brfss20$white = Recode(brfss20$racegr3, 
                      recodes="1=1; 9=NA; else=0")

brfss20$other = Recode(brfss20$racegr3,
                      recodes="3:4=1; 9=NA; else=0")

brfss20$hispanic = Recode(brfss20$racegr3,
                         recodes="5=1; 9=NA; else=0")

brfss20$race_eth = Recode(brfss20$racegr3,
recodes="1='nh white'; 2='nh black'; 3='nh other';4='nh multirace'; 5='hispanic'; else=NA",
as.factor = T)

brfss20$race_eth = relevel(brfss20$race_eth,
                          ref = "nh white")

#Education level
brfss20$educ = Recode(brfss20$educa,
                     recodes="1:2='Less than HS'; 
                     3='Some HS'; 
                     4='HS grad'; 
                     5='Some college'; 
                     6='College grad';9=NA",
                     as.factor=T)

brfss20$educ = fct_relevel(brfss20$educ,'Less than HS','Some HS','HS grad','Some college','College grad' ) 

#Insurance
brfss20$ins = Recode(brfss20$hlthpln1, recodes = "7:9= NA; 1=1; 2=0")

#Income grouping
brfss20$inc = Recode(brfss20$incomg, recodes = "9= NA; 1= '1_less15k'; 2= '2_15-25k'; 3= '3_25-35k'; 4= '4_35-50k'; 5= '5_50kmore'",
                     as.factor=T)

#Age cut into intervals
brfss20$agec = cut(brfss20$age80, 
                  breaks=c(0,24,39,59,79,99))

#Filter cases
brfss20 = brfss20 %>% 
  select(hlthcare, black, hispanic, white, other, race_eth, educ, ins, inc, agec, age80, mmsawt) %>% 
  filter(complete.cases(.))

Linear Regression

brfss_glm = gam(hlthcare ~ age80 + black + white + hispanic + other + educ + inc,
                data = brfss20,
                weights = mmsawt/mean(mmsawt, na.rm=T),
                family = binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
summary(brfss_glm)

Regression Splines

brfss_gam = gam(hlthcare ~ s(age80) + black + white + hispanic + other + educ + inc,
                data = brfss20,
                weights = mmsawt/mean(mmsawt, na.rm =T),
                family = binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
summary(brfss_gam)
plot(brfss_gam)

Using LRT to Determine Difference in Model Fit

According to the results of the Chi-Square test of difference, the GAM demonstrates a better model fit.

anova(brfss_glm,
      brfss_gam,
      test = "Chisq")