Homework Assignment 9

1. Define an outcome, this can be of any form, but use an appropriate distribution for your outcome. Also define at least 1 continuous predictor.

I will use data from BRFSS with the outcome variable is general health and the continuous predictor variable will be sleep time.

2. Using the gam() function, estimate a model with only linear terms in your model

library(mgcv)

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## This is mgcv 1.8-39. For overview type 'help("mgcv-package")'.

library(ggplot2)
library(dplyr)

dat <- data.frame(brfss2020)

gamfit <- gam(genhealth ~ sleptim1 + checkup1 + educa, data = dat, 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(gamfit)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## genhealth ~ sleptim1 + checkup1 + educa
## 
## Parametric coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.587463   0.122241  21.167  < 2e-16 ***
## sleptim1          -0.097657   0.003745 -26.076  < 2e-16 ***
## checkup11last5yrs -0.248483   0.017413 -14.270  < 2e-16 ***
## checkup12never    -0.397103   0.069727  -5.695 1.23e-08 ***
## educa2elementary   0.082819   0.124066   0.668    0.504    
## educa3somehs      -0.110699   0.121673  -0.910    0.363    
## educa4hsgrad      -0.486922   0.119734  -4.067 4.77e-05 ***
## educa5somecol     -0.572138   0.119638  -4.782 1.73e-06 ***
## educa6baorhigher  -1.021040   0.119510  -8.544  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =   0.02   Deviance explained = 1.88%
## UBRE = 0.062364  Scale est. = 1         n = 185745

3. Repeat Step 2, but include a smooth term for at least one continuous variable

dat <- data.frame(brfss2020)

gamfitsm <- gam(genhealth ~ s(sleptim1) + checkup1 + educa, data = dat, 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!

summary(gamfitsm)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## genhealth ~ s(sleptim1) + checkup1 + educa
## 
## Parametric coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.87319    0.11923  15.711  < 2e-16 ***
## checkup11last5yrs -0.25422    0.01746 -14.563  < 2e-16 ***
## checkup12never    -0.40985    0.06979  -5.873 4.29e-09 ***
## educa2elementary   0.09944    0.12406   0.802    0.423    
## educa3somehs      -0.11805    0.12168  -0.970    0.332    
## educa4hsgrad      -0.47258    0.11974  -3.947 7.92e-05 ***
## educa5somecol     -0.55438    0.11965  -4.634 3.60e-06 ***
## educa6baorhigher  -0.97341    0.11953  -8.144 3.83e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value    
## s(sleptim1) 7.686   8.27   1355  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0243   Deviance explained = 2.25%
## UBRE = 0.058399  Scale est. = 1         n = 185745

4. Test if the model in step 3 is fitting the data better than the purely linear model.

anova(gamfit, gamfitsm, test="Chisq")

5. Produce a plot of the smooth effect from the model in step 3

plot(gamfitsm)