3

Here is the code sample

cq = c(15.5193589369923,15.6082097993496,15.3048276641115,15.887944383963,15.3813224284544,14.9723432922121,14.8742353464212,15.0448020475332,15.1584221729435,15.3692219904727,15.2369219681739,15.0804950645883,15.0836397511495,14.8821059462034,14.6827696388115,14.5701385743889,14.8506248103639,14.8475325690146,14.7377458445001,14.6258765734272,15.3585770134881,15.2994209401567,15.5178103826596,15.2411668198437,15.3413307248142,15.3645926457095,15.2241340874265,15.7516405898009,15.7683360263607,15.5852354340738,14.7451372367313,14.650625258402,14.7596201108925,15.0504977144055,15.0178091754821,15.100874342289,14.5156700740607,15.0530667717479,14.754595621435,15.5879633065185,15.3449828894141,15.3112460363113,15.232600495493,15.4378070492087,15.1621663266126,16.0120124580213,16.2104534293435,16.2765899877946,17.1446379330444,17.1717364140053,17.0155350105157,15.5218922723681,15.4543443324508,15.5282690363252,15.0202919978723,15.0410524376083,15.1169661551775,15.335220483258,15.3191814464955,15.0679651604846,14.7270263787123,14.70717761566,14.7907442084919,15.8468089268423,15.6714073529734,15.5478017041242,14.6949593095613,14.7537769900696,14.830942214569,15.0820225358985,15.3454125813989,15.304399073517,15.4159319040107,16.1250033895004,15.5359407225865,15.3251900155103,15.1571914994646,15.412721442436,15.913112918988,15.9852823695227,16.0912887332562,15.4897399161851,15.0710262650299,15.3517226832146,17.0001128578501,17.1040579605654,16.9578316599788,15.8842918497549,15.7016383123704,15.8513519332371,16.9420990886101,17.0793832045434,16.9288868492911,14.9127628979216,14.7689529893246,15.0534122173222,15.3185448628303,15.5507864243439,15.3737185073511,15.4350799532271,15.2414612478027,15.361320770232,15.7401140808761,15.8582795450189,15.7634036480016,15.5797995263497,15.9126261329496,15.9256641722586,15.1308493265056,15.2450158090279,15.0699176510971,15.0368959001792,14.8828877909216,14.8852035927172,15.8253506435753,15.8938440960183,15.888311876759,15.4872886586516,15.5492199156675,15.7313291529313,16.5365758222542,16.8386981731158,16.7239280992675,15.9356391540897,15.7383049532238,15.9409000309973,15.2005952554035,15.0390142751348,15.154888655127,14.6373767323354,14.3087397097081,14.3970067065903,14.6453627024929,14.8109205614192,14.6266778290643,15.5170574352528,15.359943766027,15.5869322081508,15.5246550838727,15.4670382654415,15.4211907882731,16.9534561402918,17.4848334482537,17.3182067272327,15.7804318020053,15.86794322314,15.6532944587946,16.543432367992,16.6848617423114,16.8344939905775,15.5212254647114,15.8348559815603,15.6592827767612,15.3027400892518,15.5498124465958,15.8362202772445,14.8415823671167,15.7307379811374,14.8529575353737,16.6466266958983,16.1687733596343,16.0342973266029,14.5976161739123,14.776507726931,14.6780484406283,15.3927619991024,15.3106866267163,15.2920260038624,15.9666798099925,16.2595244266754,16.1035265916681,16.018233002759,15.8460056716414,16.0722176294152,16.2763177549617,16.364250121284,16.2995041975045,16.3975912697976,16.182759197402,16.1164022801451,16.5026752161837,16.2401540005223,16.3715573563274,18.4119769797938,18.1208386122385,18.0068316479116,17.1455993749728,17.0558275544137,16.9150038143768)
sample = c("CD4 LM","CD4 LM","CD4 JC","CD4 JC","CD4 JC","CD4 BM","CD4 BM","CD4 BM","CD4 MC","CD4 MC","CD4 MC","CD4 TM","CD4 TM","CD4 TM","CD4 MM","CD4 MM","CD4 MM","CD4 SRits","CD4 SRits","CD4 SRits","CD4 GV","CD4 GV","CD4 GV","CD4 WW","CD4 WW","CD4 WW","CD4 CH","CD4 CH","CD4 FJ","CD4 FJ","CD4 KS","CD4 KS","CD4 KS","CD4 NG","CD4 NG","CD4 NG","CD4 CG","CD4 CG","CD4 CG","CD4 CSR","CD4 CSR","CD4 CSR","CD4 JM","CD4 JM","CD4 JM","CD4 DF","CD4 DF","CD4 DF","CD4 AM","CD4 AM","CD4 AM","CD4 BP","CD4 BP","CD4 BP","CD4 ER","CD4 ER","CD4 ER","CD4 SRusse","CD4 SRusse","CD4 SRusse","CD4 DS","CD4 DS","CD4 DS","CD4 KJ","CD4 KJ","CD4 KJ","CD4 GD","CD4 GD","CD4 GD","CD4 KG","CD4 KG","CD4 KG","CD4 KR","CD4 KR","CD4 KR","CD4 FN","CD4 FN","CD4 FN","CD4 RM","CD4 RM","CD4 RM","CD4 LA","CD4 LA","CD4 LA","CD4 EC","CD4 EC","CD4 EC","CD4 KW","CD4 KW","CD4 KW","CD4 HB","CD4 HB","CD4 HB","CD8 LM","CD8 LM","CD8 LM","CD8 JC","CD8 JC","CD8 JC","CD8 BM","CD8 BM","CD8 BM","CD8 MC","CD8 MC","CD8 MC","CD8 TM","CD8 TM","CD8 TM","CD8 MM","CD8 MM","CD8 MM","CD8 SRits","CD8 SRits","CD8 SRits","CD8 GV","CD8 GV","CD8 GV","CD8 WW","CD8 WW","CD8 WW","CD8 CH","CD8 CH","CD8 CH","CD8 FJ","CD8 FJ","CD8 FJ","CD8 KS","CD8 KS","CD8 KS","CD8 NG","CD8 NG","CD8 NG","CD8 CG","CD8 CG","CD8 CG","CD8 CSR","CD8 CSR","CD8 CSR","CD8 JM","CD8 JM","CD8 JM","CD8 DF","CD8 DF","CD8 DF","CD8 AM","CD8 AM","CD8 AM","CD8 BP","CD8 BP","CD8 BP","CD8 ER","CD8 ER","CD8 ER","CD8 SRusse","CD8 SRusse","CD8 SRusse","CD8 DS","CD8 DS","CD8 DS","CD8 KJ","CD8 KJ","CD8 KJ","CD8 GD","CD8 GD","CD8 GD","CD8 KG","CD8 KG","CD8 KG","CD8 KR","CD8 KR","CD8 KR","CD8 FN","CD8 FN","CD8 FN","CD8 RM","CD8 RM","CD8 RM","CD8 LA","CD8 LA","CD8 LA","CD8 EC","CD8 EC","CD8 EC","CD8 KW","CD8 KW","CD8 KW","CD8 HB","CD8 HB","CD8 HB")
df = data.frame(cq, sample)

