-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCTOverdoseCode.R
103 lines (87 loc) · 4.23 KB
/
CTOverdoseCode.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
# Will Clifford
# ECON403 Econometrics I
library(car)
library(lmtest)
# Overdose data from Data.gov and the State of Maryland (See https://data.ct.gov/Health-and-Human-Services/Accidental-Drug-Related-Deaths-2012-2018/rybz-nyjw)
ods_csv <- read.csv("Accidental_Drug_Related_Deaths_2012-2018.csv")
### Create New Variables
ods_csv$Female <- ifelse(ods_csv$Sex == "Female", 1, 0)
ods_csv$Black <- ifelse(ods_csv$Race == "Black", 1, 0)
ods_csv$White <- ifelse(ods_csv$Race == "White", 1, 0)
ods_csv$Hispanic <- ifelse(ods_csv$Race == "Hispanic, White" |
ods_csv$Race == "Hispanic, Black", 1, 0)
ods_csv$Asian <- ifelse(ods_csv$Race == "Asian Indian" |
ods_csv$Race == "Asian, Other" |
ods_csv$Race == "Chinese", 1, 0)
ods_csv$Unknown.Other <- ifelse(ods_csv$Race == "Hawaiian" |
ods_csv$Race == "Native American, Other" |
ods_csv$Race == "Other" |
ods_csv$Race == "Unknown" |
ods_csv$Race == "", 1, 0)
ods_csv$Heroin <- ifelse(ods_csv$Heroin == "Y", 1, 0)
ods_csv$Cocaine <- ifelse(ods_csv$Cocaine == "Y", 1, 0)
ods_csv$Fentanyl <- ifelse(ods_csv$Fentanyl == "Y", 1, 0)
ods_csv$FentanylAnalogue <- ifelse(ods_csv$FentanylAnalogue == "Y", 1, 0)
ods_csv$Oxycodone <- ifelse(ods_csv$Oxycodone == "Y", 1, 0)
ods_csv$Oxymorphone <- ifelse(ods_csv$Oxymorphone == "Y", 1, 0)
ods_csv$Ethanol <- ifelse(ods_csv$Ethanol == "Y", 1, 0)
ods_csv$Hydrocodone <- ifelse(ods_csv$Hydrocodone == "Y", 1, 0)
ods_csv$Benzodiazepine <- ifelse(ods_csv$Benzodiazepine == "Y", 1, 0)
ods_csv$Methadone <- ifelse(ods_csv$Methadone == "Y", 1, 0)
ods_csv$Amphet <- ifelse(ods_csv$Amphet == "Y", 1, 0)
ods_csv$Tramad <- ifelse(ods_csv$Tramad == "Y", 1, 0)
ods_csv$Morphine_NotHeroin <- ifelse(ods_csv$Morphine_NotHeroin == "Y", 1, 0)
ods_csv$Hydromorphone <- ifelse(ods_csv$Hydromorphone == "Y", 1, 0)
ods_csv$OpiateNOS <- ifelse(ods_csv$OpiateNOS == "Y", 1, 0)
ods_csv$NumDrugs <- ods_csv$Heroin + ods_csv$Cocaine + ods_csv$Fentanyl +
ods_csv$FentanylAnalogue + ods_csv$Oxycodone + ods_csv$Oxymorphone +
ods_csv$Ethanol + ods_csv$Benzodiazepine + ods_csv$Methadone +
ods_csv$Amphet + ods_csv$Tramad + ods_csv$Morphine_NotHeroin +
ods_csv$Hydromorphone + ods_csv$OpiateNOS
ods_csv$NumDrugs <- ifelse(ods_csv$NumDrugs == 0, 1, ods_csv$NumDrugs)
### Display Tables
# Gender:
table(ods_csv$Female)
# Race:
table(ods_csv$Black)
table(ods_csv$White)
table(ods_csv$Hispanic)
table(ods_csv$Asian)
table(ods_csv$Unknown.Other)
# Drug:
table(ods_csv$Heroin)
table(ods_csv$Cocaine)
table(ods_csv$Fentanyl)
table(ods_csv$FentanylAnalogue)
table(ods_csv$Oxycodone)
table(ods_csv$Oxymorphone)
table(ods_csv$Ethanol)
table(ods_csv$Hydrocodone)
table(ods_csv$Benzodiazepine)
table(ods_csv$Methadone)
table(ods_csv$Amphet)
table(ods_csv$Tramad)
table(ods_csv$Morphine_NotHeroin)
table(ods_csv$Hydromorphone)
table(ods_csv$OpiateNOS)
# Number of drugs detected after overdose
table(ods_csv$NumDrugs)
### Create Models
# 1. Age, All Linear
ods_lm1 <- lm(Age~Female+Black+Hispanic+Asian+Unknown.Other+Cocaine+Fentanyl+FentanylAnalogue+Oxycodone+Oxymorphone
+Ethanol+Hydrocodone+Benzodiazepine+Methadone+Amphet+Tramad+Morphine_NotHeroin+Hydromorphone+OpiateNOS, data=ods_csv)
# 1a. Age, Race/Gender Interaction
ods_lm1a <- lm(Age~Female+Black+Hispanic+Asian+Unknown.Other + Black*Female+Hispanic*Female+Unknown.Other*Female
+Cocaine+Fentanyl+FentanylAnalogue+Oxycodone+Oxymorphone+Ethanol+Hydrocodone+Benzodiazepine+Methadone+Amphet+Tramad
+Morphine_NotHeroin+Hydromorphone+OpiateNOS, data=ods_csv)
# 2. Quadratic model
ods_quad <- lm(NumDrugs~Age+I(Age^2)+Female+Black+Hispanic+Asian+Unknown.Other, data=ods_csv)
### Display Model Summaries
bptest(ods_lm1)
coeftest(ods_lm1, hccm(ods_lm1))
bptest(ods_quad)
coeftest(ods_quad, hccm(ods_quad))
tp <- coefficients(ods_quad)["Age"]/(-2*coefficients(ods_quad)["I(Age^2)"])
perc_tp <- ifelse(tp > 0, sum(ifelse(is.na(ods_csv$Age), 0, ods_csv$Age > tp)) / length(ods_csv$Age), "Negative Turning Point")
tp
perc_tp