library(car)
library(stargazer)

## 
## Please cite as:

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

##  R package version 5.2.1. 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(haven)
library(Hmisc)

## Loading required package: lattice

## Loading required package: Formula

## Loading required package: ggplot2

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:questionr':
## 
##     describe, wtd.mean, wtd.table, wtd.var

## The following object is masked from 'package:survey':
## 
##     deff

## The following objects are masked from 'package:base':
## 
##     format.pval, units

dat<-read_xpt("C:/Users/Monica/Documents/CSparks_StatsII/CNTY05.xpt")

a) My outcome variable will be employment. YES or NO

b) My research questions are as follows: How does being threatened or hurt by one’s intimate partner affect one’s employment status?

c) My predictor variable is 1) intimate partner threat and hurt

d) Covariates I will use are sex, marital status, education, and health care coverage

Rename variables

nams<-names(dat)
head(nams, n=10)

##  [1] "_CNTYNAM" "_STATE"   "FMONTH"   "IDATE"    "IMONTH"   "IDAY"    
##  [7] "IYEAR"    "INTVID"   "DISPCODE" "SEQNO"

myvariables<-tolower(names(dat))
names(dat)<-myvariables

Recode for employment status

dat$myemploy<-recode(dat$employ, recodes="1:2='1YESEmployed'; 3:8='0NotEmployed'; else=NA", as.factor.result=T)
dat$myemploy<-relevel(dat$myemploy, ref= '1YESEmployed')

Recode for intimate partner threat (predictor variable)

dat$threat<-recode (dat$ipvthrat, recodes = "7:9=NA; 1='1YESthreat'; 2='0NOthreat'")

Recode for intimate partner hurt (predictor variable)

dat$hurt<-recode(dat$ipvphhrt, recodes = "7:9=NA; 1='1YEShurt'; 2='0NOhurt'")

Recode for sex

dat$sexes<-recode(dat$sex, recodes = "2= '1female'; 1= '0male'")

Recode for marital status

dat$marst<-recode(dat$marital, recodes="1='Married'; 2='Divorced'; 3='Widowed'; 4='Separated'; 5='NevrMarried';6='Cohab'; else=NA", as.factor.result=T)
dat$marst<-relevel(dat$marst, ref='Married')

Recode for health care coverage

dat$hlthcare<-recode(dat$hlthplan, recodes = "7:9=NA; 1='0YEShplan'; 2='1NOhplan'")

Recode for education level

dat$educ<-recode(dat$educa, recodes="1:3='1Lesshs'; 4='2hsgrad'; 5='3somecol'; 6='4colgrad';9=NA", as.factor.result=T)
dat$educ<-relevel(dat$educ, ref='2hsgrad')

Create a subset

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize

## 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

dat$strata<-dat$`_ststr`
dat$weight<-dat$`_cntywt`
sub<-dat%>%
  select(myemploy, threat, hurt, sexes, marst, educ, hlthcare, weight, strata) %>%
  filter( complete.cases(.))

Create my survey design

options(survey.lonely.psu = "adjust")
mydesign<-svydesign(ids=~1, strata=~strata, weights=~weight, data =sub )

Examine percent of respondents with poor mental health by intimate partner threat, perform survey-corrected chi-square test

mysurveythreat<-svyby(formula = ~myemploy, by = ~threat, design = mydesign, FUN = svymean, na.rm=T)
svychisq(~myemploy+threat, design = mydesign)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~myemploy + threat, design = mydesign)
## F = 0.16034, ndf = 1, ddf = 29826, p-value = 0.6888

mysurveythreat

##                threat myemploy1YESEmployed myemploy0NotEmployed
## 0NOthreat   0NOthreat            0.6425252            0.3574748
## 1YESthreat 1YESthreat            0.6484555            0.3515445
##            se.myemploy1YESEmployed se.myemploy0NotEmployed
## 0NOthreat               0.00619055              0.00619055
## 1YESthreat              0.01342100              0.01342100

Examine percent of respondents with poor mental health by intimate partner hurt, perform survey-corrected chi-square test and plot it

mysurveyhurt<-svyby(formula = ~myemploy, by = ~hurt, design = mydesign, FUN = svymean, na.rm=T)
svychisq(~myemploy+hurt, design = mydesign)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~myemploy + hurt, design = mydesign)
## F = 3.8976, ndf = 1, ddf = 29826, p-value = 0.04836

mysurveyhurt

