Fit models to data

This page provides tips and recommendations for fitting linear models to data. 

Here is one important difference between model fitting in R and most other packages. In R, the ANOVA table (but not the coefficients table) tests terms sequentially (in SAS jargon, R uses "type 1 sums of squares"). ANOVA table results will differ from those given in the coefficients table if the data design is not completely balanced. Also, the P-values in the ANOVA table will change if you change the order in which terms are entered in the model. 

There are good reasons for doing it this way, but there are also ways around it, if necessary, as detailed below.

Fit a linear model to data

Linear models are implemented in the lm method in R. You can pass a data frame to the lm command and indicate using the formula which variables you want to fit:

  z <- lm(response ~ explanatory, data = mydata)


The resulting object (which I've named z) is an lm object containing all the results. You use additional commands to pull out these results. Here are some of the most useful commands to extract results from the lm object:TipsT

  summary(z) # coefficients table with parameter estimates

  plot(z)    # plots of residuals, q-q, leverage

  coef(z)    # model coefficients (means, slopes, intercepts)

  confint(z) # confidence intervals for parameters

  resid(z)   # residuals

  fitted(z)  # predicted values

  abline(z)  # adds simple linear regression line to scatter plot

  predict(z,newdata=x1) # predicted values for new observations 

  anova(z)      # ANOVA table (** terms tested sequentially **)

  anova(z1,z2)  # compare fits of 2 models, "full" vs "reduced"


Examples of linear models

The simplest linear model fits a constant, the mean, to a single variable (from which one can obtain a confidence interval for a mean and a test of the null hypothesis that the mean is zero -- equivalent to a one-sample t test). 

  z <- lm(y ~ 1)


The most familiar linear model is the linear regression of y on x.  

  z <- lm(y ~ x)


To fit a linear regression model without an intercept term (regression through the origin) add "-1" to the model statement:

  z <- lm(y ~ x - 1)


If  A  is a categorical variable (factor or character) rather than a numeric variable then the following linear model is the same as a single factor ANOVA.

  z <- lm(y ~ A)


More complicated models include more variables and their interactions. For example, the following is a multiple regression:

  z1 <- lm(y ~ x1 + x2)         # no interaction between x1 and x2

  z2 <- lm(y ~ x1 + x2 + x1:x2) # interaction term present


Analysis of covariance models include both numeric and categorical variables. The method fits a separate linear regression to the data in each group of the categorical variable. If the interaction term is left out, then the regression slopes are forced to be the same in the different groups. Including an interaction term fits different slopes in different groups.

  z <- lm(y ~ x + A)        # no interation

  z <- lm(y ~ x + A + x:A)  # interaction term present


Simple linear regression

Fit the model to data:

  z <- lm(y ~ x)


Add the regression line to a scatter plot.

  plot(x,y)

  abline(z)


The following, more complicated procedure also adds a regression line to a scatter plot, but the approach is more useful in the long term. The order(x) in brackets is needed to draw the lines from left to right (otherwise the lines go back and forth between the unordered x values).

  plot(x,y)

  yhat <- predict(z)

  lines(x[order(x)],yhat[order(x)])


Add confidence bands and prediction intervals to a scatter plot. These might not fit onto the plot in their entirety unless you adjust the limits of the y-axis using the ylim option in the plot command.

  plot(x,y)

  lines(x[order(x)],yhat[order(x)])

  yhat <- as.data.frame(predict(z,interval="confidence"),

                                level=0.95)

  lines(x[order(x)],yhat$lwr[order(x)],lty=2)

  lines(x[order(x)],yhat$upr[order(x)],lty=2)

  yhat <- as.data.frame(predict(z,interval="prediction"),

                                level=0.95)

  lines(x[order(x)],yhat$lwr[order(x)],lty=3)

  lines(x[order(x)],yhat$upr[order(x)],lty=3) 


Here are some useful commands for checking assumptions:

  plot(z)        # residual plots, etc: keep hitting enter

  hist(resid(z)) # histogram of residuals


Here are some useful commands for parameter estimation

  summary(z)             # estimates of slope, intercept, SE's

  confint(z, level=0.95) # conf. intervals for slope and intercept


To obtain the ANOVA table

  anova(z)


Single-factor ANOVA

Before fitting a linear model to the data, check that the categorical variable is a factor. If not, make it a factor, as this will help later when we fit the model.

  is.factor(A)   # or 

  class(A)       # result should be "factor" if A is a factor

  A <- factor(A) # convert categorical variable A to a factor


Before fitting a linear model to the data, check how R has sorted the groups (alphabetically, by default).

  levels(A)      # result will look like: "a" "b" "c" "d"


For reasons that will become clear below, it is often useful to rearrange the order of groups so that any control group is first, followed by treatment groups*. Let's assume that there are four groups and that group "c" is the control. Change the order in R as follows:

  A <- factor(A, levels=c("c","a","b","d"))

*NB: Do not use the command "ordered" for this purpose, as this will change the type of factor you are working with and the way the linear model is fitted.

Or, equivalently, if your variable A is in a data frame, use

  mydat$A <- factor(mydat$A, levels=c("c","a","b","d")


Fit the model to data. y is the numeric response variable and A is the categorical explanatory variable (character or factor).

  z <- lm(y ~ A)

  z <- lm(y ~ A, data=mydata) # if variables are in a data frame


Visualize the fitted model in a graph that includes the data points.
  1. A fast but crude way to add the fitted values to the plot (the lines indicating fitted values may be interrupted)

     stripchart(y ~ A, vertical=TRUE, method="jitter",                             pch=16, col="red")
     stripchart(fitted(z) ~ A, vertical=TRUE, add=TRUE,                            pch="------", method="jitter")

  2. A slightly more tedious way that provides finer control.

      stripchart(y ~ A, vertical=TRUE, method="jitter",
                    pch=16, col="red")
      yhat<-tapply(fitted(z),A,mean)
      for(i in 1:length(yhat)){
        lines(c(i-.2,i+.2),rep(yhat[i],2))
        }

Checking assumptions:

  plot(z)        # residual plots, etc

  hist(resid(z)) # histogram of residuals


See how R is representing the categorical variable and the constant in the linear model by using indicator (dummy) variables.

  model.matrix(z)


Parameter estimation. By default*, R uses the intercept term in the model to estimate the mean of the first group (hopefully, you've set this to be the control group). Remaining parameters estimate the difference between the mean of each group and the first (control) group. Ignore the P-values because they are invalid for tests of the treatment effects except in special cases (i.e., planned contrasts).

  summary(z)             # parameter estimates

  confint(z, level=0.95) # conf. intervals for parameters

* Warning to users of SPlus: the default is different from that in R.

To obtain the ANOVA table

  anova(z)


You can change the way that R models the categorical variable. For example, a useful alternative approach uses a model in which the parameters represent the population means of the different groups. To accomplish this, refit the linear model with a "- 1" in the formula (this suppresses the intercept). 

  z <- lm(y ~ A - 1, data=mydata)

  model.matrix(z)    # optional: shows how R represents categories
  summary(z)         # estimates of population means with SE's*

  confint(z, level=0.95) # conf. intervals for means*
  anova(z)           # results will be the same as before


*The standard errors and confidence intervals for the fitted means are not the same as those you would obtain if you calculated them separately on the raw data for each group separately. This is because the fitted means use the MSresidual from all the data to calculate standard errors. This should generally result in smaller SE's and narrower confidence intervals. The calculation is valid if the assumption of equal variances within different groups is met.

Fit more than one variable

To show the procedures I use an example of a linear model having a numeric response variable (y) and two explanatory variables, one categorical (A) and the other numeric (x).  This combination of variables is known as ANCOVA, or analysis of covariance. Procedures are similar for multiple regression (all variables numeric).

As in the previous section, check that the categorical variable is a factor. If not, convert it to a factor, as this will help later when we fit the model.

  A <- factor(A) # convert categorical variable A to a factor


Also check how R has sorted the groups of the categorical variable (alphabetically, by default). 

  levels(A)      # result will look like: "a" "b" "c" "d"


If necessary, change the order* of the groups so that any control group comes first, followed by treatment groups:

  A <- factor(A, levels=c("c","a","b","d"))

*Do not use the command "ordered" for this purpose, as this will change the type of factor you are working with and the way the linear model is fitted.

Fit the model to data. Generally, fitting the model without an interaction term is preferred, if you can assume that any interaction is weak or absent.

  z <- lm(y ~ x + A)       # no interaction term included, or
  z <- lm(y ~ x + A + x:A) # interaction term present


Visualize the fitted model in a graph that includes the data points.

  plot(x,y,pch=as.numeric(A))

  legend(locator(1), as.character(levels(A)),     
         pch=1:length(levels(A))) # click on plot to place legend

  groups<-levels(A)               # stores group names

  for(i in 1:length(groups)){

    xi<-x[A==groups[i]]            # grabs x-values for group i

    yhati<-fitted(z)[A==groups[i]] # grabs yhat's for group i

    lines(xi[order(xi)],yhati[order(xi)]) # connects the dots

    }


Check model assumptions:

  plot(z)        # residual plots, etc

  hist(resid(z)) # histogram of residuals


Parameter estimation.

  summary(z)             # parameter estimates

  confint(z, level=0.95) # confidence intervals for parameters


The ANOVA table. Notice that the anova command fits terms in the model sequentially1 from the top down ("type I sums of squares"). P-values will therefore differ from those in the coefficients table if the design is unbalanced. The results obtained will then also depend on the order in which you entered the variables in the model formula. Try changing the order to see if this changes your results.

  anova(z)


To exert control over the desired test, fit explicit "full" and "reduced" models yourself and then compare the two model fits using

  anova(z1,z2)  # compare fits of two models, "full" vs "reduced"


1 In SAS jargon, R uses "type 1 sums of squares". Each term is tested by adding it to a model that includes only previous terms in the sequence, with interactions fitted after component variables have been entered. The default in most other statistical packages is instead to use "type III sums of squares", in which each term is tested by adding it to a model that includes all other terms (this can have weird consequences, e.g., testing an individual variable with its interactions present). Click here for a discussion of the pros and cons of different testing approaches.

AIC and related quantities

Then apply the function of interest to the fitted model object, z.

  logLik(z)       # log-likelihood of the model fit

  AIC(z)          # Akaike Information Criterion

  extractAIC(z)   # same, but differs from AIC(z) by a constant

  AIC(z,k=log(n)) # Bayesian Information Criterion (n=sample size)


The function AIC can be used to compare two or more models fitted to the same data. The model with the smallest AIC values is the "best". The same response variable must be fitted in each case, but the explanatory variables can differ.

  z1 <- lm(y ~ x1 +x 2 + x3) # first model

  z2 <- lm(y ~ x3 + x4)      # second model

  c(AIC(z1),AIC(z2))         # vector of AIC values

  AIC(z1,z2)                 # same but returns a data frame


What matters is not AIC itself but the difference in AIC between fitted models. The model with the smallest AIC will have a difference of zero.

  x <- c(AIC(z1),AIC(z2),AIC(z3))  # stores AIC values in a vector

  delta <- x - min(x)              # AIC differences


Finally, we can use these quantities to calculate the Akaike weights for each of a set of fitted models. The model with the highest weight is the "best". The magnitude of the weight of a given model indicates the weight of evidence in its favor, given that one of the models in the set being compared is the best model. The weights are used in model averaging, a technique in multimodel inference. Weights should sum to 1.

  x <- c(AIC(z1),AIC(z2),AIC(z3)) # stores AIC values in a vector

  delta <- x - min(x)             # AIC differences

  L <- exp(-0.5*delta)            # likelihoods of models

  w <- L/sum(L)                   # Akaike weights


A 95% confidence set for the “best” model can be obtained by ranking the models and summing the weights until that sum is ≥ 0.95.