
R> message("Installing replication dependencies...")

R> pks = c("bizicount", "dplyr", "tidyr", "ggplot2", 
+      "copula", "doParallel", "doRNG", "RhpcBLASctl")

R> needed = setdiff(pks, rownames(installed.packages()))

R> options(Ncpus = parallel::detectCores() - 1)

R> if (length(needed) > 0) {
+      message("Installing needed packages: ", paste(needed, collapse = ", "))
+      install.packages(needed, repos = "h ..." ... [TRUNCATED] 

R> message("Attempted installation of dependencies completed.")

R> suppressMessages(invisible(sapply(pks, library, character.only = TRUE)))

R> sink("sessioninfo.txt", split = TRUE)

R> print(sessionInfo())
R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux trixie/sid

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Vienna
tzcode source: system (glibc)

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] RhpcBLASctl_0.23-42 doRNG_1.8.6         rngtools_1.5.2     
 [4] doParallel_1.0.17   iterators_1.0.14    foreach_1.5.2      
 [7] copula_1.1-3        ggplot2_3.5.0       tidyr_1.3.1        
[10] dplyr_1.1.4         bizicount_1.3.2    

loaded via a namespace (and not attached):
 [1] utf8_1.2.4          generics_0.1.3      lattice_0.22-5     
 [4] digest_0.6.35       lme4_1.1-35.1       magrittr_2.0.3     
 [7] grid_4.3.3          mvtnorm_1.2-4       Matrix_1.6-5       
[10] Formula_1.2-5       httr_1.4.7          purrr_1.0.2        
[13] fansi_1.0.6         scales_1.3.0        stabledist_0.7-1   
[16] pbivnorm_0.6.0      codetools_0.2-19    numDeriv_2016.8-1.1
[19] cli_3.6.2           rlang_1.1.3         pspline_1.0-19     
[22] texreg_1.39.3       gsl_2.1-8           munsell_0.5.0      
[25] splines_4.3.3       withr_3.0.0         tools_4.3.3        
[28] nloptr_2.0.3        minqa_1.2.6         colorspace_2.1-0   
[31] boot_1.3-30         vctrs_0.6.5         R6_2.5.1           
[34] stats4_4.3.3        lifecycle_1.0.4     ADGofTest_0.3      
[37] MASS_7.3-60.0.1     pcaPP_2.0-4         pkgconfig_2.0.3    
[40] pillar_1.9.0        gtable_0.3.4        glue_1.7.0         
[43] Rcpp_1.0.12         tibble_3.2.1        tidyselect_1.2.1   
[46] nlme_3.1-164        DHARMa_0.4.6        compiler_4.3.3     

R> sink()

R> message("Beginning empirical replication...")

R> rm(list = ls())

R> library(bizicount)

R> suppressMessages(library(texreg))

R> suppressMessages(library(DHARMa))

R> dat = terror

R> dat[, c("xcoord", "ycoord", "pop", "mtns")] = scale(dat[, 
+      c("xcoord", "ycoord", "pop", "mtns")])

R> fmla.ful = att.ful ~ pop + mtns + xcoord * ycoord

R> fmla.bok = att.bok ~ pop + mtns + xcoord * ycoord

R> bivpois = bizicount(fmla.ful, fmla.bok, data = dat, 
+      cop = "frank", margins = c("pois", "pois"), keep = TRUE)

R> summary(bivpois)
Call:
bizicount(fmla1 = fmla.ful, fmla2 = fmla.bok, data = dat, cop = "frank", 
    margins = c("pois", "pois"), keep = TRUE)

================================================= 
Count Model: att.ful | 
---------------------- 

              Estimate Std. Err.  Z value  Pr(>|z|)    
(Intercept)   -1.26185   0.12291 -10.2661 < 2.2e-16 ***
pop           -0.44110   0.20166  -2.1873   0.02872 *  
mtns          -0.44051   0.10896  -4.0429 5.278e-05 ***
xcoord         0.72976   0.14268   5.1147 3.143e-07 ***
ycoord        -0.69300   0.12264  -5.6505 1.600e-08 ***
xcoord:ycoord -1.03695   0.17133  -6.0523 1.428e-09 ***

------------------------------------------------- 
           Estimate Std. Err. Z value Pr(>|z|)  
dependence  1.64169   0.75955  2.1614  0.03066 *
------------------------------------------------- 

Count Model: att.bok | 
---------------------- 

               Estimate Std. Err.  Z value  Pr(>|z|)    
(Intercept)   -2.378585  0.197090 -12.0685 < 2.2e-16 ***
pop            0.944876  0.062597  15.0945 < 2.2e-16 ***
mtns          -0.050484  0.056749  -0.8896    0.3737    
xcoord         2.634690  0.161948  16.2688 < 2.2e-16 ***
ycoord         0.952284  0.153673   6.1968 5.762e-10 ***
xcoord:ycoord -0.585055  0.119683  -4.8884 1.017e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

================================================= 

R> zi_test(bivpois)

====================================================
He et al. (2019)'s Test for Zero Modification
--------------------------------------

H_0:  Pr(y = 0 | x) = dpois(0 | x) 


             H_a   Z_score      p_value   n
att.ful inflated 10.809096 1.558638e-27 312
att.bok inflated  8.404112 2.155557e-17 312

====================================================

R> fmla.zi.ful = att.ful ~ pop + mtns + xcoord * ycoord | 
+      pop + mtns + xcoord * ycoord

R> fmla.zi.bok = att.bok ~ pop + mtns + xcoord * ycoord | 
+      pop + mtns + xcoord * ycoord

R> zi_bivpois = bizicount(fmla.zi.ful, fmla.zi.bok, data = dat, 
+      cop = "frank", margins = c("zip", "zip"), keep = T)

R> summary(zi_bivpois)
Call:
bizicount(fmla1 = fmla.zi.ful, fmla2 = fmla.zi.bok, data = dat, 
    cop = "frank", margins = c("zip", "zip"), keep = T)

================================================= 
Count Model: att.ful | 
---------------------- 

              Estimate Std. Err. Z value  Pr(>|z|)    
(Intercept)    0.28469   0.12039  2.3647   0.01805 *  
pop           -0.76644   0.36025 -2.1275   0.03338 *  
mtns           0.25733   0.10566  2.4354   0.01487 *  
xcoord        -2.04704   0.26001 -7.8729 3.465e-15 ***
ycoord        -1.00707   0.14800 -6.8044 1.015e-11 ***
xcoord:ycoord -1.57148   0.30833 -5.0967 3.456e-07 ***

++++++++++++++++++++++++ 
Zero Inflation: att.ful | 
------------------------- 

               Estimate Std. Err. Z value  Pr(>|z|)    
(Intercept)    -3.46119   1.07864 -3.2089 0.0013326 ** 
pop            -0.43413   0.62324 -0.6966 0.4860746    
mtns            1.73496   0.48970  3.5429 0.0003957 ***
xcoord        -10.54553   2.38677 -4.4183 9.947e-06 ***
ycoord         -4.74161   1.22946 -3.8567 0.0001150 ***
xcoord:ycoord  -6.89191   1.97450 -3.4905 0.0004822 ***

------------------------------------------------- 
           Estimate Std. Err. Z value Pr(>|z|)
dependence -0.76749   1.31529 -0.5835   0.5595
------------------------------------------------- 

Count Model: att.bok | 
---------------------- 

               Estimate Std. Err. Z value Pr(>|z|)    
(Intercept)    0.567623  0.193178  2.9383 0.003300 ** 
pop            0.232673  0.076434  3.0441 0.002334 ** 
mtns          -0.019449  0.052042 -0.3737 0.708611    
xcoord         1.217380  0.148355  8.2059 2.29e-16 ***
ycoord         0.162435  0.210258  0.7725 0.439790    
xcoord:ycoord -0.402957  0.127004 -3.1728 0.001510 ** 

++++++++++++++++++++++++ 
Zero Inflation: att.bok | 
------------------------- 

              Estimate Std. Err. Z value Pr(>|z|)    
(Intercept)    1.88502   0.27614  6.8264 8.71e-12 ***
pop           -0.75792   0.27215 -2.7849 0.005354 ** 
mtns          -0.15286   0.20591 -0.7424 0.457857    
xcoord        -0.77939   0.34073 -2.2874 0.022173 *  
ycoord        -0.52526   0.27307 -1.9236 0.054407 .  
xcoord:ycoord -0.90781   0.31491 -2.8827 0.003942 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

================================================= 

R> zip_dharmas = make_DHARMa(zi_bivpois, nsim = 5000, 
+      seed = 789443)

R> class(zip_dharmas)
[1] "list"

R> lapply(zip_dharmas, class)
[[1]]
[1] "DHARMa"

[[2]]
[1] "DHARMa"


R> lapply(zip_dharmas, testResiduals, plot = FALSE)
$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.041527, p-value = 0.6549
alternative hypothesis: two-sided


$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 4.0577, p-value = 4e-04
alternative hypothesis: two.sided


$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 2, observations = 312, p-value < 2.2e-16
alternative hypothesis: two.sided
 percent confidence interval:
 0.000000000 0.003205128
sample estimates:
outlier frequency (expected: 0.00016025641025641 ) 
                                       0.006410256 


$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.067001, p-value = 0.1214
alternative hypothesis: two-sided


$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 4.0115, p-value < 2.2e-16
alternative hypothesis: two.sided


$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 4, observations = 312, p-value < 2.2e-16
alternative hypothesis: two.sided
 percent confidence interval:
 0 0
sample estimates:
outlier frequency (expected: 0 ) 
                      0.01282051 


[[1]]
[[1]]$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.041527, p-value = 0.6549
alternative hypothesis: two-sided


[[1]]$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 4.0577, p-value = 4e-04
alternative hypothesis: two.sided


[[1]]$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 2, observations = 312, p-value < 2.2e-16
alternative hypothesis: two.sided
 percent confidence interval:
 0.000000000 0.003205128
sample estimates:
outlier frequency (expected: 0.00016025641025641 ) 
                                       0.006410256 



[[2]]
[[2]]$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.067001, p-value = 0.1214
alternative hypothesis: two-sided


[[2]]$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 4.0115, p-value < 2.2e-16
alternative hypothesis: two.sided


[[2]]$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 4, observations = 312, p-value < 2.2e-16
alternative hypothesis: two.sided
 percent confidence interval:
 0 0
sample estimates:
outlier frequency (expected: 0 ) 
                      0.01282051 




R> zi_bivnb = bizicount(fmla.zi.ful, fmla.zi.bok, data = dat, 
+      cop = "frank", margins = c("zinb", "zinb"), keep = T)

R> chi_stat_nb = -2 * (logLik(zi_bivpois) - logLik(zi_bivnb))

R> pchisq(chi_stat_nb, df = 2, lower.tail = FALSE)
[1] 1.202813e-53

R> zinb_dharmas = make_DHARMa(zi_bivnb, nsim = 5000, 
+      seed = 12473)

R> lapply(zinb_dharmas, testResiduals, plot = FALSE)
$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.021347, p-value = 0.9989
alternative hypothesis: two-sided


$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 1.4655, p-value = 0.2976
alternative hypothesis: two.sided


$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 1, observations = 312, p-value = 0.14
alternative hypothesis: two.sided
 percent confidence interval:
 0.000000000 0.003205128
sample estimates:
outlier frequency (expected: 0.000224358974358974 ) 
                                        0.003205128 


$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.039582, p-value = 0.7126
alternative hypothesis: two-sided


$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 0.53827, p-value = 0.5048
alternative hypothesis: two.sided


$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 0, observations = 312, p-value = 1
alternative hypothesis: two.sided
 percent confidence interval:
 0.000000000 0.003205128
sample estimates:
outlier frequency (expected: 0.000192307692307692 ) 
                                                  0 


[[1]]
[[1]]$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.021347, p-value = 0.9989
alternative hypothesis: two-sided


[[1]]$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 1.4655, p-value = 0.2976
alternative hypothesis: two.sided


[[1]]$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 1, observations = 312, p-value = 0.14
alternative hypothesis: two.sided
 percent confidence interval:
 0.000000000 0.003205128
sample estimates:
outlier frequency (expected: 0.000224358974358974 ) 
                                        0.003205128 



[[2]]
[[2]]$uniformity

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.039582, p-value = 0.7126
alternative hypothesis: two-sided


[[2]]$dispersion

	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
	simulated

data:  simulationOutput
dispersion = 0.53827, p-value = 0.5048
alternative hypothesis: two.sided


[[2]]$outliers

	DHARMa bootstrapped outlier test

data:  simulationOutput
outliers at both margin(s) = 0, observations = 312, p-value = 1
alternative hypothesis: two.sided
 percent confidence interval:
 0.000000000 0.003205128
sample estimates:
outlier frequency (expected: 0.000192307692307692 ) 
                                                  0 




R> mods = list(bivpois, zi_bivpois, zi_bivnb)

R> texreg(mods, groups = list(`Count Model` = 1:6, `Zero Inflation` = 7:12), 
+      reorder.coef = c(1:6, 8:13, 14, 7), digits = 2, custom.model.name .... [TRUNCATED] 

\begin{table}
\begin{center}
\begin{tabular}{l c c c c c c}
\hline
 & Poisson: Fulani & Poisson: BH & ZIP: Fulani & ZIP: BH & ZINB: Fulani & ZINB: BH \\
\hline
Count Model             &               &               &                &              &                &               \\
                        &               &               &                &              &                &               \\
\quad ct\_(Intercept)   & $-1.26^{***}$ & $-2.38^{***}$ & $0.28^{*}$     & $0.57^{**}$  & $0.20$         & $-0.17$       \\
                        & $(0.12)$      & $(0.20)$      & $(0.12)$       & $(0.19)$     & $(0.24)$       & $(0.31)$      \\
\quad ct\_pop           & $-0.44^{*}$   & $0.94^{***}$  & $-0.77^{*}$    & $0.23^{**}$  & $-0.72$        & $0.51^{**}$   \\
                        & $(0.20)$      & $(0.06)$      & $(0.36)$       & $(0.08)$     & $(0.77)$       & $(0.16)$      \\
\quad ct\_mtns          & $-0.44^{***}$ & $-0.05$       & $0.26^{*}$     & $-0.02$      & $0.17$         & $-0.17$       \\
                        & $(0.11)$      & $(0.06)$      & $(0.11)$       & $(0.05)$     & $(0.17)$       & $(0.19)$      \\
\quad ct\_xcoord        & $0.73^{***}$  & $2.63^{***}$  & $-2.05^{***}$  & $1.22^{***}$ & $-2.13^{***}$  & $2.24^{***}$  \\
                        & $(0.14)$      & $(0.16)$      & $(0.26)$       & $(0.15)$     & $(0.48)$       & $(0.29)$      \\
\quad ct\_ycoord        & $-0.69^{***}$ & $0.95^{***}$  & $-1.01^{***}$  & $0.16$       & $-1.08^{**}$   & $-1.26^{***}$ \\
                        & $(0.12)$      & $(0.15)$      & $(0.15)$       & $(0.21)$     & $(0.33)$       & $(0.37)$      \\
\quad ct\_xcoord:ycoord & $-1.04^{***}$ & $-0.59^{***}$ & $-1.57^{***}$  & $-0.40^{**}$ & $-1.12$        & $-0.22$       \\
                        & $(0.17)$      & $(0.12)$      & $(0.31)$       & $(0.13)$     & $(0.60)$       & $(0.26)$      \\
Zero Inflation          &               &               &                &              &                &               \\
                        &               &               &                &              &                &               \\
\quad zi\_(Intercept)   &               &               & $-3.46^{**}$   & $1.89^{***}$ & $-6.04^{**}$   & $-0.67$       \\
                        &               &               & $(1.08)$       & $(0.28)$     & $(2.19)$       & $(1.02)$      \\
\quad zi\_pop           &               &               & $-0.43$        & $-0.76^{**}$ & $-0.52$        & $-1.81^{*}$   \\
                        &               &               & $(0.62)$       & $(0.27)$     & $(0.91)$       & $(0.72)$      \\
\quad zi\_mtns          &               &               & $1.73^{***}$   & $-0.15$      & $2.20^{**}$    & $0.25$        \\
                        &               &               & $(0.49)$       & $(0.21)$     & $(0.82)$       & $(0.51)$      \\
\quad zi\_xcoord        &               &               & $-10.55^{***}$ & $-0.78^{*}$  & $-14.38^{***}$ & $0.47$        \\
                        &               &               & $(2.39)$       & $(0.34)$     & $(4.05)$       & $(0.71)$      \\
\quad zi\_ycoord        &               &               & $-4.74^{***}$  & $-0.53$      & $-6.58^{**}$   & $-5.33^{**}$  \\
                        &               &               & $(1.23)$       & $(0.27)$     & $(2.12)$       & $(1.72)$      \\
\quad zi\_xcoord:ycoord &               &               & $-6.89^{***}$  & $-0.91^{**}$ & $-8.79^{**}$   & $-1.85$       \\
                        &               &               & $(1.97)$       & $(0.31)$     & $(2.96)$       & $(1.17)$      \\
disp                    &               &               &                &              & $0.64^{**}$    & $0.68^{***}$  \\
                        &               &               &                &              & $(0.21)$       & $(0.16)$      \\
dependence              & $1.64^{*}$    & $1.64^{*}$    & $-0.77$        & $-0.77$      & $0.04$         & $0.04$        \\
                        & $(0.76)$      & $(0.76)$      & $(1.32)$       & $(1.32)$     & $(1.35)$       & $(1.35)$      \\
\hline
N                       & $312$         & $312$         & $312$          & $312$        & $312$          & $312$         \\
BIC                     & $1442.68$     & $1442.68$     & $1097.95$      & $1097.95$    & $855.63$       & $855.63$      \\
AIC                     & $1459.67$     & $1459.67$     & $1130.62$      & $1130.62$    & $890.91$       & $890.91$      \\
LogLik                  & $-716.84$     & $-716.84$     & $-540.31$      & $-540.31$    & $-418.46$      & $-418.46$     \\
\hline
\multicolumn{7}{l}{\scriptsize{$^{***}p<0.001$; $^{**}p<0.01$; $^{*}p<0.05$}}
\end{tabular}
\caption{Statistical models}
\label{table:coefficients}
\end{center}
\end{table}

R> screenreg(mods, groups = list(`Count Model` = 1:6, 
+      `Zero Inflation` = 7:12), reorder.coef = c(1:6, 8:13, 14, 
+      7), digits = 2)

========================================================================================================
                      Biv: att.ful  Biv: att.bok  Biv: att.ful  Biv: att.bok  Biv: att.ful  Biv: att.bok
--------------------------------------------------------------------------------------------------------
Count Model                                                                                             
                                                                                                        
    ct_(Intercept)      -1.26 ***     -2.38 ***      0.28 *        0.57 **       0.20         -0.17     
                        (0.12)        (0.20)        (0.12)        (0.19)        (0.24)        (0.31)    
    ct_pop              -0.44 *        0.94 ***     -0.77 *        0.23 **      -0.72          0.51 **  
                        (0.20)        (0.06)        (0.36)        (0.08)        (0.77)        (0.16)    
    ct_mtns             -0.44 ***     -0.05          0.26 *       -0.02          0.17         -0.17     
                        (0.11)        (0.06)        (0.11)        (0.05)        (0.17)        (0.19)    
    ct_xcoord            0.73 ***      2.63 ***     -2.05 ***      1.22 ***     -2.13 ***      2.24 *** 
                        (0.14)        (0.16)        (0.26)        (0.15)        (0.48)        (0.29)    
    ct_ycoord           -0.69 ***      0.95 ***     -1.01 ***      0.16         -1.08 **      -1.26 *** 
                        (0.12)        (0.15)        (0.15)        (0.21)        (0.33)        (0.37)    
    ct_xcoord:ycoord    -1.04 ***     -0.59 ***     -1.57 ***     -0.40 **      -1.12         -0.22     
                        (0.17)        (0.12)        (0.31)        (0.13)        (0.60)        (0.26)    
Zero Inflation                                                                                          
                                                                                                        
    zi_(Intercept)                                  -3.46 **       1.89 ***     -6.04 **      -0.67     
                                                    (1.08)        (0.28)        (2.19)        (1.02)    
    zi_pop                                          -0.43         -0.76 **      -0.52         -1.81 *   
                                                    (0.62)        (0.27)        (0.91)        (0.72)    
    zi_mtns                                          1.73 ***     -0.15          2.20 **       0.25     
                                                    (0.49)        (0.21)        (0.82)        (0.51)    
    zi_xcoord                                      -10.55 ***     -0.78 *      -14.38 ***      0.47     
                                                    (2.39)        (0.34)        (4.05)        (0.71)    
    zi_ycoord                                       -4.74 ***     -0.53         -6.58 **      -5.33 **  
                                                    (1.23)        (0.27)        (2.12)        (1.72)    
    zi_xcoord:ycoord                                -6.89 ***     -0.91 **      -8.79 **      -1.85     
                                                    (1.97)        (0.31)        (2.96)        (1.17)    
disp                                                                             0.64 **       0.68 *** 
                                                                                (0.21)        (0.16)    
dependence               1.64 *        1.64 *       -0.77         -0.77          0.04          0.04     
                        (0.76)        (0.76)        (1.32)        (1.32)        (1.35)        (1.35)    
--------------------------------------------------------------------------------------------------------
N                      312           312           312           312           312           312        
BIC                   1442.68       1442.68       1097.95       1097.95        855.63        855.63     
AIC                   1459.67       1459.67       1130.62       1130.62        890.91        890.91     
LogLik                -716.84       -716.84       -540.31       -540.31       -418.46       -418.46     
========================================================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

R> screenreg(mods, groups = list(`Count Model` = 1:6, 
+      `Zero Inflation` = 7:12), reorder.coef = c(1:6, 8:13, 14, 
+      7), digits = 2, file = .... [TRUNCATED] 

R> name = c("fulani", "bh")

R> for (i in seq_along(name)) {
+      pdf(paste0("Figures/appendix_figure_B", i, ".pdf"), width = 8, 
+          height = 4.5)
+      plot(zip_dharma .... [TRUNCATED] 

R> message("Empirical replication complete.")

R> message("Beginning copula plots (Figure 1)...")

R> library(copula)

R> cops = list(frankCopula(4), frankCopula(-2), normalCopula(-0.8), 
+      normalCopula(0.5))

R> ns = rep(c("fr", "norm"), each = 2)

R> dep = c(4, 2, 8, 5)

R> row = 1:4

R> col = 1:4

R> for (i in seq_along(cops)) {
+      iter = 1
+      a = wireframe2(cops[[i]], dCopula, delta = 0.01, col = "gray50", 
+          shade = TRUE, draw .... [TRUNCATED] 

R> message("Figure 1 Copula plots complete.")

R> message("Beginning simulations...")

R> rm(list = ls())

R> pks = c("bizicount", "doRNG", "parallel", "doParallel", 
+      "RhpcBLASctl")

R> invisible(sapply(pks, library, character.only = TRUE))

R> gen = function(n, b1, b2, psi1, psi2, dep) {
+      k1 = length(b1)
+      k2 = length(b2)
+      X1 = cbind(1, matrix(rbinom(n * (k1 - 1), 1, 0.5) .... [TRUNCATED] 

R> estimate = function(...) {
+      gen_args = list(...)
+      data = gen(...)
+      if (min(data$y1) > 0 || min(data$y2) > 0) 
+          return(N .... [TRUNCATED] 

R> n = 500

R> dep = c(0.15, 0.85)

R> b1 = c(1, 3.25, -2.3)

R> b2 = c(2, -1.75, 3.5)

R> psi1 = psi2 = c(0.1, 0.6)

R> grid = expand.grid(n = n, psi1 = psi1, psi2 = psi2, 
+      dep = dep)

R> nsims = 500

R> grid = do.call(rbind.data.frame, replicate(nsims, 
+      grid, simplify = F))

R> set.seed(987465)

R> ncores = parallel::detectCores() - 1

R> cl = makeCluster(ncores, type = if (grepl("unix", 
+      .Platform$OS.type)) "FORK" else "PSOCK")

R> registerDoParallel(cl)

R> mcoptions <- list(preschedule = T)

R> RhpcBLASctl::blas_set_num_threads(1)

R> RhpcBLASctl::omp_set_num_threads(1)

R> start = Sys.time()

R> res = foreach(i = seq_len(nrow(grid)), .packages = c("bizicount", 
+      "RhpcBLASctl"), .export = ls(), .errorhandling = "pass", 
+      .verbose .... [TRUNCATED] 

R> end = Sys.time()

R> cat("Total time:\n")
Total time:

R> end - start
Time difference of 1.352509 hours

R> stopCluster(cl)

R> stopImplicitCluster()

R> session = sessionInfo()

R> save(res, session, file = "output_montes_small.RData")

R> message("\n\n---Monte carlo simulations complete. Data available at ./output_montes_small.RData---\n\n")

R> message("Beginning simulation plots...")

R> library(dplyr)

R> library(ggplot2)

R> library(tidyr)

R> load("output_montes_small.RData")

R> res = res[!sapply(res, is.null)]

R> out = as.data.frame(do.call(rbind, res))

R> names(out) = tolower(gsub("[(|)]", "", names(out)))

R> names(out) = gsub("uni1.ct_", "uni.ct1_", names(out))

R> names(out) = gsub("uni1.zi_", "uni.zi1_", names(out))

R> names(out) = gsub("uni2.ct_", "uni.ct2_", names(out))

R> names(out) = gsub("uni2.zi_", "uni.zi2_", names(out))

R> names(out) = gsub("(.*)(\\d)[.](\\d)", "\\1\\2\\3", 
+      names(out))

R> ests = out %>% pivot_longer(cols = starts_with("est"), 
+      names_to = c("model", "param"), names_pattern = "est[.](.*)[.](.*)", 
+      values_ .... [TRUNCATED] 

R> ses = out %>% pivot_longer(cols = starts_with("se"), 
+      names_to = c("model", "param"), names_pattern = "se[.](.*)[.](.*)", 
+      values_to  .... [TRUNCATED] 

R> clean = left_join(ests, ses, by = "id")

R> clean = clean %>% dplyr::select(-ends_with("x"))

R> names(clean) = sub(".y", "", names(clean))

R> clean = clean %>% dplyr::select(psi1_tru, psi2_tru, 
+      dep_true, model, param, est, se) %>% filter(!is.na(se))

R> temp = filter(clean, param == "ct1_x11", psi2_tru == 
+      0.1)

R> n = 500

R> dep = c(0.15, 0.85)

R> b1 = c(1, 3.25, -2.3)

R> b2 = c(2, -1.75, 3.5)

R> psi1 = psi2 = c(0.1, 0.6)

R> temp$b_true = b1[2]

R> temp$est.c = temp$est - temp$b_true

R> temp = ungroup(filter(temp, between(est.c, -2.25, 
+      2.25)))

R> clrs = RColorBrewer::brewer.pal(3, "Set1")[c(3, 1, 
+      2)]

R> p = ggplot(temp, aes(x = est.c, color = model, fill = model)) + 
+      geom_density(trim = TRUE, alpha = 0.15, bw = "sj") + geom_vline(xintercept  .... [TRUNCATED] 

R> pdf("Figures/figure_2.pdf")

R> print(p)

R> dev.off()
pdf 
  2 

R> print(p)

R> sums = clean %>% filter(param == "ct1_x11", psi2_tru == 
+      0.1) %>% group_by(model, psi1_tru, dep_true) %>% mutate(b_true = b1[2]) %>% 
+      .... [TRUNCATED] 

R> p = ggplot(sums, aes(x = model, y = overconf, fill = model)) + 
+      geom_col(alpha = 0.75) + facet_grid(dep_true ~ psi1_tru, 
+      labeller =  .... [TRUNCATED] 

R> pdf("Figures/figure_3.pdf")

R> print(p)

R> dev.off()
pdf 
  2 

R> print(p)

R> clrs = clrs[-3]

R> dat = filter(clean, param == "dependence", psi2_tru == 
+      0.1)

R> dat$est.c = dat$est - dat$dep_true

R> p = ggplot(data = dat, aes(x = est.c, color = model, 
+      fill = model)) + geom_density(alpha = 0.1, trim = TRUE, bw = "sj") + 
+      geom_vlin .... [TRUNCATED] 

R> pdf("Figures/figure_4.pdf")

R> print(p)

R> dev.off()
pdf 
  2 

R> print(p)

R> message("Simulation plots complete.")
