Data and Examples from Cameron and Trivedi (1998)
CameronTrivedi1998.RdThis manual page collects a list of examples from the book. Some solutions might not be exact and the list is certainly not complete. If you have suggestions for improvement (preferably in the form of code), please contact the package maintainer.
References
Cameron, A.C. and Trivedi, P.K. (1998). Regression Analysis of Count Data. Cambridge: Cambridge University Press.
Examples
# \donttest{
library("MASS")
library("pscl")
###########################################
## Australian health service utilization ##
###########################################
## data
data("DoctorVisits", package = "AER")
## Poisson regression
dv_pois <- glm(visits ~ . + I(age^2), data = DoctorVisits, family = poisson)
dv_qpois <- glm(visits ~ . + I(age^2), data = DoctorVisits, family = quasipoisson)
## Table 3.3
round(cbind(
Coef = coef(dv_pois),
MLH = sqrt(diag(vcov(dv_pois))),
MLOP = sqrt(diag(vcovOPG(dv_pois))),
NB1 = sqrt(diag(vcov(dv_qpois))),
RS = sqrt(diag(sandwich(dv_pois)))
), digits = 3)
#> Coef MLH MLOP NB1 RS
#> (Intercept) -2.224 0.190 0.144 0.219 0.254
#> genderfemale 0.157 0.056 0.041 0.065 0.079
#> age 1.056 1.001 0.750 1.153 1.364
#> income -0.205 0.088 0.062 0.102 0.129
#> illness 0.187 0.018 0.014 0.021 0.024
#> reduced 0.127 0.005 0.004 0.006 0.008
#> health 0.030 0.010 0.007 0.012 0.014
#> privateyes 0.123 0.072 0.056 0.083 0.095
#> freepooryes -0.440 0.180 0.116 0.207 0.290
#> freerepatyes 0.080 0.092 0.070 0.106 0.126
#> nchronicyes 0.114 0.067 0.051 0.077 0.091
#> lchronicyes 0.141 0.083 0.059 0.096 0.123
#> I(age^2) -0.849 1.078 0.809 1.242 1.460
## Table 3.4
## NM2-ML
dv_nb <- glm.nb(visits ~ . + I(age^2), data = DoctorVisits)
summary(dv_nb)
#>
#> Call:
#> glm.nb(formula = visits ~ . + I(age^2), data = DoctorVisits,
#> init.theta = 0.9284725333, link = log)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -2.190007 0.233592 -9.375 < 2e-16 ***
#> genderfemale 0.216644 0.069697 3.108 0.00188 **
#> age -0.216159 1.266701 -0.171 0.86450
#> income -0.142202 0.108417 -1.312 0.18965
#> illness 0.214341 0.023579 9.090 < 2e-16 ***
#> reduced 0.143754 0.007311 19.662 < 2e-16 ***
#> health 0.038060 0.013654 2.788 0.00531 **
#> privateyes 0.118064 0.085806 1.376 0.16884
#> freepooryes -0.496611 0.210803 -2.356 0.01848 *
#> freerepatyes 0.144982 0.115970 1.250 0.21124
#> nchronicyes 0.099355 0.079303 1.253 0.21026
#> lchronicyes 0.190327 0.104357 1.824 0.06818 .
#> I(age^2) 0.609158 1.383245 0.440 0.65966
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> (Dispersion parameter for Negative Binomial(0.9285) family taken to be 1)
#>
#> Null deviance: 3928.7 on 5189 degrees of freedom
#> Residual deviance: 3028.3 on 5177 degrees of freedom
#> AIC: 6425.5
#>
#> Number of Fisher Scoring iterations: 1
#>
#>
#> Theta: 0.9285
#> Std. Err.: 0.0864
#>
#> 2 x log-likelihood: -6397.4880
## NB1-GLM = quasipoisson
summary(dv_qpois)
#>
#> Call:
#> glm(formula = visits ~ . + I(age^2), family = quasipoisson, data = DoctorVisits)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -2.223848 0.218725 -10.167 < 2e-16 ***
#> genderfemale 0.156882 0.064686 2.425 0.01533 *
#> age 1.056299 1.153198 0.916 0.35972
#> income -0.205321 0.101839 -2.016 0.04384 *
#> illness 0.186948 0.021065 8.875 < 2e-16 ***
#> reduced 0.126846 0.005801 21.868 < 2e-16 ***
#> health 0.030081 0.011637 2.585 0.00977 **
#> privateyes 0.123185 0.082551 1.492 0.13570
#> freepooryes -0.440061 0.207197 -2.124 0.03373 *
#> freerepatyes 0.079798 0.106081 0.752 0.45194
#> nchronicyes 0.114085 0.076789 1.486 0.13742
#> lchronicyes 0.141158 0.095808 1.473 0.14072
#> I(age^2) -0.848704 1.241930 -0.683 0.49440
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> (Dispersion parameter for quasipoisson family taken to be 1.327793)
#>
#> Null deviance: 5634.8 on 5189 degrees of freedom
#> Residual deviance: 4379.5 on 5177 degrees of freedom
#> AIC: NA
#>
#> Number of Fisher Scoring iterations: 6
#>
## overdispersion tests (page 79)
lrtest(dv_pois, dv_nb) ## p-value would need to be halved
#> Likelihood ratio test
#>
#> Model 1: visits ~ gender + age + income + illness + reduced + health +
#> private + freepoor + freerepat + nchronic + lchronic + I(age^2)
#> Model 2: visits ~ gender + age + income + illness + reduced + health +
#> private + freepoor + freerepat + nchronic + lchronic + I(age^2)
#> #Df LogLik Df Chisq Pr(>Chisq)
#> 1 13 -3355.5
#> 2 14 -3198.7 1 313.6 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dispersiontest(dv_pois, trafo = 1)
#>
#> Overdispersion test
#>
#> data: dv_pois
#> z = 6.5428, p-value = 3.019e-11
#> alternative hypothesis: true alpha is greater than 0
#> sample estimates:
#> alpha
#> 0.4144272
#>
dispersiontest(dv_pois, trafo = 2)
#>
#> Overdispersion test
#>
#> data: dv_pois
#> z = 7.5046, p-value = 3.08e-14
#> alternative hypothesis: true alpha is greater than 0
#> sample estimates:
#> alpha
#> 0.9574298
#>
##########################################
## Demand for medical care in NMES 1988 ##
##########################################
## select variables for analysis
data("NMES1988", package = "AER")
nmes <- NMES1988[,-(2:6)]
## dependent variable
## Table 6.1
table(cut(nmes$visits, c(0:13, 100)-0.5, labels = 0:13))
#>
#> 0 1 2 3 4 5 6 7 8 9 10 11 12 13
#> 683 481 428 420 383 338 268 217 188 171 128 115 86 500
## NegBin regression
nmes_nb <- glm.nb(visits ~ ., data = nmes)
## NegBin hurdle
nmes_h <- hurdle(visits ~ ., data = nmes, dist = "negbin")
## from Table 6.3
lrtest(nmes_nb, nmes_h)
#> Warning: original model was of class "negbin", updated model is of class "hurdle"
#> Likelihood ratio test
#>
#> Model 1: visits ~ health + chronic + adl + region + age + afam + gender +
#> married + school + income + employed + insurance + medicaid
#> Model 2: visits ~ .
#> #Df LogLik Df Chisq Pr(>Chisq)
#> 1 18 -12202
#> 2 35 -12110 17 183.35 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## from Table 6.4
AIC(nmes_nb)
#> [1] 24440.34
AIC(nmes_nb, k = log(nrow(nmes)))
#> [1] 24555.37
AIC(nmes_h)
#> [1] 24290.98
AIC(nmes_h, k = log(nrow(nmes)))
#> [1] 24514.66
## Table 6.8
coeftest(nmes_h, vcov = sandwich)
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> count_(Intercept) 1.6309834 0.2710457 6.0174 1.918e-09 ***
#> count_healthpoor 0.3325087 0.0567090 5.8634 4.869e-09 ***
#> count_healthexcellent -0.3775071 0.0875557 -4.3116 1.656e-05 ***
#> count_chronic 0.1429373 0.0135873 10.5199 < 2.2e-16 ***
#> count_adllimited 0.1290354 0.0515406 2.5036 0.012331 *
#> count_regionnortheast 0.1040669 0.0527117 1.9743 0.048414 *
#> count_regionmidwest -0.0163183 0.0475000 -0.3435 0.731207
#> count_regionwest 0.1232470 0.0504222 2.4443 0.014553 *
#> count_age -0.0753010 0.0322000 -2.3385 0.019404 *
#> count_afamyes 0.0016161 0.0700041 0.0231 0.981583
#> count_gendermale 0.0041273 0.0421279 0.0980 0.921960
#> count_marriedyes -0.0920323 0.0436135 -2.1102 0.034899 *
#> count_school 0.0216106 0.0056511 3.8242 0.000133 ***
#> count_income -0.0022357 0.0058893 -0.3796 0.704241
#> count_employedyes 0.0296559 0.0739627 0.4010 0.688471
#> count_insuranceyes 0.2271511 0.0566849 4.0073 6.245e-05 ***
#> count_medicaidyes 0.1847927 0.0665406 2.7771 0.005507 **
#> zero_(Intercept) -1.4753118 0.6463277 -2.2826 0.022501 *
#> zero_healthpoor 0.0708379 0.1687129 0.4199 0.674599
#> zero_healthexcellent -0.3285110 0.1422327 -2.3097 0.020953 *
#> zero_chronic 0.5565120 0.0527626 10.5475 < 2.2e-16 ***
#> zero_adllimited -0.1881658 0.1299284 -1.4482 0.147625
#> zero_regionnortheast 0.1292212 0.1250363 1.0335 0.301441
#> zero_regionmidwest 0.1008883 0.1146224 0.8802 0.378810
#> zero_regionwest 0.2016633 0.1336291 1.5091 0.131339
#> zero_age 0.1904976 0.0811377 2.3478 0.018927 *
#> zero_afamyes -0.3269720 0.1334512 -2.4501 0.014320 *
#> zero_gendermale -0.4644473 0.0985088 -4.7148 2.495e-06 ***
#> zero_marriedyes 0.2472641 0.1039403 2.3789 0.017407 *
#> zero_school 0.0542073 0.0131932 4.1087 4.051e-05 ***
#> zero_income 0.0067446 0.0184949 0.3647 0.715373
#> zero_employedyes -0.0123197 0.1450825 -0.0849 0.932333
#> zero_insuranceyes 0.7624604 0.1172920 6.5005 8.901e-11 ***
#> zero_medicaidyes 0.5535139 0.1812055 3.0546 0.002267 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
logLik(nmes_h)
#> 'log Lik.' -12110.49 (df=35)
1/nmes_h$theta
#> count
#> 0.7437966
###################################################
## Recreational boating trips to Lake Somerville ##
###################################################
## data
data("RecreationDemand", package = "AER")
## Poisson model:
## Cameron and Trivedi (1998), Table 6.11
## Ozuna and Gomez (1995), Table 2, col. 3
fm_pois <- glm(trips ~ ., data = RecreationDemand, family = poisson)
summary(fm_pois)
#>
#> Call:
#> glm(formula = trips ~ ., family = poisson, data = RecreationDemand)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.264993 0.093722 2.827 0.00469 **
#> quality 0.471726 0.017091 27.602 < 2e-16 ***
#> skiyes 0.418214 0.057190 7.313 2.62e-13 ***
#> income -0.111323 0.019588 -5.683 1.32e-08 ***
#> userfeeyes 0.898165 0.078985 11.371 < 2e-16 ***
#> costC -0.003430 0.003118 -1.100 0.27131
#> costS -0.042536 0.001670 -25.467 < 2e-16 ***
#> costH 0.036134 0.002710 13.335 < 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: 4849.7 on 658 degrees of freedom
#> Residual deviance: 2305.8 on 651 degrees of freedom
#> AIC: 3074.9
#>
#> Number of Fisher Scoring iterations: 7
#>
logLik(fm_pois)
#> 'log Lik.' -1529.431 (df=8)
coeftest(fm_pois, vcov = sandwich)
#>
#> z test of coefficients:
#>
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.2649934 0.4324810 0.6127 0.5400559
#> quality 0.4717259 0.0488508 9.6565 < 2.2e-16 ***
#> skiyes 0.4182137 0.1938713 2.1572 0.0309922 *
#> income -0.1113232 0.0503083 -2.2128 0.0269101 *
#> userfeeyes 0.8981653 0.2469086 3.6376 0.0002751 ***
#> costC -0.0034297 0.0146973 -0.2334 0.8154852
#> costS -0.0425364 0.0117348 -3.6248 0.0002892 ***
#> costH 0.0361336 0.0093860 3.8497 0.0001183 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
## Negbin model:
## Cameron and Trivedi (1998), Table 6.11
## Ozuna and Gomez (1995), Table 2, col. 5
library("MASS")
fm_nb <- glm.nb(trips ~ ., data = RecreationDemand)
coeftest(fm_nb, vcov = vcovOPG)
#>
#> z test of coefficients:
#>
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -1.1219363 0.1909098 -5.8768 4.183e-09 ***
#> quality 0.7219990 0.0399627 18.0668 < 2.2e-16 ***
#> skiyes 0.6121388 0.1395255 4.3873 1.148e-05 ***
#> income -0.0260588 0.0401183 -0.6495 0.516
#> userfeeyes 0.6691676 0.4488554 1.4908 0.136
#> costC 0.0480087 0.0103573 4.6353 3.565e-06 ***
#> costS -0.0926910 0.0060193 -15.3990 < 2.2e-16 ***
#> costH 0.0388357 0.0087604 4.4331 9.288e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
logLik(fm_nb)
#> 'log Lik.' -825.5576 (df=9)
## ZIP model:
## Cameron and Trivedi (1998), Table 6.11
fm_zip <- zeroinfl(trips ~ . | quality + income, data = RecreationDemand)
summary(fm_zip)
#>
#> Call:
#> zeroinfl(formula = trips ~ . | quality + income, data = RecreationDemand)
#>
#> Pearson residuals:
#> Min 1Q Median 3Q Max
#> -6.3255 -0.2714 -0.1809 -0.1646 13.3126
#>
#> Count model coefficients (poisson with log link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.099163 0.111397 18.844 < 2e-16 ***
#> quality 0.033833 0.023914 1.415 0.157
#> skiyes 0.471691 0.058187 8.106 5.21e-16 ***
#> income -0.099780 0.020779 -4.802 1.57e-06 ***
#> userfeeyes 0.610488 0.079435 7.685 1.53e-14 ***
#> costC 0.002369 0.003818 0.620 0.535
#> costS -0.037600 0.002038 -18.454 < 2e-16 ***
#> costH 0.025234 0.003355 7.522 5.40e-14 ***
#>
#> Zero-inflation model coefficients (binomial with logit link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.29191 0.51608 6.379 1.79e-10 ***
#> quality -1.91407 0.20619 -9.283 < 2e-16 ***
#> income -0.04502 0.10797 -0.417 0.677
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Number of iterations in BFGS optimization: 23
#> Log-likelihood: -1181 on 11 Df
logLik(fm_zip)
#> 'log Lik.' -1180.795 (df=11)
## Hurdle models
## Cameron and Trivedi (1998), Table 6.13
## poisson-poisson
sval <- list(count = c(2.15, 0.044, .467, -.097, .601, .002, -.036, .024),
zero = c(-1.88, 0.815, .403, .01, 2.95, 0.006, -.052, .046))
fm_hp0 <- hurdle(trips ~ ., data = RecreationDemand, dist = "poisson",
zero = "poisson", start = sval, maxit = 0)
fm_hp1 <- hurdle(trips ~ ., data = RecreationDemand, dist = "poisson",
zero = "poisson", start = sval)
fm_hp2 <- hurdle(trips ~ ., data = RecreationDemand, dist = "poisson",
zero = "poisson")
sapply(list(fm_hp0, fm_hp1, fm_hp2), logLik)
#> [1] -1209.582 -1181.612 -1181.612
## negbin-negbin
fm_hnb <- hurdle(trips ~ ., data = RecreationDemand, dist = "negbin", zero = "negbin")
summary(fm_hnb)
#>
#> Call:
#> hurdle(formula = trips ~ ., data = RecreationDemand, dist = "negbin",
#> zero.dist = "negbin")
#>
#> Pearson residuals:
#> Min 1Q Median 3Q Max
#> -1.664e+00 -2.496e-01 -1.540e-04 -5.297e-07 1.035e+01
#>
#> Count model coefficients (truncated negbin with log link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.84193 0.38278 2.200 0.02784 *
#> quality 0.17170 0.07234 2.374 0.01762 *
#> skiyes 0.62236 0.19013 3.273 0.00106 **
#> income -0.05709 0.06452 -0.885 0.37629
#> userfeeyes 0.57634 0.38508 1.497 0.13448
#> costC 0.05707 0.02169 2.632 0.00850 **
#> costS -0.07752 0.01155 -6.713 1.9e-11 ***
#> costH 0.01237 0.01490 0.830 0.40640
#> Log(theta) -0.53031 0.26114 -2.031 0.04228 *
#> Zero hurdle model coefficients (censored negbin with log link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -8.7012 11.3563 -0.766 0.443560
#> quality 65.6372 90.7017 0.724 0.469274
#> skiyes 14.1385 21.1554 0.668 0.503933
#> income -0.2729 1.5223 -0.179 0.857707
#> userfeeyes 344.0499 494.3586 0.696 0.486459
#> costC 0.7184 1.1478 0.626 0.531373
#> costS -1.4139 1.9876 -0.711 0.476873
#> costH 0.6650 1.0733 0.620 0.535567
#> Log(theta) -4.6672 1.3829 -3.375 0.000738 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Theta: count = 0.5884, zero = 0.0094
#> Number of iterations in BFGS optimization: 968
#> Log-likelihood: -718.3 on 18 Df
logLik(fm_hnb)
#> 'log Lik.' -718.2823 (df=18)
sval <- list(count = c(0.841, 0.172, .622, -.057, .576, .057, -.078, .012),
zero = c(-3.046, 4.638, -.025, .026, 16.203, 0.030, -.156, .117),
theta = c(count = 1/1.7, zero = 1/5.609))
fm_hnb2 <- try(hurdle(trips ~ ., data = RecreationDemand,
dist = "negbin", zero = "negbin", start = sval))
if(!inherits(fm_hnb2, "try-error")) {
summary(fm_hnb2)
logLik(fm_hnb2)
}
#> 'log Lik.' -718.3423 (df=18)
## geo-negbin
sval98 <- list(count = c(0.841, 0.172, .622, -.057, .576, .057, -.078, .012),
zero = c(-2.88, 1.44, .4, .03, 9.43, 0.01, -.08, .071),
theta = c(count = 1/1.7))
sval96 <- list(count = c(0.841, 0.172, .622, -.057, .576, .057, -.078, .012),
zero = c(-2.882, 1.437, .406, .026, 11.936, 0.008, -.081, .071),
theta = c(count = 1/1.7))
fm_hgnb <- hurdle(trips ~ ., data = RecreationDemand, dist = "negbin", zero = "geometric")
summary(fm_hgnb)
#>
#> Call:
#> hurdle(formula = trips ~ ., data = RecreationDemand, dist = "negbin",
#> zero.dist = "geometric")
#>
#> Pearson residuals:
#> Min 1Q Median 3Q Max
#> -1.4391 -0.2612 -0.1456 -0.1104 10.1255
#>
#> Count model coefficients (truncated negbin with log link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.84193 0.38278 2.200 0.02784 *
#> quality 0.17170 0.07234 2.374 0.01762 *
#> skiyes 0.62236 0.19013 3.273 0.00106 **
#> income -0.05709 0.06452 -0.885 0.37629
#> userfeeyes 0.57634 0.38508 1.497 0.13448
#> costC 0.05707 0.02169 2.632 0.00850 **
#> costS -0.07752 0.01155 -6.713 1.9e-11 ***
#> costH 0.01237 0.01490 0.830 0.40640
#> Log(theta) -0.53031 0.26114 -2.031 0.04228 *
#> Zero hurdle model coefficients (censored geometric with log link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -2.881e+00 4.356e-01 -6.615 3.73e-11 ***
#> quality 1.437e+00 1.058e-01 13.584 < 2e-16 ***
#> skiyes 4.055e-01 3.297e-01 1.230 0.21876
#> income 2.648e-02 8.400e-02 0.315 0.75259
#> userfeeyes 1.566e+01 1.497e+03 0.010 0.99165
#> costC 8.255e-03 2.707e-02 0.305 0.76038
#> costS -8.093e-02 1.713e-02 -4.726 2.29e-06 ***
#> costH 7.112e-02 2.168e-02 3.281 0.00104 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Theta: count = 0.5884
#> Number of iterations in BFGS optimization: 18
#> Log-likelihood: -741.3 on 17 Df
logLik(fm_hgnb)
#> 'log Lik.' -741.3096 (df=17)
## logLik with starting values from Gurmu + Trivedi 1996
fm_hgnb96 <- hurdle(trips ~ ., data = RecreationDemand, dist = "negbin", zero = "geometric",
start = sval96, maxit = 0)
logLik(fm_hgnb96)
#> 'log Lik.' -741.5969 (df=17)
## logit-negbin
fm_hgnb2 <- hurdle(trips ~ ., data = RecreationDemand, dist = "negbin")
summary(fm_hgnb2)
#>
#> Call:
#> hurdle(formula = trips ~ ., data = RecreationDemand, dist = "negbin")
#>
#> Pearson residuals:
#> Min 1Q Median 3Q Max
#> -1.4391 -0.2612 -0.1456 -0.1104 10.1255
#>
#> Count model coefficients (truncated negbin with log link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.84193 0.38278 2.200 0.02784 *
#> quality 0.17170 0.07234 2.374 0.01762 *
#> skiyes 0.62236 0.19013 3.273 0.00106 **
#> income -0.05709 0.06452 -0.885 0.37629
#> userfeeyes 0.57634 0.38508 1.497 0.13448
#> costC 0.05707 0.02169 2.632 0.00850 **
#> costS -0.07752 0.01155 -6.713 1.9e-11 ***
#> costH 0.01237 0.01490 0.830 0.40640
#> Log(theta) -0.53031 0.26114 -2.031 0.04228 *
#> Zero hurdle model coefficients (binomial with logit link):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -2.881e+00 4.356e-01 -6.615 3.73e-11 ***
#> quality 1.437e+00 1.058e-01 13.584 < 2e-16 ***
#> skiyes 4.055e-01 3.297e-01 1.230 0.21876
#> income 2.648e-02 8.400e-02 0.315 0.75259
#> userfeeyes 1.566e+01 1.497e+03 0.010 0.99165
#> costC 8.255e-03 2.707e-02 0.305 0.76038
#> costS -8.093e-02 1.713e-02 -4.726 2.29e-06 ***
#> costH 7.112e-02 2.168e-02 3.281 0.00104 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Theta: count = 0.5884
#> Number of iterations in BFGS optimization: 18
#> Log-likelihood: -741.3 on 17 Df
logLik(fm_hgnb2)
#> 'log Lik.' -741.3096 (df=17)
## Note: quasi-complete separation
with(RecreationDemand, table(trips > 0, userfee))
#> userfee
#> no yes
#> FALSE 417 0
#> TRUE 229 13
# }