library(testthat) library(sleev) test_that("logistic2ph Simulation 1", { set.seed(918) # Set sample sizes ---------------------------------------- N <- 1000 # Phase-I = N n <- 250 # Phase-II/audit size = n # Generate true values Y, Xb, Xa -------------------------- Xa <- rbinom(n = N, size = 1, prob = 0.25) Xb <- rbinom(n = N, size = 1, prob = 0.5) Y <- rbinom(n = N, size = 1,prob = (1 + exp(-(- 0.65 - 0.2 * Xb - 0.1 * Xa))) ^ (- 1)) # Generate error-prone Xb* from error model P(Xb*|Xb,Xa) -- sensX <- specX <- 0.75 delta0 <- - log(specX / (1 - specX)) delta1 <- - delta0 - log((1 - sensX) / sensX) Xbstar <- rbinom(n = N, size = 1, prob = (1 + exp(- (delta0 + delta1 * Xb + 0.5 * Xa))) ^ (- 1)) # Generate error-prone Y* from error model P(Y*|Xb*,Y,Xb,Xa) sensY <- 0.95 specY <- 0.90 theta0 <- - log(specY / (1 - specY)) theta1 <- - theta0 - log((1 - sensY) / sensY) Ystar <- rbinom(n = N, size = 1, prob = (1 + exp(- (theta0 - 0.2 * Xbstar + theta1 * Y - 0.2 * Xb - 0.1 * Xa))) ^ (- 1)) ## V is a TRUE/FALSE vector where TRUE = validated -------- V <- seq(1, N) %in% sample(x = seq(1, N), size = n, replace = FALSE) # Build dataset -------------------------------------------- sdat <- cbind(Y, Xb, Ystar, Xbstar, Xa) # Make Phase-II variables Y, Xb NA for unaudited subjects --- sdat[!V, c("Y", "Xb")] <- NA # Fit models ----------------------------------------------- ## (1) Naive model ----------------------------------------- naive <- glm(Ystar ~ Xbstar + Xa, family = "binomial", data = data.frame(sdat)) ## (4) Generalized raking ---------------------------------- ### Influence function for logistic regression ### Taken from: https://github.com/T0ngChen/multiwave/blob/master/sim.r inf.fun <- function(fit) { dm <- model.matrix(fit) Ihat <- (t(dm) %*% (dm * fit$fitted.values * (1 - fit$fitted.values))) / nrow(dm) ## influence function infl <- (dm * resid(fit, type = "response")) %*% solve(Ihat) infl } naive_infl <- inf.fun(naive) # error-prone influence functions based on naive model colnames(naive_infl) <- paste0("if", 1:3) # Add naive influence functions to sdat ----------------------------------------------- sdat <- cbind(id = 1:N, sdat, naive_infl) ### Construct B-spline basis ------------------------------- ### Since Xb* and Xa are both binary, reduces to indicators -- nsieve <- 4 B <- matrix(0, nrow = N, ncol = nsieve) B[which(Xa == 0 & Xbstar == 0), 1] <- 1 B[which(Xa == 0 & Xbstar == 1), 2] <- 1 B[which(Xa == 1 & Xbstar == 0), 3] <- 1 B[which(Xa == 1 & Xbstar == 1), 4] <- 1 colnames(B) <- paste0("bs", seq(1, nsieve)) sdat <- cbind(sdat, B) smle <- logistic2ph(Y_unval = "Ystar", Y = "Y", X_unval = "Xbstar", X = "Xb", Z = "Xa", Bspline = colnames(B), data = sdat, noSE = FALSE, MAX_ITER = 1000, TOL = 1E-4) # [["col"]] is equivalent to unname(unlist(["col"])) expect_equal(smle[["coefficients"]][["Estimate"]], unname(unlist(smle[["coefficients"]]["Estimate"]))) expect_equal(smle[["coefficients"]][["Estimate"]], unname(c(-0.587801556603343, -0.231437717735943, 0.139139424615871)), tolerance=1e-6) expect_equal(smle[["coefficients"]][["SE"]], c(0.159867158182108, 0.240522213900834, 0.203146591074255), tolerance=1e-6) expect_equal(smle[["outcome_err_coefficients"]][["Estimate"]], c(-2.50954051127845, -0.763949870590926, 5.6328754487314, 0.50008005682218, -0.345208034572803), tolerance=1e-6) # expect_equal(as.vector(smle[["Bspline_coefficients"]]), c(1,1,0.744983827508479,0.255016172491521,0.24073358078479,0.75926641921521,0.65244609543555,0.34755390456445,0.352216846922798,0.647783153077202), tolerance=1e-6) expect_equal(as.vector(smle[["covariance"]]), c(0.0255575082652231,-0.0291909056522233,-0.0106240625411377,-0.0291909056522233,0.0578509353797586,-0.000415041049293844,-0.0106240625411377,-0.000415041049293844,0.0412685374650904), tolerance=1e-6) # expect_equal(smle[["od_loglik_at_conv"]], -847.635948985586, tolerance=1e-6) # expect_equal(smle[["iterations"]], 42) # expect_equal(smle[["converged_msg"]], "Converged") expect_true(smle[["converge"]]) expect_true(smle[["converge_cov"]]) }) test_that("logistic2ph simulation 4", { skip_on_cran() skip_if(SKIP_CRAN_TESTS) set.seed(918) # Set sample sizes ---------------------------------------- N <- 1000 # Phase-I = N n <- 250 # Phase-II/audit size = n # Generate true values Y, Xb, Xa ---------------------------- Xa <- rbinom(n = N, size = 1, prob = 0.25) Xb <- rnorm(n = N, mean = 0, sd = 1) Y <- rbinom(n = N, size = 1, prob = (1 + exp(-(- 1 + Xb - 0.5 * Xa))) ^ (- 1)) # Generate error-prone Xb* = Xb + U ------------------------- ## For U ~ N(mean, var) ----------------------------------- muU <- 0 ## mean muU = 0 for unbiased, != 0 for biased ---- s2U <- 0.1 ## variance ------------------------------------ U <- rnorm(n = N, mean = muU, sd = sqrt(s2U)) Xbstar <- Xb + U # Parameters for error model P(Y*|Xb*,Y,Xb,Xa) --------------- ## Varied outcome error rates, ## assumed specificity = sensitivity - 0.05 sensY <- 0.95 # 0.95, 0.85, 0.75, 0.65, 0.55 specY <- 0.90 # 0.90, 0.80, 0.70, 0.60, 0.50 theta0 <- - log(specY / (1 - specY)) theta1 <- - theta0 - log((1 - sensY) / sensY) theta2 <- 1; theta3 <- 1; theta4 <- - 0.5 # Generate error-prone Y* from error model P(Y*|Xb*,Y,Xb,Xa) -- Ystar <- rbinom(n = N, size = 1, prob = (1 + exp(- (theta0 + theta1 * Y + theta2 * Xb + theta3 * Xbstar + theta4 * Xa))) ^ (- 1)) V <- seq(1, N) %in% sample(x = seq(1, N), size = n, replace = FALSE) # Build dataset -------------------------------------------- sdat <- data.frame(cbind(Y, Xb, Ystar, Xbstar, Xa, V)) # Make Phase-II variables Y, Xb NA for unaudited subjects --- sdat[!V, c("Y", "Xb")] <- NA ## (5) SMLE ------------------------------------------------ ### Construct B-spline basis ------------------------------- ### We chose cubic B-splines, ------------------------------ ### with 20 df for N = 1000 and 24 df for N = 2000 --------- nsieve <- 20 B <- matrix(0, nrow = N, ncol = nsieve) B[which(Xa == 0),1:(0.75 * nsieve)] <- splines::bs(x = Xbstar[which(Xa == 0)], df = 0.75 * nsieve, Boundary.knots = range(Xbstar[which(Xa == 0)]), intercept = TRUE) B[which(Xa == 1),(0.75 * nsieve + 1):nsieve] <- splines::bs(x = Xbstar[which(Xa == 1)], df = 0.25 * nsieve, Boundary.knots = range(Xbstar[which(Xa == 1)]), intercept = TRUE) colnames(B) <- paste0("bs", seq(1, nsieve)) sdat <- cbind(sdat, B) smle <- logistic2ph(Y_unval = "Ystar", Y = "Y", X_unval = "Xbstar", X = "Xb", Z = "Xa", Bspline = colnames(B), data = sdat, noSE = FALSE, MAX_ITER = 1000, TOL = 1E-4) expect_equal(smle[["coeff"]][["coeff"]], c(-1.06156561253257,0.91350196339448,-0.473970680353318), tolerance=1e-6) expect_equal(smle[["coeff"]][["se"]], c(0.140532937583599,0.161854198753782,0.312779897372571), tolerance=1e-6) expect_equal(smle[["outcome_err_coeff"]][["coeff"]], c(-2.20598715562998,2.13713609055863,4.64259829097443,-0.336402001120205,-0.180227475248351), tolerance=1e-6) expect_equal(smle[["outcome_err_coeff"]][["se"]], c(NA, NA, NA, NA, NA)) expect_equal(as.vector(smle[["Bspline_coeff"]]), c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,0.873135540339248,0,0.000665030925900439,0,0.0267758649846098,0.0364028508845507,0,0.00372704504818533,0,0.0435125769816376,0.00641615833384844,0,0.00462517789359775,0.000137273533526904,0,0,0.00181450181688563,0,0,0,0,0.000155831208924469,0,4.28236115071513e-06,0.00255782550725901,0,0,0,0,0,0,7.00401806756385e-05,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.247306942600529,0,0.0342926751781834,0,0.0941138404067509,0.104782359954657,0,0.0515839697533328,0,0.112348139439171,0.0602593707697967,0.00274190765600045,0.0552133320830148,0.0278215298253969,0.0113713679281572,0,0.0441258480841587,0,0,0,0.00903705856539449,0.0290635884169333,0,0.021855730345537,0.0493118427213797,0.00236586045132863,3.58067931523538e-05,0,0,0,0.0146852373791687,0.0271073263232447,0,0,0,0,0,0,0,0,0,0,0,0,9.83471279812554e-05,0,0,0,0.000192613711122855,0,0,2.88101275392609e-08,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.49830650764358e-06,0.000283777368972979,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0280067742504664,0,0.0527670673421414,0,0.06160688977229,0.0618020740370741,0,0.0576361539035877,0,0.0620690007187219,0.0592133602144477,0.0241790880395466,0.0585627739952524,0.0509685401647177,0.0394934781848068,0,0.0565998014226952,0,0,0,0.0368498890289141,0.0520448956989846,0.00298779971095603,0.0484358580873744,0.0585573949070004,0.0232793408620137,0.00760031640519269,0,0,0.000312638807257201,0.0434083721263479,0.051583834571141,0,0,0.00332823955319508,0,0.000183889094305212,0,0.000190876000845103,0,0.00302196827207108,0,0.00131176225709951,0,0.00942170204881729,0,0,0,0.0110361491894115,0.00420581172094894,0,0.00457031860036286,1.59028718084491e-10,0.000939053251580361,0,0,3.65179403842503e-08,0.000998623527971229,0.00106776382220891,0,2.38116828753327e-06,6.81675842329237e-06,0,0,0,0,0,9.97478136927625e-05,0,0,7.30586011948816e-06,0,0,0,0,0,0.00545905650927825,0.0123551401256788,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00382800941939119,0,0,0,6.07841268241878e-09,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00136638897644547,0,0.0224189656060695,0,0.0139277222974207,0.0128650719117854,0,0.0194789589976615,0,0.0122549971766381,0.0182147749398814,0.0343907182506732,0.019041300112376,0.024348845777896,0.0297877643211108,0,0.021066579748747,0,0,0,0.0311824118237596,0.0244088020956564,0.0279366926481005,0.0263318165818942,0.0205771030031526,0.0350860762506147,0.0332804531867199,0,0,0.0162117558927629,0.0289849973121234,0.0251806293302733,0.00308498241442266,0,0.0287043525079286,0,0.0144771357558218,0,0.0145948168278672,0.00209575470186935,0.0281303955288347,0,0.0229953328255759,0.000121078479652519,0.0344587297644937,0,0.00566926341080491,0.00165156105697891,0.0350765698090003,0.0302325545259613,0.00357143057331306,0.0307403658534765,0.00662376048025746,0.0211777702337045,0,0.00597239159772081,0.00694748501016091,0.0215171271245653,0.0218801470391064,0,0.00817925756049586,0.00891062870357474,0,0,0,0,0,0.0128953390352256,0.000360591661854556,0.00484994909812122,0.00899065607708682,0,0,0,0,0,0.0319075117571718,0.0356930988907384,0.00368635057597999,0.00270027220981559,0,0,0,0,0,0,0,0,0,0.000965746961879957,0.0021769974115905,0,0,0.00348164589707762,0,0,0.000132067745846294,0,0,0,0,0,0,0,0,0.0299108148437142,0,2.15154952217551e-05,0,0.00689209620409414,0,0,0,0,0.000179628107231725,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00147902597074552,0,0,4.860507765219e-05,0,0,0,0.00389245040616318,0,0.000150780253048681,0,0.0170552664058977,0,0,0.00183064673124963,0.00975372956879548,0,0,0.0289681223536292,1.02738138527865e-05,0,0.0293929590203449,0,0.0163882538556539,0,0.0303104148066606,0,0.0302418468803647,0.0269210043287472,0.0172071330680571,0.00312547907334812,0.0229463359151407,0.0127761747880833,0.00828778996984311,0,0.032434108663831,0.0253603286465003,0.00707904416039218,0.0146675569572162,0.0302675659648889,0.0140199910004807,0.0328457025354305,0.024840454044843,0,0.0326056095061915,0.032927955987795,0.024551839272285,0.0241981122376592,0,0.0330142313148591,0.0329317416513295,0,0,0,0.00499681037202932,0,0.0314243388396309,0.0168083809175764,0.0318541896449473,0.0329456832060305,0,0,0,0.00089703357580758,0,0.0127906082756698,0.00640950878049484,0.0304590728324605,0.0285637018151825,0,0,0.00303671161864888,0,0,0.000476822330110154,0,0,0,0.0219723877174729,0.0271647782611366,0,0,0.0301330984894077,0,0.00236292928786002,0.0130384584201281,0,0,0,0,0,0,0,0,0.0158363637830556,3.74900015275775e-05,0.00905531364628862,0,0.0330232225712495,0.00222358119895628,0,0,0,0.014043951876392,0,0,0,0,0.000117863800677663,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00119293721300889,0,0,0,0,0,0,0,0,0.0048935565690216,0,0,0.00170262176906069,0,0,0,0,0,0,0,0,0,0,7.1903696863619e-06,0,0,0,0,0,8.58585501312865e-07,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0399970095767825,0,0,0,1.21651156352069e-05,0,0,0,0,0,0,0.00206998929525444,0,0,0.0195998933031334,0,4.91819494866436e-06,0,0.00290088629878602,0,0.00284386602138292,0.0238628311002404,1.36827915224516e-05,0.038417551288571,0.000364970621732083,0.0395511204601857,0,0,0.0122825546479906,0.0262177626942847,0,3.46464352628266e-09,0.0178715615685199,0,0.0104703752630448,0.000666801517201118,0.00950306136763963,0.0116756422225772,0.00993573598187755,0.000609127613217378,0.000546060416603046,0,0.00817079831783308,0.00728591984618973,0,0,0,0.0405530465373702,0,0.0038691336510066,0.0361797509423316,0.0141840645857139,0.00721094864041019,9.80293118589735e-05,0.000189861057265643,0,0.0299991374438136,0,0,0,0.0174910680314153,0.0210798053097912,4.80282734120096e-05,0,0.0378989123179369,0,0,0.0262454141301731,0,0,0.0014246887388988,0.0305479861073035,0.0233821945982441,0,0,0.0181611090983169,0,0.0362268346009193,0.0389822035775429,0,0,0.00254745716706201,0,0,0,0,0.0020223486346111,1.47009802571576e-06,0.0169355164466343,0.0406989590833564,0,0.0100920749887959,0.0356572738226433,0.00018478709049527,0,0.00777096115919118,0.0380785985680327,3.69152640343114e-05,0,0,0,0.0199884996864984,0,0.00236168010127216,0.00639968375437168,0,0,0,0,0.00243289660994629,0.00408890887421435,0,0,0,0.00248257312294231,0,0,0,0,0,0.0312344320585682,0.00172656249039649,0,0,0,0,0,0,0.0040213516671357,0.0395807083590578,0,4.9794723455496e-06,0.0334345989729828,0,2.42411638214021e-05,0,0,0,0,0,0,0,0,0.0139840484192722,0,0,5.13229315103664e-06,0,0,0.0122449622850733,0,0,0.0002086660292609,0.003103206373286,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.016466935585271,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00014897855011931,0,0,0,0,0,0,0.00043643428260324,0,0.0184824291871774,0,0.00551299033470171,0,0,1.43996259388686e-06,0.000691928915593624,0,0,8.54784350625462e-05,0,0,0,0.0394096193000666,4.22088409045491e-07,0,0,0,0,0,0,0,0,0,0.014001494431629,0,0,0.00329825413127715,1.2761229823873e-05,0,0.0153518635532261,0.0174358857861389,0,0.0281108689877413,0,0,0,7.66929531192045e-05,0.000237953909966526,0.0136060265471663,0,0.0186703693619967,0,0,0.0313911046190791,0,0,0.0277082804975226,0.00146624275040979,0.00041149806204954,0,0,0.000100025844190929,0,0.0208870181814967,0.0053328406728967,0,0,0.031718324567599,0.00139543830186011,0.00649789538974602,0,0,0.0300644178338251,0,0.0374203873029318,0.00849400836180325,0.0018405316746274,0,0.0213550183895864,0.0172644236170244,0,0.0385330453335341,0.00471315756369577,0.0130200883479142,0.00477203207282255,0.000454401755021597,0.00068530495535718,0.0355840195266637,8.79792673456826e-05,0.0311075475310286,0.0376271591679664,0,0,0,0.000464278402219711,0.031275302278705,0.0349136946514039,0,0,0,0.0314113344374626,0.000142666154426099,0,0,0,0,0.0259868982715172,0.0288295642528031,0,0,0,0,0,0,0.0347353931347867,0.014051997484392,0,0.0103099888653332,0.0233527987206078,0.00776505151780677,0.0122051126815747,0,0,0,0,0,0,0,0,0.0380783976254594,0.00343766275179475,0,0.0103054636366236,1.47943818099533e-11,3.08237579832614e-05,0.0384839386954696,0,0,0.0175366865809253,0.0328010338136822,0.000804585441700981,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00107293384060019,0,0,0,0,3.33738662454534e-05,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.86629263378971e-06,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.7076409605323e-05,0,0,0,0,0,0,0,0,0,0,0,0,0.00949144328022396,0,0,0,0,0,0,0,0,0,0,2.39871512335081e-09,0,0,0,0,0,0.0351101627183694,0.034061344054894,0,0.00067962829490752,0,0,0,0,0,0.0357683769271465,0,3.2068058087756e-05,0,0,0.00130139773799332,0,0,0.0263031718397173,0,0,0,0.000548148293289985,0,0,9.28099221151507e-05,0,0,0,0.0220320995239715,0.0274604656669092,0.0358823645243262,0,0,0.0238116982398465,0,0.00416342956206607,0,0.0290916761675754,0,0.000113436968687608,0.033940040325323,0,0.0113265755519705,0,0.0357873404916351,0.034558882934464,0.0216403214962249,0.0235590932899383,0.0029100470035318,0.0161032020584497,0.0225445015386735,0.0132737950062302,0,0,0,0.0217062175748553,0.0222491781889421,0.0176959783893106,0,0,0,0.0220766330922666,0.0173862203598152,0,0.000981781272147509,0,0,0.000451997530340828,0.0247455871188117,0,0,0,0.000125024329929782,0,0.00487214186677564,0.0177821624297605,3.08329956931707e-08,0,0.0361865833823625,0.000230550290591108,0.0360913091636466,0.035783905791958,0,0,0,0.00929736307389768,0.000131667047630554,0,0,0.00316470606157932,0.00552595309112894,0.0326093015850325,0,0.0360874099328264,0.00932598263910052,0.0138650272904149,0.00664803107915672,3.06694829486391e-05,0,0.0333155654317732,0.0199651825820727,0.0242226147866834,0,0,0.00596805753929825,0,0,0,0.000262241577660075,0,0,0,2.29499439560346e-05,0,0,1.33974365852171e-08,0,0,0,0,0,0,0,0.0256670815010599,0,0,0,0,0.0139054137643114,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.55432279320061e-06,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00660611263440491,0.00540895255952212,0,0,0,0,0,0,0,0.00776197597730606,0,0,0,0,0,0,0,0.00159431776350375,0,0,0,0.0288023094327411,0,0,0,0,0,0,0.000753343103024758,0.030223381155883,0.0154567226962077,0,0,0.00105372647562232,0,0,0,0.0280558937225462,0,0,0.00544543740405167,0,1.59436597199148e-05,0,0.00814314062861919,0.0188249732464314,0.0366772366614051,0.0347048600508145,0,0.0413744255699168,0.000841425753813002,5.57969234623707e-05,0,0,0,0.0365861416249723,0.00080055942468341,0.000282867569814501,0,0,0,0.000775028369213648,0.0404111641995205,0,0.0323017704954065,0,0,0,0.00127033499011378,0,0,0.00124191761335219,0.022457970093197,0,0.0427676974142199,0.000292565321966089,0,0,0.0104364589178983,0,0.0134210505520769,0.00871708844990917,0,0.006500502064765,0,0.0444976516146297,0.0226562603837703,1.62188534780191e-08,0,0.0402437825992255,0,0.0222094011652159,0,0.0104122242896398,0.0445154754828146,0.0428164781226068,0,0.0188095866851886,0,0.00519013183090037,0.000513199448961176,0.0338182655733255,0,0.000140230197357977,0.0437703114215448,0,0.000670652736576829,0,0.0253558004079279,0.0112313726611699,0.0102840674989489,0,0.0182732716970042,0,0.00617943690798123,0.0136242207325346,0,0.00116935102704526,0.00403265958200389,0,0,0,0,0.0320348685496334,0,0.000709769079579144,0.000383945628603074,0.00617692161819713,0.0428231036684187,0,0,3.5799446763697e-06,0,0,0,0.00225157310263037,0,0.000704949560706211,0.00178460035681507,0,0,0,0,0,0,0,0,0,0,0.00123861707039519,0,0.00140957629221977,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0296219675540757,0,0,0,0,0,0,0,0.000701401269577025,2.21997893778815e-06,0,0,0,0,0,0,0.000472553651632636,0,0,0,0,0,0,0,3.09818840004569e-05,0.00191430660590786,0.0014321906671005,0,0.00387301057639944,0,0,0,0,0,0.00188891752931986,0,0,0,0,0,0,0.00332131363911948,0,0.0268501272648067,0,0,0,0,0,0,0.0298898487511455,0.0344572066488636,0,0.0144560595604722,0,0,0,0,0,5.8064988557526e-10,0,0,0.0384612857232883,0,0.00832625030151792,0.034433452694009,0.014019690705338,0.00814019352311245,0.018436271357623,0,0.000115414053282884,0,0,0.0082974648436309,0.00497349488252172,0,0.0367702466409151,0,0,0,0.00125128978916066,0,0.0214146999287907,0.0125294652779072,0,0.0270566924408352,0,0.0327417171388632,0.0394762687615431,0.0395397906466329,0.00276558016632707,0.0371398180989042,0.0101529927002687,0.0385051527996317,0.0390653669380191,0,0.0298295666075773,0.0365949909396078,0.00715791290377889,0.00607212051312623,0,0,0.000945307816614431,0,0.0274547580979347,0.0248597043449755,0.0386871993582739,0.00493002197689506,0,0,0.0159564224063218,0.00902386579158437,0.0131733869247079,0,0.0335760570435491,0,0.0275165801456197,0.0323291015987608,0,0,0.000358060040777005,0,0.000331876810992564,0,0,0,0,1.64721196702842e-05,0.0304475759532617,0.00706084279215634,0.031153470239683,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00190438058176495,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00121155174398122,0,0,0,0,0,0,0.0295765240629565,0.00363858871295563,0,6.06802832711326e-05,0,0,0,0,0,0,0,0,0.0154760827904824,0,1.5224055932882e-12,0.00355944232840809,0.0433163339116563,0.0447065604718688,0.000224636862011431,0,0,0,0,0,0,0,0.00501512317620315,0,0,0,0,0,0.0380178764549119,2.08465384674259e-05,0,0.0326196830143002,0,0.00272600737062381,0.00971119712998417,0.0105987425438163,0.0398816708239181,0.00523540475381655,0.0442640493420325,0.0158298753772635,0.00781123450976568,0.00668661382172917,0.0294692191894631,0.019964561280585,0.0439525237504608,0.0434836795724126,0,0,0,0,0.0320019697860759,0.0345746243538433,0.0156872274668472,0,0,0,0.0414347662017433,0.0437865441855683,0.042833710670017,0,0.0246727495286574,0,0.0317948438597374,0.0262811546904526,0.00568132190897523,0.000857528522824365,0.0285545379177703,0.00285785514740987,0.0281671882401245,0,0.000154302067938417,0,0.00519830926272572,0.0206576204274853,0.0284471514194799,0.0428701707806769,0.0276196369863432,0,0.000662232443319869,0.00527614347366823,0.00537712740508501,0.00239802546941148,2.00377346993827e-05,0,0.000932344352203854,0,1.19772935710042e-08,0,1.07306635365328e-05,0,0,0.0022270426554532,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0010137394992257,0,0,0,0,0,0,0,0,0,0,0,5.83425922793873e-05,0,0,0,0.00603967702779425,0.0108843475294714,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00295755137768937,0,0,0.00155449609041109,0,0,9.34355773335724e-07,3.25358600482078e-06,0.0203589699753612,0,0.00884381326257768,6.99960736275969e-05,0,0.046102815702102,0.00105829950539373,0.000217161890487375,0.0120666756997181,0.0135841112873082,0,0,0,0,0.00149491076149094,0.00205504550533893,7.00028746252704e-05,0,0,0,0.00504548780092842,0.00990491732754741,0.00656322867773867,0,0.000547201054045566,0,0.00149486893543342,0.000712115996892053,0.0457025859893125,0.0420321930327251,0.0315454547317867,0.0449167338154801,0.031765357325949,0,0.038208801441521,0.0256248437847315,0.045514953731705,0.0378089082771651,0.000997199212791542,0.0121110435675543,0.000891286141484695,0,0.0410550216317033,0.0451803751614118,0.0451369366452822,0.044034505794661,0.0351578856239145,0,0.0414939495150696,0,0.0322864897636622,0.00598167684353259,0.0341823067237327,0.0315285857620449,0.0309243067557558,0.0431668413610945,0,0,0.0143419744133683,0.0117945476025193,0.00774247287957804,0,0,0.00322032840506244,0.0183551686980456,0.00358696894522122,0.00100833203286771,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.60679310135829e-05,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00045956118537183,0,5.74354608185199e-06,0,0,0.0163219155445857,0,0,4.8662018448859e-05,8.81858279927007e-05,0,0,0,0,0,0,0,0,0,0,0,1.50658493168065e-05,2.70363780640805e-08,0,0,0,0,0,0.0177610279267895,0.0335688623324525,0.00224853721016676,0.0244607674218799,0.00231829942349444,0,0.0413023837285502,0.05891256289444,0.0185712295676435,0.004584364806904,0,5.31921439211783e-05,0,0,0.0347964014323832,0.0182089993536901,0.017986000334127,0.0256178716197424,0.0455235227142295,0,0.0324067430107688,0,0.0490011545275818,0.0497164206425502,0.0458500008573144,0.049878999332144,0.0506291771290424,0.0256999633644427,0,0,0.0613619916485146,0.0597265396464681,0.0537123887132609,0,0,0.038037431093995,0.0604173203955071,0.0397211462607946,0.0209614715280121,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000258830971189284,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000371824526462737,0.00330835773036715,0,0.00117791857425127,0,0,0.00655252990429387,0.0287306429545008,0.000445547088942411,0,0,0,0,0,0.00378093701118167,0.000421108303129628,0.000402716889927171,0.00141213867583087,0.00943085264600061,0,0.00307668799539557,0,0.0128383214218848,0.127795696843871,0.0100274576660831,0.0139110774353108,0.014871608619763,0.00148607106273976,0,0,0.0708966981319771,0.0854633948744668,0.112703245776243,0,0,0.146772508676228,0.0496959924101356,0.143475303748574,0.15069253006125,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000622068396815279,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.100681369047758,0,1.30687526180483e-07,1.54284429521584e-06,0,0,0,0.0160827644692476,0.0282581671268846,0.0670739420093113,0,0,0.194159612423996,0.00544915106241153,0.175253416170545,0.412417835761209,0,0,0.157827934281085,0,0.0545434235616411,0,0,0.115339871948223,0,0.147939738306762,0,0,0,0,0,0,0.10664446526079,0,0.0533226513535226,0,0.0597094847922927,0,0,0,0,0,0,0,0.0689603268057841,0.0169151813275227,0,0,0,0,0.024905099392443,0,0.0260873075865942,0,0.003026565901202,0,0,0,0,0,0,0,0.0134301749508685,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0244696377627851,0,0,0.00257959363582232,0.00397266796748703,0.0168526328907544,0,0.000482997552234244,0,0,0,0,0,0,0.0218761698916998,0,0.0120891609082774,0,0,0,0,0,0,0,0.00441540228911942,0.00220510637535283,0,0.0051306093054737,0.000916166792638535,0,0,0,0.0101378406347505,0,0,0.00401363143344554,0,0,0.0211410435021098,6.19869119281349e-05,0,0,0,0.0105878328403954,0.00505595551859188,0,0,0,0,0,0,0,0,1.39761964385531e-05,0,0,0,0,0,0,0,0,0,0,0.00200062988454317,0.000711747870738238,0.000452206252059969,0,0,0,4.94282460100616e-06,0.00017822853615301,0,0,0,0.000546224841397989,0,0.000536696462534278,0,0,0,0.000738841108725442,0.000172298712296065,0,0,0,0,0,0,0,0,0,0,0,3.54562891505981e-06,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.022981676660235,0,0.0268232614439888,0,0,0.025390971670881,0,0.0235605230242345,0,0,0,0,0,0,0.0257160689401942,0,0.0267096960576386,0,0.0267435844950841,0,0,0,0,0,0,0,0.0266599674003603,0.0239645035132647,0,0,0,0,0.0250608773200815,0,0.0251742331864314,0,0.0191136238175201,0,0,0,0,0,0,0,0.0231921300992138,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0249269067664903,0,0,0.0187033938071169,0.0197172044305255,0.0238183871820583,0,0.0158119261763477,0,0,0,0,0,0,0.0245645220821317,0,0.022783845686492,0,0,0,0,0,0.0108498326856614,0,0.0199194716331508,0.0183032405096084,0,0.0203009712441583,0.0166824582308734,0,0,0,0.022215643487656,0,0,0.019657026887322,0,0,0.0243671916110964,0.0138522032617943,0,0,0,0.0222839363229485,0.0202001691030251,0,0,0,0,0,0,0,0,0.01311159483343,0,0,0,0,0,0,0,0,0,0,0.0179777648716758,0.0161873282279517,0.0155657461218696,0,0,0,0.0127383743152357,0.0145712153745845,0.00892769695244952,0,0,0.0157925822095957,0,0.0157662980884726,0.011351946317657,0,0,0.0162012252193961,0.0145159778927478,0,0,0.0111158408267525,0,0,0,0.00893684737667822,0,0,0,0,0.0126031417807001,0,0.00912593164179596,0,0,0,0,0,0,0,0.00948755551771467,0,0,0,0,0,0.00808533763844686,0,0,0,0.0104708817041872,0,0,0.00800880102494937,0,0.0113775775233761,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00857423476462861,0.00838332817725385,0,0.00352885193438619,0,0,0,0,0.00525382934678465,0.00372371030280928,0,0,0,0.00280408240179495,0,0.00610935453240833,0,0,0,0,0,0,0,0.00314516731121214,0,0,0,0,0,0,0,0.00266612749411734,0,0,0,0,0,0.00653731995442799,0,0.00269366506914411,0,0,0,0,0,0,0.00120769985008285,0.00238630218547441,0,0,0,0.000604566856469505,0.0003002945746018,0,0,0,0,0.000112351047152053,0,0.00337754495416173,0,0.00909720152024084,0,0,0.00508369270085278,0,0.00370035034569253,0,0,0,0,0,0,0.00548202695587093,0,0.00911416766554096,0,0.00853782408520004,0,0,0,0,0,0,0,0.00776885663482147,0.0135495453094255,0,0,0,0,0.0122283056785542,0,0.0120470760719111,0,0.0170158067098059,0,0,0,0,0,0,0,0.014141331214367,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0122201657306656,0,0,0.0171439415483377,0.0165815719906081,0.0134315862904574,0,0.0183825048088504,0,0,0,0,0,0,0.0125746897483473,0,0.0143352326897266,0,0,0,0,0,0.0192461071066713,0,0.0163506943232537,0.0172408808925184,0,0.0161017739857753,0.0179692749827504,0,0,0,0.0147336956659754,0,0,0.0164700785022331,0,0,0.0126165010648553,0.0188211307060838,0,0,0,0.0145777392280281,0.0160643175770348,0,0,0,0,0,0,0,0,0.0189486007720679,0,0,0,0,0,0,0,0,0,0,0.0172054155849065,0.0179886449045883,0.0182013939366971,0,0,0,0.018933893014653,0.0185126765043981,0.018827291294371,0,0,0.0180968457988265,0,0.0181020236146064,0.0190491134377657,0,0,0.0179234034862838,0.0184920200588556,0,0,0.0190395757871202,0,0,0,0.0187891228642521,0,0,0,0,0.0188764173932001,0,0.0188125485117568,0,0,0,0,0,0,0,0.0188530658169472,0,0,0,0,0,0.0184666589802202,0,0,0,0.0189359193735647,0,0,0.0184231680060102,0,0.01891405534442,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0185708252085995,0.0185099996571974,0,0.0144933824579116,0,0,0,0,0.0166130801439809,0.0147802615621723,0,0,0,0.0131984010752198,0,0.0173170169222017,0,0,0,0,0,0,0,0.0138266547773883,0,0,0,0,0,0,0,0.0129071083026461,0,0,0,0,0,0.0175474292642543,0,0.0129612988021732,0,0,0,0,0,0,0.00889990304977381,0.0123092441226381,0,0,0,0.00617850035969796,0.00415935679659836,0,0,0,0,0.00233006632341733,0,0.000317841881829826,0,0.00185552910560009,0,0,0.000637026882037333,0,0.000369003326203604,0,0,0,0,0,0,0.000725655425202573,0,0.00185667509726184,0,0.00163550771118624,0,0,0,0,0,0,0,0.00136687579479542,0.0043093244981236,0,0,0,0,0.00341317265265583,0,0.0033041555259845,0,0.00801420415308411,0,0,0,0,0,0,0,0.00479402907893777,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00341869301751478,0,0,0.0082709926501503,0.00742376580467373,0.00424796743421419,0,0.010920051069066,0,0,0,0,0,0,0.00365203914538861,0,0.00498645769389256,0,0,0,0,0,0.0171334206709283,0,0.00716105907092348,0.00850835251353629,0,0.00684497924706557,0.00997011743865362,0,0,0,0.0053630049296972,0,0,0.00734117796728616,0,0,0.00369804766038841,0.012863250396179,0,0,0,0.00524153392438446,0.00684188990565559,0,0,0,0,0,0,0,0,0.0137205389638905,0,0,0,0,0,0,0,0,0,0,0.00860476989328295,0.0102551694031532,0.0108565632358222,0,0,0,0.0140845161156878,0.011897896947803,0.0202470601304208,0,0,0.0106027076686255,0,0.0106243754232051,0.0159934696502844,0,0,0.0101800102357252,0.0119159232818944,0,0,0.0163316206729727,0,0,0,0.0201423949535532,0,0,0,0,0.0141534245004117,0,0.019732527257331,0,0,0,0,0,0,0,0.0189926268221625,0,0,0,0,0,0.0217566286612055,0,0,0,0.0172226194017156,0,0,0.0218897583045778,0,0.0157537675206957,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0206318109964925,0.0210179839305441,0,0.0335759855058843,0,0,0,0,0.0285288601870688,0.0329662040893158,0,0,0,0.0357491075364132,0,0.0262287238011498,0,0,0,0,0,0,0,0.0347270943541951,0,0,0,0,0,0,0,0.0362093658232623,0,0,0,0,0,0.0251166947026013,0,0.0363050497991255,0,0,0,0,0,0,0.0402123251769425,0.0374672210267777,0,0,0,0.0400650638685737,0.0373708577439281,0,0,0,0,0.0323834807408035,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.6726752690855e-05,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000483762063035983,0,0,0,0,0,3.09598830736569e-06,0,0,0,0,0,0,8.10073484014526e-06,0,0,0,0.000464641660683364,0,0,0,0,0,0,0.000368555678349424,0,0,0,0,0,0,0,0.000227833951154821,0,0,0,0,0,0.00103534477295074,0,0,0,4.13464574919621e-05,0,0,0.00109750983405131,0,1.54292917071832e-06,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000619552069405691,0.000748690898093885,0,0.0215714484929796,0,0,0,0,0.00778421467127231,0.0191884258151295,0,0,0,0.0325283338955783,0,0.00446812079266808,0,0,0,0,0,0,0,0.026560923417563,0,0,0,0,0,0,0,0.0350311802463579,0,0,0,0,0,0.00317440729491338,0,0.0343509323617785,0,0,0,0,0,0,0.0896817695982358,0.0413757601263463,0,0,0,0.150679945272468,0.215844188421209,0,0,0,0,0.312633645803274), tolerance=1e-6) expect_equal(as.vector(smle[["vcov"]]), c(0.0197495065458757,-0.00599685066646061,-0.017684772949906,-0.00599685066646061,0.0261967816542286,-0.00256929509198125,-0.017684772949906,-0.00256929509198125,0.0978312642003962), tolerance=1e-6) expect_equal(smle[["od_loglik_at_conv"]], -1488.70457730476, tolerance=1e-6) expect_equal(smle[["iterations"]], 47) expect_equal(smle[["converged_msg"]], "Converged") expect_true(smle[["converged"]]) expect_true(smle[["se_converged"]]) }) test_that("missing values", { set.seed(918) # Set sample sizes ---------------------------------------- N <- 1000 # Phase-I = N n <- 250 # Phase-II/audit size = n # Generate true values Y, Xb, Xa -------------------------- Xa <- rbinom(n = N, size = 1, prob = 0.25) Xb <- rbinom(n = N, size = 1, prob = 0.5) Y <- rbinom(n = N, size = 1,prob = (1 + exp(-(- 0.65 - 0.2 * Xb - 0.1 * Xa))) ^ (- 1)) # Generate error-prone Xb* from error model P(Xb*|Xb,Xa) -- sensX <- specX <- 0.75 delta0 <- - log(specX / (1 - specX)) delta1 <- - delta0 - log((1 - sensX) / sensX) Xbstar <- rbinom(n = N, size = 1, prob = (1 + exp(- (delta0 + delta1 * Xb + 0.5 * Xa))) ^ (- 1)) # Generate error-prone Y* from error model P(Y*|Xb*,Y,Xb,Xa) sensY <- 0.95 specY <- 0.90 theta0 <- - log(specY / (1 - specY)) theta1 <- - theta0 - log((1 - sensY) / sensY) Ystar <- rbinom(n = N, size = 1, prob = (1 + exp(- (theta0 - 0.2 * Xbstar + theta1 * Y - 0.2 * Xb - 0.1 * Xa))) ^ (- 1)) ## V is a TRUE/FALSE vector where TRUE = validated -------- V <- seq(1, N) %in% sample(x = seq(1, N), size = n, replace = FALSE) # Build dataset -------------------------------------------- sdat <- cbind(Y, Xb, Ystar, Xbstar, Xa) # Make Phase-II variables Y, Xb NA for unaudited subjects --- sdat[!V, c("Y", "Xb")] <- NA # Fit models ----------------------------------------------- ## (1) Naive model ----------------------------------------- naive <- glm(Ystar ~ Xbstar + Xa, family = "binomial", data = data.frame(sdat)) ## (4) Generalized raking ---------------------------------- ### Influence function for logistic regression ### Taken from: https://github.com/T0ngChen/multiwave/blob/master/sim.r inf.fun <- function(fit) { dm <- model.matrix(fit) Ihat <- (t(dm) %*% (dm * fit$fitted.values * (1 - fit$fitted.values))) / nrow(dm) ## influence function infl <- (dm * resid(fit, type = "response")) %*% solve(Ihat) infl } naive_infl <- inf.fun(naive) # error-prone influence functions based on naive model colnames(naive_infl) <- paste0("if", 1:3) # Add naive influence functions to sdat ----------------------------------------------- sdat <- cbind(id = 1:N, sdat, naive_infl) ### Construct B-spline basis ------------------------------- ### Since Xb* and Xa are both binary, reduces to indicators -- nsieve <- 4 B <- matrix(0, nrow = N, ncol = nsieve) B[which(Xa == 0 & Xbstar == 0), 1] <- 1 B[which(Xa == 0 & Xbstar == 1), 2] <- 1 B[which(Xa == 1 & Xbstar == 0), 3] <- 1 B[which(Xa == 1 & Xbstar == 1), 4] <- 1 colnames(B) <- paste0("bs", seq(1, nsieve)) sdat <- cbind(sdat, B) ## missing Y_unval # no error # expect_error(logistic2ph(Y = "Y", # X_unval = "Xbstar", # X = "Xb", # Z = "Xa", # Bspline = colnames(B), # data = sdat), regexp = NA) ## missing Y expect_error(logistic2ph( Y_unval = "Ystar", # Y = "Y", X_unval = "Xbstar", X = "Xb", Z = "Xa", Bspline = colnames(B), data = sdat)) ## missing X_unval expect_error(logistic2ph( Y_unval = "Ystar", Y = "Y", # X_unval = "Xbstar", X = "Xb", Z = "Xa", Bspline = colnames(B), data = sdat)) ## missing X expect_error(logistic2ph( Y_unval = "Ystar", Y = "Y", X_unval = "Xbstar", # X = "Xb", Z = "Xa", Bspline = colnames(B), data = sdat)) ## missing Z # no error # expect_error(logistic2ph( # Y_unval = "Ystar", # Y = "Y", # X_unval = "Xbstar", # X = "Xb", # # Z = "Xa", # Bspline = colnames(B), # data = sdat), regexp = NA) ## missing Bspline expect_error(logistic2ph( Y_unval = "Ystar", Y = "Y", X_unval = "Xbstar", X = "Xb", Z = "Xa", # Bspline = colnames(B), data = sdat)) ## missing data expect_error(logistic2ph( Y_unval = "Ystar", Y = "Y", X_unval = "Xbstar", X = "Xb", Z = "Xa", Bspline = colnames(B))) })