Soc 712 Homework #5

Demographic Traits and Marjiuana Use

How do independent variables such as race, sex, income and education affect the chance that someone will use marijuana?

This chunk is simply reading in the csv data set into R using the read_csv function. The head function is used to see if the data set was properly imported into R. The names function was used to look at all the names of the variables to better understand the data set and to help pick which variables to observe for this study.

library(readr)
marij <- read_csv("/Users/paulkim/Downloads/balexturner-drug-use-employment-work-absence-income-race-education/data/nsduh_workforce_adults.csv")

Parsed with column specification:
cols(
  .default = col_character(),
  column_a = col_integer(),
  irpinc3 = col_integer(),
  irfamin3 = col_integer(),
  countofdrugs_ever = col_integer(),
  countofdrugs_month = col_integer(),
  countofdrugs_year = col_integer(),
  personalincome = col_integer(),
  familyincome = col_integer(),
  questid2 = col_integer(),
  employmentstatus = col_integer(),
  preemploymentdrugtest = col_integer(),
  randomdrugtest = col_integer(),
  everdrugtest = col_integer(),
  race_num = col_integer(),
  education = col_integer(),
  wouldworkfordrugtester = col_integer(),
  selectiveleave = col_integer(),
  skipsick = col_integer(),
  sex = col_integer()
)
See spec(...) for full column specifications.

head(marij)

names(marij)

 [1] "column_a"               "irpinc3"                "irfamin3"               "marij_ever"             "marij_month"           
 [6] "marij_year"             "cocaine_ever"           "cocaine_month"          "cocaine_year"           "crack_ever"            
[11] "crack_month"            "crack_year"             "heroin_ever"            "heroin_month"           "heroin_year"           
[16] "hallucinogen_ever"      "hallucinogen_month"     "hallucinogen_year"      "inhalant_ever"          "inhalant_month"        
[21] "inhalant_year"          "meth_ever"              "meth_month"             "meth_year"              "painrelieve_ever"      
[26] "painrelieve_month"      "painrelieve_year"       "tranq_ever"             "tranq_month"            "tranq_year"            
[31] "stimulant_ever"         "stimulant_month"        "stimulant_year"         "sedative_ever"          "sedative_month"        
[36] "sedative_year"          "anydrugever"            "pharmamonth"            "illicitmonth"           "illicitmonth_nomj"     
[41] "anydrugmonth"           "anydrugyear"            "anydrugever_nomj"       "anydrugmonth_nomj"      "anydrugyear_nomj"      
[46] "countofdrugs_ever"      "countofdrugs_month"     "countofdrugs_year"      "personalincome"         "familyincome"          
[51] "questid2"               "employmentstatus"       "preemploymentdrugtest"  "randomdrugtest"         "everdrugtest"          
[56] "everdrugtest2"          "race_str"               "race_num"               "education"              "wouldworkfordrugtester"
[61] "selectiveleave"         "skipsick"               "sex"                    "druglist"              

The unique function is used here to look at how the variables for marijuana use ever and sex are coded. The mutate function is used to change the way marijuana use ever, race and education are coded. This mutate function allows for the data to be analyzed more easily.

unique(marij$marij_ever)

[1] "true"  "false"

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

marij2 <- marij%>%
  mutate(marij_ever = ifelse(marij_ever == "true", 1, 0),
         race_str = factor(race_str, levels = c("White", "Hispanic", "Asian", "Black/African American",
                                                "Native American/Alaskan Native", "Hawaiian/Pacific Islander", "Mixed")),
         education = ifelse(education == 1, "<HS",
                             ifelse(education == 2, "High School",
                                    ifelse(education == 3, "Some College/Assoc Degree",
                                           ifelse(education == 4, "College Graduate", NA)))),
         sex = ifelse(sex == 1, "Male","Female"))
unique(marij2$sex)

[1] "Male"   "Female"

Analyzing Data

library(Zelig)

Loading required package: survival

m0 <- glm(marij_ever ~ race_str, family = "binomial", data = marij2)
summary(m0)

Call:
glm(formula = marij_ever ~ race_str, family = "binomial", data = marij2)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.487  -1.345   1.019   1.019   1.616  

Coefficients:
                                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)                             0.38536    0.01466  26.288  < 2e-16 ***
race_strHispanic                       -0.65539    0.03088 -21.223  < 2e-16 ***
race_strAsian                          -1.37590    0.06057 -22.714  < 2e-16 ***
race_strBlack/African American         -0.31602    0.03473  -9.099  < 2e-16 ***
race_strNative American/Alaskan Native  0.31444    0.10096   3.114  0.00184 ** 
race_strHawaiian/Pacific Islander      -0.37309    0.15734  -2.371  0.01773 *  
race_strMixed                           0.31881    0.06599   4.831 1.36e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 44149  on 32039  degrees of freedom
Residual deviance: 43109  on 32033  degrees of freedom
AIC: 43123

