library(car)

## Loading required package: carData

library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

library(survey)

## Loading required package: grid

## Loading required package: Matrix

## Loading required package: survival

## 
## Attaching package: 'survey'

## The following object is masked from 'package:graphics':
## 
##     dotchart

library(questionr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

nhis<-read.csv("C:/Users/canda/Documents/Mortality/nhis_00003.csv/nhis_00003.csv")
names(nhis)

##  [1] "YEAR"        "SERIAL"      "STRATA"      "PSU"         "NHISHID"    
##  [6] "HHWEIGHT"    "REGION"      "METRO"       "LIVINGQTR"   "PERNUM"     
## [11] "NHISPID"     "HHX"         "FMX"         "PX"          "PERWEIGHT"  
## [16] "SAMPWEIGHT"  "FWEIGHT"     "SUPP3WT"     "ASTATFLG"    "CSTATFLG"   
## [21] "AGE"         "SEX"         "SEXORIEN"    "MARSTAT"     "RACEA"      
## [26] "VETSERVWHEN" "OWNERSHIP"   "HINOTCOVE"   "DPMHMDR"     "MORTDODY"   
## [31] "MORTUCOD"    "MORTUCODLD"  "MORTWT"

#sex
nhis$male<-as.factor(ifelse(nhis$SEX==1, "Male", "Female"))

#Age cut into intervals
nhis$agec<-cut(nhis$AGE, breaks=c(0,19,29,39,49,59,69,79,89,99))

#race/ethnicity
nhis$black<-Recode(nhis$RACEA, recodes="200=1; 900:990=NA; else=0")
nhis$white<-Recode(nhis$RACEA, recodes="100=1; 900:990=NA; else=0")
nhis$other<-Recode(nhis$RACEA, recodes="300:580=1; 900:990=NA; else=0")
nhis$multi<-Recode(nhis$RACEA, recodes="600:617=1; 900:990=NA; else=0")

###Recode binary outcome variable. Veteran = 1, Non-veteran = 0

nhis$VETSERVWHEN<-Recode(nhis$VETSERVWHEN, recodes = "00= NA; 90=NA; 20:80=1; 10=0")

Research question

Did more veterans see a doctor for Depression than non-veterans? Who did not have insurance coverage?

nhis$DPMHMDR<-Recode(nhis$DPMHMDR, recodes = "1= 'No';  2='Yes'; 0=NA", as.factor=T)
nhis$DPMHMDR<-relevel(nhis$DPMHMDR, ref='No')

nhis$REGION<-Recode(nhis$REGION, recodes = "1='Northeast'; 2='Midwest'; 3='South'; 4='West'; 8:9=NA", as.factor=T)
nhis$REGION<-relevel(nhis$REGION, ref='Northeast')

nhis$HINOTCOVE<-Recode(nhis$HINOTCOVE, recodes = "1= 'Yes'; 2='No'; 7:9= NA; 0=NA", as.factor=T)
nhis$HINOTCOVE<-relevel(nhis$HINOTCOVE, ref='Yes')

You drop all your cases here

library(dplyr)
mydata<-nhis%>%
  select( agec, male, black, white, other, multi, VETSERVWHEN, HINOTCOVE, DPMHMDR, REGION, SAMPWEIGHT, STRATA, PSU) %>%
  filter( complete.cases(agec, male, SAMPWEIGHT)) 

options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~PSU, strata=~STRATA, weights=~SAMPWEIGHT, data =mydata, nest = T )

#fit.logit<-svyglm(VETSERVWHEN~DPMHMDR+agec+black+white, design= des, family=binomial)

summary(fit.logit)

str(nhis$DPMHMDR)

##  Factor w/ 2 levels "No","Yes": NA NA NA NA NA NA NA NA NA NA ...

summary(nhis$DPMHMDR)

##      No     Yes    NA's 
##     171    1512 5988471

table(nhis$DPMHMDR, nhis$VETSERVWHEN)

##      
##          0    1
##   No   149   21
##   Yes 1314  186

