## Load the oneinfl package
install.packages("oneinfl")
## Installing package into 'C:/Users/ryang/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)
## package 'oneinfl' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\ryang\AppData\Local\Temp\RtmpYjoew5\downloaded_packages
library(oneinfl)
## Section 5: Simulation studies
# For the "Poisson" part of Table 1. Reproduces the top row of the table. Increase
# "n" for additional rows. "nrep" can be reduced to save computational time, but
# won't reproduce results exactly. At "n <- 200", takes approx. 10 minutes to run.
set.seed(1)
n <- 200
b0 <- -2
b1 <- .4
b2 <- 0.2
g0 <- -21
g1 <- 2
g2 <- 0.5
b <- c(b0, b1, b2)
g <- c(g0, g1, g2)
formula <- y ~ x1 + x2 | x1 + x2
nrep <- 10000
bs <- matrix(, nrow=nrep, ncol=3)
gs <- matrix(, nrow=nrep, ncol=3)
ppbs <- matrix(, nrow=nrep, ncol=3)
for(i in 1:nrep) {
x1 <- rnorm(n,10,1)
x2 <- sample(c(0,1), n, replace=T)
X <- data.frame(rep(1,n), x1, x2)
Z <- X
y <- roipp(b, g, X, Z)
#barplot(table(y))
data <- data.frame(y, x1, x2)
mymodel <- oneinfl(formula, data, dist = "Poisson")
bs[i,] <- mymodel$b
gs[i,] <- mymodel$g
ppmodel <- truncreg(y ~ x1 + x2, data, dist = "Poisson")
ppbs[i,] <- ppmodel$b
}
100 * (colSums(bs) / nrep - b) / b
## [1] 0.16711987 0.04664234 -0.23491981
100 * (colSums(gs) / nrep - g) / g
## [1] 3.244641 3.275834 3.021585
100 * (colSums(ppbs) / nrep - b) / b
## [1] -258.07852 -141.20306 -75.61516
# For the first row of the "Negative Binomial" part of Table 1. n <- 200 must be
# increased to reproduce additional rows. Takes approx. 30 minutes.
set.seed(1)
n <- 200
b0 <- -2
b1 <- .4
b2 <- 0.2
g0 <- -21
g1 <- 2
g2 <- 0.5
alpha <- 10
b <- c(b0, b1, b2)
g <- c(g0, g1, g2)
formula <- y ~ x1 + x2 | x1 + x2
nrep <- 10000
bs <- matrix(, nrow=nrep, ncol=3)
gs <- matrix(, nrow=nrep, ncol=3)
as <- numeric(n)
ztnbbs <- matrix(, nrow=nrep, ncol=3)
ztnbas <- numeric(n)
for(i in 1:nrep) {
x1 <- rnorm(n,10,1)
x2 <- sample(c(0,1), n, replace=T)
X <- data.frame(rep(1,n), x1, x2)
Z <- X
y <- roiztnb(b, g, alpha, X, Z)
df <- data.frame(y, x1, x2)
mymodel <- oneinfl(formula, df, dist = "negbin")
bs[i,] <- mymodel$b
gs[i,] <- mymodel$g
as[i] <- mymodel$alpha
ztnbmodel <- truncreg(y ~ x1 + x2, df, dist = "negbin")
ztnbbs[i,] <- ztnbmodel$b
ztnbas[i] <- ztnbmodel$alpha
}
100 * (colSums(bs) / nrep - b) / b
## [1] -0.10047578 -0.09975944 -0.66880222
100 * (colSums(gs) / nrep - g) / g
## [1] 3.216478 3.247665 2.988251
100 * (sum(as) / nrep - alpha) / alpha
## [1] 27.76931
100 * (colSums(ztnbbs) / nrep - b) / b
## [1] -331.58534 -194.81544 -89.33777
100 * (sum(ztnbas) / nrep - alpha) / alpha
## [1] -96.02316
## Section 6.1: MedPar data
# Load the data
install.packages("msme")
## Installing package into 'C:/Users/ryang/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)
## package 'msme' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\ryang\AppData\Local\Temp\RtmpYjoew5\downloaded_packages
library(msme)
## Loading required package: MASS
## Loading required package: lattice
data(medpar)
# Table 2
df <- medpar
formula <- los ~ white + died + type2 + type3 | white + died + type2 + type3
OIZTNB <- oneinfl(formula, df, dist="negbin")
summary(OIZTNB)
## Call:
## formula: los ~ white + died + type2 + type3 | white + died + type2 + type3
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) 2.29913 0.07184 32.005 0.000e+00 ***
## bwhite -0.09708 0.07138 -1.360 1.738e-01
## bdied -0.06814 0.04492 -1.517 1.293e-01
## btype2 0.23413 0.05349 4.377 1.204e-05 ***
## btype3 0.75580 0.07884 9.586 0.000e+00 ***
##
## Coefficients (gamma):
## Estimate Std.Error z_value p.value
## g(Intercept) -4.2002 0.5060 -8.300 0.00000 ***
## gwhite 0.6594 0.4814 1.370 0.17075
## gdied 2.3349 0.2359 9.899 0.00000 ***
## gtype2 -0.5413 0.2829 -1.913 0.05569 .
## gtype3 -0.7507 0.4472 -1.679 0.09324 .
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 2.267 0.146 15.53 0 ***
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## average one-inflation: 0.0416855179433309
## average absolute one-inflation: 0.0680909723883016
## Log-likelihood: -4671.06423844655
margins(OIZTNB, df)
## Call:
## formula: los ~ white + died + type2 + type3 | white + died + type2 + type3
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## white -1.258 0.7344 -1.713 8.668e-02 .
## died -2.189 0.3964 -5.522 3.343e-08 ***
## type2 2.575 0.5875 4.382 1.174e-05 ***
## type3 10.142 1.4663 6.917 4.616e-12 ***
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
formula <- los ~ white + died + type2 + type3
ZTNB <- truncreg(formula, df, dist="negbin")
summary(ZTNB)
## Call:
## formula: los ~ white + died + type2 + type3
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) 2.3334 0.07499 31.115 0.000e+00 ***
## bwhite -0.1318 0.07469 -1.765 7.755e-02 .
## bdied -0.2512 0.04468 -5.622 1.889e-08 ***
## btype2 0.2601 0.05529 4.704 2.550e-06 ***
## btype3 0.7692 0.08259 9.314 0.000e+00 ***
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 1.881 0.1034 18.19 0 ***
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## Log-likelihood: -4736.77246562658
margins(ZTNB, df)
## Call:
## formula: los ~ white + died + type2 + type3
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## white -1.296 0.7762 -1.669 9.503e-02 .
## died -2.238 0.3868 -5.785 7.233e-09 ***
## type2 2.639 0.6163 4.283 1.846e-05 ***
## type3 10.159 1.5303 6.638 3.171e-11 ***
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
# Testing
oneWald(OIZTNB)
## $W
## [1] 469.1062
##
## $pval
## [1] 0
oneLRT(OIZTNB, ZTNB)
## $LRTstat
## [1] 131.4165
##
## $pval
## [1] 0
signifWald(OIZTNB, "white")
## $W
## [1] 3.725885
##
## $pval
## [1] 0.1552153
# Figure 1
oneplot(OIZTNB, ZTNB, df=df, maxpred=15)
## Section 6.2: Arizona Medicare data
library(COUNT)
## Loading required package: sandwich
data(azdrg112)
df <- azdrg112
# Table 3
OIZTNB.2 <- oneinfl(los ~ gender + type1 + age75 | gender + type1 + age75, df=df, dist="negbin")
summary(OIZTNB.2)
## Call:
## formula: los ~ gender + type1 + age75 | gender + type1 + age75
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) 0.7436 0.06622 11.228 0.000000 ***
## bgender -0.1778 0.04659 -3.817 0.000135 ***
## btype1 0.9226 0.05605 16.462 0.000000 ***
## bage75 0.1554 0.04975 3.123 0.001790 **
##
## Coefficients (gamma):
## Estimate Std.Error z_value p.value
## g(Intercept) -2.1308 0.1989 -10.7126 0.000e+00 ***
## ggender 0.5973 0.2054 2.9088 3.628e-03 **
## gtype1 -1.3396 0.1829 -7.3229 2.427e-13 ***
## gage75 -0.1225 0.2058 -0.5953 5.517e-01
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 2.143 0.2289 9.364 0 ***
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## average one-inflation: -0.16397156546163
## average absolute one-inflation: 0.16397156546163
## Log-likelihood: -4123.85968057303
margins(OIZTNB.2, df)
## Call:
## formula: los ~ gender + type1 + age75 | gender + type1 + age75
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## gender -0.7266 0.1654 -4.392 1.125e-05 ***
## type1 2.6973 0.1355 19.902 0.000e+00 ***
## age75 0.5520 0.1810 3.049 2.294e-03 **
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
ZTNB.2 <- truncreg(los ~ gender + type1 + age75, df, dist="negbin")
summary(ZTNB.2)
## Call:
## formula: los ~ gender + type1 + age75
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) 1.0505 0.04280 24.540 0.000e+00 ***
## bgender -0.1659 0.03534 -4.695 2.673e-06 ***
## btype1 0.7367 0.03977 18.524 0.000e+00 ***
## bage75 0.1307 0.03775 3.462 5.369e-04 ***
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 3.782 0.2804 13.49 0 ***
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## Log-likelihood: -4177.06389065033
margins(ZTNB.2, df)
## Call:
## formula: los ~ gender + type1 + age75
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## gender -0.7292 0.1592 -4.580 4.650e-06 ***
## type1 2.6991 0.1338 20.166 0.000e+00 ***
## age75 0.5795 0.1728 3.354 7.954e-04 ***
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
# Testing
oneWald(OIZTNB.2)
## $W
## [1] 690.5574
##
## $pval
## [1] 0
oneLRT(OIZTNB.2, ZTNB.2)
## $LRTstat
## [1] 106.4084
##
## $pval
## [1] 0
# Figure 2
oneplot(OIZTNB.2, ZTNB.2, df=df, maxpred=11, ylimit=390)
## Section 6.3: Recreation demand. Data is also from "COUNT" package.
data(loomis)
df <- loomis
df <- df[complete.cases(df), ]
# Table 4
formula <- anvisits ~ gender + income2 + income3 + income4 + travel2 + travel3 | gender + income2 + income3 + income4 + travel2 + travel3
OIZTNB.3 <- oneinfl(formula, df, dist="negbin")
summary(OIZTNB.3)
## Call:
## formula: anvisits ~ gender + income2 + income3 + income4 + travel2 + travel3 |
## formula: gender + income2 + income3 + income4 + travel2 + travel3
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) 3.4145 0.4297 7.9458 1.998e-15 ***
## bgender 0.4694 0.1963 2.3914 1.679e-02 *
## bincome2 0.1887 0.2848 0.6627 5.075e-01
## bincome3 -0.5263 0.3133 -1.6800 9.297e-02 .
## bincome4 -0.7804 0.2697 -2.8939 3.804e-03 **
## btravel2 -0.4987 0.2086 -2.3907 1.682e-02 *
## btravel3 -2.1424 0.2886 -7.4244 1.132e-13 ***
##
## Coefficients (gamma):
## Estimate Std.Error z_value p.value
## g(Intercept) -3.480031 0.9862 -3.52857 4.178e-04 ***
## ggender -0.530470 0.3427 -1.54771 1.217e-01
## gincome2 0.346603 0.6042 0.57363 5.662e-01
## gincome3 1.141761 0.5877 1.94261 5.206e-02 .
## gincome4 0.003421 0.5554 0.00616 9.951e-01
## gtravel2 1.982466 0.7668 2.58542 9.726e-03 **
## gtravel3 4.827657 0.7633 6.32452 2.540e-10 ***
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 0.4872 0.09618 5.065 4.081e-07 ***
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## average one-inflation: 0.242633905334975
## average absolute one-inflation: 0.264602639179767
## Log-likelihood: -1116.286418086
margins(OIZTNB.3, df)
## Call:
## formula: anvisits ~ gender + income2 + income3 + income4 + travel2 + travel3 |
## formula: gender + income2 + income3 + income4 + travel2 + travel3
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## gender 8.554 1.990 4.2989 1.717e-05 ***
## income2 4.154 7.756 0.5356 5.922e-01
## income3 -14.199 6.082 -2.3347 1.956e-02 *
## income4 -17.004 5.664 -3.0024 2.679e-03 **
## travel2 -16.782 5.991 -2.8011 5.093e-03 **
## travel3 -40.751 4.588 -8.8823 0.000e+00 ***
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
ZTNB.3 <- truncreg(anvisits ~ gender + income2 + income3 + income4 + travel2 + travel3, df, dist="negbin")
summary(ZTNB.3)
## Call:
## formula: anvisits ~ gender + income2 + income3 + income4 + travel2 + travel3
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) -0.3259 1.0241 -0.3182 0.750299
## bgender 0.8408 0.2810 2.9919 0.002773 **
## bincome2 1.1661 0.4840 2.4093 0.015983 *
## bincome3 -0.9575 0.4195 -2.2827 0.022450 *
## bincome4 -0.3259 0.4297 -0.7584 0.448204
## btravel2 -0.5822 0.3578 -1.6272 0.103700
## btravel3 -4.1408 0.3557 -11.6410 0.000000 ***
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 0.007729 0.006475 1.194 0.2326
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## Log-likelihood: -1165.82439770087
margins(ZTNB.3, df)
## Call:
## formula: anvisits ~ gender + income2 + income3 + income4 + travel2 + travel3
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## gender 14.01 3.55 3.9463 7.936e-05 ***
## income2 41.26 23.85 1.7301 8.361e-02 .
## income3 -24.20 9.56 -2.5315 1.136e-02 *
## income4 -10.28 12.19 -0.8433 3.991e-01
## travel2 -20.31 13.97 -1.4538 1.460e-01
## travel3 -54.93 12.83 -4.2803 1.867e-05 ***
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
# Testing
oneWald(OIZTNB.3)
## $W
## [1] 99.79372
##
## $pval
## [1] 0
oneLRT(OIZTNB.3, ZTNB.3)
## $LRTstat
## [1] 99.07596
##
## $pval
## [1] 0
# Figure 3
oneplot(OIZTNB.3, ZTNB.3, df=df, maxpred=9, ylim=120)
# Comparison to glm.nb()
library(MASS)
glmlmnb <- glm.nb(anvisits ~ gender + factor(income) + factor(travel), data=df)
summary(glmlmnb)
##
## Call:
## glm.nb(formula = anvisits ~ gender + factor(income) + factor(travel),
## data = df, init.theta = 0.8128559417, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5065 0.2868 12.225 < 2e-16 ***
## gender 0.4863 0.1326 3.668 0.000245 ***
## factor(income)2 0.2511 0.2022 1.242 0.214368
## factor(income)3 -0.6834 0.2144 -3.187 0.001439 **
## factor(income)4 -0.5818 0.1977 -2.942 0.003261 **
## factor(travel)2 -0.5700 0.1566 -3.639 0.000273 ***
## factor(travel)3 -2.6389 0.1716 -15.381 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.8129) family taken to be 1)
##
## Null deviance: 689.22 on 341 degrees of freedom
## Residual deviance: 363.59 on 335 degrees of freedom
## AIC: 2582.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.8129
## Std. Err.: 0.0593
##
## 2 x log-likelihood: -2566.6200
## Section 6.4: Fishing trips
install.packages("foreign")
## Installing package into 'C:/Users/ryang/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)
## package 'foreign' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\ryang\AppData\Local\Temp\RtmpYjoew5\downloaded_packages
library(foreign)
df <- read.dta("https://rtgodwin.com/data/subject.dta")
# Table 5
formula <- NTRIPS ~ AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1 | AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1
OIZTNB.4 <- oneinfl(formula, df, dist = "negbin")
summary(OIZTNB.4)
## Call:
## formula: NTRIPS ~ AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1 |
## formula: AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) 2.188192 0.103373 21.1678 0.000e+00 ***
## bAVLONG -0.283995 0.026906 -10.5553 0.000e+00 ***
## bBAD -0.147301 0.077993 -1.8887 5.894e-02 .
## bCROWD1 -0.025764 0.033047 -0.7796 4.356e-01
## bINC 0.001001 0.001032 0.9705 3.318e-01
## bMPERMILE 1.012869 0.443413 2.2843 2.236e-02 *
## bFOFT1 -0.190088 0.034755 -5.4694 4.516e-08 ***
## bWORK1 0.118410 0.033658 3.5180 4.348e-04 ***
##
## Coefficients (gamma):
## Estimate Std.Error z_value p.value
## g(Intercept) -2.32275 0.283616 -8.1898 2.220e-16 ***
## gAVLONG 0.25548 0.046180 5.5323 3.161e-08 ***
## gBAD 0.19118 0.206209 0.9271 3.539e-01
## gCROWD1 0.07463 0.090504 0.8245 4.096e-01
## gINC 0.00250 0.002629 0.9512 3.415e-01
## gMPERMILE -3.02748 1.461473 -2.0715 3.831e-02 *
## gFOFT1 0.41578 0.086906 4.7842 1.716e-06 ***
## gWORK1 0.02545 0.091523 0.2781 7.810e-01
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 1.397 0.1632 8.561 0 ***
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## average one-inflation: -0.0221051956832215
## average absolute one-inflation: 0.047671414492766
## Log-likelihood: -2851.82877040832
margins(OIZTNB.4, df)
## Call:
## formula: NTRIPS ~ AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1 |
## formula: AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## AVLONG -1.674383 0.166314 -10.0676 0.000e+00 ***
## BAD -0.795008 0.365812 -2.1733 2.976e-02 *
## CROWD1 -0.175464 0.185951 -0.9436 3.454e-01
## INC 0.004344 0.005786 0.7508 4.528e-01
## MPERMILE 6.941088 2.531932 2.7414 6.117e-03 **
## FOFT1 -1.232849 0.199962 -6.1654 7.030e-10 ***
## WORK1 0.637671 0.192160 3.3184 9.052e-04 ***
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
ZTNB.4 <- truncreg(NTRIPS ~ AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1, df, dist="negbin")
summary(ZTNB.4)
## Call:
## formula: NTRIPS ~ AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1
## distribution: negbin
##
## Coefficients (beta):
## Estimate Std.Error z_value p.value
## b(Intercept) 2.1523785 0.0905367 23.7735 0.000e+00 ***
## bAVLONG -0.2568544 0.0217481 -11.8104 0.000e+00 ***
## bBAD -0.1507364 0.0698623 -2.1576 3.096e-02 *
## bCROWD1 -0.0341243 0.0298342 -1.1438 2.527e-01
## bINC 0.0004125 0.0009202 0.4482 6.540e-01
## bMPERMILE 1.1877897 0.4073606 2.9158 3.548e-03 **
## bFOFT1 -0.2140346 0.0306926 -6.9735 3.092e-12 ***
## bWORK1 0.0967090 0.0301915 3.2032 1.359e-03 **
##
## alpha:
## Estimate Std.Error z_value p.value
## 1 1.451 0.1189 12.2 0 ***
##
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
##
## Log-likelihood: -2858.71598109464
margins(ZTNB.4, df)
## Call:
## formula: NTRIPS ~ AVLONG + BAD + CROWD1 + INC + MPERMILE + FOFT1 + WORK1
## distribution: negbin
##
## Marginal effects:
## Marginal.effects Std.Error z_value p.value
## AVLONG -1.445694 0.138637 -10.4279 0.000e+00 ***
## BAD -0.814024 0.363445 -2.2397 2.511e-02 *
## CROWD1 -0.192067 0.168140 -1.1423 2.533e-01
## INC 0.002322 0.005181 0.4481 6.541e-01
## MPERMILE 6.685423 2.316452 2.8861 3.901e-03 **
## FOFT1 -1.204684 0.181925 -6.6219 3.547e-11 ***
## WORK1 0.544322 0.172017 3.1644 1.554e-03 **
##
## Signif. codes: 0 ***' 0.001 **' 0.01 *' 0.05 .' 0.1 ' 1
# Testing
oneWald(OIZTNB.4)
## $W
## [1] 395.965
##
## $pval
## [1] 0
oneLRT(OIZTNB.4, ZTNB.4)
## $LRTstat
## [1] 13.77442
##
## $pval
## [1] 0.08783784
# Figure 4
oneplot(OIZTNB.4, ZTNB.4, df=df, maxpred=14, ylimit=170)
sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_Canada.utf8 LC_CTYPE=English_Canada.utf8 LC_MONETARY=English_Canada.utf8 LC_NUMERIC=C LC_TIME=English_Canada.utf8
##
## time zone: America/Winnipeg
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] foreign_0.8-88 COUNT_1.3.4 sandwich_3.1-0 msme_0.5.3 lattice_0.22-6 MASS_7.3-61 oneinfl_1.0.1
##
## loaded via a namespace (and not attached):
## [1] zoo_1.8-12 compiler_4.4.2 tools_4.4.2 grid_4.4.2 knitr_1.49 xfun_0.50 evaluate_0.23