Number of Fisher Scoring iterations: 4

The results of the above r chunk are showing that race plays a factor into the likelihood of people using marijuana. For example, Hispancis have a 0.65539 less log odds chance of using marijuana compared to the reference group which are Whites.

m1 <- glm(marij_ever ~ sex, family = "binomial", data = marij2)
summary(m1)

Call:
glm(formula = marij_ever ~ sex, family = "binomial", data = marij2)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.308  -1.208   1.052   1.148   1.148  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.07077    0.01561   4.532 5.84e-06 ***
sexMale      0.23191    0.02249  10.311  < 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: 44149  on 32039  degrees of freedom
Residual deviance: 44043  on 32038  degrees of freedom
AIC: 44047

Number of Fisher Scoring iterations: 3

The results of this r chunk on sex is showing that males have a 0.23191 higher log odds chance of using marijuana than compared to females.

m2 <- glm(marij_ever ~ factor(personalincome), family = "binomial", data = marij2)
summary(m2)

Call:
glm(formula = marij_ever ~ factor(personalincome), family = "binomial", 
    data = marij2)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.306  -1.263   1.054   1.094   1.157  

Coefficients:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)              0.04857    0.02311   2.102 0.035561 *  
factor(personalincome)2  0.11515    0.03428   3.359 0.000782 ***
factor(personalincome)3  0.17965    0.03781   4.751 2.02e-06 ***
factor(personalincome)4  0.15107    0.04102   3.683 0.000230 ***
factor(personalincome)5  0.19083    0.04376   4.361 1.30e-05 ***
factor(personalincome)6  0.24841    0.04047   6.138 8.34e-10 ***
factor(personalincome)7  0.21617    0.04088   5.288 1.24e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 44149  on 32039  degrees of freedom
Residual deviance: 44092  on 32033  degrees of freedom
AIC: 44106

Number of Fisher Scoring iterations: 3

Personal income is represented as follows: 1 = Less than $10,000 (Including Loss) 2 = $10,000 - $19,999 3 = $20,000 - $29,999 4 = $30,000 - $39,999 5 = $40,000 - $49,999 6 = $50,000 - $74,999 7 = $75,000 or more

The results of this analysis are showing that someone who is making $50,000 - $74,999 a year has the highest log odds chance of using marijuana than compared to people who are making less than $10,000 a year, they have a log odds chance of 0.24841.

m2 <- glm(marij_ever ~ personalincome, family = "binomial", data = marij2)
summary(m2)

Call:
glm(formula = marij_ever ~ personalincome, family = "binomial", 
    data = marij2)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.313  -1.250   1.047   1.107   1.137  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    0.058769   0.021709   2.707  0.00679 ** 
personalincome 0.036419   0.005452   6.679  2.4e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 44149  on 32039  degrees of freedom
Residual deviance: 44105  on 32038  degrees of freedom
AIC: 44109

Number of Fisher Scoring iterations: 3

m3 <- glm(marij_ever ~ education, family = "binomial", data = marij2)
summary(m3)

Call:
glm(formula = marij_ever ~ education, family = "binomial", data = marij2)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.313  -1.260   1.048   1.097   1.232  

Coefficients:
                                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)                        -0.12860    0.03349  -3.840 0.000123 ***
educationCollege Graduate           0.27086    0.03955   6.849 7.43e-12 ***
educationHigh School                0.32045    0.04009   7.994 1.31e-15 ***
educationSome College/Assoc Degree  0.44137    0.03862  11.427  < 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: 44149  on 32039  degrees of freedom
Residual deviance: 44013  on 32036  degrees of freedom
AIC: 44021

Number of Fisher Scoring iterations: 3

The results of the above code are showing that people who have some college or have an Associates degree have a 0.44137 log odds chance of using marijuana compared to people with less than a high school education. People with some college or with an Associates degree have the highest log odds chance in this analysis than compared to those who have different education levels.

Adding Multiple Variables to the Analysis

m4 <- glm(marij_ever ~ sex + race_str + education + personalincome, family = "binomial", data = marij2)
summary(m4)

