We first import the tab-delimited file, rat_pup.dat, assumed to be located in the C:\temp directory, into a data frame object in R named ratpup. The h = T argument in the read.table() function indicates that the first row (the header) in the rat_pup.dat file contains variable names. After reading the data, we “attach” the ratpup data frame to R’s working memory so that the columns (i.e., variables) in the data

frame can be easily accessed as separate objects. Note that we show the “>” prompt for each command as it would appear in R, but this prompt is not typed as part of the commands.
> ratpup <- read.table("c:\\temp\\rat_pup.dat", h = T)
> attach(ratpup)
To facilitate comparisons with the analyses performed using the other software procedures, we recode the variable SEX into SEX1, which is an indicator variable for females, so that males will be the reference group:
> ratpup$sex1[sex == "Female"] <- 1
> ratpup$sex1[sex == "Male"] <- 0Step 1: Fit a model with a “loaded” mean structure (Model 3.1). We first load the nlme package,* so that the lme() function will be available for model fitting:
> library(nlme)
We fit the initial LMM, Model 3.1, to the Rat Pup data using the lme() function:
> # Model 3.1.
> model3.1.fit <- lme(weight ~ treatment + sex1 + litsize +treatment*sex1, random = ~1 | litter, ratpup, method = "REML")

We next explain each part of the syntax used for the lme() function:

• model3.1.fit is the name of the object that will contain the results of the fitted model.
• The first argument of the function, weight ~ treatment + sex1 + litsize+ treatment*sex1, defines a model formula. The response variable, WEIGHT, and the terms that have associated fixed effects in the model (TREATMENT,SEX1, LITSIZE, and the TREATMENT×SEX1 interaction), are listed.
• The factor() function is not necessary for the categorical variable TREATMENT, because the original treatment variable has string values High, Low, and Control, and will therefore be considered as a factor automatically. We also do not need to declare SEX1 as a factor, because it is an indicator variable having only values of 0 and 1.
• The second argument, random = ~1 | litter, includes a random effect for each level of LITTER in the model. These random effects will be associated with the intercept, as indicated by ~1.
• The third argument of the function, ratpup, indicates the name of the data frame object to be used in the analysis.
• The final argument, method = "REML", specifies that the default REML estimation method is to be used.
  

By default, the lme() function treats the lowest level (alphabetically or numerically) of a categorical fixed factor as the reference category. This means that “Control” will be the reference category of TREATMENT because “Control” is the lowest level of treatment alphabetically.
We obtain estimates from the model fit by using the summary() function:
> summary(model3.1.fit)
Additional results for the fit of Model 3.1 can be obtained by using other functions in conjunction with the model3.1.fit object. For example, we can obtain F-tests for the fixed effects in the model by using the anova() function:
> anova(model3.1.fit)
The anova() function performs a series of Type I (or sequential) F-tests for the fixed effects in the model, each of which are conditional on the preceding terms in the model specification. For example, the F-test for SEX1 is conditional on the TREATMENT effects, but the F-test for TREATMENT is not conditional on the SEX1 effect. The random.effects() function can be used to display the EBLUPs for the random litter effects:
> # Display the random effects (EBLUPs) from the model.
> random.effects(model3.1.fit)