Case Study 2 comprises of analysis of 2018 data, of adolescents age 12-17 years old on heavy alcohol intake and the characteristics of alcohol use. We will be using the article Characteristics of drinking events associated with heavy episodic drinking among adolescents in the United States (Rossheim et al., 2017) as reference.

Load the packages and dataset

library(broom)
library(knitr)
library(MASS)
library(lmtest)
library(survey)
library(tidyverse)
load("NSDUH_2018.RData")
NSDUH2018 <- PUF2018_100819

Creating HED, Binge Variables and Modifying Age Variable

1. Creating a new HED variable for males and females

The dataset will be renamed to “NSDUH2018” which will be used in further coding. The table function is used here to observe the ages included in the dataset.

To create a new variables for heavy episodic of drinking (HED) in adolescents, we will use this following codes:

#find y-variable, heavy episodic drinking, cadrlast
#need to code so that females > or = 4 and males > or = 5 is 1 
#filter data for males and females
male <- filter(NSDUH2018, irsex == 1)
female <- filter(NSDUH2018, irsex == 2)
#create variable for HED and males 
male$HED[male$cadrlast>= 5] <- "1"
male$HED[male$cadrlast< 5] <- "0"
#create variable for HED and females
female$HED[female$cadrlast>= 4] <- "1"
female$HED[female$cadrlast< 4] <- "0"
#bind the two data frames
NSDUH2018 <- rbind(male, female)
#new dichotomous variable is HED 


# Another way to code the dichotomous variable for HED
NSDUH2018 <- mutate(NSDUH2018, HED = ifelse(irsex == 1 & cadrlast >= 5 & cadrlast <= 90 | irsex == 2 & cadrlast >= 4 & cadrlast <= 90, 1, 0))
table(NSDUH2018$HED)
## 
##     0     1 
## 51136  5177

2. Binge variable (dichotomous)

NSDUH2018 <- mutate(NSDUH2018, binge = ifelse(irsex == 1 & cadrlast >= 10 & cadrlast <= 90
 | irsex == 2 & cadrlast >= 8 & cadrlast <= 90, 1, 0))

3. Creating a new variable for age to add age 11 for our further analyses

NSDUH2018 <- mutate(NSDUH2018, newage = AGE2 + 11)

Sources of Alcohol Variables

### 1. Brought the alcohol by themselves at an off-premise retailer. Creating a new variable called *alcselfofprem* for adolescents age 12-17 years who bought alcohol by themselves from off-premise retailers:
NSDUH2018 <- mutate(NSDUH2018, alcselfofprem = ifelse(cagvmony == 1 & cabplace == 1 | CAFRESP2 == 13 & cabuyfre == 1 & CAFRESP2 != 6 & CAFRESP2 !=5 & CAFRESP2 != 4 & CAFRESP2 != 3 & cabuyfre != 2 & cagvmony != 2 & cabplace  != 2 & cabuywho != 1 & cagvwho != 1 &  cagvwho != 2 &   cagvwho != 3 &   cagvwho != 4 & cafrewho != 1 & cafrewho != 2 & cafrewho != 3 & cafrewho != 4 & cafrewho != 5 & cafrewho != 6 & cafrewho != 7, 1, 0))



### 2. Creating a new variable called *alcselfonprem* for adolescents age 12-17 years who bought alcohol by themselves from on-premise alcohol outlets:


NSDUH2018 <- mutate(NSDUH2018, alcselfonprem = ifelse(cagvmony == 1 | cabuyfre == 1 & cabuyfre != 2 & cagvmony != 2 & cabplace == 2 & cabplace != 1 & cagvwho != 1 &  cagvwho != 2 &   cagvwho != 3 &   cagvwho != 4 & cafrewho != 1 & cafrewho != 2 & cafrewho != 3 & cafrewho != 4 & cafrewho != 5 & cafrewho != 6 & cafrewho != 7, 1, 0))




### 3. Bought the alcohol from another person *alcelse*

NSDUH2018 <- mutate(NSDUH2018, alcelse = ifelse(cabuywho == 2 | cagvmony == 1 & cabuywho != 1 & cagvmony != 2 & cagvwho != 1 &  cagvwho != 2 &   cagvwho != 3 &   cagvwho != 4 & cafrewho != 1 & cafrewho != 2 & cafrewho != 3 & cafrewho != 4 & cafrewho != 5 & cafrewho != 6 & cafrewho != 7,  1, 0))


### 4. Was given the alcohol by their parent/guardian *alcparent*

NSDUH2018 <- mutate(NSDUH2018, alcparent = ifelse(cafrewho == 1 | cagvwho == 1 & cafrewho != 2 | cagvwho != 2 & cafrewho != 3 & cafrewho != 4 & cafrewho != 5 & cafrewho != 6 & cafrewho != 7, 1, 0))



### 5. Was given the alcohol by a family member 21 years or older *alcfam21*

NSDUH2018 <- mutate(NSDUH2018, alcfam21 = ifelse(cafrewho == 2 | cagvwho == 2 &  cafrewho != 1 | cagvwho != 1 & cafrewho != 3 & cafrewho != 4 & cafrewho != 5 & cafrewho != 6 & cafrewho != 7, 1, 0))



### 6. Was given the alcohol by someone 21 years or older who was not family member 
NSDUH2018 <- mutate(NSDUH2018, alcelse21 = ifelse(cagvwho == 3 | cafrewho == 3 & cabuyfre == 2 | CAFRESP2 == 3 & cabuyfre != 1 & cagvwho != 2 & cagvwho != 1 & cagvwho != 4 & cafrewho !=  5 & cafrewho != 6 & cafrewho != 7 & cagvmony != 1 & cabuywho != 1, 1, 0))



### 7. Was given the alcohol by someone under 21 years of age *alcunder21*

NSDUH2018 <- mutate(NSDUH2018, alcunder21 = ifelse(cagvwho == 4 | cafrewho == 4 | CAFRESP2 == 4 & cabuyfre == 2 & cabuyfre != 1 & cagvwho != 2 & cagvwho != 1 & cagvwho != 3 & cafrewho !=  5 & cafrewho != 6 & cafrewho != 7, 1, 0))



### 8. Took alcohol from their home *alchome*

NSDUH2018 <- mutate(NSDUH2018, alchome = ifelse(cafrewho == 5 | CAFRESP2 == 5 & cabuyfre == 2 & cabuyfre != 1 & cagvmony != 1 & cagvmony != 2 & cagvwho != 1 & cagvwho != 2 & cagvwho != 1 & cagvwho != 3 & cafrewho != 6 & cafrewho != 7, 1, 0))



### 9. Took alcohol from another person’s home *alcelsehome*

