Using lme to perform a nested ANCOVA in r

I want to perform a nested ANCOVA using the function lme in r. My data concerns the growth rate (sgr) of fish at different temperatures (temp). In the trial there were 3 replicate tanks each housing 8 fish at each level of temperature (15,18,21 and 24 degrees Celsius) so that a total of 12 tanks were included in the trial. I want to control for the effects that individual tanks may have had on growth and since the initial mass of the fish also effects growth rate I want to include it as a covariate.

These are the variables I plan to include in the model:

sgr= continuous dependent variable

temp= fixed main effect

mass= covariate

tank= random factor (tank coded 1-12)

Using lme I have coded my model as such:

m1=lme(sgr~temp+mass,random=~1|tank)

I want to know two things:

1. If the code I have used to include tank as a random factor has been done correctly?

2. Can a covariate (mass) be added this way using lme (just like a regular ANCOVA) when the model also includes a random effect?

Thank you in advance for any assistance provided?

• The code is for a random intercept, another option is to also include random slopes. Your code will estimate a regression for each fish, but the intercept in the regression may be different for each tank. If you include a random effect on e.g. temp that would reflect that the coefficient of temp may also change for each tank.
– user83346
Jul 22, 2016 at 6:36
• Thank you @fcop. Would changing the code to m1=lme(sgr~temp+mass,random=~1|tank/temp) add the random effect of tank on temp? Jul 22, 2016 at 7:01
• You should write m1=lme(sgr~temp+mass,random=~1+temp|tank)
– user83346
Jul 22, 2016 at 7:49
• Thank you @fcop. I am new to using lme and trying to get my head around coding the random factors. A lot of people seem to include the / in their code for the random factor. Is / used to add another random intercept for another random factor? For example if I had another random factor (call it "x") I wanted to include a random intercept for in the model in addition to tank, the code would look like m1=lme(sgr~tank+mass,random=~1|tank/x). Jul 23, 2016 at 2:12
• @fcop could you explain the difference between having (random=~1|tank) and (random=~1|tank/temp). Jul 23, 2016 at 4:07

If you use the following notation: $sgr_{ft}$ is the growth rate for fish $f$ in tank $t$, $temp_t$ temperature in tank $t$ and $mass_f$ the mass of fish $f$ then you estimate a regression of the form:

$sgr_{ft} = \beta_{0t} + \beta_1 temp_t + \beta_2 mass_f$ using the R-code lmeModel<-lme(sgr ~ mass + temp , random=~1|tank). Note that the intercept can be different for each tank ($\beta_{0t}$ has $t$ as subscript).

The $\beta_{0t}$ can be extracted from lmeModel using the extraction functions fixef and ranef. In fact, $\beta_{0t}$ has two components, i.e. $\beta_{0t}=\beta_0+b_{0t}$, where $\beta_0$ is a fixed effect and $b_{0t}$ a random effect. The function fixef gived you the maximum likelihood estimator of the fixed effects, the function ranef gived the best linear unbiased predictors (BLUP) of the random effects, you will see that ranef yields one value for each tank.

You may also try a model like $sgr_{ft} = \beta_{0t} + \beta_{1t} temp_t + \beta_2 mass_f$, where the coefficient $\beta_{1t}$ also depends on the tank, using the code lmeModel<-lme(sgr ~ mass + temp , random=~1+temp|tank)

Using the R-code lmeModel<-lme(sgr ~ mass + temp , random=~1|tank/temp) is similar to the first model, only that in the the first case you assume a different intercept for each tank while in the model with random=~1|tank/temp you assume a different intercept for each temperature within each tank. So it is something like $sgr_{ft} = \beta_{0t,T} + \beta_1 temp_t + \beta_2 mass_f$, where T in $\beta_{0t,T}$ stands for temperature. However, if you treat temperature as a factor (i.e. a categorical variable) then I think it does not make much difference.

EDIT 27-07-2016:

About the question in your comment; if temperature is categorical, having four values, then this categorical variable is replaces by three (the number of categories minus 1) dummy variables, $D_1, D_2, D_3$.

In that case lmeModel<-lme(sgr ~ mass + temp , random=~1|tank) estimates $sgr_{ft} = \beta_{0t} + \beta_{11} D_{1t} + \beta_{12} D_{2t} + \beta_{13} D_{3t}+ \beta_2 mass_f$ while lmeModel<-lme(sgr ~ mass + temp , random=~1+temp|tank) estimates $sgr_{ft} = \beta_{0t} + \beta_{11t} D_{1t} + \beta_{12t} D_{2t} + \beta_{13t} D_{3t}+ \beta_2 mass_f$

• do you have any advice on model selection between models with different random effects structures created using lme. Using 'anova(model1,model2)' I have compared to models 'model1<-lme(sgr~temp+mass,random=~1|tank' and 'model2<-sgr~temp+mass,random=~1+temp|tank)'. Model 1 has a BIC OF -6.5 and model 2 a BIC of 27.8 with a pvalue 0.81 for the test between models. Can this be interpreted as model 1 being a better fit due to lower BIC and accepting null hypothesis (p>0.05) of the test between models? Jul 26, 2016 at 23:02
• You can use a likelihood ratio test, however the likelihood ratio no longer has a chisquare distribution but a 50:50 mixture of chisquare with df 1 and 2. I am not sure whether anova takes this particularity into account.
– user83346
Jul 27, 2016 at 4:46
• Hi thanks @fcop I found a nice reference for the liklihood ratio test. One thing I can't understand is how random slopes and intercepts (=~1+temp|tank) can be included when temp is a categorical variable. Would it be true that including random slopes when treating temp as a categorical variable does not make sense? Jul 27, 2016 at 6:54
• Can you give me the reference ? For the categorical temperature; I will edit my answer, but do you know how categorical variables are treated, I mean, do you know how they are replaces by dummies ?
– user83346
Jul 27, 2016 at 7:42
• I'll read up on dummy coding but If you could try explain that would be great. The book is "Mixed Effects Models and Extensions in Ecology with R". In the chapter "Mixed Effects Modelling for Nested Data" they run through a model fitting protocol that includes the use of anova(). Jul 27, 2016 at 7:53