Smoking, Drinking & Employment

CDE 498 - Spring 2021

Joshua Geenen

compiled on April 22, 2021

Setting Up

#Library packages
library(rmdformats)
library(dplyr)
library(tidyverse)
library(GGally)
library(leaps)
library(sjPlot)
library(ggplot2)
library(corrplot)
library(RColorBrewer)
library(car)
library(ggsci)
library(jtools)
library(ppcor)
library(rmarkdown)
library(breakDown)
library(janitor)
library(nnet)


#Set working directory
setwd("C:/Users/joshu/Desktop/COLLEGE_GENERAL/ASU_GENERAL/ASU_Spring_2021/CDE498")

#Read in data set
nhis <- read.csv("nhis.csv", head=T)

#Add column for gender
nhis$gender <- ""
nhis$gender <- ifelse(nhis$male == "0", "Female",nhis$gender)
nhis$gender <- ifelse(nhis$male == "1", "Male", nhis$gender)

Summary Statistics

summary(nhis)

##       year         nhispid               age           health         
##  Min.   :2013   Min.   :2.013e+13   Min.   :18.00   Length:137548     
##  1st Qu.:2014   1st Qu.:2.014e+13   1st Qu.:33.00   Class :character  
##  Median :2015   Median :2.015e+13   Median :48.00   Mode  :character  
##  Mean   :2015   Mean   :2.015e+13   Mean   :48.21                     
##  3rd Qu.:2016   3rd Qu.:2.016e+13   3rd Qu.:62.00                     
##  Max.   :2017   Max.   :2.017e+13   Max.   :84.00                     
##    hrsleep           sleepfall            sleepcat           vet         
##  Length:137548      Length:137548      Min.   :0.0000   Min.   :0.00000  
##  Class :character   Class :character   1st Qu.:0.0000   1st Qu.:0.00000  
##  Mode  :character   Mode  :character   Median :0.0000   Median :0.00000  
##                                        Mean   :0.4818   Mean   :0.09551  
##                                        3rd Qu.:1.0000   3rd Qu.:0.00000  
##                                        Max.   :2.0000   Max.   :1.00000  
##       male            fborn        racedummy           married         
##  Min.   :0.0000   Min.   :0.000   Length:137548      Length:137548     
##  1st Qu.:0.0000   1st Qu.:0.000   Class :character   Class :character  
##  Median :0.0000   Median :0.000   Mode  :character   Mode  :character  
##  Mean   :0.4616   Mean   :0.162                                        
##  3rd Qu.:1.0000   3rd Qu.:0.000                                        
##  Max.   :1.0000   Max.   :1.000                                        
##   smokedummy           bmicat            alcohol             employed     
##  Length:137548      Length:137548      Length:137548      Min.   :0.0000  
##  Class :character   Class :character   Class :character   1st Qu.:0.0000  
##  Mode  :character   Mode  :character   Mode  :character   Median :1.0000  
##                                                           Mean   :0.6064  
##                                                           3rd Qu.:1.0000  
##                                                           Max.   :1.0000  
##   educdummy          incomecat               k6         fairpoorhealth  
##  Length:137548      Length:137548      Min.   : 0.000   Min.   :0.0000  
##  Class :character   Class :character   1st Qu.: 0.000   1st Qu.:0.0000  
##  Mode  :character   Mode  :character   Median : 1.000   Median :0.0000  
##                                        Mean   : 2.211   Mean   :0.1398  
##                                        3rd Qu.: 3.000   3rd Qu.:0.0000  
##                                        Max.   :20.000   Max.   :1.0000  
##    homeowner         health2        sleepfall2       gender         
##  Min.   :0.0000   Min.   :1.000   Min.   :0.000   Length:137548     
##  1st Qu.:0.0000   1st Qu.:1.000   1st Qu.:0.000   Class :character  
##  Median :1.0000   Median :2.000   Median :0.000   Mode  :character  
##  Mean   :0.5942   Mean   :2.313   Mean   :1.402                     
##  3rd Qu.:1.0000   3rd Qu.:3.000   3rd Qu.:2.000                     
##  Max.   :1.0000   Max.   :5.000   Max.   :7.000

summary(nhis$age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   33.00   48.00   48.21   62.00   84.00

count_gender <- nhis %>%
  count(male)
count_gender

##   male     n
## 1    0 74055
## 2    1 63493

count_race <- nhis %>%
  count(racedummy)
count_race

##   racedummy     n
## 1  hispanic 20960
## 2  nh black 17719
## 3  nh other  9622
## 4  nh white 89247

count_education <- nhis %>%
  count(educdummy)
count_education

##     educdummy     n
## 1    college+ 42465
## 2          hs 36817
## 3        lths 14918
## 4 somecollege 43348

Distribution Charts

count_smoking1 <- nhis %>%
  count(smokedummy)
count_smoking1

##                smokedummy     n
## 1 current everyday smoker 18394
## 2  current someday smoker  5612
## 3           former smoker 31857
## 4            never smoker 81685

smokingbar3 <- ggplot(count_smoking1, aes(x=reorder(smokedummy, -n), y=n)) + 
  geom_bar(stat = "identity", color = "black", fill = "purple", width = .25) + 
  geom_text(aes(label = n, vjust = -.25)) + 
  labs(x = "Type of Smoker", y = "Count") + 
  ggtitle("Smoking Distribution") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, size = 22, face = "bold"), 
        axis.text = element_text(size = 12), 
        axis.title = element_text(size = 12, face = "bold"))