NSDUH2018 <- mutate(NSDUH2018, alcelsehome = ifelse(cafrewho == 6 | CAFRESP2 == 6 & cabuyfre == 2 & cabuyfre != 1 & cagvmony != 1 & cagvmony != 2 & cagvwho != 2 & cagvwho != 1 & cagvwho != 3 & cagvwho != 4 & cafrewho !=  5 & cafrewho != 7, 1, 0))



### 10. Got the alcohol for free by some other means *alcother*

NSDUH2018 <- mutate(NSDUH2018, alcother = ifelse(cafrewho == 7 | CAFRESP2 <= 33 & CAFRESP2 != 13 & CAFRESP2 != 6 & CAFRESP2 !=5 & CAFRESP2 != 4 & CAFRESP2 != 3 & cagvwho != 2 & cagvwho != 1 & cagvwho != 3 & cagvwho == 4 & cafrewho != 1 & cafrewho != 2 & cafrewho != 3 & cafrewho !=  5 & cafrewho != 6, 1, 0))


### 11. Alcohol got for free

NSDUH2018 <- mutate(NSDUH2018, alcother = ifelse(cafrewho == 7 & cagvwho != 2 & cagvwho != 1 & cagvwho != 3 & cagvwho == 4 & cafrewho !=  5 & cafrewho != 6, 1, 0))

Drinking Locations Variables:

### 1. Car:

NSDUH2018 <- mutate(NSDUH2018, drinkloccar = ifelse(uadcar == 1 & uadhome != 1 & uadothm != 1 & uadpubl != 1 & uadbar != 1 & uadevnt != 1 & uadschl != 1 & uadotsp != 1 & uadotsp != 2
 & uadotsp != 3 & uadotsp != 4 & uadotsp != 5
& uadotsp != 6 , 1, 0))





### 2. Home:

NSDUH2018 <- mutate(NSDUH2018, drinklochome = ifelse(uadhome == 1 & uadcar != 1 & uadothm != 1  & uadpubl != 1 & uadbar != 1 & uadevnt != 1  & uadschl != 1 & uadotsp != 1 & uadotsp != 2 & uadotsp != 3 & uadotsp != 4 & uadotsp != 5 & uadotsp != 6, 1, 0))




### 3. Someone else’s home:

NSDUH2018 <- mutate(NSDUH2018, drinklocelse = ifelse(uadothm == 1 & uadhome != 1 & uadcar != 1  & uadpubl != 1 & uadbar != 1 & uadevnt != 1   & uadschl != 1 & uadotsp != 1 & uadotsp != 2
   & uadotsp != 3 & uadotsp != 4 & uadotsp != 5
 & uadotsp != 6, 1, 0))

table(NSDUH2018$drinklocelse)
## 
##    0    1 
## 1699 1225
### 4. Park, beach, or parking lot:

NSDUH2018 <- mutate(NSDUH2018, drinklocpark = ifelse(uadpubl == 1 & uadhome != 1 & uadothm != 1  & uadcar != 1 & uadbar != 1 & uadevnt != 1 & uadschl != 1 & uadotsp != 1 & uadotsp != 2 & uadotsp != 3 & uadotsp != 4 & uadotsp != 5 & uadotsp != 6, 1, 0))

table(NSDUH2018$drinklocpark)
## 
##    0    1 
## 2857   68
### 5. Restaurant, bar, or club:

NSDUH2018 <- mutate(NSDUH2018, drinklocbar = ifelse(uadbar == 1 & uadhome != 1 & uadothm != 1 & uadpubl != 1 & uadcar != 1 & uadevnt != 1  & uadschl != 1 & uadotsp != 1 & uadotsp != 2
  & uadotsp != 3 & uadotsp != 4 & uadotsp != 5
  & uadotsp != 6, 1, 0))


table(NSDUH2018$drinklocbar)
## 
##    0    1 
## 2794  131
### 6. Concert or sports event:

NSDUH2018 <- mutate(NSDUH2018, drinklocevent = ifelse(uadevnt == 1 & uadhome != 1 & uadothm != 1  & uadpubl != 1 & uadbar != 1 & uadcar != 1  & uadschl != 1 & uadotsp != 1 & uadotsp != 2& uadotsp != 3 & uadotsp != 4 & uadotsp != 5& uadotsp != 6, 1, 0))

table(NSDUH2018$drinklocevent)
## 
##    0    1 
## 2895   30
### 7. School:

NSDUH2018 <- mutate(NSDUH2018, drinklocschl = ifelse(uadschl == 1 & uadhome != 1 & uadothm != 1  & uadpubl != 1 & uadbar != 1 & uadevnt != 1 & uadcar != 1 & uadotsp != 1 & uadotsp != 2
  & uadotsp != 3 & uadotsp != 4 & uadotsp != 5
 & uadotsp != 6, 1, 0))


table(NSDUH2018$drinklocschl)
## 
##    0    1 
## 2894   31
### 8. Party, wedding, or other celebration only:

NSDUH2018 <- mutate(NSDUH2018, drinklocceleb = ifelse(uadotsp == 1, 1, 0))

table(NSDUH2018$drinklocceleb)
## 
##     0     1 
## 55771    58
### 9. Other only:

NSDUH2018 <- mutate(NSDUH2018, drinklocoth = ifelse(uadotsp >= 2 & uadotsp <=6, 1, 0))

table(NSDUH2018$drinklocoth)
## 
##     0     1 
## 55758    71

Number of People Accompanying While Drinking Variables

### 1. Drank Alone last time, *nopalone*:

# Create source of number of people accompanying cadrpeop:
NSDUH2018 <- mutate(NSDUH2018, nopalone =  ifelse (cadrpeop == 1, 1, 0))




### 2. Drank with  one other person, *nopone*:

NSDUH2018 <- mutate(NSDUH2018, nopone = ifelse (cadrpeop == 2, 1, 0))




### 3. Drank with more than one other person, *nopmo*:

NSDUH2018<- mutate(NSDUH2018, nopmo = ifelse (cadrpeop == 3, 1, 0))

Counfounding variables to be created

### 1. Race/Ethnicity

# We can factor the races according to various categories:

NSDUH2018 <- mutate(NSDUH2018, fNEWRACE2 = factor(NEWRACE2, levels = c(1, 2, 3, 4, 5, 6, 7),
                                                   labels = c("NHWhite", "NHBlack",  "NHNA", "NHHI", "NHAsian", "NHMulti", "Hispanic")))





### 2. Age of Drinking Onset

#need to remove values over 87 (missing variables)
NSDUH2018 <- mutate(NSDUH2018, drinkonset = ifelse(alctry <= 87, alctry, NA))



### 3. No. of days drank in past year

