R version 4.4.0 RC (2024-04-16 r86435 ucrt) -- "Puppy Cup" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > require(glmMisrep) Loading required package: glmMisrep > > data <- data.frame( Y = c(0, 0, 2, 0, 3, 0, 36, 0, 2, 1, 0, 2, 6, 9, 0, 0, 0, 0, 7, 1, 1, 2, 50, 4, 0, 0, 0, 1, + 0, 0, 0, 3, 0, 0, 1, 0, 1, 3, 176, 0, 0, 0, 0, 0, 2, 286, 0, 0, 0, 18), + X1 = c(0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1), + X2 = c(1.71870212, -0.55840901, 1.22589915, 0.53000107, 0.62571132, 0.02873955, 0.30989954, 1.35514993, 0.15587503, + 0.27987513, 0.48892178, 0.35218767, 0.52382778, 1.58126751, -0.07855081, -0.57128802, -0.92500953, -2.48543328, + 0.03810910, 0.39929906, -0.54854763, -0.10505694, 0.45120734, 0.32295222, -0.68595918, -0.66892486, 1.72253431, + -0.28425276, -0.67719912, -0.39644260, -0.16843500, 0.90540261, -1.38574804, 0.14456841, 0.44142810, -1.89442541, + -0.65961894, 2.13148776, 1.72410805, -1.60207312, -1.09525034, -1.31327168, -0.43378445, 1.33644956, 1.28938359, + 0.90232362, -0.94112768, -0.61851917, 0.37033085, -0.47019541), + X3 = c(0.6882029, 0.9934165, 0.9173388, 0.9406660, 0.5130041, 0.8590187, 0.4468488, 0.4186652, 0.5098278, 0.3339481, 0.6922477, + 0.6793977, 0.4983724, 0.6079911, 0.7763041, 0.8529067, 0.8287771, 0.9125900, 0.4802076, 0.8981448, 0.3570093, 0.9209584, + 0.4353817, 0.9426418, 0.9550002, 0.4869851, 0.9560156, 0.8247537, 0.1939687, 0.6103839, 0.7721900, 0.5980044, 0.8683831, + 0.7004518, 0.8577210, 0.3576712, 0.9540088, 0.9880046, 0.4304899, 0.7416618, 0.6656063, 0.8920356, 0.6097593, 0.8008748, + 0.8110432, 0.5967969, 0.6983106, 0.9471680, 0.9929186, 0.9059200), + Sex = c("Female", "Male", "Male", "Female", "Male", "Female", "Male", "Male", "Female", "Female", "Male", "Male", + "Male", "Female", "Male", "Male", "Female", "Female", "Male", "Female", "Female", "Male", "Female", "Male", + "Male", "Male", "Male", "Female", "Female", "Female", "Male", "Male", "Male", "Male", "Female", "Female", + "Female", "Female", "Male", "Male", "Male", "Male", "Male", "Male", "Female", "Male", "Male", "Female", + "Female", "Male"), + Race = c("Black", "White", "White", "White", "White", "White", "Black", "Black", "Black", "Black", "Black", "Black", "White", "Other", + "Black", "White", "Black", "White", "Other", "Black", "White", "Other", "Other", "White", "White", "White", "Black", "Other", + "White", "Other", "White", "Black", "Other", "Black", "Other", "Other", "Other", "Black", "Other", "Black", "White", "Other", + "Black", "White", "Black", "Other", "White", "Black", "Black", "Other"), + V_star = c(0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0)) > > data$Race <- as.factor(data$Race) > data$Sex <- as.factor(data$Sex) > > t1 <- tryCatch(poisRegMisrepEM(formula = y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6,0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > # The response above is inappropriately specified (y, not Y) > stopifnot( + t1$message == "object 'y' not found" + ) > > > t2 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_Star", + data = data, + lambda = c(0.6,0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > # Argument to 'v_star' is misspelled > stopifnot( + t2$message == "variable V_Star not present in dataframe" + ) > > > data$V_star <- ifelse(data$V_star == 1, yes = "yes", no = "no") > > t3 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6,0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > > # v* variable is type character (yes and no) > stopifnot( + t3$message == "v_star variable must be of class 'factor' or 'numeric'" + ) > > data$V_star <- ifelse(data$V_star == "yes", yes = 1, no = 0) > > > data$V_star[10] <- -1 > > > t4 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6,0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > # v* variable must be binary > stopifnot( + t4$message == "v_star variable must contain two unique values" + ) > > data$V_star[10] <- 0 > > > data$V_star <- ifelse(data$V_star == 1, yes = 1, no = 2) > > t5 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6,0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > # v* must be binary, but more specifically 0/1; > stopifnot( + t5$message == "v_star variable must be coded with ones and zeroes" + ) > > > > data$V_star <- ifelse(data$V_star == 1, yes = 1, no = 0) > > t6 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.49, 0.52), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > # Inappropriately specified lambda argument > stopifnot( + t6$message == "Lambda vector must sum to one" + ) > > > t7 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(1/3, 1/3, 1/3), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > # Inappropriately specified lambda argument > stopifnot( + t7$message == "Lambda vector must contain two elements" + ) > > > > data$X4 <- data$X2*0.3 > > t8 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + X4 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6, 0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > # Linearly dependent covariates/degenerate design matrix > stopifnot( + t8$message == "Linear dependencies exist in the covariates" + ) > > > t9 <- tryCatch(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race, + v_star = "V_star", + data = data, + lambda = c(0.6, 0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + error = function(x) x ) > > #V_star variable absent from formula > stopifnot( + t9$message == "v_star variable must be specified in 'formula'" + ) > > > > # EM algorithm should fail to converge within the specified number of attempts > t10 <- tryCatch( + capture.output(poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6, 0.4), + epsilon = 1e-08, + maxit = 2, + maxrestarts = 1)), + error = function(x) x + ) > > stopifnot( + t10$message == "NOT CONVERGENT! Failed to converge after 1 attempts" + ) > > > # On the first attempt, fails to converge, and restarts with new mixing props. > # Succeeds on the second attempt. > msg <- capture.output( + t11 <- poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6, 0.4), + epsilon = 1e-08, + maxit = 16, + maxrestarts = 4, verb = TRUE), + type = "message" + ) > > stopifnot( + any(msg == "Warning: Failed to converge. Restarting with new mixing proportions") + ) > > > msg <- capture.output( + t12 <- poisRegMisrepEM(formula = Y ~ X1 + X2 + X3 + Sex + Race + V_star, + v_star = "V_star", + data = data, + lambda = c(0.6, 0.4), + epsilon = 1e-08, + maxit = 10000, + maxrestarts = 20), + type = "message" + ) > > # Output validation; > > # Output should be a list > stopifnot( + is.list(t12) + ) > > # With 14 elements > stopifnot( + length(t12) == 14 + ) > > # Fisher information matrix should be symmetric > stopifnot( + isSymmetric(t12$cov.estimates) + ) > > > # The returned list should have elements with the following names > # and types > stopifnot( + all.equal(lapply(t12, class), + lapply(list(y = 0.1, lambda = 0.2, params = 0.3, loglik = 0.4, + posterior = 0.5, all.loglik = 0.6, cov.estimates = matrix(data = c(1,2,3,4), 2, 2), + std.error = 0.7, z.values = 0.8, p.values = 0.9, ICs = 1.0, ft = "*", + formula = y ~ x, v_star_name = "v*" ), class) ) + ) > > > # Verifying the function can correctly calculate things; > stopifnot( + all.equal(t12$lambda, 0.2140361, tolerance = 3e-7 ) + ) > > stopifnot( + all.equal(as.numeric(t12$params), c(-1.1486407, 1.9813423, 1.0606001, -4.0278727, 1.8886687, 2.1052795, -0.1815366, 2.2582766), tolerance = 3e-7 ) + ) > > stopifnot( + all.equal( t12$loglik, -77.17481, tolerance = 3e-7 ) + ) > > stopifnot( + all.equal( t12$posterior, c(9.131063e-01, 9.167867e-01, 9.987564e-01, 8.007038e-01, 1.000000e+00, 7.980578e-01, 1.761218e-22, 9.999998e-01, + 5.707215e-02, 4.986677e-01, 9.591568e-01, 9.998818e-01, 1.596141e-04, 1.000000e+00, 8.834283e-01, 8.269226e-01, + 8.267814e-01, 7.909421e-01, 9.995369e-01, 3.002507e-01, 7.724813e-01, 5.033734e-01, 1.413033e-29, 7.033365e-04, + 9.184187e-01, 9.118807e-01, 9.752457e-01, 4.103920e-01, 8.589374e-01, 9.278243e-01, 8.652112e-01, 1.000000e+00, + 9.998403e-01, 9.999676e-01, 5.418147e-01, 8.817133e-01, 3.273919e-01, 1.253503e-01, 2.099810e-102, 8.124995e-01, + 8.350465e-01, 9.998171e-01, 9.066149e-01, 1.000000e+00, 4.335937e-01, 2.097168e-171, 8.364673e-01, 7.911482e-01, + 7.981111e-01, 1.769426e-10), tolerance = 3e-7 ) + ) > > > stopifnot( + all.equal( t12$all.loglik, c(-91.21048, -77.54635, -77.24049, -77.18672, -77.17811, -77.17606, -77.17532, -77.17502, -77.17490, -77.17484, -77.17482, + -77.17481, -77.17481, -77.17481, -77.17481, -77.17481, -77.17481, -77.17481, -77.17481, -77.17481, -77.17481), tolerance = 3e-7 ) + ) > > > stopifnot( + all.equal(t12$cov.estimates, + matrix(data = c(0.0114492285, -0.0004087145, 0.0016883055, -0.001674646, -0.01493787, 0.004194090, 0.004620007, 0.003212934, 0.0008121915, + -0.0004087145, 0.1376023242, -0.0103938967, -0.021763262, -0.10804979, 0.001850089, -0.001523425, -0.033947804, -0.0502985311, + 0.0016883055, -0.0103938967, 0.0153809909, 0.005402047, -0.02063417, 0.003663544, 0.001754054, 0.008152043, 0.0006989913, + -0.0016746456, -0.0217632624, 0.0054020474, 0.012549504, 0.02752248, -0.006829300, -0.008499404, 0.001601437, 0.0045586028, + -0.0149378678, -0.1080497865, -0.0206341716, 0.027522475, 0.26797912, -0.040935476, -0.039438494, -0.023606980, 0.0234374153, + 0.0041940900, 0.0018500887, 0.0036635440, -0.006829300, -0.04093548, 0.024832614, 0.009305140, 0.004723165, -0.0062878481, + 0.0046200066, -0.0015234254, 0.0017540537, -0.008499404, -0.03943849, 0.009305140, 0.033156029, 0.023238618, -0.0090039233, + 0.0032129343, -0.0339478042, 0.0081520434, 0.001601437, -0.02360698, 0.004723165, 0.023238618, 0.095149535, 0.0125163473, + 0.0008121915, -0.0502985311, 0.0006989913, 0.004558603, 0.02343742, -0.006287848, -0.009003923, 0.012516347, 0.0479783934), + ncol = 9, nrow = 9, byrow = TRUE, dimnames = list( c("lambda", names(t12$params)), c("lambda", names(t12$params)) ) ), tolerance = 3e-7 + ) + ) > > stopifnot( + all.equal(as.numeric(t12$std.error), c(0.1070011, 0.3709479, 0.1240201, 0.1120246, 0.5176670, 0.1575837, 0.1820880, 0.3084632, 0.2190397), tolerance = 3e-7) + ) > > stopifnot( + all.equal(as.numeric(t12$z.values), c(-3.0965016, 15.9759747, 9.4675669, -7.7808179, 11.9851800, 11.5618810, -0.5885196, 10.3098958), tolerance = 3e-7) + ) > > > stopifnot( + all.equal(as.numeric(t12$p.values), c(1.958188e-03, 1.878943e-57, 2.864376e-21, 7.205715e-15, 4.249202e-33, 6.428496e-31, 5.561836e-01, 6.357023e-25), tolerance = 3e-7) + ) Error: as.numeric(t12$p.values) and c(0.001958188, 1.878943e-57, 2.864376e-21, 7.205715e-15, 4.249202e-33, .... are not equal: Mean relative difference: 7.240945e-07 Execution halted