##              hurt myemploy1YESEmployed myemploy0NotEmployed
## 0NOhurt   0NOhurt            0.6391226            0.3608774
## 1YEShurt 1YEShurt            0.6671559            0.3328441
##          se.myemploy1YESEmployed se.myemploy0NotEmployed
## 0NOhurt              0.006267509             0.006267509
## 1YEShurt             0.012560535             0.012560535

calculate threat, hurt, bad health cross tabulation

bothsurveys<-svyby(formula = ~myemploy, by = ~threat+hurt, design = mydesign, FUN = svymean, na.rm=T)
bothsurveys

##                         threat     hurt myemploy1YESEmployed
## 0NOthreat.0NOhurt    0NOthreat  0NOhurt            0.6391418
## 1YESthreat.0NOhurt  1YESthreat  0NOhurt            0.6383489
## 0NOthreat.1YEShurt   0NOthreat 1YEShurt            0.7289230
## 1YESthreat.1YEShurt 1YESthreat 1YEShurt            0.6502103
##                     myemploy0NotEmployed se.myemploy1YESEmployed
## 0NOthreat.0NOhurt              0.3608582              0.00636871
## 1YESthreat.0NOhurt             0.3616511              0.03372137
## 0NOthreat.1YEShurt             0.2710770              0.02255089
## 1YESthreat.1YEShurt            0.3497897              0.01461394
##                     se.myemploy0NotEmployed
## 0NOthreat.0NOhurt                0.00636871
## 1YESthreat.0NOhurt               0.03372137
## 0NOthreat.1YEShurt               0.02255089
## 1YESthreat.1YEShurt              0.01461394

descriptive statistics

describe(sub$myemploy)

## sub$myemploy 
##        n  missing distinct 
##    29975        0        2 
##                                     
## Value      1YESEmployed 0NotEmployed
## Frequency         18001        11974
## Proportion        0.601        0.399

describe(sub$threat)

## sub$threat : SEX PARTNER EVER THREATENED VIOLENCE 
##        n  missing distinct 
##    29975        0        2 
##                                 
## Value       0NOthreat 1YESthreat
## Frequency       25209       4766
## Proportion      0.841      0.159

describe(sub$hurt)

## sub$hurt : SEX PARTNER EVER VIOLENT 
##        n  missing distinct 
##    29975        0        2 
##                             
## Value       0NOhurt 1YEShurt
## Frequency     24875     5100
## Proportion     0.83     0.17

describe(sub$sexes)

## sub$sexes : SEX 
##        n  missing distinct 
##    29975        0        2 
##                           
## Value        0male 1female
## Frequency    11818   18157
## Proportion   0.394   0.606

describe(sub$marst)

## sub$marst 
##        n  missing distinct 
##    29975        0        6 
##                                                                       
## Value          Married       Cohab    Divorced NevrMarried   Separated
## Frequency        16531         876        4405        4591         610
## Proportion       0.551       0.029       0.147       0.153       0.020
##                       
## Value          Widowed
## Frequency         2962
## Proportion       0.099

describe(sub$educ)

## sub$educ 
##        n  missing distinct 
##    29975        0        4 
##                                               
## Value       2hsgrad  1Lesshs 3somecol 4colgrad
## Frequency      8273     2469     8309    10924
## Proportion    0.276    0.082    0.277    0.364

describe(sub$hlthcare)

## sub$hlthcare : HAVE HEALTH CARE COVERAGE 
##        n  missing distinct 
##    29975        0        2 
##                               
## Value      0YEShplan  1NOhplan
## Frequency      26545      3430
## Proportion     0.886     0.114

Examine percent of respondents with poor health by intimate partner hurt, perform survey-corrected chi-square test and plot it

mysurvey1<-svyby(formula = ~myemploy, by = ~hurt, design = mydesign, FUN = svymean, na.rm=T)
svychisq(~myemploy+hurt, design = mydesign)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~myemploy + hurt, design = mydesign)
## F = 3.8976, ndf = 1, ddf = 29826, p-value = 0.04836

mysurvey1

##              hurt myemploy1YESEmployed myemploy0NotEmployed
## 0NOhurt   0NOhurt            0.6391226            0.3608774
## 1YEShurt 1YEShurt            0.6671559            0.3328441
##          se.myemploy1YESEmployed se.myemploy0NotEmployed
## 0NOhurt              0.006267509             0.006267509
## 1YEShurt             0.012560535             0.012560535

