diff --git a/DESCRIPTION b/DESCRIPTION index 285d90c..9cbcf0d 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -21,5 +21,6 @@ Suggests: testthat (>= 3.0.0) Config/testthat/edition: 3 Depends: - R (>= 2.10) + R (>= 2.10), + statmod VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 2e0d785..5ff8bcb 100755 --- a/NAMESPACE +++ b/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) diff --git a/R/fit_beta_regression_random_effect.R b/R/fit_beta_regression_random_effect.R index 7eb72be..9fa9258 100755 --- a/R/fit_beta_regression_random_effect.R +++ b/R/fit_beta_regression_random_effect.R @@ -56,6 +56,7 @@ calc_beta_loglik <- function(para, Z.aug, Y, subject.n, time.n, #' } #' #' @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) { @@ -79,7 +80,7 @@ fit_beta_random_effect <- function(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 diff --git a/R/fit_logistic_regression_random_effect.R b/R/fit_logistic_regression_random_effect.R index 35999db..dd8eb11 100755 --- a/R/fit_logistic_regression_random_effect.R +++ b/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 diff --git a/R/fit_zero_inflated_beta_regression_random_effect.R b/R/fit_zero_inflated_beta_regression_random_effect.R index 58ecd38..e52b8d3 100755 --- a/R/fit_zero_inflated_beta_regression_random_effect.R +++ b/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 diff --git a/R/generate_gaussian_quad_points.R b/R/generate_gaussian_quad_points.R deleted file mode 100755 index e14cec8..0000000 --- a/R/generate_gaussian_quad_points.R +++ /dev/null @@ -1,151 +0,0 @@ -generate_gaussian_quad_points = function(quad.n = 30) { - ## library(statmod) - ## gherm <- gauss.quad(100, kind = "hermite") - if (quad.n == 10) { - gherm <- - structure(list( - nodes = c( - -3.43615911883774,-2.53273167423279,-1.75668364929988,-1.03661082978951,-0.342901327223705, 0.342901327223705, - 1.03661082978951, 1.75668364929988, - 2.53273167423279, 3.43615911883774 - ), - weights = c( - 7.64043285523267e-06, 0.00134364574678123, - 0.0338743944554811, 0.240138611082315, - 0.610862633735326, 0.610862633735326, - 0.240138611082315, 0.0338743944554811, - 0.00134364574678123, 7.64043285523262e-06 - ) - ), - .Names = c("nodes", "weights")) - } - else if (quad.n == 20) { - gherm <- structure(list( - nodes = c( - -5.38748089001124,-4.60368244955074,-3.94476404011562,-3.34785456738321,-2.78880605842813,-2.25497400208927,-1.73853771211659,-1.23407621539532,-0.737473728545395,-0.245340708300901,0.245340708300901, 0.737473728545395, 1.23407621539532, 1.73853771211659,2.25497400208927, 2.78880605842813, 3.34785456738321, 3.94476404011562,4.60368244955074, 5.38748089001124 - ), - weights = c( - 2.22939364553406e-13, - 4.39934099227314e-10, 1.08606937076928e-07, 7.80255647853211e-06, - 0.000228338636016354, 0.00324377334223786, 0.0248105208874636, - 0.109017206020023, 0.286675505362834, 0.46224366960061, 0.46224366960061, - 0.286675505362834, 0.109017206020023, 0.0248105208874636, 0.00324377334223786, - 0.000228338636016356, 7.80255647853206e-06, 1.08606937076928e-07, - 4.39934099227317e-10, 2.22939364553411e-13 - ) - ), .Names = c("nodes","weights")) - } - else if (quad.n == 30) { - gherm <- - structure(list( - nodes = c( - -6.86334529352989,-6.13827922012394,-5.53314715156749,-4.98891896858994,-4.48305535709252,-4.00390860386123,-3.54444387315535,-3.09997052958644,-2.66713212453562,-2.2433914677615,-1.82674114360369,-1.41552780019819,-1.00833827104672,-0.603921058625552,-0.201128576548872, 0.201128576548871, 0.603921058625552, 1.00833827104672, - 1.41552780019819, 1.82674114360369, 2.2433914677615, 2.66713212453562, - 3.09997052958644, 3.54444387315535, 4.00390860386123, 4.48305535709252, - 4.98891896858994, 5.5331471515675, 6.13827922012394, 6.86334529352989 - ), weights = c( - 2.90825470013124e-21, 2.81033360275083e-17, 2.87860708054873e-14, - 8.10618629746283e-12, 9.17858042437863e-10, 5.10852245077593e-08, - 1.57909488732472e-06, 2.93872522892296e-05, 0.000348310124318684, - 0.00273792247306766, 0.0147038297048267, 0.0551441768702344, - 0.14673584754089, 0.280130930839213, 0.386394889541813, 0.386394889541814, - 0.280130930839212, 0.14673584754089, 0.0551441768702342, 0.0147038297048267, - 0.00273792247306767, 0.000348310124318684, 2.93872522892297e-05, - 1.57909488732472e-06, 5.10852245077599e-08, 9.1785804243787e-10, - 8.10618629746318e-12, 2.87860708054873e-14, 2.81033360275085e-17, - 2.90825470013124e-21 - ) - ), .Names = c("nodes", "weights")) - } - else if (quad.n == 50) { - gherm <- structure(list( - nodes = c( - -9.18240695812932,-8.5227710309178,-7.97562236820564,-7.4864094298642,-7.03432350977062,-6.60864797385536,-6.20295251927467,-5.8129946754204,-5.43578608722496,-5.06911758491724,-4.71129366616904,-4.36097316045458,-4.01706817285814,-3.67867706251527,-3.34503831393789,-3.01549776957452,-2.68948470226775,-2.36649390429866,-2.04607196868641,-1.7278065475159,-1.4113177548983,-1.09625112895768,-0.782271729554607,-0.469059056678236,-0.156302546889468, 0.156302546889468, - 0.469059056678236, 0.782271729554606, 1.09625112895768, 1.4113177548983, - 1.7278065475159, 2.04607196868641, 2.36649390429866, 2.68948470226775, - 3.01549776957452, 3.34503831393789, 3.67867706251527, 4.01706817285813, - 4.36097316045457, 4.71129366616904, 5.06911758491724, 5.43578608722495, - 5.8129946754204, 6.20295251927468, 6.60864797385537, 7.03432350977062, - 7.4864094298642, 7.97562236820564, 8.5227710309178, 9.18240695812932 - ), weights = c( - 1.83379404857344e-37, 1.67380166790783e-32, 1.21524412340454e-28, - 2.13765830836019e-25, 1.41709359957333e-22, 4.47098436540803e-20, - 7.74238295704311e-18, 8.09426189346557e-16, 5.46594403181537e-14, - 2.50665552389962e-12, 8.11187736493057e-11, 1.90904054381193e-09, - 3.34679340402142e-08, 4.45702996681781e-07, 4.58168270795558e-06, - 3.68401905378075e-05, 0.000234269892109254, 0.00118901178174964, - 0.00485326382617197, 0.0160319410684122, 0.0430791591567656, - 0.0945489354770863, 0.170032455677164, 0.251130856332002, 0.305085129204399, - 0.305085129204398, 0.251130856332004, 0.170032455677164, 0.0945489354770862, - 0.0430791591567658, 0.0160319410684122, 0.00485326382617194, - 0.00118901178174965, 0.000234269892109256, 3.68401905378074e-05, - 4.58168270795551e-06, 4.45702996681783e-07, 3.34679340402146e-08, - 1.90904054381202e-09, 8.11187736492997e-11, 2.50665552389959e-12, - 5.46594403181576e-14, 8.09426189346522e-16, 7.74238295704339e-18, - 4.47098436540784e-20, 1.4170935995734e-22, 2.13765830836017e-25, - 1.21524412340456e-28, 1.67380166790793e-32, 1.83379404857339e-37 - ) - ), .Names = c("nodes", "weights")) - } - else if (quad.n == 100) { - gherm <- - structure(list( - nodes = c( - -13.4064873381449,-12.8237997494878,-12.3429642228597,-11.9150619431142,-11.521415400787,-11.1524043855851,-10.8022607536847,-10.4671854213428,-10.1445099412928,-9.83226980777795,-9.5289658233901,-9.23342089021916,-8.94468921732547,-8.66199616813451,-8.38469694041627,-8.11224731116279,-7.84418238446082,-7.58010080785749,-7.31965282230454,-7.06253106024886,-6.80846335285879,-6.55720703192154,-6.30854436111214,-6.0622788326143,-5.81823213520352,-5.57624164932992,-5.33615836013836,-5.09784510508914,-4.86117509179121,-4.62603063578716,-4.39230207868269,-4.15988685513103,-3.92868868342767,-3.69861685931849,-3.46958563641859,-3.24151367963101,-3.01432358033115,-2.78794142398199,-2.56229640237261,-2.33732046390688,-2.11294799637119,-1.88911553742701,-1.66576150874151,-1.44282597021593,-1.22025039121895,-0.997977436098106,-0.775950761540146,-0.554114823591618,-0.332414692342232,-0.11079587242244, - 0.110795872422439, 0.332414692342232, 0.554114823591617, 0.775950761540145, - 0.997977436098105, 1.22025039121895, 1.44282597021593, 1.66576150874151, - 1.88911553742701, 2.11294799637119, 2.33732046390688, 2.5622964023726, - 2.78794142398199, 3.01432358033115, 3.24151367963101, 3.46958563641859, - 3.69861685931849, 3.92868868342767, 4.15988685513103, 4.39230207868269, - 4.62603063578716, 4.86117509179121, 5.09784510508914, 5.33615836013836, - 5.57624164932993, 5.81823213520351, 6.0622788326143, 6.30854436111214, - 6.55720703192153, 6.80846335285879, 7.06253106024886, 7.31965282230453, - 7.58010080785749, 7.84418238446082, 8.11224731116279, 8.38469694041626, - 8.66199616813451, 8.94468921732548, 9.23342089021915, 9.52896582339012, - 9.83226980777795, 10.1445099412928, 10.4671854213428, 10.8022607536847, - 11.1524043855851, 11.521415400787, 11.9150619431142, 12.3429642228597, - 12.8237997494878, 13.4064873381449 - ), weights = c( - 5.90806786503149e-79, - 1.97286057487953e-72, 3.08302899000321e-67, 9.01922230369242e-63, - 8.51888308176111e-59, 3.45947793647603e-55, 7.19152946346349e-52, - 8.59756395482676e-49, 6.42072520534849e-46, 3.18521787783596e-43, - 1.10047068271428e-40, 2.74878488435709e-38, 5.11623260438594e-36, - 7.27457259688812e-34, 8.06743427870884e-32, 7.10181222638517e-30, - 5.03779116621273e-28, 2.91735007262926e-26, 1.39484152606877e-24, - 5.56102696165936e-23, 1.86499767513029e-21, 5.30231618313167e-20, - 1.28683292112113e-18, 2.68249216476057e-17, 4.82983532170314e-16, - 7.5488968779154e-15, 1.02887493735098e-13, 1.22787851441009e-12, - 1.28790382573158e-11, 1.19130063492903e-10, 9.74792125387112e-10, - 7.07585728388942e-09, 4.568127508485e-08, 2.62909748375372e-07, - 1.35179715911036e-06, 6.22152481777778e-06, 2.56761593845487e-05, - 9.51716277855096e-05, 0.000317291971043304, 0.000952692188548621, - 0.00257927326005907, 0.00630300028560806, 0.0139156652202317, - 0.0277791273859335, 0.0501758126774289, 0.0820518273912242, 0.121537986844105, - 0.163130030502782, 0.198462850254188, 0.218892629587438, 0.21889262958744, - 0.198462850254186, 0.163130030502783, 0.121537986844104, 0.082051827391225, - 0.0501758126774289, 0.0277791273859336, 0.0139156652202318, 0.00630300028560809, - 0.00257927326005912, 0.000952692188548612, 0.000317291971043303, - 9.51716277855086e-05, 2.5676159384549e-05, 6.22152481777782e-06, - 1.35179715911039e-06, 2.62909748375376e-07, 4.56812750848495e-08, - 7.07585728388942e-09, 9.74792125387167e-10, 1.19130063492907e-10, - 1.28790382573154e-11, 1.22787851441012e-12, 1.02887493735101e-13, - 7.5488968779154e-15, 4.82983532170362e-16, 2.68249216476036e-17, - 1.28683292112121e-18, 5.30231618313197e-20, 1.86499767513026e-21, - 5.56102696165912e-23, 1.39484152606877e-24, 2.91735007262916e-26, - 5.03779116621305e-28, 7.10181222638506e-30, 8.06743427870919e-32, - 7.2745725968875e-34, 5.1162326043855e-36, 2.74878488435732e-38, - 1.10047068271418e-40, 3.18521787783605e-43, 6.42072520534922e-46, - 8.59756395482676e-49, 7.1915294634638e-52, 3.45947793647628e-55, - 8.51888308176039e-59, 9.01922230369063e-63, 3.08302899000303e-67, - 1.97286057487992e-72, 5.90806786503182e-79 - ) - ), .Names = c("nodes", - "weights")) - } - else { - cat("quad.n should be 10,20,30,50 or 100. Use 30 instead.\n") - gherm <- generate_gaussian_quad_points(quad.n = 30) - } - return(gherm) -}