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(sjPlot)
## Warning: package 'sjPlot' was built under R version 4.0.4
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
library(ggplot2)


load("~/ipumstor3/data.RData")
summary(data)
##      STRATA          PSU          PERWEIGHT       SAMPWEIGHT    
##  Min.   :6001   Min.   :1.000   Min.   :    0   Min.   :     0  
##  1st Qu.:6081   1st Qu.:1.000   1st Qu.: 2093   1st Qu.:     0  
##  Median :6158   Median :2.000   Median : 3356   Median :     0  
##  Mean   :6156   Mean   :1.507   Mean   : 3380   Mean   :  3380  
##  3rd Qu.:6238   3rd Qu.:2.000   3rd Qu.: 4394   3rd Qu.:  5466  
##  Max.   :6300   Max.   :2.000   Max.   :28756   Max.   :104188  
##     FWEIGHT           SEX           MARSTAT         ALCAMT       
##  Min.   :  897   Min.   :1.000   Min.   : 0.0   Min.   :  0.000  
##  1st Qu.: 1917   1st Qu.:1.000   1st Qu.:10.0   1st Qu.:  0.000  
##  Median : 3123   Median :2.000   Median :10.0   Median :  0.000  
##  Mean   : 3184   Mean   :1.516   Mean   :20.8   Mean   :  5.678  
##  3rd Qu.: 4184   3rd Qu.:2.000   3rd Qu.:50.0   3rd Qu.:  0.000  
##  Max.   :26845   Max.   :2.000   Max.   :99.0   Max.   :999.000  
##     CIGSDAY     
##  Min.   : 1.00  
##  1st Qu.:96.00  
##  Median :96.00  
##  Mean   :91.27  
##  3rd Qu.:96.00  
##  Max.   :99.00
#Using IPUMs data we use the count variable cigarettes per day since we also used alcohol consumed per day we did not need a offset variable
#research question, does marital status relate to smoking cigarettes
#results show that the best model is the poisson model fit2 because it has the lowest AIC
data$gender<-recode(data$SEX,
                    recodes="1=1; 2=0; else=NA")
#1=married, 2=widowed, 3=divorced, 4=seperated, 5=NM

data$mar<-recode(data$MARSTAT,
                 recodes="10:13 =1; 20=2; 30=2; 40=4; 50=5; else=NA")
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
sub<-data%>%
  select(gender, mar, MARSTAT,
         SEX,CIGSDAY, ALCAMT, FWEIGHT, PERWEIGHT, PSU,
         SAMPWEIGHT, STRATA) %>%
  filter( complete.cases(.))

#First we tell R our survey design
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1, strata=~STRATA, weights=~SAMPWEIGHT, data =sub )
#OR THE BRFSS, R GAVE ME A WARNING AND I NEEDED TO ADD:
#YOU MAY NOT NEED TO DO THIS!!!!
#First we tell R our survey design
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1, strata=~STRATA, 
               weights=~SAMPWEIGHT,
               data = data[is.na(data$SAMPWEIGHT)==F,])
svyhist(~CIGSDAY, des)

