R Under development (unstable) (2025-06-30 r88369 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 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. > if(require("lme4")){ + library(distfreereg) + all.equal.distfreereg <- distfreereg:::all.equal.distfreereg + test_dfr_functions <- distfreereg:::test_dfr_functions + + n <- 1e2 + set.seed(20250516) + df_lmerMod <- data.frame(x = rnorm(n), + z = rnorm(n), + state = factor(sample(c("CA", "OR", "WA"), size = n, + replace = TRUE, prob = c(2,1,1))), + party = factor(sample(c("D", "R"), size = n, replace = TRUE)), + g = rep(1:10, 10)) + state_num <- ifelse(df_lmerMod$state == "CA", 7, ifelse(df_lmerMod$state == "OR", 2, 3)) + party_num <- ifelse(df_lmerMod$party == "D", 1, -1) + theta <- c(2,1) + df_lmerMod$y <- theta[1]*df_lmerMod$x + theta[2]*df_lmerMod$z + state_num + rnorm(n) + + m <- lmer(y ~ x + (1|state), data = df_lmerMod) + + set.seed(20250516) + dfr_lmerMod <- distfreereg(m, verbose = FALSE, control = list(return_on_error = FALSE)) + print(dfr_lmerMod) + + set.seed(20250516) + dfr_lmerMod_J <- distfreereg(test_mean = m, + override = list(J = dfr_lmerMod[["J"]]), + control = list(return_on_error = FALSE)) + + set.seed(20250516) + dfr_lmerMod_fitted <- distfreereg(test_mean = m, + override = list(fitted_values = dfr_lmerMod[["fitted_values"]]), + control = list(return_on_error = FALSE)) + + stopifnot(all.equal(dfr_lmerMod, dfr_lmerMod_J)) + stopifnot(all.equal(dfr_lmerMod, dfr_lmerMod_fitted)) + + set.seed(20250516) + newdata_lmerMod <- data.frame(y = rnorm(10), x = rnorm(10), z = rnorm(10)) + test_dfr_functions(dfr_lmerMod, newdata = newdata_lmerMod) + print(signif(predict(dfr_lmerMod, newdata = newdata_lmerMod, re.form = NA), 3)) + + set.seed(20250516) + dfr_form_lmerMod <- distfreereg(y ~ x + (1|state), data = df_lmerMod, + method = "lmer", verbose = FALSE) + + print(dfr_form_lmerMod) + + set.seed(20250516) + dfr_form_lmerMod_verbose <- distfreereg(y ~ x + (1|state), data = df_lmerMod, + method = "lmer") + + stopifnot(all.equal(dfr_lmerMod, dfr_form_lmerMod)) + + cdfr_form_lmerMod <- asymptotics(dfr_form_lmerMod, reps = 5) + cdfr_lmerMod <- asymptotics(dfr_lmerMod, reps = 5) + + print(signif(rejection(cdfr_form_lmerMod, alpha = c(0.1, 0.5))[,2:3], digits = 3)) + print(signif(rejection(cdfr_lmerMod, alpha = c(0.1, 0.5))[,2:3], digits = 3)) + + # Orderings + + set.seed(20250516) + dfr_lmerMod_asis <- update(dfr_lmerMod, ordering = "asis") + set.seed(20250516) + dfr_form_lmerMod_asis <- update(dfr_form_lmerMod, ordering = "asis") + stopifnot(all.equal(dfr_lmerMod_asis, dfr_form_lmerMod_asis)) + + set.seed(20250516) + dfr_lmerMod_optimal <- update(dfr_lmerMod, ordering = "optimal") + set.seed(20250516) + dfr_form_lmerMod_optimal <- update(dfr_form_lmerMod, ordering = "optimal") + stopifnot(all.equal(dfr_lmerMod_optimal, dfr_form_lmerMod_optimal)) + + set.seed(20250516) + dfr_lmerMod_natural <- update(dfr_lmerMod, ordering = "natural") + set.seed(20250516) + dfr_form_lmerMod_natural <- update(dfr_form_lmerMod, ordering = "natural") + stopifnot(all.equal(dfr_lmerMod_natural, dfr_form_lmerMod_natural)) + + set.seed(20250516) + dfr_lmerMod_g <- update(dfr_lmerMod, ordering = list("g")) + set.seed(20250516) + dfr_form_lmerMod_g <- update(dfr_form_lmerMod, ordering = list("g")) + stopifnot(all.equal(dfr_lmerMod_g, dfr_form_lmerMod_g)) + + print(df_lmerMod[dfr_lmerMod_g[["res_order"]],][["g"]]) + + set.seed(20250516) + dfr_lmerMod_g_grouped <- update(dfr_lmerMod_g, group = TRUE) + set.seed(20250516) + dfr_form_lmerMod_g_grouped <- update(dfr_form_lmerMod_g, group = TRUE) + stopifnot(all.equal(dfr_lmerMod_g_grouped, dfr_form_lmerMod_g_grouped)) + + + + + ### Partial output + + dfr_lmerMod_partial <- distfreereg(test_mean = m, verbose = FALSE, + control = list(orth_tol = 1e-100)) + names(dfr_lmerMod_partial) + } Loading required package: lme4 Loading required package: Matrix Number of observations: 100 Length of EPSP: 100 Monte Carlo simulations: 10000 Estimated parameter values: (Intercept) x 3.825e+00 2.027e+00 Observed statistics: Stat Value Pr(>Value) MCSE KS 6.846e-01 6.372e-01 4.808e-03 CvM 1.036e-01 5.613e-01 4.962e-03 --- `MCSE' is the Monte Carlo standard error of the estimated p-value. Extracting covariance structure from lmerMod object... Using supplied Jacobian... Retrieving fitted values from 'lmerMod' object... Calculating the inverse square root of the covariance matrix... Calculating mu... Ordering observations by simplex method... Calculating transformation anchors... Calculating r_tilde... Calculating residuals... All covariate observations are unique; no grouping done... Calculating empirical partial sum process... Calculating observed statistic(s)... Running Monte Carlo simulation... ...1000 of 10000 ...2000 of 10000 ...3000 of 10000 ...4000 of 10000 ...5000 of 10000 ...6000 of 10000 ...7000 of 10000 ...8000 of 10000 ...9000 of 10000 ...10000 of 10000 Extracting covariance structure from lmerMod object... Calculating Jacobian from 'lmerMod' object... Calculating Jacobian from lmerMod object... Using supplied fitted values... Calculating the inverse square root of the covariance matrix... Calculating mu... Ordering observations by simplex method... Calculating transformation anchors... Calculating r_tilde... Calculating residuals... All covariate observations are unique; no grouping done... Calculating empirical partial sum process... Calculating observed statistic(s)... Running Monte Carlo simulation... ...1000 of 10000 ...2000 of 10000 ...3000 of 10000 ...4000 of 10000 ...5000 of 10000 ...6000 of 10000 ...7000 of 10000 ...8000 of 10000 ...9000 of 10000 ...10000 of 10000 (Intercept) x 3.83 2.03 Computing profile confidence intervals ... 2.5 % 97.5 % .sig01 1.05 5.97 .sigma 1.16 1.54 (Intercept) 0.56 7.09 x 1.74 2.31 1 2 3 4 5 6 7 8 9 10 11 4.580 0.780 2.400 4.330 0.358 2.740 3.870 5.580 3.300 3.810 4.560 12 13 14 15 16 17 18 19 20 21 22 3.490 3.650 8.320 3.500 2.270 4.980 5.700 6.350 2.150 5.540 2.290 23 24 25 26 27 28 29 30 31 32 33 3.730 4.070 -1.820 4.830 5.580 1.840 2.310 3.660 2.310 1.480 5.720 34 35 36 37 38 39 40 41 42 43 44 3.390 4.590 2.470 2.550 0.773 0.226 3.050 3.630 2.640 5.440 4.340 45 46 47 48 49 50 51 52 53 54 55 5.100 4.220 3.940 2.920 3.030 3.420 3.900 2.020 0.902 3.430 2.030 56 57 58 59 60 61 62 63 64 65 66 5.260 4.350 5.600 3.290 4.720 7.100 2.750 4.390 4.730 2.810 3.080 67 68 69 70 71 72 73 74 75 76 77 3.940 0.499 1.560 5.460 3.100 5.600 3.480 4.740 4.520 7.180 3.950 78 79 80 81 82 83 84 85 86 87 88 3.000 4.530 -0.565 7.250 4.270 4.170 6.540 -1.700 5.480 1.700 2.620 89 90 91 92 93 94 95 96 97 98 99 2.410 4.530 7.940 1.780 5.970 3.130 5.570 2.560 3.930 2.480 4.260 100 2.370 1 2 3 4 5 6 7 8 9 10 7.3500 3.5500 1.5600 3.4900 3.1200 1.9000 1.9400 4.7400 6.0600 6.5800 11 12 13 14 15 16 17 18 19 20 7.3300 2.6500 6.4200 11.1000 2.6600 5.0300 7.7500 3.7800 4.4200 4.9200 21 22 23 24 25 26 27 28 29 30 3.6100 1.4500 6.5000 2.1500 -2.6600 3.9900 8.3500 1.0000 0.3870 2.8200 31 32 33 34 35 36 37 38 39 40 5.0700 4.2400 8.4800 1.4700 7.3600 5.2400 5.3200 -0.0671 -1.7000 5.8200 41 42 43 44 45 46 47 48 49 50 6.3900 5.4000 4.6000 2.4200 7.8700 6.9800 2.0200 2.0800 2.1900 6.1900 51 52 53 54 55 56 57 58 59 60 1.9700 4.7900 3.6700 2.5900 1.1900 3.3300 7.1100 8.3600 1.3600 7.4800 61 62 63 64 65 66 67 68 69 70 5.1800 0.8260 7.1600 2.8000 1.9700 1.1500 6.7100 3.2700 -0.3640 3.5300 71 72 73 74 75 76 77 78 79 80 2.2600 3.6800 6.2400 3.9000 2.5900 9.9500 2.0200 5.7700 3.6900 -2.4900 81 82 83 84 85 86 87 88 89 90 6.4100 2.3400 6.9400 5.7000 -3.6200 4.6400 -0.2310 0.6980 5.1800 7.2900 91 92 93 94 95 96 97 98 99 100 7.1000 -0.1490 8.7400 5.8900 3.6500 5.3300 3.0900 0.5500 7.0300 5.1400 NULL 2 x 2 Matrix of class "dsyMatrix" (Intercept) x (Intercept) 2.05000 0.00201 x 0.00201 0.02120 1 2 3 4 5 6 7 8 9 10 4.56 3.49 3.65 8.32 3.50 2.27 4.98 5.70 6.35 2.15 Number of observations: 100 Length of EPSP: 100 Monte Carlo simulations: 10000 Estimated parameter values: (Intercept) x 3.825e+00 2.027e+00 Observed statistics: Stat Value Pr(>Value) MCSE KS 6.846e-01 6.372e-01 4.808e-03 CvM 1.036e-01 5.613e-01 4.962e-03 --- `MCSE' is the Monte Carlo standard error of the estimated p-value. Extracting covariance structure from lmerMod object... Calculating Jacobian from 'lmerMod' object... Calculating Jacobian from lmerMod object... Retrieving fitted values from 'lmerMod' object... Calculating the inverse square root of the covariance matrix... Calculating mu... Ordering observations by simplex method... Calculating transformation anchors... Calculating r_tilde... Calculating residuals... All covariate observations are unique; no grouping done... Calculating empirical partial sum process... Calculating observed statistic(s)... Running Monte Carlo simulation... ...1000 of 10000 ...2000 of 10000 ...3000 of 10000 ...4000 of 10000 ...5000 of 10000 ...6000 of 10000 ...7000 of 10000 ...8000 of 10000 ...9000 of 10000 ...10000 of 10000 Running simulation... Repetition 1 of 5 Repetition 2 of 5 Repetition 3 of 5 Repetition 4 of 5 Repetition 5 of 5 Running simulation... Repetition 1 of 5 Repetition 2 of 5 Repetition 3 of 5 Repetition 4 of 5 Repetition 5 of 5 alpha rate 1 0.1 0.0 2 0.5 0.2 3 0.1 0.0 4 0.5 0.2 alpha rate 1 0.1 0.2 2 0.5 0.4 3 0.1 0.0 4 0.5 0.4 [1] 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 [26] 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 [51] 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 [76] 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 [1] "call" "data" "test_mean" [4] "covariance" "theta_hat" "optimization_output" [7] "fitted_values" "J" "" Warning message: Error encountered, partial results returned: Error in calc_mu(J = J_for_mu, solve_tol = solve_tol, orth_tol = orth_tol): crossprod(mu) is not equal to the identity matrix > > proc.time() user system elapsed 15.90 0.45 16.37