mysurvey2<-svyby(formula = ~myemploy, by = ~sexes, design = mydesign, FUN = svymean, na.rm=T)
svychisq(~myemploy+sexes, design = mydesign)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~myemploy + sexes, design = mydesign)
## F = 154.44, ndf = 1, ddf = 29826, p-value < 2.2e-16

mysurvey2

##           sexes myemploy1YESEmployed myemploy0NotEmployed
## 0male     0male            0.7160036            0.2839964
## 1female 1female            0.5733063            0.4266937
##         se.myemploy1YESEmployed se.myemploy0NotEmployed
## 0male               0.008274118             0.008274118
## 1female             0.007476035             0.007476035

mysurvey3<-svyby(formula = ~myemploy, by = ~marst, design = mydesign, FUN = svymean, na.rm=T)
svychisq(~myemploy+marst, design = mydesign)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~myemploy + marst, design = mydesign)
## F = 50.145, ndf = 4.2892e+00, ddf = 1.2793e+05, p-value < 2.2e-16

mysurvey3

##                   marst myemploy1YESEmployed myemploy0NotEmployed
## Married         Married            0.6766474            0.3233526
## Cohab             Cohab            0.7060758            0.2939242
## Divorced       Divorced            0.6539658            0.3460342
## NevrMarried NevrMarried            0.6318396            0.3681604
## Separated     Separated            0.6079811            0.3920189
## Widowed         Widowed            0.2247798            0.7752202
##             se.myemploy1YESEmployed se.myemploy0NotEmployed
## Married                  0.00680245              0.00680245
## Cohab                    0.03243026              0.03243026
## Divorced                 0.01625499              0.01625499
## NevrMarried              0.01500992              0.01500992
## Separated                0.05545342              0.05545342
## Widowed                  0.01564661              0.01564661

mysurvey4<-svyby(formula = ~myemploy, by = ~educ, design = mydesign, FUN = svymean, na.rm=T)
svychisq(~myemploy+educ, design = mydesign)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~myemploy + educ, design = mydesign)
## F = 74.727, ndf = 2.856, ddf = 85183.000, p-value < 2.2e-16

mysurvey4

##              educ myemploy1YESEmployed myemploy0NotEmployed
## 2hsgrad   2hsgrad            0.5935981            0.4064019
## 1Lesshs   1Lesshs            0.4336176            0.5663824
## 3somecol 3somecol            0.6330242            0.3669758
## 4colgrad 4colgrad            0.7391183            0.2608817
##          se.myemploy1YESEmployed se.myemploy0NotEmployed
## 2hsgrad              0.011199000             0.011199000
## 1Lesshs              0.024021582             0.024021582
## 3somecol             0.010632294             0.010632294
## 4colgrad             0.007707763             0.007707763

mysurvey5<-svyby(formula = ~myemploy, by = ~hlthcare, design = mydesign, FUN = svymean, na.rm=T)
svychisq(~myemploy+hlthcare, design = mydesign)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~myemploy + hlthcare, design = mydesign)
## F = 5.1098, ndf = 1, ddf = 29826, p-value = 0.0238

mysurvey5

##            hlthcare myemploy1YESEmployed myemploy0NotEmployed
## 0YEShplan 0YEShplan            0.6499711            0.3500289
## 1NOhplan   1NOhplan            0.6047895            0.3952105
##           se.myemploy1YESEmployed se.myemploy0NotEmployed
## 0YEShplan             0.005688612             0.005688612
## 1NOhplan              0.019582796             0.019582796

Create Logit Mode on employment

fit.mylogit<-svyglm(myemploy~hurt+sexes+marst+educ+hlthcare,design= mydesign, family=binomial)

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

fit.mylogit

## Stratified Independent Sampling design (with replacement)
## svydesign(ids = ~1, strata = ~strata, weights = ~weight, data = sub)
## 
## Call:  svyglm(formula = myemploy ~ hurt + sexes + marst + educ + hlthcare, 
##     design = mydesign, family = binomial)
## 
## Coefficients:
##      (Intercept)      hurt1YEShurt      sexes1female        marstCohab  
##         -0.80879          -0.26270           0.61087          -0.40297  
##    marstDivorced  marstNevrMarried    marstSeparated      marstWidowed  
##          0.02445           0.12588           0.06301           1.72885  
##      educ1Lesshs      educ3somecol      educ4colgrad  hlthcare1NOhplan  
##          0.67181          -0.13378          -0.59664           0.09265  
## 
## Degrees of Freedom: 29974 Total (i.e. Null);  29815 Residual
## Null Deviance:       39060 
## Residual Deviance: 36390     AIC: 33510