df.res <- lm(cq~sample, data = df)
require(lsmeans)
t<- pairs(lsmeans(df.res, "sample"))
system.time(tc <- confint(t, level=0.95))

The time taken by the confint function is 

   user  system elapsed 
  10.58    0.00   10.60

I have tried using confint.default but I get an error

tc <- confint.default(t, level=0.95)

Error: $ operator not defined for this S4 class

cheedep
  • 937
  • 6
  • 19
  • I'd suspect that some kind of bootstrapping approach is used. – Roland Jul 29 '16 at 21:07
  • if I use glm instead of lm the confint executes in 0.27 seconds but the confidence intervals are not exactly what is expected. – cheedep Jul 29 '16 at 21:22
  • the other question : do you need this to be faster? Is it a bottleneck in your workflow? – Ben Bolker Jul 29 '16 at 21:32
  • Yes, I have multiple models to be computed and at 10 seconds a piece it sums up. I can use the TukeyHSD function with aov to get the confidence intervals and it takes only a second. However I do two way anova with type III sum of squares for unbalanced designs using car package and cannot use the TukeyHSD function anymore. I am looking at your answer below. Looks like the "mvt" option does not calculate the CI's – cheedep Jul 29 '16 at 21:53

2 Answers2

4

It's a bit buried in the documentation, but what's slowing you down is the multiple-comparisons correction computations. There's wide variation in the elapsed time for the available methods. See the Confidence-limit and P-value adjustments section of ?summary.ref.grid for details, and decide which method meets your criteria of being both fast enough and reliable enough for your use case ...

