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.


 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
  • $\begingroup$ pt(-32.506, 59346) = 4.61944e-230 $\endgroup$
    – kmm
    Jul 13 '17 at 1:13
  • $\begingroup$ @kmm lol how did you do that? $\endgroup$
    – Mark White
    Jul 13 '17 at 1:17
  • $\begingroup$ @kmm Thanks!! 2*pt(-abs(t),df=n-1) stats.stackexchange.com/questions/45153/… $\endgroup$
    – TM-F
    Jul 13 '17 at 18:58
  • $\begingroup$ Above will produce UNADJUSTED P values, not the mvt-adjusted ones. $\endgroup$
    – Russ Lenth
    Jul 13 '17 at 22:14
  • $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Jul 13 '17 at 23:02

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.


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 + 2.22e-17 == 1
[1] TRUE

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

See relevant discussion here.

  • 1
    $\begingroup$ There's an earlier, more extensive discussion on our own site here $\endgroup$
    – Glen_b
    Jul 13 '17 at 1:39
  • $\begingroup$ @Glen_b nice and comprehensive! $\endgroup$
    – Mark White
    Jul 13 '17 at 3:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.