1

So I ran a Two-way Ordinal Regression with CLM with the specific columns: gra, location, species and it worked: 

gra.f <- clm(gra ~ location + species + location:species, data = gra)

now I'm tring to run the pairwise using the lsmeans function:

pwise <- lsmeans(gra.f, 
                   pairwise ~ location + species, 
                   adjust="tukey")

And I am getting: Error in eval(predvars, data, env) : object 'location' not found

I'm not sure I know what the problem is

PS: If I remove location and try to just do species, the error still read 'location' not found

kjgrad
  • 21
  • 2
  • Comments cannot be attached to deleted answers but I thought that using the `maintainer` function would be more economical than the earlier (misspelled) suggestion: `maintainer("lsmeans")` – IRTFM Dec 21 '18 at 22:29

1 Answers1

1

lsmeans has to reconstruct the dataset (including all predictor values) in order to determine the reference grid. The error message you report results from it’s being unable to do that. The only thing I can guess is that the data frame used to fit the model is no longer in the workspace or on the search path.

If you restore the dataset or make it visible, lsmeans should work again. Alternatively, add , data = gra to the lsmeans() call.

Update 1

It turns out that the problem occurs in the code to handle the fact that this model is rank-deficient. Ironically, I think it relates to this entry in the NEWS file for the ordinal package: 2014-11-12: - Reimplementation of formula, model.frame and design matrix processing motivated by a bug in model.matrix.clm and predict.clm reported by Russell Lenth 2014-11-07 when implementing lsmeans support for clm::ordinal.

I had put-in some code to work around that, and apparently,

  1. my workaround only worked until the aforementioned update
  2. not many users over the past 4 years are encountering rank-deficient ordinal models; or at least I didn't hear from them if they did.

I am now trying to remember which lines of code are my workaround, and which are still needed... When I have this resolved, I'll push it to the github repository for emmeans, and in a month or so, emmeans will be updated on CRAN (it'll be a version greater than 1.3.1). Updating emmeans will make lsmeans work too, as now it is just a front-end for emmeans.

Update 2

It seems to be fixed now:

> lsmeans(gra.f, 
+         pairwise ~ location + species, 
+         adjust="tukey")

$`lsmeans`
 location species  lsmean    SE  df asymp.LCL asymp.UCL
 B10      HD       nonEst    NA  NA        NA        NA
 B30      HD       nonEst    NA  NA        NA        NA
 B50      HD       nonEst    NA  NA        NA        NA
 B70      HD      -0.2802 0.191 Inf  -0.65519    0.0949
 Black Pt HD      -0.1298 0.325 Inf  -0.76730    0.5077
 Bolongo  HD       nonEst    NA  NA        NA        NA
      ... Several rows of output omitted ...
 SMA      TT       nonEst    NA  NA        NA        NA

Confidence level used: 0.95 

$contrasts
 contrast                  estimate    SE  df z.ratio p.value
      ... MANY rows of output omitted ...
 B70,HD - Magens,HD          nonEst    NA  NA      NA      NA
 B70,HD - SMA,HD             nonEst    NA  NA      NA      NA
 B70,HD - B10,HS           -0.71283 0.225 Inf  -3.162  0.4869
 B70,HD - B30,HS            0.07174 0.227 Inf   0.316  1.0000
 B70,HD - B50,HS           -0.74197 0.228 Inf  -3.253  0.4093
 B70,HD - B70,HS           -0.72863 0.226 Inf  -3.229  0.4291
 B70,HD - Black Pt,HS      -0.10537 0.220 Inf  -0.478  1.0000
 B70,HD - Bolongo,HS         nonEst    NA  NA      NA      NA
 B70,HD - Fortuna,HS       -0.64301 0.240 Inf  -2.684  0.8656
 B70,HD - Lindberg,HS      -0.10132 0.220 Inf  -0.460  1.0000
 B70,HD - Magens,HS         0.02351 0.226 Inf   0.104  1.0000
 B70,HD - SMA,HS           -0.05100 0.219 Inf  -0.232  1.0000
 B70,HD - B10,HW             nonEst    NA  NA      NA      NA
 B70,HD - B30,HW             nonEst    NA  NA      NA      NA
 B70,HD - B50,HW             nonEst    NA  NA      NA      NA
 B70,HD - B70,HW             nonEst    NA  NA      NA      NA
 B70,HD - Black Pt,HW        nonEst    NA  NA      NA      NA
 B70,HD - Bolongo,HW         nonEst    NA  NA      NA      NA
 [ reached getOption("max.print") -- omitted 1059 rows ]

P value adjustment: tukey method for comparing a family of 50 estimates

However, your problems are just beginning. You have a large number of non-estimable combinations due to missing data. Good luck...

I will push the updated version to https://github.com/rvlenth/emmeans.

Russ Lenth
  • 5,922
  • 2
  • 13
  • 21