lsmeans R package: how to obtain exact pvalue question when doing a pairwise comparision

Question

I'm using the lsmeans R package, and I'm doing a simple pairwise comparison between two groups. Is there a way to get the exact p.value in scientific notation rather than the estimate? I've tried setting options(scipen=100000000000), without luck. It just shows 0.000000e+00 when I do so.

lsmeans(fit,pairwise~type+age,at=list(age=c("40","50","60")),adjust="mvt")

$contrasts
 contrast     estimate          SE    df  t.ratio p.value
 1 - 2    -0.084728198 0.002606520 59346  -32.506  <.0001
 1 - 3    -0.269773787 0.002146597 59346 -125.675  <.0001
 1 - 4    -0.002048158 0.001896132 59346   -1.080  0.6824
 2 - 3    -0.185045589 0.002406580 59346  -76.892  <.0001
 2 - 4     0.082680040 0.002211575 59346   37.385  <.0001
 3 - 4     0.267725629 0.001064637 59346  251.471  <.0001

@kmm Thanks!! 2*pt(-abs(t),df=n-1) https://stats.stackexchange.com/questions/45153/manually-calculating-p-value-from-t-value-in-t-test — TM-F, Jul 13 '17 at 18:58
Above will produce UNADJUSTED P values, not the mvt-adjusted ones. — Russ Lenth, Jul 13 '17 at 22:14
Looking at this a second time, I don't believe the output shown could have come from that `lsmeans` call, because there are two factors involved. I think you did not have `age` on the rhs of the formula. — Russ Lenth, Jul 13 '17 at 23:02

score 2 · Answer 1 · answered Jul 13 '17 at 22:28

What do you mean by exact? The "mvt" method uses simulation-based routines from the mvtnorm package, and they are only good to a few digits. If you re-run the commands, you'll get slightly different results.

All that said, the way the adjusted P values are printed is a result of the print method for summary.ref.grid objects. Such objects inherit from data.frame, so you can just extract the column you need. A technicality: the lsmeans call you show actually produces a list of two ref.grids, so you need to pick out just the $contrasts member. Here's how it goes:

lsm <- lsmeans(fit, pairwise ~ type + age,
               at = list(age=c("40", "50", "60")), adjust = "mvt")
pval <- summary(lsm$contrasts)$p.value

Then format and print pval however you like. You can extract other results in a similar way.

But again, some thought was put into the way the P values are presented, and I think presenting them with exaggerated precision is misleading.

score 1 · Answer 2 · answered Jul 13 '17 at 01:17

First, is it really of practical importance to report a p-value to that level of precision?

Second, I believe this is probably due to your machine's precision ability. Before you run options(scipen=100000000000), were you getting a p-value of 2.2e-16? This is what it is for most machines I've worked on. This is because the computer cannot go to more precise estimate than 16 digits.

To see what your computer is capable of, check out the documentation ?.Machine and what the .Machine function returns.

For floats/doubles (i.e., numeric non-integers), you can enter .Machine$double.eps to see how precise your machine can get. For instance, see what happens when run this code:

> .Machine$double.eps
[1] 2.220446e-16

> 1 + 2.22e-15 == 1
[1] FALSE

> 1 + 2.22e-17 == 1
[1] TRUE

My machine couldn't tell 1.000000000000000222 was different from 1.

See relevant discussion here.

There's an earlier, more extensive discussion on our own site [here](https://stats.stackexchange.com/questions/78839/how-should-tiny-p-values-be-reported-and-why-does-r-put-a-minimum-on-2-22e-1/78840#78840) — Glen_b, Jul 13 '17 at 01:39

lsmeans R package: how to obtain exact pvalue question when doing a pairwise comparision

2 Answers2