print(smokingbar3)

ggsave(plot = smokingbar3,"smoking_distribution.png", dpi = 300, width = 15, height = 6)


count_alcohol1 <- nhis %>%
  count(alcohol)
count_alcohol1

##           alcohol     n
## 1       abstainer 25818
## 2 current drinker 91259
## 3  former drinker 20471

alcoholbar <- ggplot(count_alcohol1, aes(x=reorder(alcohol, -n), y=n)) + 
  geom_bar(stat = "identity", color = "black", fill = "orange", width = .25) + 
  geom_text(aes(label = n, vjust = -.25)) + 
  labs(x = "Type of Drinker", y = "Count") + 
  ggtitle("Alcohol Distribution") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, size = 22, face = "bold"), 
        axis.text = element_text(size = 12), 
        axis.title = element_text(size = 12, face = "bold"))
print(alcoholbar)

ggsave(plot = alcoholbar,"alcohol_distribution.png", dpi = 300, width = 15, height = 6)

Logistic Regression for Employment

Smoking

nhis$smokedummy <- as.factor(nhis$smokedummy)

logitsmoking <- glm(nhis$employed ~ relevel(nhis$smokedummy, ref = "never smoker"), 
                    family = binomial("logit"))
summary(logitsmoking)

## 
## Call:
## glm(formula = nhis$employed ~ relevel(nhis$smokedummy, ref = "never smoker"), 
##     family = binomial("logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4437  -1.2747   0.9327   0.9327   1.1354  
## 
## Coefficients:
##                                                                        Estimate
## (Intercept)                                                            0.607101
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker -0.381290
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker  -0.038443
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           -0.507386
##                                                                       Std. Error
## (Intercept)                                                             0.007323
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker   0.016549
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker    0.028733
## relevel(nhis$smokedummy, ref = "never smoker")former smoker             0.013398
##                                                                       z value
## (Intercept)                                                            82.907
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker -23.040
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker   -1.338
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           -37.872
##                                                                       Pr(>|z|)
## (Intercept)                                                             <2e-16
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker   <2e-16
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker     0.181
## relevel(nhis$smokedummy, ref = "never smoker")former smoker             <2e-16
##                                                                          
## (Intercept)                                                           ***
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker ***
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker     
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184406  on 137547  degrees of freedom
## Residual deviance: 182740  on 137544  degrees of freedom
## AIC: 182748
## 
## Number of Fisher Scoring iterations: 4

exp(coef(logitsmoking))

##                                                           (Intercept) 
##                                                             1.8351034 
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker 
##                                                             0.6829796 
##  relevel(nhis$smokedummy, ref = "never smoker")current someday smoker 
##                                                             0.9622861 
##           relevel(nhis$smokedummy, ref = "never smoker")former smoker 
##                                                             0.6020676

nhis$smokedummy <- as.factor(nhis$smokedummy)
nhis$educdummy <- as.factor(nhis$educdummy)
nhis$racedummy <- as.factor(nhis$racedummy)
nhis$gender <- as.factor(nhis$gender)
nhis$health <- as.factor(nhis$health)

logitcontrolledsmoking <- glm(nhis$employed ~ relevel(nhis$smokedummy, ref = "never smoker") + 
                                relevel(nhis$health, ref = "Excellent") + 
                                relevel(nhis$educdummy, ref = "lths") + 
                                relevel(nhis$racedummy, ref = "nh white") + 
                                relevel(nhis$gender, ref = "Female") + 
                                nhis$age, family = binomial("logit"))
summary(logitcontrolledsmoking)