NSDUH2018 <- mutate(NSDUH2018, daysdrank = ifelse(alcyrtot <= 365, alcyrtot, NA))




### 4. Past 30-day cigarette use

NSDUH2018 <- mutate(NSDUH2018, ciguse = ifelse(cigrec == 1, 1,0))



### 5. “How do you feel about someone your age having one or two drinks of an alcoholic beverage nearly every day?”

#First Mutate 2s and 3s to 0
NSDUH2018 <- mutate(NSDUH2018, drinkdailyop = ifelse(yegaldly == 2 | yegaldly ==3, 0, yegaldly))

#Then recode missing variables as NA
NSDUH2018 <- mutate(NSDUH2018, drinkdailyop = ifelse(drinkdailyop <= 1, drinkdailyop, NA))




### 6. “How do you think your parents would feel about you having one or two drinks of an alcoholic beverage nearly every day?”

# First mutate 2s and 3s to 0
NSDUH2018 <- mutate(NSDUH2018, parentdrinkop = ifelse(yepaldly == 2 | yepaldly ==3, 0, yepaldly))

#Then recode missing variables as NA 
NSDUH2018 <- mutate(NSDUH2018, parentdrinkop = ifelse(parentdrinkop <= 1, parentdrinkop, NA))





### 7. “How do you think your close friends would feel about you having one or two drinks of an alcoholic beverage nearly every day?”

#recode 2s and 3s to 0
NSDUH2018 <- mutate(NSDUH2018, frienddrinkop = ifelse(yefaldly == 2 | yefaldly ==3, 0, yefaldly))
#recode missing variables as NA 
NSDUH2018 <- mutate(NSDUH2018, frienddrinkop = ifelse(frienddrinkop <= 1, frienddrinkop, NA))




### 8. “How many of the students in your grade at school would you say drink alcoholic beverages?”

#recode 3s and 4s to 1 and 1s and 2s to 0
NSDUH2018 <- mutate(NSDUH2018, studentdrinkop = ifelse(yestsalc == 1 | yestsalc == 2, 0, yestsalc))
NSDUH2018 <- mutate(NSDUH2018, studentdrinkop = ifelse(studentdrinkop == 3 | studentdrinkop == 4, 1, studentdrinkop))
#recode missing variables as NA 
NSDUH2018 <- mutate(NSDUH2018, studentdrinkop = ifelse(studentdrinkop <= 1, studentdrinkop, NA))


### 9. Non-metropolitan vs. Metropolitan

NSDUH2018 <- mutate(NSDUH2018, nonmetro = ifelse(COUTYP4 == 3, 1, 0))

Selecting our Target Population (Inclusion Criteria)

A survey design is set up with the code:

desg <-svydesign(id= ~verep, strata= ~vestr, weights= ~ANALWT_C, nest = TRUE, data = NSDUH2018, na.rm = TRUE)

To select our target group of 12 - 17 years of age who had consumed alcohol in the last 30 days, within the survey design:

desg2018 <- subset(desg, catage == 1 & alcrec == 1)

Exploration of the data

To check the mean weight of the whole sample, and the mean weight of the subset.

#Whole Sample
svymean(~verep, design = desg, na.rm = TRUE)
##         mean     SE
## verep 1.5058 0.0729
## The mean weight for the whole sample is 1.5

#Age greater than 12 and less than 17
svymean(~verep, design = desg2018, na.rm = TRUE)
##        mean     SE
## verep 1.527 0.0752
## The mean weight for subset desg2018 is also 1.5

# To check the number of people in different age groups.
table(NSDUH2018$catage)
## 
##     1     2     3     4 
## 13287 13637  8794 20595
## The number of participants aged 12 - 17 (catage = 1) years is 13287.


table(NSDUH2018$HED)
## 
##     0     1 
## 51136  5177
## The number of participants who was categorized as HED is 5177.

# The number of participants who had a drink in the last 30 days
table(NSDUH2018$alcrec)
## 
##     1     2     3     8     9    11    85    91    97    98 
## 25243  8263  6310   357   253    13     4 15849     8    13
## The total number for people who had a drink in the last 30 days (alcrec = 1) is 25243.

Analyzing the data for HED in the population

svyglm(HED ~ alcselfonprem + alcselfofprem + alcelse + alcparent + alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + alcother, design=desg,
       family = "binomial") 
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (100) clusters.
## svydesign(id = ~verep, strata = ~vestr, weights = ~ANALWT_C, 
##     nest = TRUE, data = NSDUH2018, na.rm = TRUE)
## 
## Call:  svyglm(formula = HED ~ alcselfonprem + alcselfofprem + alcelse + 
##     alcparent + alcfam21 + alcelse21 + alcunder21 + alchome + 
##     alcelsehome + alcother, design = desg, family = "binomial")
## 
## Coefficients:
##   (Intercept)  alcselfonprem  alcselfofprem        alcelse      alcparent  
##       -2.1893         2.4632        -0.2077        -0.1034        -1.4379  
##      alcfam21      alcelse21     alcunder21        alchome    alcelsehome  
##        1.1364         1.7872         1.8505         0.8398         1.3016  
##      alcother  
##            NA  
## 
## Degrees of Freedom: 56312 Total (i.e. Null);  41 Residual
## Null Deviance:       31850 
## Residual Deviance: 31170     AIC: 29740

Table 1 & Data Exploration

This table shows shows the mean for different demographic variables like age, sex and various races/ethnicities. These values will be compared with that of Dr. Rossheim’s paper.

#Mean age 
svymean(~newage, design = desg2018, na.rm = T)
##          mean     SE
## newage 15.944 0.0437
#Percentage of males
svytable(~irsex, design = desg2018, Ntotal = 100)
## irsex
##        1        2 
## 48.51159 51.48841
#Mean of males = 48.5%

#Percentage of races
svytable(~fNEWRACE2, design = desg2018, Ntotal = 100)
## fNEWRACE2
##    NHWhite    NHBlack       NHNA       NHHI    NHAsian    NHMulti   Hispanic 
## 65.7697320  7.2897253  0.3807356  0.1483516  3.0548808  3.6704080 19.6861666
#Average age of drinking onset
svymean(~drinkonset, design = desg2018, na.rm = T)
##              mean    SE
## drinkonset 13.839 0.086
#Average no of days drank
svymean(~daysdrank, design = desg2018, na.rm = T)
##           mean     SE
## daysdrank 40.6 2.6139
# No. of participants that had HED event
svytable(~HED, design = desg2018, Ntotal = 100)
## HED
##        0        1 
## 70.46581 29.53419

Table 2 Analyses: Compute table with Number of people (this is weighted up to the population of interest)

The purpose of this table is to compare the proportions between people who drank alcohol in different locatiotions and had gotten alcohol from various sources. This comparison is made between one location with a source.

