Use statmod package instead of hardcoded structs

DESCRIPTION

@@ -21,5 +21,6 @@ Suggests:
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
 Depends: 
-    R (>= 2.10)
+    R (>= 2.10),
+    statmod
 VignetteBuilder: knitr

NAMESPACE

@@ -2,6 +2,7 @@
 
 export(simulate_zero_inflated_beta_random_effect_data)
 export(zibr)
+importFrom(statmod,gauss.quad)
 importFrom(stats,nlminb)
 importFrom(stats,pchisq)
 importFrom(stats,rbeta)

R/fit_beta_regression_random_effect.R

@@ -1,7 +1,7 @@
 calc_beta_loglik <- function(para, Z.aug, Y, subject.n, time.n,
                             prod.mat = prod.mat,
                             Z.test.coeff.index,
                             gh.weights, gh.nodes, quad.n) {
   ### Y must be a vector here
   Y <- as.vector(Y)
   s2 <- para[1]
@@ -30,7 +30,7 @@
   #### logY==0 -> infinite, we want to exclude Y==0 in the likelihood calculation
   #### replace with infinite with 0, it will be ignored, e^0*e^log(p2)= p2
   log.i[is.infinite(log.i)] <- 0
   logL <- sum(
     log(rowSums(gh.weights / sqrt(pi) * exp(prod.mat %*% log.i),
       na.rm = TRUE
     )),
@@ -56,8 +56,9 @@
 #' }
 #'
 #' @importFrom stats nlminb pchisq
+#' @importFrom statmod gauss.quad
 fit_beta_random_effect <- function(Z = Z, Y = Y,
                                    subject.ind = subject.ind, time.ind = time.ind,
                                    quad.n = 30, verbose = FALSE) {
   Z <- as.matrix(Z)
   Y <- as.matrix(Y)
@@ -79,7 +80,7 @@
                      ncol = subject.n * time.n)
 
   #### generate quad points
-  gherm <- generate_gaussian_quad_points(quad.n)
+  gherm <- gauss.quad(quad.n, kind = "hermite")
   gh.weights <- matrix(rep(gherm$weights, subject.n), nrow = subject.n, byrow = TRUE)
   gh.nodes <- matrix(rep(gherm$nodes, subject.n * time.n),
     nrow = subject.n * time.n, byrow = TRUE

R/fit_logistic_regression_random_effect.R

@@ -45,6 +45,7 @@ calc_logistic_loglik <- function(para, X.aug, Y, subject.n, time.n,
 #' }
 #'
 #' @importFrom stats nlminb pchisq
+#' @importFrom statmod gauss.quad
 fit_logistic_random_effect <- function(X = X, Y = Y,
                                        subject.ind = subject.ind, time.ind = time.ind,
                                        quad.n = 30, verbose = FALSE) {
@@ -68,7 +69,7 @@ fit_logistic_random_effect <- function(X = X, Y = Y,
                      ncol = subject.n * time.n)
 
   #### generate quad points
-  gherm <- generate_gaussian_quad_points(quad.n)
+  gherm <- gauss.quad(quad.n, kind = "hermite")
   gh.weights <- matrix(rep(gherm$weights, subject.n), nrow = subject.n, byrow = TRUE)
   gh.nodes <- matrix(rep(gherm$nodes, subject.n * time.n),
     nrow = subject.n * time.n, byrow = TRUE

R/fit_zero_inflated_beta_regression_random_effect.R

@@ -56,6 +56,7 @@ calc_zibeta_loglik <- function(para,
 #' }
 #'
 #' @importFrom stats nlminb pchisq
+#' @importFrom statmod gauss.quad
 fit_zero_inflated_beta_random_effect <- function(X = X, Z = Z, Y = Y,
                                                  subject_ind = subject_ind, time_ind = time_ind,
                                                  component_wise_test = TRUE,
@@ -81,7 +82,7 @@ fit_zero_inflated_beta_random_effect <- function(X = X, Z = Z, Y = Y,
                      ncol = subject.n * time.n)
 
   #### generate quad points
-  gherm <- generate_gaussian_quad_points(quad_n)
+  gherm <- gauss.quad(quad_n, kind = "hermite")
   gh.weights <- matrix(rep(gherm$weights, subject.n), nrow = subject.n, byrow = TRUE)
   gh.nodes <- matrix(rep(gherm$nodes, subject.n * time.n),
     nrow = subject.n * time.n, byrow = TRUE