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

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  • pt(-32.506, 59346) = 4.61944e-230 – kmm Jul 13 '17 at 1:13 • @kmm lol how did you do that? Jul 13 '17 at 1:17 • @kmm Thanks!! 2*pt(-abs(t),df=n-1) stats.stackexchange.com/questions/45153/… – TM-F Jul 13 '17 at 18:58 • Above will produce UNADJUSTED P values, not the mvt-adjusted ones. 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. Jul 13 '17 at 23:02 ## 2 Answers 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] 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 Jul 13 '17 at 1:39
• @Glen_b nice and comprehensive! Jul 13 '17 at 3:32