#percentages of drinking location
svytable(~drinkloccar, design = desg2018, Ntotal = 100)
## drinkloccar
##         0         1 
## 98.169819  1.830181
svytable(~drinklochome, design = desg2018, Ntotal = 100)
## drinklochome
##        0        1 
## 67.65391 32.34609
svytable(~drinklocelse, design = desg2018, Ntotal = 100)
## drinklocelse
##        0        1 
## 54.85728 45.14272
svytable(~drinklocpark, design = desg2018, Ntotal = 100)
## drinklocpark
##         0         1 
## 97.318091  2.681909
svytable(~drinklocbar, design = desg2018, Ntotal = 100)
## drinklocbar
##         0         1 
## 98.602484  1.397516
svytable(~drinklocevent, design = desg2018, Ntotal = 100)
## drinklocevent
##          0          1 
## 99.3714833  0.6285167
svytable(~drinklocschl, design = desg2018, Ntotal = 100)
## drinklocschl
##          0          1 
## 99.7554377  0.2445623
svytable(~drinklocceleb, design = desg2018, Ntotal = 100)
## drinklocceleb
##         0         1 
## 97.042574  2.957426
svytable(~drinklocoth, design = desg2018, Ntotal = 100)
## drinklocoth
##         0         1 
## 97.147845  2.852155
#two or more locations were still analyzed even though we know they have a value of 0, for consistency 

#percentages of source of alcohol 
svytable(~alcselfofprem, design = desg2018, Ntotal = 100)
## alcselfofprem
##         0         1 
## 98.190293  1.809707
svytable(~alcselfonprem, design = desg2018, Ntotal = 100)
## alcselfonprem
##         0         1 
## 96.148869  3.851131
svytable(~alcelse, design = desg2018, Ntotal = 100)
## alcelse
##         0         1 
## 98.884263  1.115737
svytable(~alcparent, design = desg2018, Ntotal = 100)
## alcparent
##        0        1 
## 58.01304 41.98696
svytable(~alcfam21, design = desg2018, Ntotal = 100)
## alcfam21
##        0        1 
## 56.59452 43.40548
svytable(~alcelse21, design = desg2018, Ntotal = 100)
## alcelse21
##        0        1 
## 79.67402 20.32598
svytable(~alcunder21, design = desg2018, Ntotal = 100)
## alcunder21
##        0        1 
## 77.86929 22.13071
svytable(~alchome, design = desg2018, Ntotal = 100)
## alchome
##        0        1 
## 86.41833 13.58167
svytable(~alcelsehome, design = desg2018, Ntotal = 100)
## alcelsehome
##         0         1 
## 97.331881  2.668119
svytable(~alcother, design = desg2018, Ntotal = 100)
## alcother
##   0 
## 100

Next step is to understand what the proportion of each source of alcohol was for each drinking locations. First, new variables for drinking locations and sources. Then the proportions will be computed.

#create new drinking location variable

#Drank at car
NSDUH2018$drinkloc[NSDUH2018$drinkloccar==1] <- "1"
#Drank at home
NSDUH2018$drinkloc[NSDUH2018$drinklochome==1] <- "2"
#Drank at someone else's home
NSDUH2018$drinkloc[NSDUH2018$drinklocelse==1] <- "3"
#Drank at a park, beach, or parking lot
NSDUH2018$drinkloc[NSDUH2018$drinklocpark==1] <- "4"
#Drank at restaurant, bar, or club
NSDUH2018$drinkloc[NSDUH2018$drinklocbar==1] <- "5"
#Drank at concert or sports game 
NSDUH2018$drinkloc[NSDUH2018$drinklocevent==1] <- "6"
#Drank at school
NSDUH2018$drinkloc[NSDUH2018$drinklocschl==1] <- "7"
#Drank at party, wedding, or celebration
NSDUH2018$drinkloc[NSDUH2018$drinklocceleb==1] <- "8"
#Drank at other location only 
NSDUH2018$drinkloc[NSDUH2018$drinklocoth==1] <- "9"
#Drank in two locations
NSDUH2018$drinkloc[NSDUH2018$drinktwoloc==1] <- "10"
#Drank in three or more locations?
NSDUH2018$drinkloc[NSDUH2018$drinkmoreloc==1] <- "11"
#View new variable
table(NSDUH2018$drinkloc)
## 
##    1    2    3    4    5    6    7    8    9 
##   49  948 1225   68  131   30   31   58   71
#To create new source variable 

#Self off-premises
NSDUH2018$drinksource[NSDUH2018$alcselfofprem==1] <- "1"
#Self on-premises
NSDUH2018$drinksource[NSDUH2018$alcselfonprem==1] <- "2"
#Other person
NSDUH2018$drinksource[NSDUH2018$alcelse==1] <- "3"
#Parent
NSDUH2018$drinksource[NSDUH2018$alcparent==1] <- "4"
#Family over 21
NSDUH2018$drinksource[NSDUH2018$alcfam21==1] <- "5"
#Someone else over 21
NSDUH2018$drinksource[NSDUH2018$alcelse21==1] <- "6"
#Under 21
NSDUH2018$drinksource[NSDUH2018$alcunder21==1] <- "7"
#Home
NSDUH2018$drinksource[NSDUH2018$alchome==1] <- "8"
#Someone else's home
NSDUH2018$drinksource[NSDUH2018$alcelsehome==1] <- "9"
#Other
NSDUH2018$drinksource[NSDUH2018$alcother==1] <- "10"
#To view new variable
table(NSDUH2018$drinksource)
## 
##     4     5     6     7     8     9 
##    24 54394  1008   444   190    69

The survey design and the subset will need to be re-run after the above variables are created.

#Survey design & subset codes
desg <-svydesign(id= ~verep, strata= ~vestr, weights= ~ANALWT_C, nest = TRUE, data = NSDUH2018, na.rm = TRUE)

desg2018 <- subset(desg, catage == 1 & alcrec == 1)

To create the table, the following code must be run. This codes for each categorical variable against each other.

Table2 <- svytable(~drinksource + drinkloc, design = desg2018)
prop.table(Table2, margin = 2)
##            drinkloc
## drinksource           1           2           3           4           5
##           4 0.000000000 0.006514432 0.000000000 0.000000000 0.000000000
##           5 0.373135029 0.501806463 0.203191716 0.204169006 0.856669697
##           6 0.118910141 0.075586458 0.327204686 0.250582859 0.143330303
##           7 0.268779889 0.048124950 0.390433154 0.262725913 0.000000000
##           8 0.239174941 0.361287268 0.025889560 0.230629423 0.000000000
##           9 0.000000000 0.006680430 0.053280884 0.051892798 0.000000000
##            drinkloc
## drinksource           6           7           8           9
##           4 0.000000000 0.000000000 0.000000000 0.000000000
##           5 0.086942232 0.765995666 0.180783094 0.595596020
##           6 0.172947921 0.000000000 0.385222807 0.256992755
##           7 0.740109847 0.187732974 0.375404735 0.140238726
##           8 0.000000000 0.046271361 0.058589364 0.007172499
##           9 0.000000000 0.000000000 0.000000000 0.000000000