fit.logit<-svyglm(HINOTCOVE~agec+male+black+white,
                  design= des,
                  family=binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(fit.logit)

## 
## Call:
## svyglm(formula = HINOTCOVE ~ agec + male + black + white, design = des, 
##     family = binomial)
## 
## Survey design:
## svydesign(ids = ~PSU, strata = ~STRATA, weights = ~SAMPWEIGHT, 
##     data = mydata, nest = T)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.105452   0.027054 -77.824   <2e-16 ***
## agec(19,29]  1.157177   0.013800  83.851   <2e-16 ***
## agec(29,39]  0.766432   0.011971  64.023   <2e-16 ***
## agec(39,49]  0.518687   0.012848  40.370   <2e-16 ***
## agec(49,59]  0.211149   0.015119  13.966   <2e-16 ***
## agec(59,69] -0.501136   0.020348 -24.628   <2e-16 ***
## agec(69,79] -2.945905   0.066619 -44.221   <2e-16 ***
## agec(79,89] -3.383211   0.111926 -30.227   <2e-16 ***
## maleMale     0.200461   0.008088  24.785   <2e-16 ***
## black       -0.030147   0.027259  -1.106    0.269    
## white       -0.222900   0.024701  -9.024   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1.299363)
## 
## Number of Fisher Scoring iterations: 7

#probit model
fit.probit<-svyglm(HINOTCOVE~agec+male+black+white,
                   design=des,
                   family=binomial(link= "probit"))

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(fit.probit)

## 
## Call:
## svyglm(formula = HINOTCOVE ~ agec + male + black + white, design = des, 
##     family = binomial(link = "probit"))
## 
## Survey design:
## svydesign(ids = ~PSU, strata = ~STRATA, weights = ~SAMPWEIGHT, 
##     data = mydata, nest = T)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.220534   0.014865 -82.109   <2e-16 ***
## agec(19,29]  0.641922   0.007702  83.344   <2e-16 ***
## agec(29,39]  0.414677   0.006503  63.769   <2e-16 ***
## agec(39,49]  0.275733   0.006883  40.058   <2e-16 ***
## agec(49,59]  0.109754   0.007920  13.858   <2e-16 ***
## agec(59,69] -0.248048   0.009912 -25.026   <2e-16 ***
## agec(69,79] -1.242575   0.023633 -52.577   <2e-16 ***
## agec(79,89] -1.389820   0.037539 -37.024   <2e-16 ***
## maleMale     0.103471   0.004381  23.617   <2e-16 ***
## black       -0.024169   0.015255  -1.584    0.113    
## white       -0.129490   0.013829  -9.364   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1.299748)
## 
## Number of Fisher Scoring iterations: 7

knitr::kable(data.frame(OR = exp(coef(fit.logit)), ci = exp(confint(fit.logit))))

	OR	ci.2.5..	ci.97.5..
(Intercept)	0.1217906	0.1155009	0.1284228
agec(19,29]	3.1809412	3.0960551	3.2681547
agec(29,39]	2.1520737	2.1021674	2.2031649
agec(39,49]	1.6798200	1.6380462	1.7226591
agec(49,59]	1.2350968	1.1990352	1.2722429
agec(59,69]	0.6058419	0.5821555	0.6304920
agec(69,79]	0.0525545	0.0461215	0.0598846
agec(79,89]	0.0339383	0.0272533	0.0422631
maleMale	1.2219658	1.2027477	1.2414910
black	0.9703030	0.9198234	1.0235529
white	0.8001945	0.7623780	0.8398870

myexp<-function(x){exp(x)}

stargazer(fit.logit, fit.probit,
          type = "text",
          style="demography",
          covariate.labels=c("Age 19-29, Age 30-39","Age 40-49","Age 50-59", "Age 60-69", "Age 70-79", "Age 80+", "Male", "Black","White"),
          ci=T, column.sep.width = "10pt")

HINOTCOVE survey-weighted survey-weighted logistic probit Model 1 Model 2

Age 19-29, Age 30-39 1.157*** 0.642 (1.130, 1.184) (0.627, 0.657) Age 40-49 0.766 0.415 (0.743, 0.790) (0.402, 0.427) Age 50-59 0.519 0.276 (0.494, 0.544) (0.262, 0.289) Age 60-69 0.211 0.110 (0.182, 0.241) (0.094, 0.125) Age 70-79 -0.501 -0.248 (-0.541, -0.461) (-0.267, -0.229) Age 80+ -2.946 -1.243 (-3.076, -2.815) (-1.289, -1.196) Male -3.383 -1.390 (-3.603, -3.164) (-1.463, -1.316) Black 0.200 0.103 (0.185, 0.216) (0.095, 0.112) White -0.030 -0.024 (-0.084, 0.023) (-0.054, 0.006) white -0.223 -0.129 (-0.271, -0.174) (-0.157, -0.102) Constant -2.105 -1.221*** (-2.158, -2.052) (-1.250, -1.191) N 2,027,180 2,027,180 Log Likelihood -975,834.400 -975,938.300 AIC 1,951,691.000 1,951,899.000

p < .05; p < .01; p < .001