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