svyby(~CIGSDAY, ~mar+gender, des, svymean, na.rm=T)
##     mar gender  CIGSDAY        se
## 1.0   1      0 84.78112 0.4192275
## 2.0   2      0 79.48747 0.6204091
## 4.0   4      0 74.34465 1.7711725
## 5.0   5      0 81.80322 0.5862962
## 1.1   1      1 83.22563 0.4470258
## 2.1   2      1 71.93251 0.9536398
## 4.1   4      1 64.18817 2.5132948
## 5.1   5      1 78.05711 0.6367085
svyby(~CIGSDAY, ~ALCAMT, des, svymean, na.rm=T)
##     ALCAMT  CIGSDAY         se
## 0        0 90.55462  0.1769934
## 10      10 85.47080  0.4300409
## 20      20 79.78102  0.5350991
## 30      30 70.97386  0.9690191
## 40      40 67.08518  1.3840793
## 50      50 61.28640  2.0790308
## 60      60 56.09178  1.9782904
## 70      70 59.32912  4.7994932
## 80      80 55.55627  4.1955653
## 90      90 51.58078  7.7145131
## 100    100 57.75151  4.3975446
## 110    110 62.69831 28.8588833
## 120    120 58.67694  4.3756821
## 130    130 96.00000  0.0000000
## 140    140 69.92419 19.7054872
## 150    150 37.56220  9.6902262
## 160    160 96.00000  0.0000000
## 180    180 25.32328 12.7048950
## 200    200 52.81000 17.5101385
## 240    240 84.97549  7.5211834
## 250    250 40.00000  0.0000000
## 300    300 82.81882 11.4038980
## 350    350 20.00000  0.0000000
## 600    600  1.00000  0.0000000
## 870    870 10.00000  0.0000000
## 950    950 65.49671 27.0041024
## 997    997 96.00000  0.0000000
## 998    998 64.06579 28.3919526
## 999    999 72.37294  4.7130137
fit1<-svyglm(CIGSDAY~factor(mar)+factor(gender), design=des, family=poisson)
summary(fit1)
## 
## Call:
## svyglm(formula = CIGSDAY ~ factor(mar) + factor(gender), design = des, 
##     family = poisson)
## 
## Survey design:
## svydesign(ids = ~1, strata = ~STRATA, weights = ~SAMPWEIGHT, 
##     data = data[is.na(data$SAMPWEIGHT) == F, ])
## 
## Coefficients:
##                  Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      4.451252   0.004442 1002.098  < 2e-16 ***
## factor(mar)2    -0.095191   0.007793  -12.215  < 2e-16 ***
## factor(mar)4    -0.182648   0.021326   -8.565  < 2e-16 ***
## factor(mar)5    -0.050006   0.006567   -7.615 2.71e-14 ***
## factor(gender)1 -0.040979   0.005618   -7.294 3.08e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 12.88573)
## 
## Number of Fisher Scoring iterations: 4
round(exp(summary(fit1)$coef[-1,1]), 3)
##    factor(mar)2    factor(mar)4    factor(mar)5 factor(gender)1 
##           0.909           0.833           0.951           0.960
fit2<-glm(CIGSDAY~factor(mar)+factor(gender), data=data, family=poisson)
summary(fit2)
## 
## Call:
## glm(formula = CIGSDAY ~ factor(mar) + factor(gender), family = poisson, 
##     data = data)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -13.2115    0.3301    0.4611    0.7552    1.5290  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      4.5304653  0.0006775 6687.40   <2e-16 ***
## factor(mar)2    -0.0795048  0.0011940  -66.59   <2e-16 ***
## factor(mar)4    -0.1063431  0.0026095  -40.75   <2e-16 ***
## factor(mar)5    -0.0306387  0.0008925  -34.33   <2e-16 ***
## factor(gender)1 -0.0135545  0.0007986  -16.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 599581  on 71017  degrees of freedom
## Residual deviance: 593674  on 71013  degrees of freedom
##   (18958 observations deleted due to missingness)
## AIC: 1036483
## 
## Number of Fisher Scoring iterations: 4
fit3<-glm(CIGSDAY~factor(mar)+factor(gender)+factor(ALCAMT), data=data, family=quasipoisson)
summary(fit3)
## 
## Call:
## glm(formula = CIGSDAY ~ factor(mar) + factor(gender) + factor(ALCAMT), 
##     family = quasipoisson, data = data)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -13.4191    0.0413    0.2313    0.5941    9.4359  
## 
## Coefficients:
##                    Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)        4.558095   0.001584 2877.938  < 2e-16 ***
## factor(mar)2      -0.055006   0.002704  -20.341  < 2e-16 ***
## factor(mar)4      -0.072830   0.005890  -12.364  < 2e-16 ***
## factor(mar)5      -0.019484   0.002019   -9.651  < 2e-16 ***
## factor(gender)1    0.002034   0.001811    1.123 0.261315    
## factor(ALCAMT)10  -0.095045   0.003355  -28.328  < 2e-16 ***
## factor(ALCAMT)20  -0.175898   0.003689  -47.685  < 2e-16 ***
## factor(ALCAMT)30  -0.271161   0.005579  -48.600  < 2e-16 ***
## factor(ALCAMT)40  -0.353632   0.008154  -43.367  < 2e-16 ***
## factor(ALCAMT)50  -0.397382   0.011363  -34.972  < 2e-16 ***
## factor(ALCAMT)60  -0.473496   0.011744  -40.319  < 2e-16 ***
## factor(ALCAMT)70  -0.447303   0.028606  -15.637  < 2e-16 ***
## factor(ALCAMT)80  -0.518704   0.025462  -20.372  < 2e-16 ***
## factor(ALCAMT)90  -0.511140   0.049415  -10.344  < 2e-16 ***
## factor(ALCAMT)100 -0.507966   0.025511  -19.912  < 2e-16 ***
## factor(ALCAMT)110 -0.568319   0.218891   -2.596 0.009424 ** 
## factor(ALCAMT)120 -0.437940   0.026016  -16.833  < 2e-16 ***
## factor(ALCAMT)130  0.004219   0.230007    0.018 0.985365    
## factor(ALCAMT)140 -0.340621   0.193246   -1.763 0.077967 .  
## factor(ALCAMT)150 -1.068695   0.076996  -13.880  < 2e-16 ***
## factor(ALCAMT)160  0.023704   0.230009    0.103 0.917920    
## factor(ALCAMT)180 -1.603230   0.231211   -6.934 4.12e-12 ***
## factor(ALCAMT)200 -0.524177   0.091553   -5.725 1.04e-08 ***
## factor(ALCAMT)240 -0.292540   0.085559   -3.419 0.000628 ***
## factor(ALCAMT)250 -0.798419   0.356365   -2.240 0.025065 *  
## factor(ALCAMT)300 -0.389409   0.141973   -2.743 0.006093 ** 
## factor(ALCAMT)350 -1.544912   0.503911   -3.066 0.002171 ** 
## factor(ALCAMT)600 -4.540645   2.253543   -2.015 0.043920 *  
## factor(ALCAMT)870 -2.236026   0.712635   -3.138 0.001704 ** 
## factor(ALCAMT)950 -0.367679   0.160974   -2.284 0.022369 *  
## factor(ALCAMT)997  0.041824   0.115012    0.364 0.716123    
## factor(ALCAMT)998 -0.589837   0.218889   -2.695 0.007047 ** 
## factor(ALCAMT)999 -0.283721   0.028811   -9.848  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.078453)
## 
##     Null deviance: 599581  on 71017  degrees of freedom
## Residual deviance: 534882  on 70985  degrees of freedom
##   (18958 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.0.4
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich)
## Warning: package 'sandwich' was built under R version 4.0.4
coeftest(fit2, vcov=vcovHC(fit2, type="HC1",cluster="STRATA"))
## 
## z test of coefficients:
## 
##                   Estimate Std. Error   z value  Pr(>|z|)    
## (Intercept)      4.5304653  0.0013317 3401.9212 < 2.2e-16 ***
## factor(mar)2    -0.0795048  0.0033507  -23.7277 < 2.2e-16 ***
## factor(mar)4    -0.1063431  0.0085903  -12.3794 < 2.2e-16 ***
## factor(mar)5    -0.0306387  0.0019825  -15.4548 < 2.2e-16 ***
## factor(gender)1 -0.0135545  0.0018260   -7.4231 1.144e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
fit.nb1<-glm.nb(CIGSDAY~factor(mar),
              data=sub,
              weights=SAMPWEIGHT/mean(SAMPWEIGHT, na.rm=T))
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in glm.nb(CIGSDAY ~ factor(mar), data = sub, weights = SAMPWEIGHT/
## mean(SAMPWEIGHT, : alternation limit reached
fit.nb2<-glm.nb(CIGSDAY~factor(mar)+factor(gender)+factor(ALCAMT),
              data=sub,
              weights=SAMPWEIGHT/mean(SAMPWEIGHT, na.rm=T))
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in glm.nb(CIGSDAY ~ factor(mar) + factor(gender) + factor(ALCAMT), :
## alternation limit reached
#clx2(fit.nb2,cluster =sub$ststr)
tests1<-coeftest(fit.nb1, vcov=vcovHC(fit.nb2, type="HC1",cluster="STRATA"))
tests<-coeftest(fit.nb2, vcov=vcovHC(fit.nb2, type="HC1",cluster="STRATA"))
library(stargazer)

stargazer(fit.nb1, fit.nb2,style="demography", type = "text", t.auto=F,p.auto=F,coef=list(tests1[, 1],tests[,1]),  se =list(tests1[, 2], tests[, 2]), p=list(tests1[,4],tests[, 4])   )
## 
## ---------------------------------------------
##                             CIGSDAY          
##                      Model 1       Model 2   
## ---------------------------------------------
## factor(mar)2        -0.089***     -0.094***  
##                      (0.003)       (0.003)   
## factor(mar)4        -0.178***     -0.161***  
##                      (0.009)       (0.009)   
## factor(mar)5        -0.051***     -0.037***  
##                      (0.003)       (0.003)   
## factor(gender)1                     0.001    
##                                    (0.002)   
## factor(ALCAMT)10                  -0.029***  
##                                    (0.003)   
## factor(ALCAMT)20                  -0.098***  
##                                    (0.003)   
## factor(ALCAMT)30                  -0.209***  
##                                    (0.006)   
## factor(ALCAMT)40                  -0.265***  
##                                    (0.009)   
## factor(ALCAMT)50                  -0.355***  
##                                    (0.014)   
## factor(ALCAMT)60                  -0.445***  
##                                    (0.015)   
## factor(ALCAMT)70                  -0.386***  
##                                    (0.034)   
## factor(ALCAMT)80                  -0.447***  
##                                    (0.032)   
## factor(ALCAMT)90                  -0.528***  
##                                    (0.063)   
## factor(ALCAMT)100                 -0.411***  
##                                    (0.032)   
## factor(ALCAMT)110                  -0.327    
##                                    (0.194)   
## factor(ALCAMT)120                 -0.395***  
##                                    (0.031)   
## factor(ALCAMT)130                 0.061***   
##                                    (0.002)   
## factor(ALCAMT)140                  -0.256*   
##                                    (0.119)   
## factor(ALCAMT)150                 -0.826***  
##                                    (0.107)   
## factor(ALCAMT)160                 0.098***   
##                                    (0.002)   
## factor(ALCAMT)180                 -1.253***  
##                                    (0.208)   
## factor(ALCAMT)200                 -0.504***  
##                                    (0.142)   
## factor(ALCAMT)240                  -0.011    
##                                    (0.035)   
## factor(ALCAMT)250                 -0.654***  
##                                    (0.009)   
## factor(ALCAMT)300                  -0.054    
##                                    (0.057)   
## factor(ALCAMT)350                 -1.471***  
##                                    (0.002)   
## factor(ALCAMT)600                 -4.467***  
##                                    (0.002)   
## factor(ALCAMT)870                 -2.163***  
##                                    (0.002)   
## factor(ALCAMT)950                  -0.306    
##                                    (0.169)   
## factor(ALCAMT)997                 0.130***   
##                                    (0.007)   
## factor(ALCAMT)998                  -0.343    
##                                    (0.187)   
## factor(ALCAMT)999                 -0.186***  
##                                    (0.027)   
## Constant            4.431***      4.502***   
##                      (0.002)       (0.002)   
## N                    71,018        71,018    
## Log Likelihood    -888,686.100  -853,615.000 
## theta              878,673.400  1,503,626.000
## AIC               1,777,380.000 1,707,296.000
## ---------------------------------------------
## *p < .05; **p < .01; ***p < .001
AIC(fit.nb1)
## [1] 1777380
AIC(fit.nb2)
## [1] 1707296
AIC(fit1)
##        eff.p          AIC     deltabar 
## 1.898724e+02 1.379723e+06 4.746811e+01
AIC(fit2)
## [1] 1036483
AIC(fit3)
## [1] NA