1
$\begingroup$

I'm trying to plot confidence intervals for linear mixed effects models trained with lme4 and lmerTest in R. I am using this data file, which I've shared via Google Drive.

Here is my trained model. The data consists of 8 subjects (SID) and 580 items per subject (UID).

> df <- readRDS(file="data.Rda")
> summary(df)
      UID            SID             Y                  X          
 U1     :   8   H1     : 580   Min.   :-1.75000   Min.   :0.00000  
 U10    :   8   H2     : 580   1st Qu.:-0.13330   1st Qu.:0.00000  
 U100   :   8   H3     : 580   Median :-0.02470   Median :0.08054  
 U101   :   8   H4     : 580   Mean   :-0.08563   Mean   :0.14070  
 U102   :   8   H5     : 580   3rd Qu.: 0.00000   3rd Qu.:0.21053  
 U103   :   8   H6     : 580   Max.   : 0.50000   Max.   :1.20000  
 (Other):4592   (Other):1160 
> my.model <- lmer(Y ~ X + (1|UID) + (1|SID), data=df)
> summary(my.model)
Linear mixed model fit by REML 
t-tests use  Satterthwaite approximations to degrees of freedom ['merModLmerTest']
Formula: Y ~ X + (1 | UID) + (1 | SID)
   Data: df

REML criterion at convergence: -10980.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-16.4681  -0.2699   0.0042   0.3194   6.7467 

Random effects:
 Groups   Name        Variance  Std.Dev.
 UID      (Intercept) 4.573e-03 0.067624
 SID      (Intercept) 2.185e-06 0.001478
 Residual             4.109e-03 0.064099
Number of obs: 4640, groups:  UID, 580; SID, 8

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)  8.501e-04  3.255e-03  3.056e+02   0.261    0.794    
X           -6.146e-01  8.861e-03  1.644e+03 -69.362   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
  (Intr)
X -0.383

I've tried generating confidence intervals for X using a number of approaches, to no success. With lsmeans, I don't get any output.

> lsmeans(my.model)
Least Squares Means table:
     Estimate Standard Error DF t-value Lower CI Upper CI p-value

I can generate confidence intervals using confint with Wald statistics, but using the default method runs indefinitely.

> confint(my.model, method="Wald")
                   2.5 %       97.5 %
(Intercept) -0.005530357  0.007230525
X           -0.631989857 -0.597255116
> confint(my.model) # This runs indefinitely
Computing profile confidence intervals ...

>
> effect(c("X"), my.model) # This also runs indefinitely

I have no problem getting confidence intervals on the example datasets.

> m1 <- lmer(Informed.liking ~ Gender*Information +(1|Consumer), data=ham)
> lsmeans(m1)
Least Squares Means table:
                        Gender Information Estimate Standard Error  DF t-value Lower CI Upper CI p-value    
Gender  1                    1          NA    5.854          0.183  79    32.0     5.49     6.22  <2e-16 ***
Gender  2                    2          NA    5.609          0.185  79    30.3     5.24     5.98  <2e-16 ***
Information  1              NA           1    5.632          0.155 154    36.4     5.33     5.94  <2e-16 ***
Information  2              NA           2    5.831          0.155 154    37.7     5.53     6.14  <2e-16 ***
Gender:Information  1 1      1           1    5.707          0.218 154    26.2     5.28     6.14  <2e-16 ***
Gender:Information  2 1      2           1    5.556          0.220 154    25.2     5.12     5.99  <2e-16 ***
Gender:Information  1 2      1           2    6.000          0.218 154    27.6     5.57     6.43  <2e-16 ***
Gender:Information  2 2      2           2    5.662          0.220 154    25.7     5.23     6.10  <2e-16 ***

Here is my session info:

> sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] effects_3.0-3   pbkrtest_0.4-2  lmerTest_2.0-20 lme4_1.1-7      Rcpp_0.11.6     Matrix_1.2-0   

loaded via a namespace (and not attached):
 [1] acepack_1.3-3.3     bitops_1.0-6        caTools_1.17.1      cluster_2.0.1       colorspace_1.2-6    digest_0.6.8        foreign_0.8-63     
 [8] Formula_1.2-1       gdata_2.16.1        ggplot2_1.0.1       gplots_2.17.0       grid_3.1.1          gridExtra_0.9.1     gtable_0.1.2       
[15] gtools_3.4.2        Hmisc_3.16-0        KernSmooth_2.23-14  lattice_0.20-29     latticeExtra_0.6-26 magrittr_1.5        MASS_7.3-40        
[22] minqa_1.2.4         munsell_0.4.2       nlme_3.1-120        nloptr_1.0.4        nnet_7.3-9          numDeriv_2014.2-1   parallel_3.1.1     
[29] plyr_1.8.2          proto_0.3-10        RColorBrewer_1.1-2  reshape2_1.4.1      rpart_4.1-9         scales_0.2.4        splines_3.1.1      
[36] stringi_0.4-1       stringr_1.0.0       survival_2.38-1     tools_3.1.1

Can someone explain why lsmeans won't work on my model? Thanks in advance.

Update: Thanks to @aosmith, I now understand that lsmeans only displays confidence intervals on factors. So here's a related question.

I also tried computing confidence intervals on the fixed effects using the effects package. However, this seems to run indefinitely.

lvls <- c(1:10) / 10
Effect(c("X"), my.model, xlevels=lvls)

I don't think it's related to the fact that X is numeric. I tried the following example, and I got

str(mtcars)
m <- lmer(mpg ~ 1 + wt + hp + (1 + wt |gear), data=mtcars)
str(m)
Effect("wt", m)
> str(mtcars)
'data.frame':   32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
     $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
     $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
     $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
     $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
     $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...
> m <- lmer(mpg ~ 1 + wt + hp + (1 + wt |gear), data=mtcars)
> Effect("wt", m)

 wt effect
wt
       2        3        4        5 
24.59905 20.20403 15.80902 11.41400

Any suggestions why it won't converge in my dataset?

$\endgroup$
  • 3
    $\begingroup$ The description of lsmeans in lmerTest is: Calculates Least Squares Means and Confidence Intervals for the factors of a fixed part of mixed effects model of lmer object. You don't have any factors in the fixed effects part of your model, which seems a likely reason that lsmeans isn't returning anything. $\endgroup$ – aosmith May 27 '15 at 16:09
  • $\begingroup$ Good catch. That's probably the reason. $\endgroup$ – Nick Ruiz May 27 '15 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.