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, unionlibrary(INLA)## Loading required package: Matrix## Loading required package: sp## This is INLA_18.07.12 built 2018-07-12 11:05:18 UTC.
## See www.r-inla.org/contact-us for how to get help.Reading in the birth data:
natl2012 <- read.csv("C:/Users/tobik/OneDrive/Documents/DEM7083 Fertility/Term Paper/Data/natl2012.csv",
header = TRUE)
natl2013 <- read.csv("C:/Users/tobik/OneDrive/Documents/DEM7083 Fertility/Term Paper/Data/natl2013.csv",
header = TRUE)
natl2015 <- read.csv("C:/Users/tobik/OneDrive/Documents/DEM7083 Fertility/Term Paper/Data/natl2015.csv",
header = TRUE)
natl2016 <- read.csv("C:/Users/tobik/OneDrive/Documents/DEM7083 Fertility/Term Paper/Data/natl2016.csv",
header = TRUE)Drawing out the needed variables:
b12<-natl2012%>%
select(yob="dob_yy", sex, bwt="dbwt", mager9, mracerec, umhisp, meduc, lbo_rec, precare_rec, previs_rec, pay_rec, mager)
b13<-natl2013%>%
select(yob="dob_yy", sex, bwt="dbwt",mager9, mracerec, umhisp, meduc, lbo_rec, precare_rec, previs_rec, pay_rec, mager)
#note that variables in 2015/2016 needed to be recoded to match the same spelling as in the previous databases - the definitions are the same
b15<-natl2015%>%
select(yob="dob_yy", sex, bwt="dbwt",mager9, mracerec="mbrace", umhisp="mhisp_r", meduc, lbo_rec, precare_rec="precare5", previs_rec, pay_rec, mager)
b16<-natl2016%>%
select(yob="dob_yy", sex, bwt="dbwt",mager9, mracerec="mbrace", umhisp="mhisp_r", meduc, lbo_rec, precare_rec="precare5", previs_rec, pay_rec, mager)#deleting the read-in data for RAM space
rm(natl2012,natl2013,natl2015,natl2016)Variables:
yob = Birth Year U,R 2012 Year of birth
bwt = Birth Weight – Detail in Grams U,R 0227-8165 Number of grams
mager9 = Mother’s Age Recode 9 1 Under 15 years 2 15-19 years 3 20-24 years 4 25-29 years 5 30-34 years 6 35-39 years 7 40-44 years 8 45-49 years 9 50-64 years
mracerec = Mother’s Race Recode 1 White 2 Black 3 American Indian / Alaskan Native 4 Asian / Pacific Islander
umhisp = Mother’s Hispanic Origin 0 Non-Hispanic 1 Mexican 2 Puerto Rican 3 Cuban 4 Central or South American 5 Other and Unknown Hispanic 9 Origin unknown or not stated
meduc = Mother’s Education 1 8th grade or less 2 9th through 12th grade with no diploma 3 High school graduate or GED completed 4 Some college credit, but not a degree 5 Associate degree (AA, AS) 6 Bachelor’s degree (BA, AB, BS) 7 Master’s degree (MA, MS) 8 Doctorate (PHD, EdD) or Professional Degree (MD, DDS, DVM, LLB, JD) 9 Unknown Blank Not on certificate
lbo_rec = Live Birth Order Recode 1-7 Live birth order 8 Live birth order of 8 or more 9 Unknown or not stated
precare_rec = Month Prenatal Care Began Recode 1 1st to 3rd month 2 4th to 6th month 3 7th to final month 4 No prenatal care 5 Unknown or not stated Blank Not on certificate
previs_rec = Number of Prenatal Visits Recode 01 No visits 02 1 to 2 visits 03 3 to 4 visits 04 5 to 6 visits 05 7 to 8 visits 06 9 to 10 visits 07 11 to 12 visits 08 13 to 14 visits 09 15 to 16 visits 10 17 to 18 visits 11 19 or more visits 12 Unknown or not stated
pay_rec = Payment Recode 1 Medicade 2 Private Insurance 3 Self Pay 4 Other 9 Unknown
mager = Mother’s Single Year of Age (numeric 12-50)
Combining the data into one database and cleaning RAM
birthdata<-rbind(b12,b13, b15, b16)
rm(b12,b13, b15, b16) #not needed anymore - kick!Subsetting for missing values as Complete Cases (CC):
birthdataCC<-subset(birthdata, birthdata$umhisp!=9) #removed missing Hispanic staus reporting 127,593 cases, .8%
birthdataCC<-subset(birthdataCC, birthdata$meduc!='NA') #removed missing values for mother's eductation, additional 782.902 cases
birthdataCC<-subset(birthdataCC, birthdata$precare_rec!='NA')
birthdataCC<-subset(birthdataCC, birthdata$precare_rec!=5)
birthdataCC<-subset(birthdataCC, birthdata$previs_rec!=12)
birthdataCC<-subset(birthdataCC, birthdata$pay_rec!=9)rm(birthdata)Mutating variables:
birthdataCC<-birthdataCC%>%
mutate(y12=ifelse(yob==2012, 1, 0))%>% #dummy-variables for year
mutate(y13=ifelse(yob==2013, 1, 0))%>%
mutate(y15=ifelse(yob==2015, 1, 0))%>%
mutate(y16=ifelse(yob==2016, 1, 0))%>%
mutate(lbw=ifelse(bwt<=2500, 1, 0))%>% #dicotomizing low birth weight
mutate(kid1=ifelse(lbo_rec==1, 1, 0))%>% #dicotomizing for first child
mutate(umhisp=ifelse(umhisp>=1, 1, 0))%>% #combining all Hispanics into one variable
mutate(White=ifelse(mracerec==1, 1, 0))%>%
mutate(Black=ifelse(mracerec==2, 1, 0))%>%
mutate(Native=ifelse(mracerec==3, 1, 0))%>%
mutate(Asian=ifelse(mracerec==4, 1, 0))%>%
mutate(Medicare=ifelse(pay_rec==1, 1, 0))%>%
mutate(Private=ifelse(pay_rec==2, 1, 0))%>%
mutate(Self_Pay=ifelse(pay_rec==3, 1, 0))%>%
mutate(Other=ifelse(pay_rec==4, 1, 0))m1<-glm(previs_rec~lbw+sex+mager+Medicare+Private+Self_Pay+Other, data=birthdataCC, na.action = na.omit, family=poisson)
m2<-glm(previs_rec~lbw+sex+mager+Medicare+Private+Self_Pay+Other+as.factor(yob), data=birthdataCC, na.action = na.omit, family=poisson)
m3<-glm(previs_rec~lbw+sex+mager+Medicare+Private+Self_Pay+Other+as.factor(yob)+kid1, data=birthdataCC, na.action = na.omit, family=poisson)
m4<-glm(previs_rec~lbw+sex+mager+Medicare+Private+Self_Pay+Other+as.factor(yob)+kid1+umhisp+as.factor(mracerec), data=birthdataCC, na.action = na.omit, family=poisson)
m5<-glm(previs_rec~lbw+sex+mager+Medicare+Private+Self_Pay+Other+as.factor(yob)+kid1+umhisp+as.factor(mracerec)+meduc, data=birthdataCC, na.action = na.omit, family=poisson)
stargazer::stargazer(m1,m2,m3,m4,m5, type="html", out="C:/Users/tobik/OneDrive/Documents/DEM7083 Fertility/Term Paper/models.htm")##
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="5"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="5" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="5">previs_rec</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">lbw</td><td>-0.100<sup>***</sup></td><td>-0.100<sup>***</sup></td><td>-0.101<sup>***</sup></td><td>-0.098<sup>***</sup></td><td>-0.097<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0004)</td><td>(0.0004)</td><td>(0.0004)</td><td>(0.0004)</td><td>(0.0004)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">sexM</td><td>-0.005<sup>***</sup></td><td>-0.005<sup>***</sup></td><td>-0.005<sup>***</sup></td><td>-0.005<sup>***</sup></td><td>-0.005<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0002)</td><td>(0.0002)</td><td>(0.0002)</td><td>(0.0002)</td><td>(0.0002)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">mager</td><td>0.003<sup>***</sup></td><td>0.003<sup>***</sup></td><td>0.004<sup>***</sup></td><td>0.004<sup>***</sup></td><td>0.003<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00002)</td><td>(0.00002)</td><td>(0.00002)</td><td>(0.00002)</td><td>(0.00002)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Medicare</td><td>-0.063<sup>***</sup></td><td>-0.061<sup>***</sup></td><td>-0.058<sup>***</sup></td><td>-0.053<sup>***</sup></td><td>-0.046<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Private</td><td>0.004<sup>***</sup></td><td>0.005<sup>***</sup></td><td>0.001</td><td>-0.005<sup>***</sup></td><td>-0.013<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Self_Pay</td><td>-0.204<sup>***</sup></td><td>-0.203<sup>***</sup></td><td>-0.202<sup>***</sup></td><td>-0.198<sup>***</sup></td><td>-0.192<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Other</td><td>-0.062<sup>***</sup></td><td>-0.060<sup>***</sup></td><td>-0.060<sup>***</sup></td><td>-0.058<sup>***</sup></td><td>-0.057<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">as.factor(yob)2013</td><td></td><td>0.005<sup>***</sup></td><td>0.005<sup>***</sup></td><td>0.006<sup>***</sup></td><td>0.006<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(0.0004)</td><td>(0.0004)</td><td>(0.0004)</td><td>(0.0004)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">as.factor(yob)2015</td><td></td><td>-0.005<sup>***</sup></td><td>-0.005<sup>***</sup></td><td>-0.003<sup>***</sup></td><td>-0.004<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(0.0003)</td><td>(0.0003)</td><td>(0.0003)</td><td>(0.0003)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">as.factor(yob)2016</td><td></td><td>-0.009<sup>***</sup></td><td>-0.009<sup>***</sup></td><td>-0.007<sup>***</sup></td><td>-0.008<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(0.0003)</td><td>(0.0003)</td><td>(0.0003)</td><td>(0.0003)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">kid1</td><td></td><td></td><td>0.036<sup>***</sup></td><td>0.036<sup>***</sup></td><td>0.030<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td>(0.0003)</td><td>(0.0003)</td><td>(0.0003)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">umhisp</td><td></td><td></td><td></td><td>-0.027<sup>***</sup></td><td>-0.019<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td>(0.0003)</td><td>(0.0003)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">as.factor(mracerec)2</td><td></td><td></td><td></td><td>-0.039<sup>***</sup></td><td>-0.037<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td>(0.0003)</td><td>(0.0003)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">as.factor(mracerec)3</td><td></td><td></td><td></td><td>-0.092<sup>***</sup></td><td>-0.088<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td>(0.001)</td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">as.factor(mracerec)4</td><td></td><td></td><td></td><td>-0.030<sup>***</sup></td><td>-0.033<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td>(0.0005)</td><td>(0.0005)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">meduc</td><td></td><td></td><td></td><td></td><td>0.011<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td>(0.0001)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>1.913<sup>***</sup></td><td>1.913<sup>***</sup></td><td>1.871<sup>***</sup></td><td>1.887<sup>***</sup></td><td>1.869<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>10,868,813</td><td>10,868,813</td><td>10,868,813</td><td>10,868,813</td><td>10,868,813</td></tr>
## <tr><td style="text-align:left">Log Likelihood</td><td>-23,863,076.000</td><td>-23,862,062.000</td><td>-23,851,510.000</td><td>-23,839,030.000</td><td>-23,830,225.000</td></tr>
## <tr><td style="text-align:left">Akaike Inf. Crit.</td><td>47,726,168.000</td><td>47,724,145.000</td><td>47,703,044.000</td><td>47,678,092.000</td><td>47,660,484.000</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="5" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>AIC(m1, m2,m3,m4,m5)## df AIC
## m1 8 47726168
## m2 11 47724145
## m3 12 47703044
## m4 16 47678092
## m5 17 47660484INLA does not work right now with this data #{r} fit_est<- inla(previs_rec~lbw+sex+mager+Medicare+Private+Self_Pay+Other+y12+y13+y15+y16+kid1+umhisp+White+Black+Native+Asian+meduc, data=birthdataCC, family = "poisson",Ntrials = 1, num.threads = 22,verbose=TRUE, control.predictor = list(link=1))