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.

> library(distfreereg)
> 
> ### Grouping
> 
> set.seed(20250610)
> # Define a data frame with repeated covariate values.
> df <- as.data.frame(rbind(
+   c(1,3),# 1:  1, 1
+   c(2,3),# 2:  3, 2
+   c(2,4),# 3:  4, 2
+   c(1,4),# 4:  2, 1
+   c(1,4),# 5:  2, 1
+   c(1,4),# 6:  2, 1
+   c(1,3),# 7:  1, 1
+   c(2,4),# 8:  4, 2
+   c(2,4),# 9:  4, 2
+   c(2,3) # 10: 3, 2
+ ))
> 
> n <- nrow(df)
> df <- cbind(data.frame(Y = rnorm(n)), df)
> names(df) <- c("y", "x1", "x2")
> 
> form <- y ~ x1 + x2
> m <- lm(formula = form, data = df)
> 
> dfr_form_lm_natural <- distfreereg(test_mean = m, ordering = "natural",
+                                    group = TRUE, B = 10)
'x' and/or 'y' not found in lm object; refitting linear model...
Extracting covariance structure from lm object...
Calculating Jacobian from 'lm' object...
Calculating Jacobian from lm object...
Retrieving fitted values from 'lm' object...
Calculating the inverse square root of the covariance matrix...
Calculating mu...
Ordering observations by natural order...
Calculating transformation anchors...
Calculating r_tilde...
Calculating residuals...
Determining grouping...
Calculating empirical partial sum process...
Calculating observed statistic(s)...
Running Monte Carlo simulation...
...1 of 10
...2 of 10
...3 of 10
...4 of 10
...5 of 10
...6 of 10
...7 of 10
...8 of 10
...9 of 10
...10 of 10
> 
> dfr_form_lm_natural
Number of observations: 10
Length of EPSP: 4 
Monte Carlo simulations: 10 

Estimated parameter values:
(Intercept)          x1          x2 
 -1.881e+00   1.894e-01   3.694e-01 

Observed statistics:
Stat      Value  Pr(>Value)       MCSE 
  KS  6.214e-01   4.000e-01  1.549e-01 
 CvM  1.328e-01   4.000e-01  1.549e-01 
---
`MCSE' is the Monte Carlo standard error of the estimated p-value.

> stopifnot(identical(dfr_form_lm_natural[["res_order"]],
+                     as.integer(c(1, 3, 4, 2, 2, 2, 1, 4, 4, 3))))
> 
> # Verify that epsp from grouped dfr is correctly obtained from transformed residuals
> groupit_natural_1 <- ((sum(dfr_form_lm_natural$residuals$transformed[c(1,7)])/sqrt(2))/sqrt(4))
> stopifnot(all.equal(groupit_natural_1, dfr_form_lm_natural$epsp[1]))
> 
> groupit_natural_2 <- groupit_natural_1 + ((sum(dfr_form_lm_natural$residuals$transformed[4:6])/sqrt(3))/sqrt(4))
> stopifnot(all.equal(groupit_natural_2, dfr_form_lm_natural$epsp[2]))
> 
> groupit_natural_3 <- groupit_natural_2 + ((sum(dfr_form_lm_natural$residuals$transformed[c(2,10)])/sqrt(2))/sqrt(4))
> stopifnot(all.equal(groupit_natural_3, dfr_form_lm_natural$epsp[3]))
> 
> groupit_natural_4 <- groupit_natural_3 + ((sum(dfr_form_lm_natural$residuals$transformed[c(3,8,9)])/sqrt(3))/sqrt(4))
> stopifnot(all.equal(groupit_natural_4, dfr_form_lm_natural$epsp[4]))
> 
> 
> # Now test it when grouping by a subset of columns.
> 
> dfr_form_lm_x1 <- distfreereg(test_mean = m, ordering = list("x1"),
+                               group = TRUE, B = 10)
'x' and/or 'y' not found in lm object; refitting linear model...
Extracting covariance structure from lm object...
Calculating Jacobian from 'lm' object...
Calculating Jacobian from lm object...
Retrieving fitted values from 'lm' object...
Calculating the inverse square root of the covariance matrix...
Calculating mu...
Ordering observations by specified columns...
Calculating transformation anchors...
Calculating r_tilde...
Calculating residuals...
Determining grouping...
Calculating empirical partial sum process...
Calculating observed statistic(s)...
Running Monte Carlo simulation...
...1 of 10
...2 of 10
...3 of 10
...4 of 10
...5 of 10
...6 of 10
...7 of 10
...8 of 10
...9 of 10
...10 of 10
> 
> dfr_form_lm_x1
Number of observations: 10
Length of EPSP: 2 
Monte Carlo simulations: 10 

Estimated parameter values:
(Intercept)          x1          x2 
 -1.881e+00   1.894e-01   3.694e-01 

Observed statistics:
Stat      Value  Pr(>Value)       MCSE 
  KS  5.964e-02   9.000e-01  9.487e-02 
 CvM  1.779e-03   9.000e-01  9.487e-02 
---
`MCSE' is the Monte Carlo standard error of the estimated p-value.

> stopifnot(identical(dfr_form_lm_x1[["res_order"]],
+                     as.integer(c(1, 2, 2, 1, 1, 1, 1, 2, 2, 2))))
> 
> groupit_x1_1 <- ((sum(dfr_form_lm_x1$residuals$transformed[c(1,4,5,6,7)])/sqrt(5))/sqrt(2))
> stopifnot(all.equal(groupit_x1_1, dfr_form_lm_x1$epsp[1]))
> 
> groupit_x1_2 <- groupit_x1_1 + ((sum(dfr_form_lm_x1$residuals$transformed[c(2,3,8,9,10)])/sqrt(5))/sqrt(2))
> stopifnot(all.equal(groupit_x1_2, dfr_form_lm_x1$epsp[2]))
> 
> proc.time()
   user  system elapsed 
   1.90    0.17    2.06