Present model coefficients

stargazer(fit.mylogit, type= "text", style= "demography", covariate.labels=c("Intercept","1YEShurt", "1female", "Divorced", "Widowed", "Separated", "Not married", "Cohabitating","Lesshs", "SomeHS","SomeColl", "CollGrad", "1NOhplan"))

## 
## ---------------------------------
##                       myemploy   
## ---------------------------------
## Intercept            -0.263***   
##                       (0.070)    
## 1YEShurt              0.611***   
##                       (0.052)    
## 1female               -0.403*    
##                       (0.159)    
## Divorced               0.024     
##                       (0.078)    
## Widowed                0.126     
##                       (0.073)    
## Separated              0.063     
##                       (0.219)    
## Not married           1.729***   
##                       (0.099)    
## Cohabitating          0.672***   
##                       (0.113)    
## Lesshs                -0.134*    
##                       (0.066)    
## SomeHS               -0.597***   
##                       (0.062)    
## SomeColl               0.093     
##                       (0.087)    
## CollGrad             -0.809***   
##                       (0.058)    
## N                      29,975    
## Log Likelihood      -16,744.780  
## AIC                  33,513.560  
## ---------------------------------
## *p < .05; **p < .01; ***p < .001

(REWRITE) For section above, those who were threatened and hurt were more likely to experience poor mental health than those who were not, my reference group for marital status is “married” so those who are divorced, widowed, separated are more likely to experience poor mental health than those who are married (my reference group). Those who are College graduates are less likely to experience poor mental health than those who have are HS graduates (my reference group).

Present model marginal effects in data frame

log.marg<-coef(fit.mylogit)*mean(dlogis(predict(fit.mylogit)), na.rm=T)
data.frame(m.logit=log.marg)

##                       m.logit
## (Intercept)      -0.169067674
## hurt1YEShurt     -0.054914726
## sexes1female      0.127696376
## marstCohab       -0.084235645
## marstDivorced     0.005111989
## marstNevrMarried  0.026313555
## marstSeparated    0.013172153
## marstWidowed      0.361396426
## educ1Lesshs       0.140433764
## educ3somecol     -0.027966249
## educ4colgrad     -0.124720699
## hlthcare1NOhplan  0.019367674

Nested model comparisons

Fit logit 1 for hurt and sexes

fit.logit1<-svyglm(myemploy~hurt+sexes, design=mydesign, family=binomial)

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

summary(fit.logit1)

## 
## Call:
## svyglm(formula = myemploy ~ hurt + sexes, design = mydesign, 
##     family = binomial)
## 
## Survey design:
## svydesign(ids = ~1, strata = ~strata, weights = ~weight, data = sub)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.90067    0.04134  -21.79  < 2e-16 ***
## hurt1YEShurt -0.24938    0.06443   -3.87 0.000109 ***
## sexes1female  0.65390    0.05157   12.68  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1.000003)
## 
## Number of Fisher Scoring iterations: 4

Fit logit 2 for threat, hurt , sex, and marital status

fit.logit2<-svyglm(myemploy~hurt+sexes+marst, design=mydesign, family=binomial)

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

summary(fit.logit2)

## 
## Call:
## svyglm(formula = myemploy ~ hurt + sexes + marst, design = mydesign, 
##     family = binomial)
## 
## Survey design:
## svydesign(ids = ~1, strata = ~strata, weights = ~weight, data = sub)
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.02377    0.04316 -23.722  < 2e-16 ***
## hurt1YEShurt     -0.22965    0.06890  -3.333 0.000860 ***
## sexes1female      0.59600    0.05226  11.405  < 2e-16 ***
## marstCohab       -0.16145    0.15626  -1.033 0.301508    
## marstDivorced     0.10372    0.08516   1.218 0.223248    
## marstNevrMarried  0.24895    0.07183   3.466 0.000529 ***
## marstSeparated    0.33886    0.23364   1.450 0.146968    
## marstWidowed      1.85869    0.09693  19.175  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 0.9995411)
## 
## Number of Fisher Scoring iterations: 4

Fit logit 3 for hurt, sex, marital status, education, and employment

fit.logit3<-svyglm(myemploy~hurt+sexes+marst+educ+hlthcare, design=mydesign, family= binomial)

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

summary(fit.logit3)