adj <- c("tukey","scheffe","sidak","bonferroni","dunnettx","mvt") 
sapply(adj,function(a) system.time(confint(t,adjust=a))["elapsed"])

##      tukey.elapsed    scheffe.elapsed      sidak.elapsed bonferroni.elapsed 
         9.256              0.195              0.168              0.173 
##   dunnettx.elapsed        mvt.elapsed 
        14.370              1.821
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Thanks Ben for the answer. I wanted to use tukey correction. I could step through into the code for the confint and it goes to summary.ref.grid. Thanks to your answer for pointing me in this direction. I found that most of the time is wasted in the qtukey function and that can be optimized. I will look if there is a way to override an internal library method. If not, I will have to copy that function and make my own. – cheedep Jul 29 '16 at 23:05
  • it will probably be easier to copy the function. I'm surprised that you can optimize qtukey -- already coded in C ... ? – Ben Bolker Jul 29 '16 at 23:12
  • No not the qtukey function :). There is an internal function that is calling qtukey with qtukey(0.95, 64, rep(125, 2016))/sqrt(2). I think I can optimize it by calling qtukey for only unique values of degrees of freedom and fill the array. – cheedep Jul 29 '16 at 23:15
  • this function in lsmeans getAnywhere('.adj.critval') – cheedep Jul 29 '16 at 23:20
  • oh, OK, that makes sense – Ben Bolker Jul 30 '16 at 00:17
  • if you want to hack that function, try dumping `lsmeans:::.adj.critval` to a file (e.g. via `dump()`); editing it, reading it back in; then use `assignInNamespace` to replace ... – Ben Bolker Jul 30 '16 at 00:18
  • 1
    Hmmmm. I know that `qtukey` is among the slowest built-in functions in R. But it surprises the heck out of me that the `"mvt"` method, which uses a simulation algorithm in the **mvtnorm** package, is *faster*. If the speed for `"mvt"` is acceptable, then use it! (Because it is "exact" for the correlation structure that exists among these results.) If it's still too slow, then you can do `smry <- summary(t, infer = FALSE)` which produces an object that inherits from `data.frame`, having the estimates and SEs; then manually compute the CIs you need using manually-provided multipliers of the SEs. – Russ Lenth Jul 31 '16 at 01:57
  • 3
    PS -- Do you *really* want to simultaneously test 2016 pairwise comparisons? Really? Are you sure? Are you really sure? – Russ Lenth Jul 31 '16 at 02:06
  • PPS - I think `"mvt"` has a limit of fewer comparisons than that, and that's why its timing is so much better. – Russ Lenth Jul 31 '16 at 02:06
  • Thank you rvl. The limit is probably why the mvt was not returning the confidence intervals. – cheedep Aug 03 '16 at 16:50
0

Here is a fast way:

confint(t,method="Wald")

NJArap
  • 31
  • 3
  • In this question, `t` is of class "emmGrid", for which there is a `confint` method. Maybe some other `confint` methods have a `method` option, but `confint.emmGrid` does not. So this answer will not speed up anything, as the argument will be ignored. – Russ Lenth Jun 14 '22 at 22:30