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(foreign)
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
brfsss=read.xport("~/Desktop/DataFolder/brfss.XPT ")
brfsss$marital<-Recode(brfsss$MARITAL, recodes="1=1;5=0;else=NA",as.factor=T)
brfsss$sxorient<-Recode(brfsss$SXORIENT,recodes="1=0;2=1;else=NA",as.factor=T)
brfsss$income<-Recode(brfsss$INCOME2,recodes="1:5='<$35k';5:7='$35k<$75k';8='75k+';else=NA",as.factor=T)
brfsss$child<-Recode(brfsss$X_CHLDCNT,recodes="1='0';2='1';3='2';4='3';5:6='4+';else=NA",as.factor=T)
brfsss <- brfsss %>%
select(sxorient,
marital,
child,
income,
X_LLCPWT,
X_STRWT)
#Create a survey design
sub<-brfsss%>%
select(sxorient,marital,income,child,X_LLCPWT,X_STRWT) %>%
filter( complete.cases(.))
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1, strata=~X_STRWT, weights=~X_LLCPWT, data =sub )
#Nested model comparison. Sexual orientation status by marital status.
fit.logit1<-svyglm(sxorient~marital,design=des,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit1)
##
## Call:
## svyglm(formula = sxorient ~ marital, design = des, family = binomial)
##
## Survey design:
## svydesign(ids = ~1, strata = ~X_STRWT, weights = ~X_LLCPWT, data = sub)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.03153 0.05788 -52.38 <2e-16 ***
## marital1 -1.64089 0.08910 -18.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.000009)
##
## Number of Fisher Scoring iterations: 7
#Let's look at sexual orientation, marital status, and having children
fit.logit2<-svyglm(sxorient~marital+child,design=des,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit2)
##
## Call:
## svyglm(formula = sxorient ~ marital + child, design = des, family = binomial)
##
## Survey design:
## svydesign(ids = ~1, strata = ~X_STRWT, weights = ~X_LLCPWT, data = sub)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.83098 0.06219 -45.522 < 2e-16 ***
## marital1 -1.53098 0.08865 -17.269 < 2e-16 ***
## child1 -0.60214 0.15860 -3.797 0.000147 ***
## child2 -1.18927 0.17574 -6.767 1.32e-11 ***
## child3 -1.25039 0.25008 -5.000 5.74e-07 ***
## child4+ -1.48468 0.49350 -3.008 0.002626 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.015576)
##
## Number of Fisher Scoring iterations: 7
#Let's look at sexual orientation, marital status, children, and income
fit.logit3<-svyglm(sxorient~marital+child+income,design=des,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit3)
##
## Call:
## svyglm(formula = sxorient ~ marital + child + income, design = des,
## family = binomial)
##
## Survey design:
## svydesign(ids = ~1, strata = ~X_STRWT, weights = ~X_LLCPWT, data = sub)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.871261 0.087511 -32.810 < 2e-16 ***
## marital1 -1.576642 0.087463 -18.026 < 2e-16 ***
## child1 -0.604811 0.158565 -3.814 0.000137 ***
## child2 -1.192681 0.176233 -6.768 1.32e-11 ***
## child3 -1.238906 0.250251 -4.951 7.41e-07 ***
## child4+ -1.465625 0.494670 -2.963 0.003049 **
## income$35k<$75k -0.004003 0.115166 -0.035 0.972272
## income75k+ 0.171105 0.105344 1.624 0.104326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.006669)
##
## Number of Fisher Scoring iterations: 7
stargazer(fit.logit1,fit.logit2,fit.logit3,type = "text",style="demography",covariate.labels = c("Gay","1 kid","2 kids","3 kids","35-75","75+"),ci=T)
##
## -----------------------------------------------------------------
## sxorient
## Model 1 Model 2 Model 3
## -----------------------------------------------------------------
## Gay -1.641*** -1.531*** -1.577***
## (-1.816, -1.466) (-1.705, -1.357) (-1.748, -1.405)
## 1 kid -0.602*** -0.605***
## (-0.913, -0.291) (-0.916, -0.294)
## 2 kids -1.189*** -1.193***
## (-1.534, -0.845) (-1.538, -0.847)
## 3 kids -1.250*** -1.239***
## (-1.741, -0.760) (-1.729, -0.748)
## 35-75 -1.485** -1.466**
## (-2.452, -0.517) (-2.435, -0.496)
## 75+ -0.004
## (-0.230, 0.222)
## income75k+ 0.171
## (-0.035, 0.378)
## Constant -3.032*** -2.831*** -2.871***
## (-3.145, -2.918) (-2.953, -2.709) (-3.043, -2.700)
## N 112,656 112,656 112,656
## Log Likelihood -9,715.172 -9,543.806 -9,539.155
## AIC 19,434.340 19,099.610 19,094.310
## -----------------------------------------------------------------
## *p < .05; **p < .01; ***p < .001
#res<-data.frame(model1=c(round(coef(fit.logit1),3),rep("-",7)) ,
#model2=c(round(coef(fit.logit2),3), rep("-",3)),
#model3=round(coef(fit.logit3),3))
#row.names(res)<-c("Intercept", "Hispanic", "Black","MultRace" ,"Other", "Some HS","HS Graduate", "Some College", "College #Grad", "Insurance", "Current Smoker", "Former Smoker")
#knitr::kable(res, caption="Model Coefficients in nested Logit Models", align="c")
#stratified
fit.logitint<-svyglm(sxorient~marital*child+income,design=des,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
regTermTest(fit.logitint,test.terms = ~marital:child,method = "Wald",df=NULL)
## Wald test for marital:child
## in svyglm(formula = sxorient ~ marital * child + income, design = des,
## family = binomial)
## F = 1.421753 on 4 and 111673 df: p= 0.22378
fit.unrestricted<-svyglm(sxorient~marital+income,design=des,family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
fit.logit.0kids<-svyglm(sxorient~(marital+income),design=subset(des,child=='0'),family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit.0kids)
##
## Call:
## svyglm(formula = sxorient ~ (marital + income), design = subset(des,
## child == "0"), family = binomial)
##
## Survey design:
## subset(des, child == "0")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.84811 0.09286 -30.672 <2e-16 ***
## marital1 -1.62036 0.09598 -16.882 <2e-16 ***
## income$35k<$75k -0.03184 0.12849 -0.248 0.804
## income75k+ 0.16602 0.11631 1.427 0.153
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 0.9928583)
##
## Number of Fisher Scoring iterations: 7
fit.logit.1kid<-svyglm(sxorient~(marital+income),design=subset(des,child=='1'),family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit.1kid)
##
## Call:
## svyglm(formula = sxorient ~ (marital + income), design = subset(des,
## child == "1"), family = binomial)
##
## Survey design:
## subset(des, child == "1")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.49364 0.26682 -13.094 < 2e-16 ***
## marital1 -1.27271 0.28103 -4.529 5.98e-06 ***
## income$35k<$75k -0.07802 0.34769 -0.224 0.822
## income75k+ -0.05511 0.33036 -0.167 0.868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.004135)
##
## Number of Fisher Scoring iterations: 7
fit.logit.2kids<-svyglm(sxorient~(marital+income),design=subset(des,child=='2'),family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit.2kids)
##
## Call:
## svyglm(formula = sxorient ~ (marital + income), design = subset(des,
## child == "2"), family = binomial)
##
## Survey design:
## subset(des, child == "2")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.3761 0.3264 -13.407 < 2e-16 ***
## marital1 -2.2270 0.3719 -5.989 2.18e-09 ***
## income$35k<$75k 0.8520 0.4809 1.772 0.0765 .
## income75k+ 1.0855 0.4611 2.354 0.0186 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.032731)
##
## Number of Fisher Scoring iterations: 8
fit.logit.3kids<-svyglm(sxorient~(marital+income),design=subset(des,child=='3'),family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit.3kids)
##
## Call:
## svyglm(formula = sxorient ~ (marital + income), design = subset(des,
## child == "3"), family = binomial)
##
## Survey design:
## subset(des, child == "3")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.3102 0.4014 -10.738 <2e-16 ***
## marital1 -0.7065 0.5032 -1.404 0.160
## income$35k<$75k -0.1508 0.5242 -0.288 0.774
## income75k+ -0.7365 0.7513 -0.980 0.327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.007884)
##
## Number of Fisher Scoring iterations: 8
fit.logit.4kids<-svyglm(sxorient~(marital+income),design=subset(des,child=='4+'),family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit.4kids)
##
## Call:
## svyglm(formula = sxorient ~ (marital + income), design = subset(des,
## child == "4+"), family = binomial)
##
## Survey design:
## subset(des, child == "4+")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.59176 0.71368 -7.835 7.41e-15 ***
## marital1 -0.15241 1.03463 -0.147 0.883
## income$35k<$75k -0.05871 1.18339 -0.050 0.960
## income75k+ 0.96901 1.04981 0.923 0.356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 0.9917282)
##
## Number of Fisher Scoring iterations: 8
fit.logit.4kids<-svyglm(sxorient~(marital+income),design=subset(des,child=='4+'),family=binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
summary(fit.logit.4kids)
##
## Call:
## svyglm(formula = sxorient ~ (marital + income), design = subset(des,
## child == "4+"), family = binomial)
##
## Survey design:
## subset(des, child == "4+")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.59176 0.71368 -7.835 7.41e-15 ***
## marital1 -0.15241 1.03463 -0.147 0.883
## income$35k<$75k -0.05871 1.18339 -0.050 0.960
## income75k+ 0.96901 1.04981 0.923 0.356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 0.9917282)
##
## Number of Fisher Scoring iterations: 8
#beta.test<-function(model1, model2, betaname){
#s1<-summary(model1)$coef
#s2<-summary(model2)$coef
#db <- ((s2[rownames(s2)==betaname,1]-s1[rownames(s1)==betaname,1]))^2
#sd <-s2[rownames(s2)==betaname,2]^2+s1[rownames(s1)==betaname,2]^2
#td <- db/sd
#beta1=s1[rownames(s1)==betaname,1]
#beta2=s2[rownames(s2)==betaname,1]
#pv<-1- pchisq(td, df = 1)
#print(list(beta=betaname,beta1=beta1, beta2=beta2, x2=td, pvalue=pv))
#}
#beta.test(fit.logit.0kids,fit.logit.1kid,fit.logit.2kids,fit.logit.3kids,,fit.logit.4kids,'marital')
#construct a test of whether the betas are the same for each race group
#See Allison p 217 for this
#deviance from total model
d1<-fit.unrestricted$deviance
#sum of deviances from cause-specific models
otherds<- (fit.logit.0kids$deviance+fit.logit.1kid$deviance+ fit.logit.2kids$deviance+fit.logit.3kids$deviance+fit.logit.4kids$deviance)
#Chow test
test<- d1-otherds
df<-(length(coef(fit.logit.0kids))*5)-length(coef(fit.unrestricted))
#print the test results
print(list(chowtest=test, df=df,pval= pchisq(test, df=df, lower=F)))
## $chowtest
## [1] -1019.542
##
## $df
## [1] 16
##
## $pval
## [1] 1
log.marg1<-coef(fit.logit1)*mean(dlogis(predict(fit.logit1)), na.rm=T)
log.marg2<-coef(fit.logit2)*mean(dlogis(predict(fit.logit2)), na.rm=T)
log.marg3<-coef(fit.logit3)*mean(dlogis(predict(fit.logit3)), na.rm=T)
log.marg1
## (Intercept) marital1
## -0.05154990 -0.02790266
log.marg2
## (Intercept) marital1 child1 child2 child3 child4+
## -0.05112369 -0.02764749 -0.01087382 -0.02147657 -0.02258043 -0.02681128
log.marg3
## (Intercept) marital1 child1 child2
## -5.165119e-02 -2.836226e-02 -1.087996e-02 -2.145517e-02
## child3 child4+ income$35k<$75k income75k+
## -2.228672e-02 -2.636517e-02 -7.201144e-05 3.078017e-03