I initially posted this question to one of the StackExchange sites, and they suggested that I repost my problem here. After using ezANOVA as my primary way of specifying mixed ANOVAs, I've hit a stumbling block when it come to adding a covariate to the model. I am using an ANCOVA in order to determine if there is a developmental trajectory in my data; namely, I need to be able to see the F-statistic and p-values for interactions with the covariate (see p.466 onwards here [http://www.psyc.bbk.ac.uk/research/DNL/personalpages/annaz_etal_2009.pdf] if you want an example). Using ezANOVA, I can include covariates but the output does not show the F-statistic and p-values for interactions with the covariate - the main effect of the covariate is also not tested using this method. My ezANOVA model is as follows: aov.model<-ezANOVA( data=textureView.child.outliersRemoved , dv=.(x) , wid=.(ID) , within=.(Texture,View) , between=.(TNOGroup) , between_covariates=.(Age) , type=3 , return_aov=TRUE ) Another option is to use lm or Anova, but I don't know how to specify the error terms properly for either and I'm limited because I want to use Type-III sums of squares (drop1 doesn't work in the cases where I've tried to use the aov wrapper for lm; it fails while reporting 'Error in formula.default(object, env = baseenv()) : invalid formula'). Finally, I've heard about using the nlme package to specify my ANCOVA as a mixed model instead, but I don't know where to begin here (despite spending a while reading about it). To give a summary, I'm trying to do a 2 (between; TNOGroup) x 2 (within, Texture) x2 (within, View) mixed ANCOVA, with age as a covariate. I want to use Type-III sums of squares, and see the F-statistic and p-values for interactions with the covariate, as well as for the main effect of the covariate. Any advice on the best way to do this would be much appreciated. Thanks, Danielle