## 
## Call:
## glm(formula = nhis$employed ~ relevel(nhis$smokedummy, ref = "never smoker") + 
##     relevel(nhis$health, ref = "Excellent") + relevel(nhis$educdummy, 
##     ref = "lths") + relevel(nhis$racedummy, ref = "nh white") + 
##     relevel(nhis$gender, ref = "Female") + nhis$age, family = binomial("logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4321  -0.9275   0.5491   0.8570   2.8590  
## 
## Coefficients:
##                                                                         Estimate
## (Intercept)                                                            2.0045887
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker -0.0188019
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker   0.0270432
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           -0.0459380
## relevel(nhis$health, ref = "Excellent")Fair                           -1.0831554
## relevel(nhis$health, ref = "Excellent")Good                           -0.2614911
## relevel(nhis$health, ref = "Excellent")Poor                           -2.4356481
## relevel(nhis$health, ref = "Excellent")Very Good                       0.0270078
## relevel(nhis$educdummy, ref = "lths")college+                          1.1997940
## relevel(nhis$educdummy, ref = "lths")hs                                0.4974686
## relevel(nhis$educdummy, ref = "lths")somecollege                       0.7066896
## relevel(nhis$racedummy, ref = "nh white")hispanic                      0.1143729
## relevel(nhis$racedummy, ref = "nh white")nh black                     -0.1403065
## relevel(nhis$racedummy, ref = "nh white")nh other                     -0.2873158
## relevel(nhis$gender, ref = "Female")Male                               0.5000385
## nhis$age                                                              -0.0449112
##                                                                       Std. Error
## (Intercept)                                                            0.0327208
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker  0.0195836
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker   0.0331921
## relevel(nhis$smokedummy, ref = "never smoker")former smoker            0.0157962
## relevel(nhis$health, ref = "Excellent")Fair                            0.0234461
## relevel(nhis$health, ref = "Excellent")Good                            0.0175732
## relevel(nhis$health, ref = "Excellent")Poor                            0.0526874
## relevel(nhis$health, ref = "Excellent")Very Good                       0.0169299
## relevel(nhis$educdummy, ref = "lths")college+                          0.0244800
## relevel(nhis$educdummy, ref = "lths")hs                                0.0233591
## relevel(nhis$educdummy, ref = "lths")somecollege                       0.0233616
## relevel(nhis$racedummy, ref = "nh white")hispanic                      0.0197382
## relevel(nhis$racedummy, ref = "nh white")nh black                      0.0194250
## relevel(nhis$racedummy, ref = "nh white")nh other                      0.0254998
## relevel(nhis$gender, ref = "Female")Male                               0.0128093
## nhis$age                                                               0.0004068
##                                                                        z value
## (Intercept)                                                             61.263
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker   -0.960
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker     0.815
## relevel(nhis$smokedummy, ref = "never smoker")former smoker             -2.908
## relevel(nhis$health, ref = "Excellent")Fair                            -46.198
## relevel(nhis$health, ref = "Excellent")Good                            -14.880
## relevel(nhis$health, ref = "Excellent")Poor                            -46.228
## relevel(nhis$health, ref = "Excellent")Very Good                         1.595
## relevel(nhis$educdummy, ref = "lths")college+                           49.011
## relevel(nhis$educdummy, ref = "lths")hs                                 21.297
## relevel(nhis$educdummy, ref = "lths")somecollege                        30.250
## relevel(nhis$racedummy, ref = "nh white")hispanic                        5.794
## relevel(nhis$racedummy, ref = "nh white")nh black                       -7.223
## relevel(nhis$racedummy, ref = "nh white")nh other                      -11.267
## relevel(nhis$gender, ref = "Female")Male                                39.037
## nhis$age                                                              -110.406
##                                                                       Pr(>|z|)
## (Intercept)                                                            < 2e-16
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker  0.33701
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker   0.41522
## relevel(nhis$smokedummy, ref = "never smoker")former smoker            0.00364
## relevel(nhis$health, ref = "Excellent")Fair                            < 2e-16
## relevel(nhis$health, ref = "Excellent")Good                            < 2e-16
## relevel(nhis$health, ref = "Excellent")Poor                            < 2e-16
## relevel(nhis$health, ref = "Excellent")Very Good                       0.11065
## relevel(nhis$educdummy, ref = "lths")college+                          < 2e-16
## relevel(nhis$educdummy, ref = "lths")hs                                < 2e-16
## relevel(nhis$educdummy, ref = "lths")somecollege                       < 2e-16
## relevel(nhis$racedummy, ref = "nh white")hispanic                     6.85e-09
## relevel(nhis$racedummy, ref = "nh white")nh black                     5.09e-13
## relevel(nhis$racedummy, ref = "nh white")nh other                      < 2e-16
## relevel(nhis$gender, ref = "Female")Male                               < 2e-16
## nhis$age                                                               < 2e-16
##                                                                          
## (Intercept)                                                           ***
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker    
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker     
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           ** 
## relevel(nhis$health, ref = "Excellent")Fair                           ***
## relevel(nhis$health, ref = "Excellent")Good                           ***
## relevel(nhis$health, ref = "Excellent")Poor                           ***
## relevel(nhis$health, ref = "Excellent")Very Good                         
## relevel(nhis$educdummy, ref = "lths")college+                         ***
## relevel(nhis$educdummy, ref = "lths")hs                               ***
## relevel(nhis$educdummy, ref = "lths")somecollege                      ***
## relevel(nhis$racedummy, ref = "nh white")hispanic                     ***
## relevel(nhis$racedummy, ref = "nh white")nh black                     ***
## relevel(nhis$racedummy, ref = "nh white")nh other                     ***
## relevel(nhis$gender, ref = "Female")Male                              ***
## nhis$age                                                              ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184406  on 137547  degrees of freedom
## Residual deviance: 150776  on 137532  degrees of freedom
## AIC: 150808
## 
## Number of Fisher Scoring iterations: 4

exp(coef(logitcontrolledsmoking))

##                                                           (Intercept) 
##                                                             7.4230402 
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker 
##                                                             0.9813738 
##  relevel(nhis$smokedummy, ref = "never smoker")current someday smoker 
##                                                             1.0274121 
##           relevel(nhis$smokedummy, ref = "never smoker")former smoker 
##                                                             0.9551012 
##                           relevel(nhis$health, ref = "Excellent")Fair 
##                                                             0.3385257 
##                           relevel(nhis$health, ref = "Excellent")Good 
##                                                             0.7699027 
##                           relevel(nhis$health, ref = "Excellent")Poor 
##                                                             0.0875410 
##                      relevel(nhis$health, ref = "Excellent")Very Good 
##                                                             1.0273758 
##                         relevel(nhis$educdummy, ref = "lths")college+ 
##                                                             3.3194332 
##                               relevel(nhis$educdummy, ref = "lths")hs 
##                                                             1.6445530 
##                      relevel(nhis$educdummy, ref = "lths")somecollege 
##                                                             2.0272690 
##                     relevel(nhis$racedummy, ref = "nh white")hispanic 
##                                                             1.1211701 
##                     relevel(nhis$racedummy, ref = "nh white")nh black 
##                                                             0.8690918 
##                     relevel(nhis$racedummy, ref = "nh white")nh other 
##                                                             0.7502747 
##                              relevel(nhis$gender, ref = "Female")Male 
##                                                             1.6487847 
##                                                              nhis$age 
##                                                             0.9560824