## 
## Call:
## svyglm(formula = myemploy ~ hurt + sexes + marst + educ + hlthcare, 
##     design = mydesign, family = binomial)
## 
## Survey design:
## svydesign(ids = ~1, strata = ~strata, weights = ~weight, data = sub)
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.80879    0.05842 -13.844  < 2e-16 ***
## hurt1YEShurt     -0.26270    0.06999  -3.753 0.000175 ***
## sexes1female      0.61087    0.05209  11.728  < 2e-16 ***
## marstCohab       -0.40297    0.15939  -2.528 0.011469 *  
## marstDivorced     0.02445    0.07775   0.315 0.753131    
## marstNevrMarried  0.12588    0.07329   1.717 0.085902 .  
## marstSeparated    0.06301    0.21910   0.288 0.773659    
## marstWidowed      1.72885    0.09900  17.463  < 2e-16 ***
## educ1Lesshs       0.67181    0.11346   5.921 3.23e-09 ***
## educ3somecol     -0.13378    0.06644  -2.013 0.044073 *  
## educ4colgrad     -0.59664    0.06216  -9.598  < 2e-16 ***
## hlthcare1NOhplan  0.09265    0.08731   1.061 0.288595    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1.002419)
## 
## Number of Fisher Scoring iterations: 4

Odds ratios and confidence values

stargazer(fit.logit1, fit.logit2, fit.logit3,type= "text", style= "demography", coviariate.labels=c("Intercept", "1Yeshurt","1female", "Divorced", "Widowed", "Separated", "Not married", "Cohabitating","PrimarySchool", "SomeHS","SomeColl", "CollGrad", "1NOhplan"))

## 
## ----------------------------------------------------
##                               myemploy              
##                    Model 1     Model 2     Model 3  
## ----------------------------------------------------
## hurt1YEShurt      -0.249***   -0.230***   -0.263*** 
##                    (0.064)     (0.069)     (0.070)  
## sexes1female      0.654***    0.596***    0.611***  
##                    (0.052)     (0.052)     (0.052)  
## marstCohab                     -0.161      -0.403*  
##                                (0.156)     (0.159)  
## marstDivorced                   0.104       0.024   
##                                (0.085)     (0.078)  
## marstNevrMarried              0.249***      0.126   
##                                (0.072)     (0.073)  
## marstSeparated                  0.339       0.063   
##                                (0.234)     (0.219)  
## marstWidowed                  1.859***    1.729***  
##                                (0.097)     (0.099)  
## educ1Lesshs                               0.672***  
##                                            (0.113)  
## educ3somecol                               -0.134*  
##                                            (0.066)  
## educ4colgrad                              -0.597*** 
##                                            (0.062)  
## hlthcare1NOhplan                            0.093   
##                                            (0.087)  
## Constant          -0.901***   -1.024***   -0.809*** 
##                    (0.041)     (0.043)     (0.058)  
## N                  29,975      29,975      29,975   
## Log Likelihood   -17,597.450 -17,144.230 -16,744.780
## AIC              35,200.910  34,304.460  33,513.560 
## ----------------------------------------------------
## *p < .05; **p < .01; ***p < .001                    
## 
## ------------------------------------------------------------------------------------------------------------------------------
## Intercept 1Yeshurt 1female Divorced Widowed Separated Not married Cohabitating PrimarySchool SomeHS SomeColl CollGrad 1NOhplan
## ------------------------------------------------------------------------------------------------------------------------------

Research Draft_CSparks

Monica Cruz

May 6, 2018

a) My outcome variable will be employment. YES or NO

b) My research questions are as follows: How does being threatened or hurt by one’s intimate partner affect one’s employment status?

c) My predictor variable is 1) intimate partner threat and hurt

d) Covariates I will use are sex, marital status, education, and health care coverage

Rename variables

Recode for employment status

Recode for intimate partner threat (predictor variable)

Recode for intimate partner hurt (predictor variable)

Recode for sex

Recode for marital status

Recode for health care coverage

Recode for education level

Create a subset

Create my survey design

Examine percent of respondents with poor mental health by intimate partner threat, perform survey-corrected chi-square test

Examine percent of respondents with poor mental health by intimate partner hurt, perform survey-corrected chi-square test and plot it

calculate threat, hurt, bad health cross tabulation

descriptive statistics

Examine percent of respondents with poor health by intimate partner hurt, perform survey-corrected chi-square test and plot it

Create Logit Mode on employment

Present model coefficients

Present model marginal effects in data frame

Nested model comparisons

Fit logit 1 for hurt and sexes

Fit logit 2 for threat, hurt , sex, and marital status

Fit logit 3 for hurt, sex, marital status, education, and employment

Odds ratios and confidence values