To create the percentages, each cell value was multiplied by a 100.

Table 3 & Logistic Regression Model

Since our variable of interest, HED, is binary, a logistic regression model was created to include all the new variables for sources of alcohol, locations and no. of people accompaying (except the reference variables) and also confounding variables like demographic factors and other drinking variables.

svyglm <- svyglm(HED ~ alcselfofprem + alcselfonprem + alcelse + alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + alcother + drinkloccar + drinklocelse + drinklocpark + drinklocbar + drinklocevent + drinklocschl + drinklocceleb + nopone + nopmo + newage + irsex + fNEWRACE2 + nonmetro + drinkonset + daysdrank + ciguse + drinkdailyop + parentdrinkop + frienddrinkop + studentdrinkop, design = desg2018, family = 'binomial')
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
#age, sex and race/ethnicity and other variables are added to the model to control for confounding

svyglm %>% summary()
## 
## Call:
## svyglm(formula = HED ~ alcselfofprem + alcselfonprem + alcelse + 
##     alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + 
##     alcother + drinkloccar + drinklocelse + drinklocpark + drinklocbar + 
##     drinklocevent + drinklocschl + drinklocceleb + nopone + nopmo + 
##     newage + irsex + fNEWRACE2 + nonmetro + drinkonset + daysdrank + 
##     ciguse + drinkdailyop + parentdrinkop + frienddrinkop + studentdrinkop, 
##     design = desg2018, family = "binomial")
## 
## Survey design:
## subset(desg, catage == 1 & alcrec == 1)
## 
## Coefficients: (1 not defined because of singularities)
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -6.898116   1.698758  -4.061 0.000813 ***
## alcselfofprem      -2.309206   1.180545  -1.956 0.067096 .  
## alcselfonprem       2.501223   0.929017   2.692 0.015420 *  
## alcelse            -0.862445   1.853736  -0.465 0.647659    
## alcfam21            0.442398   0.271543   1.629 0.121659    
## alcelse21           0.357392   0.328656   1.087 0.292022    
## alcunder21          0.473309   0.375530   1.260 0.224561    
## alchome             1.124934   0.452966   2.483 0.023736 *  
## alcelsehome         0.214695   0.563193   0.381 0.707770    
## drinkloccar         0.108083   0.845636   0.128 0.899797    
## drinklocelse        0.308802   0.331128   0.933 0.364096    
## drinklocpark       -0.878072   0.739099  -1.188 0.251155    
## drinklocbar         0.002654   0.742226   0.004 0.997188    
## drinklocevent      -3.555068   1.235021  -2.879 0.010426 *  
## drinklocschl      -11.133472   1.144683  -9.726 2.32e-08 ***
## drinklocceleb       0.641251   0.604412   1.061 0.303554    
## nopone              0.426558   0.494497   0.863 0.400358    
## nopmo               1.435991   0.463118   3.101 0.006495 ** 
## newage              0.230980   0.113977   2.027 0.058693 .  
## irsex               0.332426   0.206242   1.612 0.125409    
## fNEWRACE2NHBlack   -0.457042   0.518770  -0.881 0.390598    
## fNEWRACE2NHNA      -0.812090   1.157714  -0.701 0.492508    
## fNEWRACE2NHHI       0.578757   1.274712   0.454 0.655551    
## fNEWRACE2NHAsian   -0.854107   0.762931  -1.120 0.278493    
## fNEWRACE2NHMulti    0.167976   0.505298   0.332 0.743628    
## fNEWRACE2Hispanic  -0.128652   0.311338  -0.413 0.684609    
## nonmetro            0.136671   0.290257   0.471 0.643726    
## drinkonset         -0.046836   0.050578  -0.926 0.367397    
## daysdrank           0.001129   0.001677   0.673 0.510088    
## ciguse              0.460367   0.309478   1.488 0.155181    
## drinkdailyop        0.043116   0.346954   0.124 0.902558    
## parentdrinkop      -0.185360   0.281585  -0.658 0.519179    
## frienddrinkop       0.168535   0.309559   0.544 0.593213    
## studentdrinkop      0.688684   0.201245   3.422 0.003248 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 0.9845113)
## 
## Number of Fisher Scoring iterations: 12
## To get the odds ratios from the coefficients
svyglm %>% coef() %>% exp()
## To get the 95% confidence intervals (CIs):
svyglm %>% confint() %>% exp()
# To get the tabulated form of ORs for drinking characteristics 95% CIs, and p-values:
Table3 <- tidy(svyglm, conf.int = T) %>% mutate(., OR = exp(estimate),
    LB = exp(conf.low),
    UB = exp(conf.high)) %>%
  select(., term, OR, LB, UB, p.value) %>%
  filter(term != "(Intercept)")

kable(Table3, caption = ("HED and Drinking Characteristics OR, CIs and p-values"), digits = 2)
HED and Drinking Characteristics OR, CIs and p-values
term OR LB UB p.value
alcselfofprem 0.10 0.01 1.00 0.07
alcselfonprem 12.20 1.97 75.34 0.02
alcelse 0.42 0.01 15.97 0.65
alcfam21 1.56 0.91 2.65 0.12
alcelse21 1.43 0.75 2.72 0.29
alcunder21 1.61 0.77 3.35 0.22
alchome 3.08 1.27 7.48 0.02
alcelsehome 1.24 0.41 3.74 0.71
drinkloccar 1.11 0.21 5.84 0.90
drinklocelse 1.36 0.71 2.61 0.36
drinklocpark 0.42 0.10 1.77 0.25
drinklocbar 1.00 0.23 4.29 1.00
drinklocevent 0.03 0.00 0.32 0.01
drinklocschl 0.00 0.00 0.00 0.00
drinklocceleb 1.90 0.58 6.21 0.30
nopone 1.53 0.58 4.04 0.40
nopmo 4.20 1.70 10.42 0.01
newage 1.26 1.01 1.58 0.06
irsex 1.39 0.93 2.09 0.13
fNEWRACE2NHBlack 0.63 0.23 1.75 0.39
fNEWRACE2NHNA 0.44 0.05 4.29 0.49
fNEWRACE2NHHI 1.78 0.15 21.70 0.66
fNEWRACE2NHAsian 0.43 0.10 1.90 0.28
fNEWRACE2NHMulti 1.18 0.44 3.18 0.74
fNEWRACE2Hispanic 0.88 0.48 1.62 0.68
nonmetro 1.15 0.65 2.02 0.64
drinkonset 0.95 0.86 1.05 0.37
daysdrank 1.00 1.00 1.00 0.51
ciguse 1.58 0.86 2.91 0.16
drinkdailyop 1.04 0.53 2.06 0.90
parentdrinkop 0.83 0.48 1.44 0.52
frienddrinkop 1.18 0.65 2.17 0.59
studentdrinkop 1.99 1.34 2.95 0.00