Alcohol

nhis$alcohol <- as.factor(nhis$alcohol)

logitalcohol <- glm(nhis$employed ~ relevel(nhis$alcohol, ref = "abstainer"), 
                    family = binomial("logit"))
summary(logitalcohol)

## 
## Call:
## glm(formula = nhis$employed ~ relevel(nhis$alcohol, ref = "abstainer"), 
##     family = binomial("logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5180  -1.1393   0.8715   0.8715   1.3119  
## 
## Coefficients:
##                                                         Estimate Std. Error
## (Intercept)                                             -0.09039    0.01246
## relevel(nhis$alcohol, ref = "abstainer")current drinker  0.86280    0.01435
## relevel(nhis$alcohol, ref = "abstainer")former drinker  -0.22033    0.01885
##                                                         z value Pr(>|z|)    
## (Intercept)                                              -7.254 4.04e-13 ***
## relevel(nhis$alcohol, ref = "abstainer")current drinker  60.122  < 2e-16 ***
## relevel(nhis$alcohol, ref = "abstainer")former drinker  -11.687  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184406  on 137547  degrees of freedom
## Residual deviance: 177481  on 137545  degrees of freedom
## AIC: 177487
## 
## Number of Fisher Scoring iterations: 4

exp(coef(logitalcohol))

##                                             (Intercept) 
##                                               0.9135784 
## relevel(nhis$alcohol, ref = "abstainer")current drinker 
##                                               2.3697789 
##  relevel(nhis$alcohol, ref = "abstainer")former drinker 
##                                               0.8022534

nhis$alcohol <- as.factor(nhis$alcohol)
nhis$educdummy <- as.factor(nhis$educdummy)
nhis$racedummy <- as.factor(nhis$racedummy)
nhis$gender <- as.factor(nhis$gender)
nhis$health <- as.factor(nhis$health)

logitcontrolledalcohol <- glm(nhis$employed ~ relevel(nhis$alcohol, ref = "abstainer") + 
                                relevel(nhis$health, ref = "Excellent") + 
                                relevel(nhis$educdummy, ref = "lths") + 
                                relevel(nhis$racedummy, ref = "nh white") + 
                                relevel(nhis$gender, ref = "Female") + 
                                nhis$age, family = binomial("logit"))
summary(logitcontrolledalcohol)

## 
## Call:
## glm(formula = nhis$employed ~ relevel(nhis$alcohol, ref = "abstainer") + 
##     relevel(nhis$health, ref = "Excellent") + relevel(nhis$educdummy, 
##     ref = "lths") + relevel(nhis$racedummy, ref = "nh white") + 
##     relevel(nhis$gender, ref = "Female") + nhis$age, family = binomial("logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4722  -0.9133   0.5350   0.8423   2.8811  
## 
## Coefficients:
##                                                           Estimate Std. Error
## (Intercept)                                              1.5531309  0.0339174
## relevel(nhis$alcohol, ref = "abstainer")current drinker  0.6740393  0.0170427
## relevel(nhis$alcohol, ref = "abstainer")former drinker   0.2342090  0.0220920
## relevel(nhis$health, ref = "Excellent")Fair             -1.0607259  0.0235054
## relevel(nhis$health, ref = "Excellent")Good             -0.2611929  0.0176134
## relevel(nhis$health, ref = "Excellent")Poor             -2.3917962  0.0528091
## relevel(nhis$health, ref = "Excellent")Very Good         0.0089289  0.0170135
## relevel(nhis$educdummy, ref = "lths")college+            1.0726183  0.0245737
## relevel(nhis$educdummy, ref = "lths")hs                  0.4456939  0.0235625
## relevel(nhis$educdummy, ref = "lths")somecollege         0.6086513  0.0235693
## relevel(nhis$racedummy, ref = "nh white")hispanic        0.1956456  0.0197311
## relevel(nhis$racedummy, ref = "nh white")nh black       -0.0576845  0.0195600
## relevel(nhis$racedummy, ref = "nh white")nh other       -0.1369143  0.0259415
## relevel(nhis$gender, ref = "Female")Male                 0.4389568  0.0128870
## nhis$age                                                -0.0441247  0.0004021
##                                                          z value Pr(>|z|)    
## (Intercept)                                               45.792  < 2e-16 ***
## relevel(nhis$alcohol, ref = "abstainer")current drinker   39.550  < 2e-16 ***
## relevel(nhis$alcohol, ref = "abstainer")former drinker    10.602  < 2e-16 ***
## relevel(nhis$health, ref = "Excellent")Fair              -45.127  < 2e-16 ***
## relevel(nhis$health, ref = "Excellent")Good              -14.829  < 2e-16 ***
## relevel(nhis$health, ref = "Excellent")Poor              -45.291  < 2e-16 ***
## relevel(nhis$health, ref = "Excellent")Very Good           0.525  0.59971    
## relevel(nhis$educdummy, ref = "lths")college+             43.649  < 2e-16 ***
## relevel(nhis$educdummy, ref = "lths")hs                   18.915  < 2e-16 ***
## relevel(nhis$educdummy, ref = "lths")somecollege          25.824  < 2e-16 ***
## relevel(nhis$racedummy, ref = "nh white")hispanic          9.916  < 2e-16 ***
## relevel(nhis$racedummy, ref = "nh white")nh black         -2.949  0.00319 ** 
## relevel(nhis$racedummy, ref = "nh white")nh other         -5.278 1.31e-07 ***
## relevel(nhis$gender, ref = "Female")Male                  34.062  < 2e-16 ***
## nhis$age                                                -109.734  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184406  on 137547  degrees of freedom
## Residual deviance: 148962  on 137533  degrees of freedom
## AIC: 148992
## 
## Number of Fisher Scoring iterations: 4

