Search code examples
remmeans

Error with emmeans function in R . How to fix it?

I have a file like this :

enter image description here

I am using this data set to predict a linear mixed model and the I want to use the function emmeans in order to calculate the estimated means for my conditions. The code that I am using is here : 

newtab5 <- read.csv(file="sliding_LMM_test.csv", header=T, sep=",")
head(newtab5)

library(lmerTest)

model <- lmer(formula = data1 ~ flabel + (1 | whichFragments), data = newtab5, control=lmerControl(check.nlev.gtr.1 = "ignore"))
sfit <-summary(model)
sfit$coefficients[,2][2]
sfit$coefficients[,5][2]

library(emmeans)

means <- emmeans(model, "flabel")
test<-summary(means)

And I have this error: 

>  Error in t(ZZ) %*% EE :Cholmod error 'X and/or Y have wrong dimensions' at file ../MatrixOps/cholmod_sdmult.c, line 90

It looks like the emmeans function does not work when the which fragments column has the same values.

Do you know any way to skip this or to make emmeans to work for that case ?

Thank you in advance

Solution

  • It is an error that comes up in the pbkrtest package, in the computation of the adjusted covariance matrix needed for the Kenward-Roger degrees of freedom (which is the default degrees-of-freedom method for lmerMod objects). Almost certainly this is related to the undisclosed warning message encountered when fitting the model:

    > model <- lmer(formula = data1 ~ flabel + (1 | whichFragments), 
+    data = data, control=lmerControl(check.nlev.gtr.1 = "ignore")) 
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined

    Surprisingly (to me), you can still get results if you use some other d.f. method, e.g. Satterthwaite:

    > emmeans(model, "flabel", mode = "satterth")
 flabel emmean    SE   df lower.CL upper.CL
 con1    0.599 0.147 0.26    -3879     3880
 con2    0.107 0.147 0.26    -3880     3880

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95

    Having 0.26 degrees of freedom isn't much: look at the CIs.

    Just as a general observation, you are always asking for trouble when you try to estimate a variance with only one observation. That is what this model tries to do, in terms of the between-whichFragments variance. I sure hope you have more data than this in your real research.

    By the way, if all you want is the means, you can just compute the means:

    > with(data, tapply(data1, flabel, mean))

  con1   con2 
0.5991 0.1072