Assumption Check for Logistic Regression Model

#newage vs HED

ggplot(NSDUH2018, aes(x=newage, y=HED)) +geom_point() + geom_smooth()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

# number of day/year reported drinking vs HED
ggplot(NSDUH2018, aes(x=daysdrank, y=HED)) +geom_point() + geom_smooth()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 22904 rows containing non-finite values (stat_smooth).
## Warning: Removed 22904 rows containing missing values (geom_point).

# age of drinking onset vs HED
ggplot(NSDUH2018, aes(x=drinkonset, y=HED)) +geom_point() + geom_smooth()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 16438 rows containing non-finite values (stat_smooth).
## Warning: Removed 16438 rows containing missing values (geom_point).

All 3 plots show some degree of linearity and we do not see a prominent curve. The scatter of 1s and 0s are not clustered together.

Sensitivity Analysis

1. Using a linear regression model for extreme binge drinking

a) Linear regression model with the variable for ‘number of drinks consumed’:

# Run glm model: a linear regression was considered for this variable since it is non-binary 

svyglm2a <- svyglm(cadrlast ~ alcselfofprem + alcselfonprem + alcelse + alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + alcother + drinkloccar + drinklocelse + drinklocpark + drinklocbar + drinklocevent + drinklocschl + drinklocceleb + nopone + nopmo + newage + irsex + fNEWRACE2 + nonmetro + drinkonset + daysdrank + ciguse + drinkdailyop + parentdrinkop + frienddrinkop + studentdrinkop + nonmetro, design = desg2018)

# Create tidy output
tidy_glm <- tidy(svyglm2a, conf.int = T) %>%
select(., term, estimate, conf.low, conf.high) %>%
filter(term != "(Intercept)")
tidy_glm
## # A tibble: 33 x 4
##    term          estimate conf.low conf.high
##    <chr>            <dbl>    <dbl>     <dbl>
##  1 alcselfofprem   90.6     -84.4      266. 
##  2 alcselfonprem    5.99    -31.0       43.0
##  3 alcelse         -9.90    -59.5       39.7
##  4 alcfam21         2.36     -9.19      13.9
##  5 alcelse21       -2.36    -17.4       12.7
##  6 alcunder21       2.58    -11.5       16.6
##  7 alchome         -7.21    -31.8       17.4
##  8 alcelsehome      0.596   -14.8       16.0
##  9 drinkloccar    -23.9     -58.0       10.1
## 10 drinklocelse    -7.09    -35.0       20.8
## # ... with 23 more rows
## To get the odds ratios from the coefficients
svyglm2a %>% coef() %>% exp()
## To get the 95% confidence intervals (CIs):
svyglm2a %>% confint() %>% exp()
## Tabulation of new HED variables with ORs, 95% CIs, p-values
Table4a <- tidy(svyglm2a, conf.int = T) %>% mutate(., OR = exp(estimate),
    LB = exp(conf.low),
    UB = exp(conf.high)) %>%
  select(., term, OR, LB, UB, p.value) %>%
  filter(term != "(Intercept)")

kable(Table4a, caption = ("Sensitivity Analysis Using Bring Drinking Variable OR, CIs and p-values"), digits = 2)
Sensitivity Analysis Using Bring Drinking Variable OR, CIs and p-values
term OR LB UB p.value
alcselfofprem 2.325636e+39 0.00 2.562785e+115 0.32
alcselfonprem 4.011900e+02 0.00 4.799285e+18 0.75
alcelse 0.000000e+00 0.00 1.763892e+17 0.70
alcfam21 1.055000e+01 0.00 1.092313e+06 0.69
alcelse21 9.000000e-02 0.00 3.124457e+05 0.76
alcunder21 1.321000e+01 0.00 1.648370e+07 0.72
alchome 0.000000e+00 0.00 3.567087e+07 0.57
alcelsehome 1.820000e+00 0.00 8.561685e+06 0.94
drinkloccar 0.000000e+00 0.00 2.554470e+04 0.19
drinklocelse 0.000000e+00 0.00 1.063815e+09 0.62
drinklocpark 0.000000e+00 0.00 7.170305e+11 0.70
drinklocbar 0.000000e+00 0.00 2.187251e+11 0.51
drinklocevent 0.000000e+00 0.00 4.054941e+09 0.42
drinklocschl 0.000000e+00 0.00 4.390636e+15 0.21
drinklocceleb 9.631578e+14 0.00 6.982500e+49 0.41
nopone 5.200000e-01 0.00 5.637719e+06 0.94
nopmo 7.772420e+03 0.00 1.973692e+16 0.55
newage 0.000000e+00 0.00 1.780500e+02 0.19
irsex 1.000000e-02 0.00 2.924265e+04 0.55
fNEWRACE2NHBlack 7.300000e-01 0.00 1.798367e+12 0.98
fNEWRACE2NHNA 0.000000e+00 0.00 2.535400e+03 0.13
fNEWRACE2NHHI 0.000000e+00 0.00 4.475243e+10 0.69
fNEWRACE2NHAsian 1.243665e+18 0.00 1.778448e+61 0.42
fNEWRACE2NHMulti 6.426066e+04 0.00 4.932788e+26 0.67
fNEWRACE2Hispanic 0.000000e+00 0.00 5.606800e+02 0.25
nonmetro 3.647000e+01 0.00 3.221931e+13 0.80
drinkonset 6.195000e+01 0.39 9.771550e+03 0.13
daysdrank 9.400000e-01 0.82 1.080000e+00 0.39
ciguse 2.000000e-02 0.00 2.744895e+05 0.64
drinkdailyop 1.000000e-02 0.00 2.800228e+04 0.56
parentdrinkop 1.428521e+21 0.00 6.807087e+55 0.25
frienddrinkop 9.712100e+02 0.00 6.941856e+10 0.47
studentdrinkop 0.000000e+00 0.00 3.577561e+06 0.35

b) Using the binge variable in the logistic regression model:

svyglm2b <- svyglm(binge ~ alcselfofprem + alcselfonprem + alcelse + alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + alcother + drinkloccar + drinklocelse + drinklocpark + drinklocbar + drinklocevent + drinklocschl + drinklocceleb + nopone + nopmo + newage + irsex + fNEWRACE2 + nonmetro + drinkonset + daysdrank + ciguse + drinkdailyop + parentdrinkop + frienddrinkop + studentdrinkop + nonmetro, design = desg2018, family = 'binomial')
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# The binge variable is added in this model

svyglm2b %>% summary()
## 
## Call:
## svyglm(formula = binge ~ alcselfofprem + alcselfonprem + alcelse + 
##     alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + 
##     alcother + drinkloccar + drinklocelse + drinklocpark + drinklocbar + 
##     drinklocevent + drinklocschl + drinklocceleb + nopone + nopmo + 
##     newage + irsex + fNEWRACE2 + nonmetro + drinkonset + daysdrank + 
##     ciguse + drinkdailyop + parentdrinkop + frienddrinkop + studentdrinkop + 
##     nonmetro, design = desg2018, family = "binomial")
## 
## Survey design:
## subset(desg, catage == 1 & alcrec == 1)
## 
## Coefficients: (1 not defined because of singularities)
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -5.624e+00  3.420e+00  -1.644  0.11844    
## alcselfofprem     -1.654e-01  1.529e+00  -0.108  0.91514    
## alcselfonprem      1.799e+00  1.061e+00   1.696  0.10821    
## alcelse           -1.666e+01  1.265e+00 -13.171 2.39e-10 ***
## alcfam21           3.966e-01  4.752e-01   0.834  0.41557    
## alcelse21         -4.110e-01  6.779e-01  -0.606  0.55234    
## alcunder21         3.252e-01  5.919e-01   0.549  0.58992    
## alchome           -1.004e+00  1.097e+00  -0.916  0.37257    
## alcelsehome       -3.141e+00  1.193e+00  -2.633  0.01746 *  
## drinkloccar       -1.527e+01  1.221e+00 -12.515 5.28e-10 ***
## drinklocelse       7.354e-01  5.350e-01   1.374  0.18715    
## drinklocpark      -3.904e-01  1.247e+00  -0.313  0.75809    
## drinklocbar       -1.195e+00  1.192e+00  -1.002  0.33030    
## drinklocevent     -1.462e+01  8.467e-01 -17.263 3.27e-12 ***
## drinklocschl      -1.531e+01  1.859e+00  -8.234 2.46e-07 ***
## drinklocceleb      1.230e+00  8.173e-01   1.506  0.15053    
## nopone            -1.600e+00  1.392e+00  -1.149  0.26632    
## nopmo             -4.247e-02  1.075e+00  -0.040  0.96894    
## newage             4.761e-02  2.237e-01   0.213  0.83403    
## irsex             -2.739e-01  3.520e-01  -0.778  0.44724    
## fNEWRACE2NHBlack  -4.113e-01  7.117e-01  -0.578  0.57085    
## fNEWRACE2NHNA     -1.528e+01  9.950e-01 -15.358 2.13e-11 ***
## fNEWRACE2NHHI      3.188e+00  1.986e+00   1.605  0.12686    
## fNEWRACE2NHAsian  -2.704e+00  1.180e+00  -2.292  0.03495 *  
## fNEWRACE2NHMulti  -3.017e+00  8.654e-01  -3.486  0.00283 ** 
## fNEWRACE2Hispanic -8.998e-02  6.646e-01  -0.135  0.89390    
## nonmetro           7.311e-01  4.929e-01   1.483  0.15628    
## drinkonset         4.674e-02  1.174e-01   0.398  0.69545    
## daysdrank         -3.928e-04  2.153e-03  -0.182  0.85742    
## ciguse             1.185e+00  4.324e-01   2.740  0.01396 *  
## drinkdailyop       8.892e-01  7.030e-01   1.265  0.22298    
## parentdrinkop      7.265e-01  5.151e-01   1.410  0.17648    
## frienddrinkop     -1.094e+00  6.171e-01  -1.773  0.09406 .  
## studentdrinkop     1.519e+00  4.865e-01   3.122  0.00621 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 0.9294412)
## 
## Number of Fisher Scoring iterations: 16
## To get the odds ratios from the coefficients
svyglm2b %>% coef() %>% exp()
## To get the 95% confidence intervals (CIs):
svyglm2b %>% confint() %>% exp()
## Tabulation of new HED variables with ORs, 95% CIs, p-values
Table4 <- tidy(svyglm2b, conf.int = T) %>% mutate(., OR = exp(estimate),
    LB = exp(conf.low),
    UB = exp(conf.high)) %>%
  select(., term, OR, LB, UB, p.value) %>%
  filter(term != "(Intercept)")

kable(Table4, caption = ("Sensitivity Analysis Using Bring Drinking Variable OR, CIs and p-values"), digits = 2)
Sensitivity Analysis Using Bring Drinking Variable OR, CIs and p-values
term OR LB UB p.value
alcselfofprem 0.85 0.04 16.97 0.92
alcselfonprem 6.04 0.76 48.34 0.11
alcelse 0.00 0.00 0.00 0.00
alcfam21 1.49 0.59 3.77 0.42
alcelse21 0.66 0.18 2.50 0.55
alcunder21 1.38 0.43 4.42 0.59
alchome 0.37 0.04 3.14 0.37
alcelsehome 0.04 0.00 0.45 0.02
drinkloccar 0.00 0.00 0.00 0.00
drinklocelse 2.09 0.73 5.95 0.19
drinklocpark 0.68 0.06 7.80 0.76
drinklocbar 0.30 0.03 3.13 0.33
drinklocevent 0.00 0.00 0.00 0.00
drinklocschl 0.00 0.00 0.00 0.00
drinklocceleb 3.42 0.69 16.98 0.15
nopone 0.20 0.01 3.09 0.27
nopmo 0.96 0.12 7.87 0.97
newage 1.05 0.68 1.63 0.83
irsex 0.76 0.38 1.52 0.45
fNEWRACE2NHBlack 0.66 0.16 2.67 0.57
fNEWRACE2NHNA 0.00 0.00 0.00 0.00
fNEWRACE2NHHI 24.25 0.49 1189.20 0.13
fNEWRACE2NHAsian 0.07 0.01 0.68 0.03
fNEWRACE2NHMulti 0.05 0.01 0.27 0.00
fNEWRACE2Hispanic 0.91 0.25 3.36 0.89
nonmetro 2.08 0.79 5.46 0.16
drinkonset 1.05 0.83 1.32 0.70
daysdrank 1.00 1.00 1.00 0.86
ciguse 3.27 1.40 7.63 0.01
drinkdailyop 2.43 0.61 9.65 0.22
parentdrinkop 2.07 0.75 5.68 0.18
frienddrinkop 0.33 0.10 1.12 0.09
studentdrinkop 4.57 1.76 11.85 0.01