exp(coef(logitcontrolledalcohol))

##                                             (Intercept) 
##                                              4.72624446 
## relevel(nhis$alcohol, ref = "abstainer")current drinker 
##                                              1.96214709 
##  relevel(nhis$alcohol, ref = "abstainer")former drinker 
##                                              1.26390864 
##             relevel(nhis$health, ref = "Excellent")Fair 
##                                              0.34620440 
##             relevel(nhis$health, ref = "Excellent")Good 
##                                              0.77013231 
##             relevel(nhis$health, ref = "Excellent")Poor 
##                                              0.09146524 
##        relevel(nhis$health, ref = "Excellent")Very Good 
##                                              1.00896884 
##           relevel(nhis$educdummy, ref = "lths")college+ 
##                                              2.92302284 
##                 relevel(nhis$educdummy, ref = "lths")hs 
##                                              1.56157335 
##        relevel(nhis$educdummy, ref = "lths")somecollege 
##                                              1.83795093 
##       relevel(nhis$racedummy, ref = "nh white")hispanic 
##                                              1.21609584 
##       relevel(nhis$racedummy, ref = "nh white")nh black 
##                                              0.94394771 
##       relevel(nhis$racedummy, ref = "nh white")nh other 
##                                              0.87204498 
##                relevel(nhis$gender, ref = "Female")Male 
##                                              1.55108827 
##                                                nhis$age 
##                                              0.95683465

Smoking & Alcohol

nhis$smokedummy <- as.factor(nhis$smokedummy)
nhis$alcohol <- as.factor(nhis$alcohol)

logitsmokingalcohol <- glm(nhis$employed ~ 
                      relevel(nhis$smokedummy, ref = "never smoker") + 
                      relevel(nhis$alcohol, ref = "abstainer"), 
                    family = binomial("logit"))
summary(logitsmokingalcohol)

## 
## Call:
## glm(formula = nhis$employed ~ relevel(nhis$smokedummy, ref = "never smoker") + 
##     relevel(nhis$alcohol, ref = "abstainer"), family = binomial("logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6330  -1.1748   0.7823   1.0123   1.4547  
## 
## Coefficients:
##                                                                        Estimate
## (Intercept)                                                           -0.006101
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker -0.513619
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker  -0.243428
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           -0.625555
## relevel(nhis$alcohol, ref = "abstainer")current drinker                1.033387
## relevel(nhis$alcohol, ref = "abstainer")former drinker                 0.008729
##                                                                       Std. Error
## (Intercept)                                                             0.012650
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker   0.017323
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker    0.029564
## relevel(nhis$smokedummy, ref = "never smoker")former smoker             0.014255
## relevel(nhis$alcohol, ref = "abstainer")current drinker                 0.015011
## relevel(nhis$alcohol, ref = "abstainer")former drinker                  0.019603
##                                                                       z value
## (Intercept)                                                            -0.482
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker -29.649
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker   -8.234
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           -43.884
## relevel(nhis$alcohol, ref = "abstainer")current drinker                68.842
## relevel(nhis$alcohol, ref = "abstainer")former drinker                  0.445
##                                                                       Pr(>|z|)
## (Intercept)                                                              0.630
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker   <2e-16
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker    <2e-16
## relevel(nhis$smokedummy, ref = "never smoker")former smoker             <2e-16
## relevel(nhis$alcohol, ref = "abstainer")current drinker                 <2e-16
## relevel(nhis$alcohol, ref = "abstainer")former drinker                   0.656
##                                                                          
## (Intercept)                                                              
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker ***
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker  ***
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           ***
## relevel(nhis$alcohol, ref = "abstainer")current drinker               ***
## relevel(nhis$alcohol, ref = "abstainer")former drinker                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184406  on 137547  degrees of freedom
## Residual deviance: 175212  on 137542  degrees of freedom
## AIC: 175224
## 
## Number of Fisher Scoring iterations: 4

exp(coef(logitsmokingalcohol))

##                                                           (Intercept) 
##                                                             0.9939180 
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker 
##                                                             0.5983262 
##  relevel(nhis$smokedummy, ref = "never smoker")current someday smoker 
##                                                             0.7839358 
##           relevel(nhis$smokedummy, ref = "never smoker")former smoker 
##                                                             0.5349643 
##               relevel(nhis$alcohol, ref = "abstainer")current drinker 
##                                                             2.8105677 
##                relevel(nhis$alcohol, ref = "abstainer")former drinker 
##                                                             1.0087671

nhis$smokedummy <- as.factor(nhis$smokedummy)
nhis$alcohol <- as.factor(nhis$alcohol)
nhis$educdummy <- as.factor(nhis$educdummy)
nhis$racedummy <- as.factor(nhis$racedummy)
nhis$gender <- as.factor(nhis$gender)
nhis$health <- as.factor(nhis$health)