Call:
glm(formula = marij_ever ~ sex + race_str + education + personalincome, 
    family = "binomial", data = marij2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6256  -1.2606   0.9166   1.0711   1.7868  

Coefficients:
                                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            -0.019141   0.042004  -0.456  0.64861    
sexMale                                 0.247366   0.023627  10.470  < 2e-16 ***
race_strHispanic                       -0.617721   0.032098 -19.245  < 2e-16 ***
race_strAsian                          -1.368094   0.060993 -22.431  < 2e-16 ***
race_strBlack/African American         -0.285489   0.035390  -8.067 7.21e-16 ***
race_strNative American/Alaskan Native  0.319887   0.101495   3.152  0.00162 ** 
race_strHawaiian/Pacific Islander      -0.359422   0.158095  -2.273  0.02300 *  
race_strMixed                           0.325981   0.066378   4.911 9.06e-07 ***
educationCollege Graduate               0.145903   0.044787   3.258  0.00112 ** 
educationHigh School                    0.221431   0.041158   5.380 7.45e-08 ***
educationSome College/Assoc Degree      0.334497   0.040400   8.280  < 2e-16 ***
personalincome                          0.017454   0.006459   2.702  0.00689 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 44149  on 32039  degrees of freedom
Residual deviance: 42902  on 32028  degrees of freedom
AIC: 42926

Number of Fisher Scoring iterations: 4

When multiple variables are added together in the analysis, the log odds for the variables are affect. For example, the log odds for people who had some college or an Associates degree decreased from 0.44137 to 0.33076. In this example controlling for race and education the log odds of personal income is 0.017454 which is a decrease from 0.036419 when income was analyzed by itself.

m5 <- glm(marij_ever ~ sex*race_str + education, family = "binomial", data = marij2)
summary(m5)

Call:
glm(formula = marij_ever ~ sex * race_str + education, family = "binomial", 
    data = marij2)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.604  -1.283   0.939   1.060   1.752  

Coefficients:
                                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                                     0.04986    0.04055   1.230 0.218839    
sexMale                                         0.18387    0.02948   6.236 4.48e-10 ***
race_strHispanic                               -0.71632    0.04489 -15.956  < 2e-16 ***
race_strAsian                                  -1.34247    0.08795 -15.264  < 2e-16 ***
race_strBlack/African American                 -0.45236    0.04673  -9.681  < 2e-16 ***
race_strNative American/Alaskan Native          0.25503    0.14108   1.808 0.070658 .  
race_strHawaiian/Pacific Islander              -0.44601    0.22284  -2.001 0.045346 *  
race_strMixed                                   0.32995    0.08988   3.671 0.000242 ***
educationCollege Graduate                       0.19399    0.04185   4.635 3.57e-06 ***
educationHigh School                            0.23371    0.04114   5.681 1.34e-08 ***
educationSome College/Assoc Degree              0.35682    0.04008   8.903  < 2e-16 ***
sexMale:race_strHispanic                        0.18161    0.06214   2.923 0.003472 ** 
sexMale:race_strAsian                          -0.04725    0.12146  -0.389 0.697258    
sexMale:race_strBlack/African American          0.35939    0.07083   5.074 3.89e-07 ***
sexMale:race_strNative American/Alaskan Native  0.11764    0.20278   0.580 0.561814    
sexMale:race_strHawaiian/Pacific Islander       0.15767    0.31618   0.499 0.618011    
sexMale:race_strMixed                          -0.03500    0.13269  -0.264 0.791949    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 44149  on 32039  degrees of freedom
Residual deviance: 42877  on 32023  degrees of freedom
AIC: 42911

Number of Fisher Scoring iterations: 4

There are only two statistically significant findings from this analysis which are for Hispanic men and for Black men. When race was analyzed by itself the log odds for Hispanics and Blacks were negative, -0.71632 and -0.45236 respectively. But when race and sex were interacted the log odds for Hispanic and Black men went up, 0.18161 and 0.35939 respectively. This means that men have a higher log odds chance of using marijuana. This also means that the affect race has on marijuana use affects men and women differently.

m6 <- glm(marij_ever ~ sex*education + race_str, family = "binomial", data = marij2)
summary(m6)

Call:
glm(formula = marij_ever ~ sex * education + race_str, family = "binomial", 
    data = marij2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5930  -1.2421   0.9263   1.0493   1.8242  

Coefficients:
                                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                                -0.092254   0.055424  -1.665  0.09601 .  
sexMale                                     0.449225   0.069492   6.464 1.02e-10 ***
educationCollege Graduate                   0.401246   0.061528   6.521 6.97e-11 ***
educationHigh School                        0.243123   0.062851   3.868  0.00011 ***
educationSome College/Assoc Degree          0.463389   0.059726   7.759 8.59e-15 ***
race_strHispanic                           -0.629108   0.032008 -19.655  < 2e-16 ***
race_strAsian                              -1.361680   0.060959 -22.338  < 2e-16 ***
race_strBlack/African American             -0.293912   0.035221  -8.345  < 2e-16 ***
race_strNative American/Alaskan Native      0.312628   0.101602   3.077  0.00209 ** 
race_strHawaiian/Pacific Islander          -0.366919   0.158367  -2.317  0.02051 *  
race_strMixed                               0.314778   0.066337   4.745 2.08e-06 ***
sexMale:educationCollege Graduate          -0.417327   0.081803  -5.102 3.37e-07 ***
sexMale:educationHigh School               -0.008978   0.082778  -0.108  0.91363    
sexMale:educationSome College/Assoc Degree -0.196378   0.079896  -2.458  0.01397 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 44149  on 32039  degrees of freedom
Residual deviance: 42857  on 32026  degrees of freedom
AIC: 42885

Number of Fisher Scoring iterations: 4

This analysis is showing the interaction between sex and education. All the log odds for education decreased when education and sex interacted together. This means the affect that educational levels have on marijuana use are different for men and women. The men who have graduated college have the greatest log odds chance of using marijuana.

library(texreg)

Version: 1.36.23 Date: 2017-03-03 Author: Philip Leifeld (University of Glasgow)

Please cite the JSS article in your publications – see citation(“texreg”).

htmlreg(list(m4, m5, m6))

Statistical models

	Model 1	Model 2	Model 3
(Intercept)	-0.02	0.05	-0.09
	(0.04)	(0.04)	(0.06)
sexMale	0.25***	0.18***	0.45***
	(0.02)	(0.03)	(0.07)
race_strHispanic	-0.62***	-0.72***	-0.63***
	(0.03)	(0.04)	(0.03)
race_strAsian	-1.37***	-1.34***	-1.36***
	(0.06)	(0.09)	(0.06)
race_strBlack/African American	-0.29***	-0.45***	-0.29***
	(0.04)	(0.05)	(0.04)
race_strNative American/Alaskan Native	0.32**	0.26	0.31**
	(0.10)	(0.14)	(0.10)
race_strHawaiian/Pacific Islander	-0.36*	-0.45*	-0.37*
	(0.16)	(0.22)	(0.16)
race_strMixed	0.33***	0.33***	0.31***
	(0.07)	(0.09)	(0.07)
educationCollege Graduate	0.15**	0.19***	0.40***
	(0.04)	(0.04)	(0.06)
educationHigh School	0.22***	0.23***	0.24***
	(0.04)	(0.04)	(0.06)
educationSome College/Assoc Degree	0.33***	0.36***	0.46***
	(0.04)	(0.04)	(0.06)
personalincome	0.02**
	(0.01)
sexMale:race_strHispanic		0.18**
		(0.06)
sexMale:race_strAsian		-0.05
		(0.12)
sexMale:race_strBlack/African American		0.36***
		(0.07)
sexMale:race_strNative American/Alaskan Native		0.12
		(0.20)
sexMale:race_strHawaiian/Pacific Islander		0.16
		(0.32)
sexMale:race_strMixed		-0.04
		(0.13)
sexMale:educationCollege Graduate			-0.42***
			(0.08)
sexMale:educationHigh School			-0.01
			(0.08)
sexMale:educationSome College/Assoc Degree			-0.20*
			(0.08)
AIC	42925.61	42911.13	42885.48
BIC	43026.10	43053.50	43002.72
Log Likelihood	-21450.80	-21438.56	-21428.74
Deviance	42901.61	42877.13	42857.48
Num. obs.	32040	32040	32040
p < 0.001, p < 0.01, p < 0.05

lmtest::lrtest(m4, m5, m6)

Likelihood ratio test

Model 1: marij_ever ~ sex + race_str + education + personalincome Model 2: marij_ever ~ sex * race_str + education Model 3: marij_ever ~ sex * education + race_str #Df LogLik Df Chisq Pr(>Chisq)
1 12 -21451
2 17 -21439 5 24.480 0.0001755 3 14 -21429 -3 19.652 0.0002004 — Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘’ 1

According to the comparisons of the models, model 3 or variable m6 on the list fits best out of the 3 models. Model 3 is the best fitting model because it has the lowest AIC, the lowest BIC and the lowest log likelihood. This means that model 3 is the most likely to produce the observed data.

library(visreg)
visreg(m5, "sex", by = "race_str", scale = "response")

This figure is displaying the interaction between sex and race. The results of this figure are showing that males in all races have a higher chance of using marijuana than compared to females. This also means that race is affecting marijuana use on people differently for males and females.

visreg(m6, "sex", by = "education", scale = "response")

This figure is displaying the interaction between sex and education. The results here are also displaying that males have a higher chance of using marjiuana than compared to women. But this figure is showing that as level of education increases the difference in log odds chance between males and females decreases. The less than high school education portion of this population is displaying the greatest difference in log odds chance between males and females but the college graduate portion shows a very small difference in log odds between males and females.

Findings

The major finding that was displayed in m5 and m6 was that males had the greater log odds chance of using marijuana than females. When sex and race were interacted together it showed that race affects marijuana use differently for males and females. Males in any race category had the higher chance of using marijuana. One interesting finding that was displayed in m6 was that as education levels increased, the difference between the log odds of males and females decreased. So as male and female education levels increase, the log odds chance of them using marijuana becomes similar.