Discrepancies in robust ancova calculated in R with the WRS2 package vs by hand

Question

I've conducted a study where participants in 2 groups (Control and Experimental) completed three scales (NJ, A, and D) at two time points (T1 and T2). I'm running ANCOVAs with the T2 responses as my dependent variables, T1 responses as my covariates, and Group as my independent variables (i.e. NJ2 ~ NJ1 + Condition, A2 ~ A1 + Condition, and D2 ~ D1 + Condition). Since my slopes are heterogeneous, I used the ancova function from the WRS2 package to perform a robust ANCOVA for each.

After realizing that the ancova function doesn't provide degrees of freedom or effect sizes, I decided to try running the Yuen comparisons by hand by selecting subsets of my data that are close to my comparison points of interest, calculating the trimmed means of these subsets, and comparing those trimmed means using the yuen function also from the WRS2 package. On page 28 of the package vignette (https://cran.r-project.org/web/packages/WRS2/vignettes/WRS2.pdf), the authors give the following formula for identifying data points that are close to a point of interest: |Xi − x| ≤ f × MADN, where Xi are the values of X, x is the point of interest, f is a smoothing parameter (which can be dropped in my case because I've set it to 1), and MADN is the median absolute deviation of X divided by the .75 quantile of the standard normal distribution. Using this formula, however, has produced different results for each of my three outcome variables.

1) For my NJ variable, I was able to replicate the results from the ancova function by rounding the MADN from 3.71 up to 4, which seems reasonable.

2) For the A variable, I was able to replicate the results from the ancova function but had to round the MADN from 5.93 down to 4.

3) For the D variable, rounding the MADN from 5.93 down to 4 gave me the same n value as ancova for one of my groups but I had to round the MADN down to 3 to get the same n value for my other group. Combining the two approaches to get the right n values for both groups though, produced a different mean difference value than the ancova function.

Does anyone have any idea why there would be a discrepancy between the ancova function and my by-hand calculations in terms of how data subsets are being selected? Is there something in the function "black-box" that I'm missing?

My data file is available here and samples of my code are available below:

For the NJ ANCOVA using a trim level of .2, smoothing values of 1, and comparing at points NJ1 = 11, 13, 15, 17, and 19

NJANCOVA <- ancova(data = Data, NJ2 ~ NJ1 + Condition, tr = .2, fr1 = 1, fr2 = 1, pts = c(11, 13, 15, 17, 19))
NJANCOVA

And here is my code for calculating the tests by hand:

#Grand median:

NJMedian <- median(Data$NJ1)

# Calculate median absolute deviation:

Data$NJDevs <- abs(Data$NJ1 - NJMedian)
NJMAD <- median(Data$NJDevs)

# Divide by .75 quantile:

NJMADN <- NJMAD/qnorm(.75)

# Determine differences from the point of comparison:

NJDiffs11 <- abs(Data$NJ1 - 11)

# Compare to MADN:

NJ11 <- subset(Data, NJDiffs11 <= round(NJMADN, 0))

# Calculate Yuen's test:

yuen(data = NJ11, NJ2 ~ Condition, tr = .2)

At NJ1 = 11, for example, both chunks of code give me n1 = 12, n2 = 16, mean difference = -2.95, test statistic = 1.8276, and p = 0.0864.

For the D variable, I used D1 = 2, 4, and 6 as comparison points.

At D1 = 2, for example, the ancova function gives me n1 = 21, n2 = 13, mean difference = .7692, test statistic = .6739, and p = .5081.

My "by hand" calculations, however, give me n1 = 26, n2 = 13, mean difference = 1.5, test statistic = 1.3966, and p = .17692.