[R] [FORGED] Error occurring in "emmeans" package for the two data sets I used. Please help.

[Rolf Turner]

On 07/01/18 02:19, Akhilesh Singh wrote:

> I am a Professor of Statistics at Indira Gandhi Krishi Vishwavidyalaya,
> Raipur, India. While teaching in class about analysis of variance using R,
> I was doing a one-way analysis for the two data-sets given below in the
> R-class. I got a typical error in "emmeans" package, please help:
> 
> Data-set-1:
> --------------
> Medley and Clements (1998) investigated the impact of zinc contamination
> (and other heavy metals) on the diversity of diatom species in the USA
> Rocky Mountains. The diversity of diatoms (number of species) and degree of
> zinc contamination (categorized as either of high, medium, low or natural
> background level) were recorded from between four and six sampling stations
> within each of six streams known to be polluted, as given below:
> 
> stream=c("Eagle", "Eagle", "Eagle", "Eagle", "Blue", "Blue",
>           "Blue", "Blue", "Blue", "Blue", "Blue", "Snake", "Snake",
>           "Snake", "Snake", "Snake", "Arkan", "Arkan", "Arkan",
>           "Arkan", "Arkan", "Arkan", "Arkan", "Chalk", "Chalk",
>           "Chalk", "Chalk", "Chalk", "Splat", "Splat", "Splat",
>           "Splat", "Splat", "Splat")
> 
> zinc=c("BACK", "HIGH", "HIGH", "MED", "BACK", "HIGH", "BACK", "BACK",
>         "HIGH", "MED", "MED", "BACK", "MED", "HIGH", "HIGH", "HIGH",
>         "LOW", "LOW", "LOW", "LOW", "MED", "MED", "LOW", "LOW",
>         "HIGH", "HIGH", "MED", "LOW", "BACK", "BACK", "MED", "LOW",
>         "MED", "BACK")
> 
> diversity=c(2.27, 1.25, 1.15, 1.62, 1.7, 0.63, 2.05, 1.98, 1.04,
>              2.19, 2.1, 2.2, 2.06, 1.9, 1.88, 0.85, 1.4, 2.18, 1.83,
>              1.88, 2.02, 1.94, 2.1, 2.38, 1.43, 1.37, 1.75, 2.83,
>              1.53, 0.76, 0.8, 1.66, 0.98, 1.89)
> 
> medley.clementis=data.frame(stream,zinc,diversity)
> 
> I did the one-way anova:
> -------------------------------
> 
> medley.clementis.aov=with(medley.clementis, aov(diversity ~ zinc))
> 
> anova(medley.clementis)
> 
> Then, I tried to do post hoc analysis using "emmeans" package following
> command:
> -----------------------------------------------------------------------------------------------
> 
> emmeans::emmeans(medley.clementis.aov, "zinc")
> 
> 
> This gives following error:
> ----------------------------------
> Error in recover_data.call(fcall, delete.response(terms(object)),
> object$na.action,  :
>    object 'possibly.random' not found
> Error in ref_grid(object, ...) :
>    Perhaps a 'data' or 'params' argument is needed
> 
> 
> 
> Data-set-2:
> ---------------
> Keough and Raimondi (1995) examined the effects of four biofilm types (SL:
> sterile unfilmed substrate, NL: netted laboratory biofilms, UL: unnetted
> laboratory biofilms and F: netted field biofilms) on the recruitment of
> serpulid larvae. Substrates treated with one of the four biofilm types were
> left in shallow marine waters for one week after which the number of newly
> recruited serpulid worms were counted, as given below:
> 
> biofilm=c("SL", "SL", "SL", "SL", "SL", "SL", "SL", "UL", "UL", "UL",
>            "UL", "UL", "UL", "UL", "NL", "NL", "NL", "NL", "NL", "NL",
>            "NL", "F", "F", "F", "F", "F", "F", "F")
> 
> serpulid=c(61, 113, 123, 75, 75, 83, 95, 143, 81, 101, 155, 156, 193,
>             163, 203, 159, 139, 161, 179, 97, 157, 128.5, 204.5,
>             108.5, 116.5, 140.5, 160.5, 87.5)
> 
> keough.raimondi=data.frame(biofilm,serpulid)
> 
> Applied log-transformation:
> -------------------------------------------
> keough.raimondi.ln=transform(keough.raimondi, serpulid.ln=log(serpulid))
> 
> I did the one-way anova, with contrasts defined below:
> ------------------------------------------------------------------------
> contrasts(keough.raimondi.ln$biofilm) <- cbind(c(0, 1, 0, -1),
>                              c(2, -1, 0, -1), c(-1, -1, 3, -1))
> keough.raimondi.ln$biofilm
> 
> keough.contr.list <- list(biofilm = list('NL vs UL' = 1,
>                'F vs (NL & UL)' = 2, 'SL vs (F & NL & UL)' = 3))
> keough.contr.list
> 
> One-way anova:
> ----------------------
> keough.raimondi.ln.aov=with(keough.raimondi.ln, aov(serpulid.ln ~ biofilm))
> 
> summary(keough.raimondi.ln.aov,split=keough.contr.list)
> 
> 
> Then, I tried to do post hoc analysis using "emmeans" package following
> command:
> -----------------------------------------------------------------------------------------------
> 
> emmeans(keough.raimondi.ln.aov, ~ biofilm)
> 
> 
> This gives following error:
> ----------------------------------
> Error in recover_data.call(fcall, delete.response(terms(object)),
> object$na.action,  :
>    object 'possibly.random' not found
> Error in ref_grid(object, ...) :
>    Perhaps a 'data' or 'params' argument is needed
> 
> 
> Help Needed:
> ------------------
> On many other data sets and data frame I successfully used "emmeans"
> package using the help available in R.
> 
> But, for the above two data-sets, I consistently got the same error as
> described above.
> 
> I do not know what is amiss. Where I am missing or whatever is wrong, I
> request the entire R-team to help me to solve above problem.

Well, you don't need the *entire* R-team!!! It probably (in some sense) 
includes millions of people. :-)

> Thanking in advance.

Thanks for your thorough and well set out description of the problem.
Your reproducible examples were flawless.

I am not *completely* certain, but this looks to me like a bug in emmeans.

I have therefore taken the liberty of cc-ing this reply to Russell Lenth 
(the maintainer of emmeans) to get his take on the issue.

cheers,

Rolf Turner

-- 
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276