How to pass paired argument to summaryM function test in Hmisc [r]

Question

I have reproducible example here

library(Hmisc)
set.seed(173)
sex <- factor(sample(c("m","f"), 500, rep=TRUE))
country <- factor(sample(c('US', 'Canada'), 500, rep=TRUE))
age <- rnorm(500, 50, 5)
sbp <- rnorm(500, 120, 12)
label(sbp) <- 'Systolic BP'
units(sbp) <- 'mmHg'
treatment <- factor(rep(c("PreTretment","PostTretment"), 250))

f <- summaryM(age + sex + sbp ~ treatment, test=TRUE)

SummaryM function from Hmisc package has test argument which as default applies Wilcoxon test to continues variable assuming they are independent. Now I would like to pass paired=TRUE to Wilcoxon. How can I do it? Thanks

Answer 1

My current effort at making a paired Wilcoxon test with Hmisc summary methods:

 conTestWP <- 
function (group, x) 
{
    st <- wilcox.test( x[as.numeric(group)==1], x[as.numeric(group)==2], paired=TRUE)
    list(P = st["p.value"], stat = st["statistic"], df = st[c("df1", "df2")], 
        testname = "Wilcoxon-paired",
        statname = "V", latexstat = "V", plotmathstat = "F[df]")}

The summaryM method splits its grouping variables and is therefore not appropriate for paired tests. The summary.formula set of methods does allow a "reverse" method where the continuous variable is on the RHS of the formula:

f <- summary.formula(treatment ~ sbp, data=dat, method="reverse", 
                      test=TRUE,conTest=conTestWP)

Trying to print f throws an error (claiming falsely that the p-value is not numeric) but you can look inside to see that the paired wilcox.test results were passed into the object and the are the same as if you had done them "by hand":

 str(f) # but did snip some of the output:

 $ testresults:List of 1
  ..$ sbp:List of 7
  .. ..$ P           :List of 1
  .. .. ..$ p.value: num 0.589
  .. ..$ stat        :List of 1
  .. .. ..$ NA: NULL
  .. ..$ df          :List of 2
  .. .. ..$ NA: NULL
  .. .. ..$ NA: NULL
  .. ..$ testname    : chr "Wilcoxon-paired"
  .. ..$ statname    : chr "V"
  .. ..$ latexstat   : chr "V"
  .. ..$ plotmathstat: chr "F[df]"

Efforts to fix the error that gets thrown by putting in hard-coded numbers for the "df" values are failing. I have not succeeding in following the traceback which I paste in only the top of:

> f
Error in log10(ifelse(pv > 0, pv, 1e-50)) : 
  non-numeric argument to mathematical function
> traceback()
5: format.pval(pval, digits = pdig, eps = eps)
4: formatTestStats(tr, prtest = prtest, latex = latex, testUsed = testUsed, 
       pdig = pdig, eps = eps, footnoteTest = footnoteTest)
3: formatCons(stats[[i]], nam, tr, x$group.freq, prmsd, sep, formatArgs, 
       round, prtest, pdig = pdig, eps = eps)
2: print.summary.formula.reverse(list(stats = list(sbp = c(97.9191028465814, 
   94.9839938500064, 100.014783572809, 97.2881910017715, 107.288034416825, 
   105.746587052709, 111.864782483651, 112.689945667021, 116.229748640414, 
   115.604190135259, 119.743427097173, 119.276780441804, 123.380695706571, 
   122.111672516175, 128.138778071723, 126.592782661592, 133.726823015259, 
   132.141219449201, 140.847941698775, 136.762891898923, 145.175812916341, 
   141.635692905295, 120.013464038065, 118.994407752318, 12.1494617994813, 
   11.9252958974706)), type = 2, group.name = "treatment", group.label = "treatment", 

Answer 2

that works for me:

conTestWP <- 
function (group, x) 
{
    st <- wilcox.test( x[as.numeric(group)==1], x[as.numeric(group)==2], paired=TRUE)
    list(P=st$p.value, stat=st$statistic, df = "NA", 
    testname = "Wilcoxon-paired",
    statname = "V", namefun = "fstat", latexstat = "V", plotmathstat = "F[df]")}