0

I have the following model

require(effects)

fit<-lme(x ~ y, data, random= ~1|item)
plot(allEffects(fit)

fit2<-lme(x ~ y, data2, random = ~1|item)
plot(allEffects(fit2)

How can I plot fit and fit2 overlaying? I have tried the par(new=T), but it does not work. The graphs plot fine individually.

zx8754
  • 52,746
  • 12
  • 114
  • 209
Katie Tetzloff
  • 55
  • 1
  • 1
  • 6
  • 1
    Please read the info about [how to ask a good question](http://stackoverflow.com/help/how-to-ask) and how to give a [reproducible example](http://stackoverflow.com/questions/5963269). This will make it much easier for others to help you. – zx8754 May 19 '16 at 20:35
  • provide some data, and what do you mean by overlay? At the very least show us an example of what you want your plot to look like – Cyrus Mohammadian May 19 '16 at 20:42

1 Answers1

2

I'm not sure there's a very nice way to do this. I usually extract the information from the effects structure and plot it with ggplot (lattice would be possible too).

Here's an example:

library(effects)
library(nlme)
library(plyr)  ## utilities

Fit a model to the first and second half of one of the standard example data sets:

fm1 <- lme(distance ~ age, random = ~1|Subject,
          data = Orthodont[1:54,])
fm2 <- update(fm1, data = Orthodont[55:108,])
a1 <- allEffects(fm1)
a2 <- allEffects(fm2)

Extract the information from the efflist object. This is the part that isn't completely general ... the hard part is getting out the predictor variable.

as.data.frame.efflist <- function(x) {
   ldply(x,
         function(z) {
           r <- with(z,data.frame(fit,
                                  var=variables[[1]]$levels,
                                  lower,upper))
           return(plyr::rename(r,setNames(z$variables[[1]]$name,"var")))
         })
}

For convenience, use ldply to put the results of both models together:

comb <- ldply(list(fm1=a1,fm2=a2),as.data.frame,.id="model")

Now plot:

library(ggplot2); theme_set(theme_bw())
ggplot(comb,aes(age,fit,
               ymin=lower,ymax=upper,
               colour=model,fill=model))+
     geom_line()+
     geom_ribbon(alpha=0.2,colour=NA)+
     geom_rug(sides="b")

The rug plot component is a little silly here.

enter image description here

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453