logitcontrolledsmokingalcohol <- glm(nhis$employed ~ relevel(nhis$smokedummy, ref = "never smoker") + 
                                relevel(nhis$alcohol, ref = "abstainer") + 
                                relevel(nhis$health, ref = "Excellent") + 
                                relevel(nhis$educdummy, ref = "lths") + 
                                relevel(nhis$racedummy, ref = "nh white") + 
                                relevel(nhis$gender, ref = "Female") + 
                                nhis$age, family = binomial("logit"))
summary(logitcontrolledsmokingalcohol)

## 
## Call:
## glm(formula = nhis$employed ~ relevel(nhis$smokedummy, ref = "never smoker") + 
##     relevel(nhis$alcohol, ref = "abstainer") + relevel(nhis$health, 
##     ref = "Excellent") + relevel(nhis$educdummy, ref = "lths") + 
##     relevel(nhis$racedummy, ref = "nh white") + relevel(nhis$gender, 
##     ref = "Female") + nhis$age, family = binomial("logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4786  -0.9109   0.5331   0.8418   2.8898  
## 
## Coefficients:
##                                                                         Estimate
## (Intercept)                                                            1.5645123
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker -0.1467149
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker  -0.1260925
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           -0.1533034
## relevel(nhis$alcohol, ref = "abstainer")current drinker                0.7182075
## relevel(nhis$alcohol, ref = "abstainer")former drinker                 0.2820788
## relevel(nhis$health, ref = "Excellent")Fair                           -1.0363449
## relevel(nhis$health, ref = "Excellent")Good                           -0.2461370
## relevel(nhis$health, ref = "Excellent")Poor                           -2.3592105
## relevel(nhis$health, ref = "Excellent")Very Good                       0.0156159
## relevel(nhis$educdummy, ref = "lths")college+                          1.0386286
## relevel(nhis$educdummy, ref = "lths")hs                                0.4396209
## relevel(nhis$educdummy, ref = "lths")somecollege                       0.5960784
## relevel(nhis$racedummy, ref = "nh white")hispanic                      0.1686115
## relevel(nhis$racedummy, ref = "nh white")nh black                     -0.0730625
## relevel(nhis$racedummy, ref = "nh white")nh other                     -0.1486506
## relevel(nhis$gender, ref = "Female")Male                               0.4513928
## nhis$age                                                              -0.0436754
##                                                                       Std. Error
## (Intercept)                                                            0.0342975
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker  0.0200312
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker   0.0335608
## relevel(nhis$smokedummy, ref = "never smoker")former smoker            0.0162846
## relevel(nhis$alcohol, ref = "abstainer")current drinker                0.0175654
## relevel(nhis$alcohol, ref = "abstainer")former drinker                 0.0225456
## relevel(nhis$health, ref = "Excellent")Fair                            0.0237086
## relevel(nhis$health, ref = "Excellent")Good                            0.0177212
## relevel(nhis$health, ref = "Excellent")Poor                            0.0529929
## relevel(nhis$health, ref = "Excellent")Very Good                       0.0170443
## relevel(nhis$educdummy, ref = "lths")college+                          0.0249058
## relevel(nhis$educdummy, ref = "lths")hs                                0.0235725
## relevel(nhis$educdummy, ref = "lths")somecollege                       0.0236574
## relevel(nhis$racedummy, ref = "nh white")hispanic                      0.0199557
## relevel(nhis$racedummy, ref = "nh white")nh black                      0.0196475
## relevel(nhis$racedummy, ref = "nh white")nh other                      0.0259709
## relevel(nhis$gender, ref = "Female")Male                               0.0129533
## nhis$age                                                               0.0004119
##                                                                        z value
## (Intercept)                                                             45.616
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker   -7.324
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker    -3.757
## relevel(nhis$smokedummy, ref = "never smoker")former smoker             -9.414
## relevel(nhis$alcohol, ref = "abstainer")current drinker                 40.888
## relevel(nhis$alcohol, ref = "abstainer")former drinker                  12.511
## relevel(nhis$health, ref = "Excellent")Fair                            -43.712
## relevel(nhis$health, ref = "Excellent")Good                            -13.889
## relevel(nhis$health, ref = "Excellent")Poor                            -44.519
## relevel(nhis$health, ref = "Excellent")Very Good                         0.916
## relevel(nhis$educdummy, ref = "lths")college+                           41.702
## relevel(nhis$educdummy, ref = "lths")hs                                 18.650
## relevel(nhis$educdummy, ref = "lths")somecollege                        25.196
## relevel(nhis$racedummy, ref = "nh white")hispanic                        8.449
## relevel(nhis$racedummy, ref = "nh white")nh black                       -3.719
## relevel(nhis$racedummy, ref = "nh white")nh other                       -5.724
## relevel(nhis$gender, ref = "Female")Male                                34.848
## nhis$age                                                              -106.028
##                                                                       Pr(>|z|)
## (Intercept)                                                            < 2e-16
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker 2.40e-13
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker  0.000172
## relevel(nhis$smokedummy, ref = "never smoker")former smoker            < 2e-16
## relevel(nhis$alcohol, ref = "abstainer")current drinker                < 2e-16
## relevel(nhis$alcohol, ref = "abstainer")former drinker                 < 2e-16
## relevel(nhis$health, ref = "Excellent")Fair                            < 2e-16
## relevel(nhis$health, ref = "Excellent")Good                            < 2e-16
## relevel(nhis$health, ref = "Excellent")Poor                            < 2e-16
## relevel(nhis$health, ref = "Excellent")Very Good                      0.359566
## relevel(nhis$educdummy, ref = "lths")college+                          < 2e-16
## relevel(nhis$educdummy, ref = "lths")hs                                < 2e-16
## relevel(nhis$educdummy, ref = "lths")somecollege                       < 2e-16
## relevel(nhis$racedummy, ref = "nh white")hispanic                      < 2e-16
## relevel(nhis$racedummy, ref = "nh white")nh black                     0.000200
## relevel(nhis$racedummy, ref = "nh white")nh other                     1.04e-08
## relevel(nhis$gender, ref = "Female")Male                               < 2e-16
## nhis$age                                                               < 2e-16
##                                                                          
## (Intercept)                                                           ***
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker ***
## relevel(nhis$smokedummy, ref = "never smoker")current someday smoker  ***
## relevel(nhis$smokedummy, ref = "never smoker")former smoker           ***
## relevel(nhis$alcohol, ref = "abstainer")current drinker               ***
## relevel(nhis$alcohol, ref = "abstainer")former drinker                ***
## relevel(nhis$health, ref = "Excellent")Fair                           ***
## relevel(nhis$health, ref = "Excellent")Good                           ***
## relevel(nhis$health, ref = "Excellent")Poor                           ***
## relevel(nhis$health, ref = "Excellent")Very Good                         
## relevel(nhis$educdummy, ref = "lths")college+                         ***
## relevel(nhis$educdummy, ref = "lths")hs                               ***
## relevel(nhis$educdummy, ref = "lths")somecollege                      ***
## relevel(nhis$racedummy, ref = "nh white")hispanic                     ***
## relevel(nhis$racedummy, ref = "nh white")nh black                     ***
## relevel(nhis$racedummy, ref = "nh white")nh other                     ***
## relevel(nhis$gender, ref = "Female")Male                              ***
## nhis$age                                                              ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184406  on 137547  degrees of freedom
## Residual deviance: 148847  on 137530  degrees of freedom
## AIC: 148883
## 
## Number of Fisher Scoring iterations: 4

