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)
dat<-read_xpt("C:/Users/Monica/Documents/CSparks_StatsII/CNTY05.xpt")
a) My ordinal outcome variable will be as the self-reporting of number of poor mental health days in the past month. I will recode the variable to define number of poor mental days into three categories based on the number of self-reported poor mental health days.
b) My research question is as follows: How do threats of violence and acts of violence from one’s intimate partner affect the number of self-reported poor mental health days? And how does marital status affect number of self-reported mental health days?
c) My predictor variables are 1) the threat of violence and 2) the act of violence and 3) marital status
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 poor mental health (outcome variable), with 1= low number of self-reported poor mental health days, 2= moderate number of self-reported poor mental health days, and 3= high number of self-reported poor mental health days
dat$poormenhlth<-recode (dat$menthlth,recodes= "1:10=1; 11:20= 2; 21:30= 3; else=NA", as.factor.result=T)
dat$poormenhlth<-relevel (dat$poormenhlth, ref = "1")
dat$menhlthdays<-car::recode(dat$menthlth, recodes="1:10=1; 11:20=2; 21:30= 3; else=NA", as.factor.result=F)
Recode for intimate partner threat (predictor variable)
dat$threat<-recode (dat$ipvthrat, recodes = "7:9=NA; 1='YESthreat'; 2='NOthreat'")
Recode for intimate partner hurt (predictor variable)
dat$hurt<-recode(dat$ipvphhrt, recodes = "7:9=NA; 1='YEShurt'; 2='NOhurt'")
Recode for marital status (predictor variable)
dat$marst<-recode(dat$marital, recodes="1='married'; 2='divorced'; 3='widowed'; 4='separated'; 5='nm';6='cohab'; else=NA", as.factor.result=T)
dat$marst<-relevel(dat$marst, ref='married')
Create survey design with subset
options (survey.lonely.psu= "adjust")
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
dat$strata<-dat$`_ststr`
dat$weight<-dat$`_cntywt`
mysub<-dat%>%
select(threat, hurt, poormenhlth, menhlthdays, marst, weight, strata)%>%
filter(complete.cases(.))
options (survey.lonely.psu= "adjust")
mydesign<-svydesign(ids = ~1, strata=~strata, weights=~weight, data = mysub)
Fit three nested models for poor mental health
fit.solr1<-svyolr(poormenhlth~threat + hurt + marst, mydesign)
summary (fit.solr1)
## Call:
## svyolr(poormenhlth ~ threat + hurt + marst, mydesign)
##
## Coefficients:
## Value Std. Error t value
## threatYESthreat 0.20196308 0.1528328 1.321464
## hurtYEShurt 0.23069595 0.1492775 1.545417
## marstcohab 0.31964765 0.2058000 1.553195
## marstdivorced 0.47110792 0.1126606 4.181656
## marstnm 0.01343992 0.1155695 0.116293
## marstseparated 1.33467781 0.2928758 4.557146
## marstwidowed 0.71328302 0.1506855 4.733588
##
## Intercepts:
## Value Std. Error t value
## 1|2 1.1894 0.0663 17.9444
## 2|3 2.0328 0.0762 26.6656
As seen above, individuals who reported YES to intimate partner threat have higher odds of self-reporting poor mental health than those who reported not experiencing intimate partner threat. Individuals who reported YES to intimate partner hurt have higher odds of self-reported poor mental health than those who reported not experiencing being hurt by their intimate partner. In addition, those who reported cohabitating, divorced, never married, separated, or widowed have higher odds of self-reporting poor mental health than those who were married, with those who are separated having even higher odds than the rest.
Calculate the AIC
fit.solr1$deviance+2*length(fit.solr1$coefficients)
## [1] 15959.43
Examine proportional odds assumption by fitting logits for each change
ex1<-svyglm(I(menhlthdays>1)~threat + hurt + marst, mydesign, family="binomial")
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
ex2<-svyglm(I(menhlthdays>2)~threat + hurt + marst, mydesign, family="binomial")
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
Plot the coefficients for each model
plot(coef(ex1)[-1], ylim=c(-3, 3), type="l",xaxt="n",
ylab="Beta", main=c("Comparison of betas for", " proportional odds assumption"))
lines(coef(ex2)[-1], col=1, lty=2)
axis(side=1, at=1:12, labels=F)
text(x=1:12, y=-4, srt = 90, pos = 1, xpd = TRUE,cex=.8,
labels = c( "Partner Threat", "Partner Hurt", "Married", "Divorced", "Widowed","Separated", "Never Married","Cohabitating"))
legend("bottomright",col=c(1,1),lty=c(1,2), legend=c(">1", ">2"))
lines(coef(fit.solr1)[c(-13:-16)], col=4, lwd=3)
Based on the results, the effects appear to be consistent across transitions.
Print the odds ratios
round(exp(rbind(coef(ex1)[-1], coef(ex2)[-1])),3)
## threatYESthreat hurtYEShurt marstcohab marstdivorced marstnm
## [1,] 1.198 1.295 1.381 1.619 1.040
## [2,] 1.382 1.108 1.314 1.511 0.874
## marstseparated marstwidowed
## [1,] 3.350 1.978
## [2,] 4.373 2.164
As provided in the odds ratios table above, those who experience intimate partner threat are approximately 20% more likely to report poor mental health days than those who have not experienced intimate partner threat. Those who have been hurt by their intimate partner are 30% more likely to report poor mental health days than those who have not been hurt by their intimate partner. Regarding marital status, those who reported cohabitation, divorced, never married, separated or widowed are more likely than married individuals to report poor mental health days, with those who are separated being 3 times more likely than those who are married to report poor mental health days.
Proportional Odds Model
library (VGAM)
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following object is masked from 'package:survey':
##
## calibrate
## The following object is masked from 'package:car':
##
## logit
fit.vgam<-vglm(as.ordered(poormenhlth)~threat + hurt + marst, dat, weights =weight/mean(weight, na.rm=T),
family=cumulative(parallel = T, reverse = T))
summary(fit.vgam)
##
## Call:
## vglm(formula = as.ordered(poormenhlth) ~ threat + hurt + marst,
## family = cumulative(parallel = T, reverse = T), data = dat,
## weights = weight/mean(weight, na.rm = T))
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logit(P[Y>=2]) -2.834 -0.3958 -0.2082 0.25383 9.40
## logit(P[Y>=3]) -4.408 -0.2535 -0.1290 -0.06641 10.53
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.18965 0.04044 -29.417 < 2e-16 ***
## (Intercept):2 -2.03262 0.04641 -43.798 < 2e-16 ***
## threatYESthreat 0.20190 0.10354 1.950 0.05117 .
## hurtYEShurt 0.23058 0.10219 2.256 0.02406 *
## marstcohab 0.31972 0.11213 2.851 0.00435 **
## marstdivorced 0.47151 0.08443 5.585 2.34e-08 ***
## marstnm 0.01353 0.06762 0.200 0.84142
## marstseparated 1.33528 0.14871 8.979 < 2e-16 ***
## marstwidowed 0.71342 0.12689 5.622 1.89e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: logit(P[Y>=2]), logit(P[Y>=3])
##
## Residual deviance: 10707.3 on 20609 degrees of freedom
##
## Log-likelihood: -5353.651 on 20609 degrees of freedom
##
## Number of iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
## Exponentiated coefficients:
## threatYESthreat hurtYEShurt marstcohab marstdivorced
## 1.223724 1.259324 1.376745 1.602419
## marstnm marstseparated marstwidowed
## 1.013621 3.801042 2.040955
Non-proportional odds
fit.vgam2<-vglm(as.ordered(poormenhlth)~threat + hurt + marst, dat, weights =weight/mean(weight, na.rm=T),
family=cumulative(parallel = F, reverse = T))
## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals
## = residuals, : fitted values close to 0 or 1
summary(fit.vgam2)
##
## Call:
## vglm(formula = as.ordered(poormenhlth) ~ threat + hurt + marst,
## family = cumulative(parallel = F, reverse = T), data = dat,
## weights = weight/mean(weight, na.rm = T))
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logit(P[Y>=2]) -2.891 -0.3962 -0.2082 0.25535 9.593
## logit(P[Y>=3]) -4.621 -0.2515 -0.1283 -0.06651 11.286
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.19682 0.04078 -29.350 < 2e-16 ***
## (Intercept):2 -1.99540 0.05263 -37.914 < 2e-16 ***
## threatYESthreat:1 0.16844 0.10555 1.596 0.110532
## threatYESthreat:2 0.32794 0.13012 2.520 0.011723 *
## hurtYEShurt:1 0.26907 0.10402 2.587 0.009691 **
## hurtYEShurt:2 0.09440 0.12949 0.729 0.465989
## marstcohab:1 0.33395 0.11407 2.928 0.003416 **
## marstcohab:2 0.26920 0.14308 1.881 0.059907 .
## marstdivorced:1 0.48472 0.08629 5.617 1.94e-08 ***
## marstdivorced:2 0.41227 0.10619 3.882 0.000103 ***
## marstnm:1 0.03987 0.06815 0.585 0.558518
## marstnm:2 -0.14306 0.09163 -1.561 0.118432
## marstseparated:1 1.19492 0.15934 7.499 6.42e-14 ***
## marstseparated:2 1.46078 0.16467 8.871 < 2e-16 ***
## marstwidowed:1 0.68641 0.13125 5.230 1.70e-07 ***
## marstwidowed:2 0.76785 0.15159 5.065 4.08e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: logit(P[Y>=2]), logit(P[Y>=3])
##
## Residual deviance: 10689.92 on 20602 degrees of freedom
##
## Log-likelihood: -5344.962 on 20602 degrees of freedom
##
## Number of iterations: 8
##
## No Hauck-Donner effect found in any of the estimates
##
## Exponentiated coefficients:
## threatYESthreat:1 threatYESthreat:2 hurtYEShurt:1 hurtYEShurt:2
## 1.1834589 1.3881037 1.3087494 1.0990025
## marstcohab:1 marstcohab:2 marstdivorced:1 marstdivorced:2
## 1.3964797 1.3089180 1.6237265 1.5102490
## marstnm:1 marstnm:2 marstseparated:1 marstseparated:2
## 1.0406738 0.8666994 3.3033086 4.3093195
## marstwidowed:1 marstwidowed:2
## 1.9865689 2.1551235
fit.null2<-vglm(as.ordered(poormenhlth)~1,
dat, weights =weight/mean(weight, na.rm=T),
family=cumulative(parallel = F, reverse = T))
(1-exp((fit.vgam2@criterion$deviance - fit.null2@criterion$deviance)/20609))/(1-exp(-fit.null2@criterion$deviance/20609))
## [1] 0.9926289
AIC(fit.vgam)
## [1] 10725.3
AIC(fit.vgam2)
## [1] 10721.92
Adjacent categories
fit.acat<-vglm(as.ordered(poormenhlth)~threat + hurt + marst, dat, weights = weight/mean(weight, na.rm=T),
family=acat(parallel = F, reverse = T))
summary(fit.acat)
##
## Call:
## vglm(formula = as.ordered(poormenhlth) ~ threat + hurt + marst,
## family = acat(parallel = F, reverse = T), data = dat, weights = weight/mean(weight,
## na.rm = T))
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## loge(P[Y=1]/P[Y=2]) -9.054 -0.31253 0.2150 0.4086 3.005
## loge(P[Y=2]/P[Y=3]) -10.842 0.05175 0.1012 0.2050 5.529
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 1.91530 0.05473 34.993 < 2e-16 ***
## (Intercept):2 -0.05090 0.07083 -0.719 0.4724
## threatYESthreat:1 -0.01971 0.13961 -0.141 0.8877
## threatYESthreat:2 -0.30657 0.17381 -1.764 0.0778 .
## hurtYEShurt:1 -0.35576 0.13620 -2.612 0.0090 **
## hurtYEShurt:2 0.18838 0.17133 1.100 0.2715
## marstcohab:1 -0.32037 0.15058 -2.128 0.0334 *
## marstcohab:2 -0.00439 0.18864 -0.023 0.9814
## marstdivorced:1 -0.46575 0.11351 -4.103 4.07e-05 ***
## marstdivorced:2 -0.03033 0.13981 -0.217 0.8283
## marstnm:1 -0.17720 0.08795 -2.015 0.0439 *
## marstnm:2 0.28543 0.11802 2.418 0.0156 *
## marstseparated:1 -0.54630 0.24470 -2.233 0.0256 *
## marstseparated:2 -1.02888 0.25120 -4.096 4.21e-05 ***
## marstwidowed:1 -0.46367 0.18278 -2.537 0.0112 *
## marstwidowed:2 -0.38707 0.21085 -1.836 0.0664 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: loge(P[Y=1]/P[Y=2]), loge(P[Y=2]/P[Y=3])
##
## Residual deviance: 10690.22 on 20602 degrees of freedom
##
## Log-likelihood: -5345.11 on 20602 degrees of freedom
##
## Number of iterations: 4
Multinomial Model
mfit<-vglm(poormenhlth~threat + hurt + marst, multinomial(refLevel = 1),dat,
weights=weight/mean(weight, na.rm=T))
summary(mfit)
##
## Call:
## vglm(formula = poormenhlth ~ threat + hurt + marst, family = multinomial(refLevel = 1),
## data = dat, weights = weight/mean(weight, na.rm = T))
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,2]/mu[,1]) -2.324 -0.3052 -0.1738 -0.09541 10.58
## log(mu[,3]/mu[,1]) -2.084 -0.3167 -0.1795 -0.09648 11.67
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.91530 0.05473 -34.993 < 2e-16 ***
## (Intercept):2 -1.86440 0.05331 -34.970 < 2e-16 ***
## threatYESthreat:1 0.01971 0.13961 0.141 0.8877
## threatYESthreat:2 0.32628 0.13413 2.433 0.0150 *
## hurtYEShurt:1 0.35576 0.13620 2.612 0.0090 **
## hurtYEShurt:2 0.16738 0.13345 1.254 0.2098
## marstcohab:1 0.32037 0.15058 2.128 0.0334 *
## marstcohab:2 0.32476 0.14586 2.227 0.0260 *
## marstdivorced:1 0.46575 0.11351 4.103 4.07e-05 ***
## marstdivorced:2 0.49608 0.10876 4.561 5.08e-06 ***
## marstnm:1 0.17720 0.08795 2.015 0.0439 *
## marstnm:2 -0.10822 0.09259 -1.169 0.2424
## marstseparated:1 0.54630 0.24470 2.233 0.0256 *
## marstseparated:2 1.57517 0.17359 9.074 < 2e-16 ***
## marstwidowed:1 0.46367 0.18278 2.537 0.0112 *
## marstwidowed:2 0.85074 0.15598 5.454 4.92e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
##
## Residual deviance: 10690.22 on 20602 degrees of freedom
##
## Log-likelihood: -5345.11 on 20602 degrees of freedom
##
## Number of iterations: 4
##
## Reference group is level 1 of the response
Calculate the odds rations and confidence intervals
round(exp(coef(mfit)), 3)
## (Intercept):1 (Intercept):2 threatYESthreat:1 threatYESthreat:2
## 0.147 0.155 1.020 1.386
## hurtYEShurt:1 hurtYEShurt:2 marstcohab:1 marstcohab:2
## 1.427 1.182 1.378 1.384
## marstdivorced:1 marstdivorced:2 marstnm:1 marstnm:2
## 1.593 1.642 1.194 0.897
## marstseparated:1 marstseparated:2 marstwidowed:1 marstwidowed:2
## 1.727 4.832 1.590 2.341
round(exp(confint(mfit)), 3)
## 2.5 % 97.5 %
## (Intercept):1 0.132 0.164
## (Intercept):2 0.140 0.172
## threatYESthreat:1 0.776 1.341
## threatYESthreat:2 1.065 1.802
## hurtYEShurt:1 1.093 1.864
## hurtYEShurt:2 0.910 1.536
## marstcohab:1 1.026 1.851
## marstcohab:2 1.040 1.842
## marstdivorced:1 1.275 1.990
## marstdivorced:2 1.327 2.032
## marstnm:1 1.005 1.418
## marstnm:2 0.748 1.076
## marstseparated:1 1.069 2.790
## marstseparated:2 3.438 6.790
## marstwidowed:1 1.111 2.275
## marstwidowed:2 1.725 3.179