2. Using the model without the ‘age of onset’ variable

svyglm3 <- svyglm(HED ~ alcselfofprem + alcselfonprem + alcelse + alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + alcother + drinkloccar + drinklocelse + drinklocpark + drinklocbar + drinklocevent + drinklocschl + drinklocceleb + nopone + nopmo + newage + irsex + fNEWRACE2 + nonmetro + daysdrank + ciguse + drinkdailyop + parentdrinkop + frienddrinkop + studentdrinkop + nonmetro, design = desg2018, family = 'binomial')
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
svyglm3 %>% summary()
## 
## Call:
## svyglm(formula = HED ~ alcselfofprem + alcselfonprem + alcelse + 
##     alcfam21 + alcelse21 + alcunder21 + alchome + alcelsehome + 
##     alcother + drinkloccar + drinklocelse + drinklocpark + drinklocbar + 
##     drinklocevent + drinklocschl + drinklocceleb + nopone + nopmo + 
##     newage + irsex + fNEWRACE2 + nonmetro + daysdrank + ciguse + 
##     drinkdailyop + parentdrinkop + frienddrinkop + studentdrinkop + 
##     nonmetro, design = desg2018, family = "binomial")
## 
## Survey design:
## subset(desg, catage == 1 & alcrec == 1)
## 
## Coefficients: (1 not defined because of singularities)
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -6.860344   1.629855  -4.209 0.000527 ***
## alcselfofprem      -2.266183   1.171331  -1.935 0.068901 .  
## alcselfonprem       2.526578   0.932065   2.711 0.014321 *  
## alcelse            -0.875868   1.858030  -0.471 0.643021    
## alcfam21            0.476006   0.273754   1.739 0.099143 .  
## alcelse21           0.377059   0.329423   1.145 0.267366    
## alcunder21          0.472490   0.373592   1.265 0.222102    
## alchome             1.151537   0.452858   2.543 0.020404 *  
## alcelsehome         0.257265   0.582003   0.442 0.663729    
## drinkloccar         0.113664   0.830932   0.137 0.892714    
## drinklocelse        0.287934   0.321081   0.897 0.381683    
## drinklocpark       -0.940300   0.734802  -1.280 0.216912    
## drinklocbar        -0.034102   0.738719  -0.046 0.963688    
## drinklocevent      -3.540753   1.223670  -2.894 0.009679 ** 
## drinklocschl      -11.278171   1.145905  -9.842 1.14e-08 ***
## drinklocceleb       0.661395   0.614369   1.077 0.295902    
## nopone              0.406890   0.503037   0.809 0.429153    
## nopmo               1.425857   0.467548   3.050 0.006898 ** 
## newage              0.186432   0.094653   1.970 0.064467 .  
## irsex               0.314196   0.195983   1.603 0.126297    
## fNEWRACE2NHBlack   -0.445830   0.521574  -0.855 0.403908    
## fNEWRACE2NHNA      -0.824195   1.148696  -0.718 0.482271    
## fNEWRACE2NHHI       0.579254   1.240798   0.467 0.646215    
## fNEWRACE2NHAsian   -0.870579   0.773942  -1.125 0.275424    
## fNEWRACE2NHMulti    0.182932   0.500381   0.366 0.718935    
## fNEWRACE2Hispanic  -0.110876   0.312096  -0.355 0.726522    
## nonmetro            0.150475   0.285803   0.527 0.604968    
## daysdrank           0.001395   0.001668   0.837 0.413802    
## ciguse              0.510880   0.314908   1.622 0.122120    
## drinkdailyop        0.018231   0.344408   0.053 0.958368    
## parentdrinkop      -0.140439   0.282837  -0.497 0.625524    
## frienddrinkop       0.160107   0.314500   0.509 0.616877    
## studentdrinkop      0.723819   0.202611   3.572 0.002177 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 0.9810257)
## 
## Number of Fisher Scoring iterations: 12
## To get the odds ratios from the coefficients
svyglm3 %>% coef() %>% exp()
## To get the 95% confidence intervals (CIs):
svyglm3 %>% confint() %>% exp()
## Tabulation of the new HED variables with ORs, 95% CIs, p-values within the sensitivity analysis:
Table5 <- tidy(svyglm3, conf.int = T) %>% mutate(., OR = exp(estimate),
    LB = exp(conf.low),
    UB = exp(conf.high)) %>%
  select(., term, OR, LB, UB, p.value) %>%
  filter(term != "(Intercept)")

kable(Table5, caption = ("Sensitivity analysis without age of drinking onset variable"), digits = 2)
Sensitivity analysis without age of drinking onset variable
term OR LB UB p.value
alcselfofprem 0.10 0.01 1.03 0.07
alcselfonprem 12.51 2.01 77.74 0.01
alcelse 0.42 0.01 15.89 0.64
alcfam21 1.61 0.94 2.75 0.10
alcelse21 1.46 0.76 2.78 0.27
alcunder21 1.60 0.77 3.34 0.22
alchome 3.16 1.30 7.68 0.02
alcelsehome 1.29 0.41 4.05 0.66
drinkloccar 1.12 0.22 5.71 0.89
drinklocelse 1.33 0.71 2.50 0.38
drinklocpark 0.39 0.09 1.65 0.22
drinklocbar 0.97 0.23 4.11 0.96
drinklocevent 0.03 0.00 0.32 0.01
drinklocschl 0.00 0.00 0.00 0.00
drinklocceleb 1.94 0.58 6.46 0.30
nopone 1.50 0.56 4.03 0.43
nopmo 4.16 1.66 10.40 0.01
newage 1.20 1.00 1.45 0.06
irsex 1.37 0.93 2.01 0.13
fNEWRACE2NHBlack 0.64 0.23 1.78 0.40
fNEWRACE2NHNA 0.44 0.05 4.17 0.48
fNEWRACE2NHHI 1.78 0.16 20.31 0.65
fNEWRACE2NHAsian 0.42 0.09 1.91 0.28
fNEWRACE2NHMulti 1.20 0.45 3.20 0.72
fNEWRACE2Hispanic 0.90 0.49 1.65 0.73
nonmetro 1.16 0.66 2.04 0.60
daysdrank 1.00 1.00 1.00 0.41
ciguse 1.67 0.90 3.09 0.12
drinkdailyop 1.02 0.52 2.00 0.96
parentdrinkop 0.87 0.50 1.51 0.63
frienddrinkop 1.17 0.63 2.17 0.62
studentdrinkop 2.06 1.39 3.07 0.00

The sensitivity analyses show that there are no significant differences beteen our main model and when binge, last drank were added and when age of onset was taken out.