exp(coef(logitcontrolledsmokingalcohol))

##                                                           (Intercept) 
##                                                             4.7803432 
## relevel(nhis$smokedummy, ref = "never smoker")current everyday smoker 
##                                                             0.8635401 
##  relevel(nhis$smokedummy, ref = "never smoker")current someday smoker 
##                                                             0.8815333 
##           relevel(nhis$smokedummy, ref = "never smoker")former smoker 
##                                                             0.8578694 
##               relevel(nhis$alcohol, ref = "abstainer")current drinker 
##                                                             2.0507540 
##                relevel(nhis$alcohol, ref = "abstainer")former drinker 
##                                                             1.3258832 
##                           relevel(nhis$health, ref = "Excellent")Fair 
##                                                             0.3547489 
##                           relevel(nhis$health, ref = "Excellent")Good 
##                                                             0.7818151 
##                           relevel(nhis$health, ref = "Excellent")Poor 
##                                                             0.0944948 
##                      relevel(nhis$health, ref = "Excellent")Very Good 
##                                                             1.0157384 
##                         relevel(nhis$educdummy, ref = "lths")college+ 
##                                                             2.8253396 
##                               relevel(nhis$educdummy, ref = "lths")hs 
##                                                             1.5521187 
##                      relevel(nhis$educdummy, ref = "lths")somecollege 
##                                                             1.8149871 
##                     relevel(nhis$racedummy, ref = "nh white")hispanic 
##                                                             1.1836602 
##                     relevel(nhis$racedummy, ref = "nh white")nh black 
##                                                             0.9295427 
##                     relevel(nhis$racedummy, ref = "nh white")nh other 
##                                                             0.8618702 
##                              relevel(nhis$gender, ref = "Female")Male 
##                                                             1.5704980 
##                                                              nhis$age 
##                                                             0.9572646

Multinomial Regression for Income

Smoking

nhis$smokedummy <- as.factor(nhis$smokedummy)

multismoking <- multinom(incomecat ~ relevel(smokedummy, ref = "never smoker"), data=nhis)

## # weights:  15 (8 variable)
## initial  value 151111.923082 
## iter  10 value 148391.570903
## final  value 148384.079815 
## converged

summary(multismoking)

## Call:
## multinom(formula = incomecat ~ relevel(smokedummy, ref = "never smoker"), 
##     data = nhis)
## 
## Coefficients:
##                 (Intercept)
## $35,000-$74,999 -0.17535207
## $75k++          -0.02497886
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                       -0.4690947
## $75k++                                                                -1.2032323
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                      -0.3231668
## $75k++                                                               -0.8045236
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                             0.12277374
## $75k++                                                     -0.02692905
## 
## Std. Errors:
##                 (Intercept)
## $35,000-$74,999 0.008689925
## $75k++          0.008353584
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                       0.01904185
## $75k++                                                                0.02249670
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                      0.03223660
## $75k++                                                               0.03561853
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                             0.01619804
## $75k++                                                      0.01601811
## 
## Residual Deviance: 296768.2 
## AIC: 296784.2

exp(coef(multismoking))

##                 (Intercept)
## $35,000-$74,999   0.8391615
## $75k++            0.9753305
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                        0.6255684
## $75k++                                                                 0.3002222
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                       0.7238531
## $75k++                                                                0.4473010
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                              1.1306286
## $75k++                                                       0.9734303

zmultismoking <- summary(multismoking)$coefficients/summary(multismoking)$standard.errors
pmultismoking <- (1 - pnorm(abs(zmultismoking), 0, 1)) * 2 
pmultismoking

##                 (Intercept)
## $35,000-$74,999 0.000000000
## $75k++          0.002787982
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                                0
## $75k++                                                                         0
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                               0
## $75k++                                                                        0
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                           3.463896e-14
## $75k++                                                    9.273130e-02

Alcohol

nhis$alcohol <- as.factor(nhis$alcohol)

multialcohol <- multinom(incomecat ~ relevel(alcohol, ref = "abstainer"), data=nhis)

## # weights:  12 (6 variable)
## initial  value 151111.923082 
## iter  10 value 147442.344411
## final  value 147442.339961 
## converged

summary(multialcohol)

## Call:
## multinom(formula = incomecat ~ relevel(alcohol, ref = "abstainer"), 
##     data = nhis)
## 
## Coefficients:
##                 (Intercept) relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999  -0.6330296                                          0.6065967
## $75k++           -0.9372108                                          1.0738120
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                        0.12660287
## $75k++                                                 0.02922135
## 
## Std. Errors:
##                 (Intercept) relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999  0.01465345                                         0.01685297
## $75k++           0.01626608                                         0.01812740
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                        0.02180176
## $75k++                                                 0.02460751
## 
## Residual Deviance: 294884.7 
## AIC: 294896.7

exp(coef(multialcohol))

##                 (Intercept) relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999   0.5309807                                           1.834179
## $75k++            0.3917189                                           2.926514
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                          1.134966
## $75k++                                                   1.029652

zmultialcohol <- summary(multialcohol)$coefficients/summary(multialcohol)$standard.errors
pmultialcohol <- (1 - pnorm(abs(zmultialcohol), 0, 1)) * 2 
pmultialcohol

##                 (Intercept) relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999           0                                                  0
## $75k++                    0                                                  0
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                      6.360152e-09
## $75k++                                               2.350316e-01

Smoking & Alcohol

nhis$smokedummy <- as.factor(nhis$smokedummy)
nhis$alcohol <- as.factor(nhis$alcohol)

multismokingalcohol <- multinom(incomecat ~ 
                                  relevel(smokedummy, ref = "never smoker") + 
                                  relevel(alcohol, ref = "abstainer"), data=nhis)

## # weights:  21 (12 variable)
## initial  value 151111.923082 
## iter  10 value 145820.960011
## final  value 145003.704163 
## converged

summary(multismokingalcohol)

## Call:
## multinom(formula = incomecat ~ relevel(smokedummy, ref = "never smoker") + 
##     relevel(alcohol, ref = "abstainer"), data = nhis)
## 
## Coefficients:
##                 (Intercept)
## $35,000-$74,999  -0.5906729
## $75k++           -0.8467637
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                       -0.5917017
## $75k++                                                                -1.3936411
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                      -0.4852569
## $75k++                                                               -1.0649401
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                             0.03097213
## $75k++                                                     -0.15874441
##                 relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999                                          0.6942142
## $75k++                                                   1.2712391
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                         0.1991360
## $75k++                                                  0.2254263
## 
## Std. Errors:
##                 (Intercept)
## $35,000-$74,999  0.01485191
## $75k++           0.01645518
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                       0.01957194
## $75k++                                                                0.02312318
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                      0.03268264
## $75k++                                                               0.03628516
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                             0.01685923
## $75k++                                                      0.01683641
##                 relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999                                         0.01750936
## $75k++                                                  0.01880875
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                        0.02267515
## $75k++                                                 0.02551349
## 
## Residual Deviance: 290007.4 
## AIC: 290031.4

exp(coef(multismokingalcohol))

##                 (Intercept)
## $35,000-$74,999   0.5539544
## $75k++            0.4288004
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                        0.5533848
## $75k++                                                                 0.2481700
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                       0.6155391
## $75k++                                                                0.3447485
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                              1.0314568
## $75k++                                                       0.8532144
##                 relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999                                           2.002135
## $75k++                                                    3.565268
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                          1.220348
## $75k++                                                   1.252857

zmultismokingalcohol <- summary(multismokingalcohol)$coefficients/summary(multismokingalcohol)$standard.errors
pmultismokingalcohol <- (1 - pnorm(abs(zmultismokingalcohol), 0, 1)) * 2 
pmultismokingalcohol

##                 (Intercept)
## $35,000-$74,999           0
## $75k++                    0
##                 relevel(smokedummy, ref = "never smoker")current everyday smoker
## $35,000-$74,999                                                                0
## $75k++                                                                         0
##                 relevel(smokedummy, ref = "never smoker")current someday smoker
## $35,000-$74,999                                                               0
## $75k++                                                                        0
##                 relevel(smokedummy, ref = "never smoker")former smoker
## $35,000-$74,999                                             0.06619485
## $75k++                                                      0.00000000
##                 relevel(alcohol, ref = "abstainer")current drinker
## $35,000-$74,999                                                  0
## $75k++                                                           0
##                 relevel(alcohol, ref = "abstainer")former drinker
## $35,000-$74,999                                                 0
## $75k++                                                          0