Fit model

Fit models to data

This page provides tips and recommendations for fitting linear and nonlinear models to data. Updated and revised frequently (click the reload button on your browser to make sure you are seeing the most recent version).

The main model-fitting commands covered on this page are:

  • lm (linear models for fixed effects)
  • lme and lmer (linear models for mixed effects)
  • glm (generalized linear models for fixed effects)
  • nls (nonlinear least squares)
  • gam (generalized additive models)

Also covered is the use of

  • visreg (to visualize scatter of points around model fits; works well for fixed effects, but see remarks about random effects here)
  • lsmeans (to estimate predicted means)

Note: The default anova in R fits terms in the model sequentially (“type I sum of squares”). Other statistical packages such as SAS and JMP use marginal testing of terms in anova (“type III sum of squares”) instead. The question of “which set of sums of squares is the Right Thing provokes low-level holy wars on R-help from time to time” (see here).

Under sequential fitting, the sum of squares for each term or factor in the anova table is the improvement in the error sum of squares when that term is added to the linear model, compared with a model including only the terms listed above it (but not those listed below it) in the table. This means that the order in which you list the variables in the lm formula affects the anova table of results, including the P-values. The formula y~A+B+A:B will lead to different sums of squares than the formula y~B+A+A:B when the design is unbalanced. Note that anova also respects hierarchy: the intercept is fitted first, before any other terms, main effects are fitted next, and interactions are fitted last. An interaction is never tested without its corresponding main effects included in the model.

Under marginal testing of terms (“type III sum of squares”), order of appearance of terms in the formula doesn’t matter, and neither does hierarchy. The contribution of each model term is measured by the improvement in the error sum of squares when that term is entered last into the model. Main effects are tested with their interactions already in the model. The Anova function in the car package, combined with a change in the contrasts used to calculate sums of squares, can be used to fit models using type III sum of squares (and also Type II, an in-between solution). Instructions are given below.

A discussion of the pros and cons of type I and type III approaches can be found here, with further explanation provided here.


Fit a linear model to data

Linear models for fixed effects are implemented in the lm method in R. You can pass a data frame to the lm command, using a formula to indicate the model you want to fit.

z <- lm(response ~ explanatory, data = mydata)
z <- lm(response ~ explanatory, data = mydata, na.action = na.exclude)

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, including residuals and predicted values. The argument na.action = na.exclude is optional --- it tells R to keep track of cases having missing values, in which case the residuals and predicted values will have NA's inserted for those cases. Otherwise R just drops missing cases.

Here are some of the most useful functions to extract results from the lm object:

summary(z)     # parameter estimates and overall model fit
plot(z)        # plots of residuals, normal quantiles, leverage
coef(z)        # model coefficients (means, slopes, intercepts)
confint(z)     # confidence intervals for parameters
resid(z)       # residuals
predict(z)     # predicted values
predict(z, newdata = mynewdata) # see below
fitted(z)      # predicted values
anova(z1, z2)  # compare fits of 2 models, "full" vs "reduced"
anova(z)       # ANOVA table (** terms tested sequentially **)

The command predict(z, newdata = mynewdata) will used the model to predict values for new observations contained in the data frame mynewdata. The explanatory variables (the X variables) in mynewdata must have exactly the same names as in mydata.

In the case of simple linear regression you can add a line to a scatter plot of the data using the abline function.

plot(response ~ explanatory, data = mydata)
abline(z)

Basic types of linear models

In the examples below, I'm assuming that all the variables needed are in a data frame mydata.

The simplest linear model fits a constant, the mean, to a single variable. This is useful if you want to obtain an estimate of the mean with a standard error and confidence interval, or wanted to test the null hypothesis that the mean is zero.

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

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

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

To fit a linear regression model without an intercept (i.e., fit a regression through the origin) use one of the following modifications:

z <- lm(y ~ x - 1, data = mydata)
z <- lm(y ~ 0 + x, data = mydata)

If A is a categorical variable (factor or character) rather than a numeric variable then the model conforms to a single factor ANOVA instead.

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

More complicated models include more variables and their interactions

z1 <- lm(y ~ x1 + x2, data = mydata)          # no interaction between x1 and x2
z2 <- lm(y ~ x1 + x2 + x1:x2, data = mydata)  # interaction term present
z2 <- lm(y ~ x1 * x2, data = mydata)          # interaction term present

Analysis of covariance models include both numeric and categorical variables. The linear model in this case is a separate linear regression for each group of the categorical variable. Interaction terms between these two types of variables, if included in the model, fit different linear regression slopes; otherwise the same slope is forced upon every group.

z <- lm(y ~ x + A + x:A, data = mydata)

Simple linear regression

Fit the model to data. Note: avoid calling the variables mydata$y and mydata$x in the lm formula, so that R knows to find the variables in mydata in later commands.

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

Add the regression line to a scatter plot.

plot(y ~ x, data = mydata)
abline(z)

To obtain the regression coefficients (parameter estimates) with confidence intervals, as well as R2 values,

summary(z)               # estimates of slope, intercept, SE's
confint(z, level = 0.95) # 95% conf. intervals for slope and intercept

To test the null hypothesis of zero slope with the ANOVA table,

anova(z)

The following plots a tidier regression line that doesn't extend beyond the data points. The method calculates the predicted values for the smallest and largest values of the explanatory variable and connects them with a line. Variables my.x and my.y are in the data frame mydata.

plot(my.y ~ my.x, data = mydata)
z <- lm(my.y ~ my.x, data = mydata)
xpoints <- range(my.x)
yhat <- predict(z, newdata = data.frame(my.x = xpoints))
lines(yhat ~ xpoints)

Add confidence bands to your previous scatter plot.
xpoints <- seq(min(my.x), max(my.x), length.out = 100) ypoints <- data.frame(predict(z, newdata = data.frame(my.x = xpoints), interval = "confidence", level = 0.95)) lines(ypoints$lwr ~ xpoints, lty = 2) lines(ypoints$upr ~ xpoints, lty = 2) Adding prediction intervals instead is similar:

xpoints <- seq(min(my.x), max(my.x), length.out = 100)
ypoints <- data.frame(predict(z, newdata = data.frame(my.x = xpoints), 
                   interval = "prediction", level = 0.95))
lines(ypoints$lwr ~ xpoints, lty = 2)
lines(ypoints$upr ~ xpoints, lty = 2)

The intervals might not all squeeze onto the plot unless you adjust the limits of the y-axis using the ylim = option in the plot command.

The visreg package provides tools for plotting the fitted regression line along with confidence intervals of predicted mean values (confidence bands). You'll need to install the package if you haven't already done so. Then load it before using the first time in an R session. The following commands work for simple linear regression.

library(visreg)
z <- lm(y ~ x, data = mydata)                        # regression model
visreg(z)                                            # basic plot
visreg(z, points.par = list(cex = 1.2, col = "red")) # with points options

Here are some useful diagnostic plots for checking assumptions. Most of the plots produced by plot(z) will be self-explanatory, except perhaps the last one. "Leverage" calculates the influence that each data point has on the estimated parameters. "Cook's distance" measures the effect of each data point on the predicted values for all the other data points. A value greater than 1 is said to be worrisome.

plot(z)        # residual plots, etc: keep hitting enter
hist(resid(z)) # histogram of residuals

Finally, here are a couple of useful tricks to constrain parameter values for simple linear regression. x and y are numeric variables in the data frame mydata.

z <- lm(y ~ x, data = mydata)               # Ordinary linear regression, for comparison
z <- lm(y ~ 1, data = mydata)               # Force the regression line through the origin
z <- lm(y ~ 1 + offset(1*x), data = mydata) # Force a slope of 1 through the data
z <- lm(y ~ 1 + offset(b*x), data = mydata) # Force a slope of "b" through the data 
                                            #   (b must be a number)

Forcing a line through the origin (forcing an intercept of 0) results in only the slope being estimated. Forcing a slope of 1 or b through the data results in only the intercept being estimated. Unfortunately, visreg doesn't work well for these models.

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(mydata$A)          # or
class(mydata$A)              # result should be "factor" if A is a factor
mydata$A <- factor(mydata$A) # convert categorical variable A to a factor

Check how R has sorted the factor groups (levels). You will see that R orders them alphabetically, by default, which isn't necessarily the most useful order.

levels(mydata$A)

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, "a" to "d", and that group "c" is the control. Change the order in R as follows:

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

*NB: Do not use the command "ordered" for this purpose, which will mess things up.

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

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

Visualize the fitted model in a graph that includes the data points. The second command below is a fast but crude way to add the fitted values to the plot.

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

The visreg package provides a fast tool to visualize the model fit with a strip chart, including optional 95% confidence intervals for the predicted means.

library(visreg)                                     # loads package, if not done yet
z <- lm(y ~ A, data = mydata)                       # fit model
visreg(z)                                           # basic plot
visreg(z, points.par = list(cex = 1.2, col = "red") # with points options
visreg(z, whitespace = 0.4)                         # adjust space between groups

Or, draw the predicted values yourself.

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

The lsmeans package provides a convenient way to estimate the group means fitted by the a single factor ANOVA model, including SE's and confidence intervals. You'll have to install and load the package first, if you haven't already done so.

library(lsmeans)
lsmeans(z, "A", data=mydata)   # use the real name of your factor in place of "A"

Note that these are not the SE's and confidence intervals you would calculate yourself on each group separately, because they use the residual mean square from the model fitted to all the data. 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.

Check assumptions:

plot(z)        # residual plots, etc
hist(resid(z)) # histogram of residuals

Estimate parameters. By default in R, the intercept term estimates the mean of the first group (set this to be the control group or the most convenience reference 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 hypothesis tests of the differences except in the special case of a planned contrast.

summary(z)               # coefficients table with estimates
confint(z, level = 0.95) # conf. intervals for parameters

Test the null hypothesis of equal group means with the ANOVA table.

anova(z)

See how R is representing the categorical variable behind the scenes by using indicator (dummy) variables.

model.matrix(z)

Fit multiple factors

When two or more categorical factors are fitted, the possibility of an interaction might also be evaluated. Remember that in R the order in which you enter the variables in the formula affects the anova results. By default, anova fits the variables sequentially ("type I sum of squares"), while at the same time respecting hierarchy. Some other statistics programs such as SAS and JMP use marginal fitting of terms instead ("type III sum of squares") instead. See notes on this at the top of the current page.

In the case of two factors, A and B, the linear model is as follows.

z <- lm(y ~ A + B, data = mydata) # additive model, no interaction fitted
z <- lm(y ~ A * B, data = mydata) # main effects and interaction terms fitted
summary(z)                        # coefficients table
confint(z, level = 0.95)          # confidence intervals for parameters
anova(z)                          # A is tested before B; interaction tested last

We can use the Anova command in the car package to give us the anova table based on marginal fitting of terms ("type III sum of squares"). To accomplish this we first need to override R's default coding of categorical variables, as follows.

z <- lm(y ~ A * B, contrasts = c("contr.sum", "contr.poly"), data = mydata)
library(car)
Anova(z, type = 3)

Use the visreg package to visualize the scatter of data around the model fitted values, and the contribution of each variable. When a subset of multiple variables is plotted, the method visualizes the predicted values for the plotted variable(s) while conditioning on the value of the other variable (adjusted to the most common category for other factor). In this case you are seeing a slice of the model but not the overall predicted means or fitted values. Note that this conditioning makes sense only when no interaction term is fitted. When all factors are plotted (e.g., both factors A and B), then the predicted values are indeed shown for every combination of groups (levels). You can choose to overlay the fits or plot in separate panels. Try different options to see which method is most effective.

library(visreg)
visreg(z, xvar = "B", whitespace = 0.4, 
    points.par = list(cex = 1.1, col = "red"))
visreg(z, xvar = "A", by = "B", whitespace = 0.4, 
    points.par = list(cex = 1.1, col = "red"))
visreg(z, xvar = "A", by = "B", whitespace = 0.5, overlay = TRUE, 
    band = FALSE, points.par = list(cex = 1.1))

Use the lsmeans package to obtain fitted group means. In the formulas below, use the real names of your factors in place of "A" and "B".

library(lsmeans)
z <- lm(y ~ A + B, data = mydata) # additive model, no interaction
lsmeans(z, "A", data = mydata)         # group means for A variable, averaged over levels of B
lsmeans(z, c("A", "B"), data = mydata) # model fitted means* for all combinations of A and B groups

z <- lm(y ~ A * B, data = mydata) # main effects and interaction terms fitted
lsmeans(z, c("A", "B"), data = mydata) # means for all combinations of A and B groups

These will not usually be the actual group means. Instead, they are the predicted means for group combinations under the constraint that no interaction between A and B is present.

Fit a factor and a continuous covariate

Here, 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). (Fitting a model with both numerical and categorical variables is also known as ancova, or analysis of covariance.) Note that the order of terms in your formula affects the sums of squares and anova tests. See the notes at the top of this page on sequential vs marginal fitting of terms. To fit terms sequentially, use the following methods.

z <- lm(y ~ x + A, data = mydata) # no interaction term included, or
z <- lm(y ~ x * A, data = mydata) # interaction term present
summary(z)                        # coefficients table
confint(z, level = 0.95)          # confidence intervals for parameters
anova(z)                          # sequential fitting of terms

Check model assumptions:

plot(z)        # residual plots, etc
hist(resid(z)) # histogram of residuals

Use Anova command in the car package to carry out marginal fitting of terms instead ("type III sum of squares"). First, we need to override R's default contrasts.

z <- lm(y ~ x * A, contrasts = c("contr.sum", "contr.poly"), data = mydata)
library(car)
Anova(z, type = 3)

Use the visreg package to help visualize the scatter of the data around the model fits and the contributions of each variable. When you plot a subset of the variables in the model, visreg shows you their contributions when the data values are adjusted for the other factors, conditioning on the most common factor level (for categorical factors) or the median value (for numeric variables). In this case you aren't looking at the overall model fitted values. If you plot more than one variable, you can choose to overlay the visuals for the different levels of the factor, or you can plot in separate panels.

library(visreg)
visreg(z, xvar = "A", whitespace = 0.4, 
    points.par = list(cex = 1.1, col = "red"))
visreg(z, xvar = "x", by = "A", 
    points.par = list(cex = 1.1, col = "red"))
visreg(z, xvar = "x", by = "A", whitespace = 0.4, overlay = TRUE, 
    band = FALSE, points.par = list(cex = 1.1))

Or, add the separate regression lines for each group to the scatter plot yourself.

plot(y ~ x, pch = as.numeric(A), data = mydata)
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(yhati[order(xi)] ~ xi[order(xi)]) # connects the dots
    }

Fit a linear mixed-effects model to data

Linear mixed-effects models are implemented in the lme method of the nlme package in R, and in the lmer method of the lme4 package. These method use restricted maximum likelihood (REML) to produce unbiased estimates of model parameters and to test hypotheses.

The basic structure of the lme command includes two formulas. One is for the fixed effects and the second is for the random effects. For example, a model for a single random effect B would looks like the following.

library(nlme)
z <- lme(y ~ 1, random = ~ 1|B, data = mydata) # avoid mydata$y and mydata$B notation

The formula for the fixed effects comes first and includes the name of the response variable (here, an intercept is fitted). The formula for the random effects part leaves out the name of the response variable (it would be redundant) but it must start with "random=". The formula for the random effects indicates the random groups, indicated here with the variable B in mydata

The resulting object (here named z) is an lme object containing all the results. Here are some of the most useful commands to extract results from the lme object:

summary(z)                # parameter estimates, fit, tests of fixed effects
plot(z)                   # plot of residuals against predicted values
intervals(z)              # conf. intervals for fixed effects and variances

resid(z)                  # residuals
fitted(z)                 # best linear unbiased predictors (BLUPs)
VarCorr(z)                # variance components for random effects  

anova(z)                  # Tests of fixed effects (fitted sequentially)
anova(z, type="marginal") # Test of fixed effects using "drop-one" testing

A frequent error when analyzing data with lme is using the same labels to identify different sampling units. For example, the same number codes 1 through 5 might be used to label the 5 subplots of every plot. This will confuse lme and lead to much user frustration. Instead, use the codes a1 through a5 to label the subplots from plot "a", b1 through b5 to label the subplots taken from plot "b", and so on.

To use lmer instead, you will need to install the lmerTest and lme4 packages if you haven't already done so. The lmer command requires only one formula, with random effects specified within parentheses, as follows (compare with the example lme formula above).

library(lmerTest)
z <- lmer(y ~ 1 + (1|B), data = mydata)

Most of the above functions such as summary work the same way with lmerTest objects as with lme objects except for intervals and anova, whose use is as follows.

confint(z)                # use instead of intervals(z)
anova(z, type = 1)        # use type=2 or type=3 for alternative sums of squares

Example: One random factor -- repeated measurements of individuals

The simplest example of a model with random effects is an observational study in which multiple measurements are taken on a random sample of individuals (or units, plots, etc) from a population. In this case "individual" is the random factor or group.

A frequent application is estimation of the repeatability of a trait or measurement. Lessels and Boag (1987, Auk 104: 116-121) define repeatability as "the proportion of variance in a character that occurs among rather than within individuals". It can be thought of as the correlation between repeated measurements made on the same individuals.

In the example below, a numeric variable y is measured two (or more) times on randomly sampled individuals whose identities are recorded in the variable B. Both variables are in the data frame mydata.

z <- lme(y ~ 1, random = ~ 1|B, data = mydata) # nlme package
z <- lmer(y ~ 1 + (1|B), data = mydata)        # lmerTest package

Two formulas in one command can take getting used to. The first formula is the fixed part, and in this case it fits an intercept, the grand mean of the fitted values of the random groups. lme will estimate this mean and provide a standard error in the output. The second formula is the random part and it fits an intercept to each of the random "groups" of measurements (here, a group consists of the repeated measurements of a single individual). lme won't automatically show you these intercepts because they are not interesting by themselves. Instead, the method uses the data to estimate the variance (and standard deviation) of intercepts in the population, which is what we are interested in.

Use the following commands to obtain estimates of the parameters (grand mean, for the fixed part; standard deviation and variances among groups for the random part, as well as the variance of the residuals). The output also includes numbers to indicate how well the model fits the data (AIC, BIC, logLik). We will learn more about these measures of fit after we cover likelihood methods.

summary(z)
intervals(z) # 95% confidence intervals, where z is lme object
confint(z)   # use instead of intervals(z) if z is lmerTest object

Repeatability of a measurement is calculated as
r = varamong / (varamong + varwithin)

where varamong is the variance among the means of groups (individuals, in this case), and varwithin is the variance among repeat measurements within groups. Both these quantities are estimated in the lme fit. A quick way to get them is

VarCorr(z)

In the output of the VarCorr command, the estimated variance among groups is shown as the variance associated with the random intercepts. The variance associated with the residual is the estimate of the within-group variance. Plugging both numbers into the equation above yields the estimated repeatability.

To see a plot of the residual values against fitted values,

plot(z)

You might notice a positive trend in the residuals, depending on the amount of variation between repeat measurements and the number of repeat measurements per individual in the data. To confirm what is happening, plot the original data and superimpose the lme fitted (predicted) values:

stripchart(y ~ B, vertical = TRUE, pch = 1)
stripchart(fitted(z) ~ B, vertical = TRUE, add = TRUE, pch = "---")

Notice that the lme predicted values are smaller than the mean for individuals having large values of y, whereas the lme predicted values are larger than the mean for individuals having small values of y. This is not a mistake. lme obtains "best linear unbiased predictors" for the random effects, and these tend to lie closer to the grand mean than the mean of the actual measurements of each individual. Basically, lme is predicting that the "true" trait value for individuals having extreme measurements is likely to be closer, on average, to the grand mean than is the average of the repeat measurements. This effect diminishes the more measurements you have taken of each individual.

Unfortunately, visreg doesn't work in this simple case because the method requires at least one fixed factor in the model (not just an intercept).

Example: Two nested random factors -- nested sampling designs

The above procedures can be extended to study designs in which sampled units are nested within other units, and there are no fixed experimental treatments.

For example, a study might measure a variable y in quadrats, where quadrats are randomly sampled within transects, transects are randomly placed within woodlots, and woodlots are randomly sampled from a population of woodlots in a region. Let y be the measurement obtained from each quadrat. Variable A identifies the transect, and B indicates the woodlot. The multiple values of y for each combination of A and B are the quadrat measurements. This sampling design is modeled as

z <- lme(y ~ 1, random = ~ 1|B/A, data = mydata) # if using nlme package
z <- lmer(y ~ 1 + (1|B/A), data = mydata)        # if using lmerTest package

where "B/A" expresses the nesting: "woodlot, and transect within woodlot".

This same formulation could be used to model the results of a half-sib breeding design in quantitative genetics, where multiple offspring are measured per dam (A), and multiple dams are mated to each sire (B).

Take special care to label sampling units (e.g., transects) uniquely. For example, if you have 5 transects in each woodlot, don't give them the same labels 1 through 5 in each woodlot. Instead, label those from woodlot 1 as w1.1 through w1.5, and those from woodlot 2 as w2.1 through w2.5. Otherwise R won't compute the correct quantities.

Example: One fixed and one random factor -- randomized block, and subjects-by-treatment repeated measures designs

A randomized complete block (RCB) design is like a paired design but for more than two treatments. In the typical RCB design, every treatment is applied exactly once within every block. This yields a single value for the response variable from every treatment by block combination. In this case it is not possible to fit an interaction term between treatment and block.

In field experiments, block are usually experimental units that are grouped by location or in time. In lab experiments blocks may be separate environment chambers, or separate shelves within a single chamber. Used this way, blocks aren't true random samples from a population of blocks. Nevertheless blocks are usually modeled as random effects when the data are analyzed.

The data from a RCB experiment include the response variable y, a categorical variable A indicating the treatment, and a categorical variable B indicating block. If you plan to use visreg to visualize the scatter of the data around the model fits, convert A and B to factors before using lme.

z <- lme(y ~ A, random = ~ 1|B, data = mydata)
z <- lmer(y ~ A + (1|B), data = mydata)        # equivalent command if using lmerTest

The only difference from the nested sampling example above is that now the fixed part of the model includes the treatment variable A as an explanatory variable. The random part of the model again fits an intercept to each of the random groups defined in B (subjects or blocks).

"Subjects by treatment" repeated measures designs are often analyzed in the same way as RCB designs. Here, each treatment is assigned to every subject in random order. The analysis assumes that there is no carry-over between treatments: earlier measurements made on a subject under one treatment should not influence later measurements taken on the same subject in another treatment. The categorical variable B indicates the identities of the subjects.

You can use interaction.plot to view how the mean y changes between treatment levels of A separately for each level of B. Parallel lines would indicate a lack of an interaction between A and B. The drawback of this function is that it does not show the data (not a problem if you have just one data point for each combination of categories of A and B).

interaction.plot(A, B, y)

The visreg function seems to work in this case to visualize lme model fits to data (although it won't include the confidence intervals). If A is the fixed (treatment) effect and B is the random effect, try one of the following.

visreg(z, xvar = "A")
visreg(z, xvar = "A", by = "B", scales=list(rot = 90))

To obtain estimates of the fixed effects and variance components,

summary(z)
intervals(z) # 95% confidence intervals
confint(z)   # use instead of intervals(z) if using lmerTest package

Use lsmeans to obtain model-based estimates of the treatment (A) means.

library(lsmeans)
lsmeans(z, "A", data = mydata)

The anova command will test the fixed effects (the grand mean and the treatment A). In nlme the anova command tests the fixed effects sequentially (see explanation at top of this page). To use marginal testing ("drop-one" or "type III sums of squares") instead, set the option type="marginal".

anova(z)                    # for sequential testing
anova(z, type = "marginal") # for drop-one testing

In lmerTest the anova command tests the fixed effects using marginal testing ("drop-one", or "type III sums of squares") by default. To use sequential testing, use the option type=1 instead.

anova(z)           # for drop-one (type III) testing
anova(z, type = 1) # for sequential testing

Example: One fixed and one random factor -- factorial design

The approach to a factorial design with one fixed and one random factor is similar to the randomized complete block design except that each combination of levels of the fixed and random effect may be replicated.

Two-factor ANOVA is more complex when one of the factors is random because random sampling of groups adds extra sampling error to the design. This sampling error adds noise to the measurement of differences between group means for the fixed factors that interact with the random factor.

If A is the fixed factor and B is the random factor, use the following if you want to include the interaction between A and B (which is also a random effect)

z <- lme(y ~ A, random = ~ 1|B/A, data = mydata) # nlme package
z <- lmer(y ~ A + (1|B/A), data = mydata)        # lmerTest package

The corresponding formulas without an interaction term are

z <- lme(y ~ A, random = ~ 1|B, data = mydata) # nlme package
z <- lmer(y ~ A + (1|B), data = mydata)        # lmerTest package

To plot the residuals against the fitted values,

plot(z)

To test the fixed effects, use

anova(z)           # if using nlme package
anova(z, type = 1) # if using lmerTest package

Example: One numeric variable and one random factor -- linear regression in multiple random groups

In the examples below, y is the response variable, x is the numeric explanatory variable, and B identifies the random groups, with all variables in a data frame mydata. x is always modeled as a fixed effect.

The simplest model for linear regression in multiple random groups assumes equal slopes within each group. In this case, only a random intercept is modeled in each group.

z <- lme(y ~ x, random = ~ 1|B, data = mydata) # nlme package
z <- lmer(y ~ x + (1|B), data = mydata)        # lmerTest package

If slopes and intercepts both vary among random groups, then the model includes an interaction between x and B, which is a random effect, as follows.

z <- lme(y ~ x, random = ~ x|B, data = mydata) # nlme package
z <- lmer(y ~ x + (x|B), data = mydata)        # lmerTest package

Example: Two random factors -- factorial design

Use lmer in the lmerTest package if fitting models to a factorial design with two crossed random factors, A and B.

z <- lmer(y ~ (1|A) + (1|B), data = mydata)           # without interaction
z <- lmer(y ~ (1|A) + (1|B) + (1|A:B), data = mydata) # with interaction

Accomplishing the same goal using lme in the nlme package is more challenging, and is explained here).

To obtain estimates of the fixed effects and variance components, use

summary(z)

In this model there are no fixed effects to test, so the anova command will not accomplish anything.

Example: Two fixed and one random factor -- split plot design

A basic split plot design has two treatment variables, which are fixed factors, A and B. One of these is applied to whole plots (ponds, subjects, etc) as in a completely randomized design. All levels of the second factor B are applied to subplots within every plots. For example, if there are two treatment levels of factor B (e.g., a treatment and a control), then each plot is split in two; one of the treatments is applied to one side (randomly chosen) and the second is applied to the other side. In this case, each side of a given plot yields a single value for the response variable. Because levels of B are repeated over multiple plots within each level of A, it is possible to investigate the effects of each treatment variable as well as their interaction.

Let "plot" refer to the variable identifying the random plot (pond, subject, etc). The response variable is modeled as

z <- lme(y ~ A * B, random = ~ 1|plot, data = mydata) # with nlme package
z <- lmer(y ~ A * B + (1|plot), data = mydata)        # with lmerTest package

Notice that "A * B" is shorthand for "A + B + A:B", where the last term is the interaction between A and B.

Plot the residuals against the fitted values in a partial check of assumptions.

plot(z)

To visualize the model fit, visreg seems to work in this case, since there is only one random factor.

visreg(z, xvar = "A", by = "B", scales=list(rot = 90))

Use lsmeans to generate model-based fitted means for treatments.

lsmeans(z, c("A", "B"), data = mydata) # means of each treatment combination
lsmeans(z, c("A"), data = mydata)      # means of "A" treatment levels averaged over levels of "B"

To get an overall test of the fixed effects and their interactions, use the following if using the nlme package.

anova(z)                    # for sequential testing
anova(z, type = "marginal") # for drop-one testing

Use the following instead if using the lmerTest package.

anova(z, type = 1) # for sequential testing
anova(z)           # for drop-one testing

Fit a generalized linear model

Generalized linear models for fixed effects are implemented using glm in R. The glm command is similar to the lm command, used for fitting linear models, except that the error distribution and link function must be specified using the "family" argument.

For example, to model a binary response variable, specify a binomial error distribution and the logit link function:

z <- glm(response ~ explanatory, family = binomial(link="logit"), data = mydata)

For count data, specify the Poisson error distribution and the log link function.

z <- glm(response ~ explanatory, family = poisson(link="log", data = mydata)

If the variance of the error distribution in the data is greater than that expected under the binomial and poisson distributions ("overdispersion"), try the following instead,

family = quasibinomial(link = "logit") # or
family = quasipoisson(link = "log"))

These error distributions model overdispersion for binary and count data, and additionally provide an estimate of the dispersion parameter (a value greater than one indicates overdispersion).

Type help(family) to see other error distributions and link functions that can be modeled using glm.

The resulting object (which I've named z) is a glm object containing all the results. Here are some of the most useful commands used to extract results from the glm object:

summary(z)                    # parameter estimates and overall model fit
coef(z)                       # model coefficients
resid(z)                      # deviance residuals
predict(z)                    # predicted values on the transformed scale
predict(z, se.fit = TRUE)     # Includes SE's of predicted values
fitted(z)                     # predicted values on the original scale
anova(z, test = "Chisq")      # Analysis of deviance - sequential
anova(z1, z2, test = "Chisq") # compare fits of 2 models, "reduced" vs "full"
anova(z, test = "F")          # Use F test for gaussian, quasibinomial or quasipoisson link

"Deviance residuals" are not the same as ordinary residuals in linear models. They are a measure of the goodness of fit of the model to each data point. Each residual can be thought of as the contribution of the corresponding data point to the residual deviance (given in the analysis of deviance table).

There is no plot(z) method for glm model objects. If you enter plot(z), the resulting plots will look strange and won't be easy to interpret. R is using the plot function for lm model objects, and they aren't valid for glm objects. Instead, it is more useful to visualize the data by adding the fitted curves or means onto scatter plots or other graphs that show the data.

A method for calculating likelihood based confidence intervals for parameters is available in the MASS library. A function to estimate LD50's is also found here, for the case of a dose-response curve with a binary response variable.

library(MASS)
confint(z, level = 0.95) # approximate 95% confidence intervals
dose.p(z, p = 0.50)      # LD50 for a dose-response curve

Logistic and log-linear regression

These methods are analogous to linear regression. We assume that x is a continuous numeric explanatory variable and that the response variable y is either binary data (logistic regression) or count data (log-linear regression). The only difference between logistic and log-linear regression is the error distribution and link function, which are specified by the family argument.

z <- glm(y ~ x, family = binomial(link="logit"), data = mydata)
z <- glm(y ~ x, family = poisson(link="log"), data = mydata)

You can use visreg to visualize the model on the transformed scale (the function uses predict(z) to generate the result). Replace x in the command below with the name of your x-variable (in quotes). The glm method of visreg fits a linear model on the transformed scale by default. This plot will show points, but they are not the transformed data, however. They are "working values" obtained by transforming the residuals of the fitted model on the original scale. glm repeatedly recalculates the working values and the fitted model as it converges on the maximum likelihood estimate. visreg shows you the results from the final iteration.

Use the scale = "response" argument to see the fitted relationship on the original scale. You can then add the data points, but make sure you have enough room on your y-axis to plot all the points.

library(visreg)
visreg(z, xvar = "x") # fit on the transformed scale

visreg(z, xvar = "x", ylim = range(y), scale = "response")  # fit on the original scale
points(y ~ x)                                               # add the data points

Use fitted(z) to generate the predicted values corresponding to your data points on the original scale of the data. Since the response variable is discrete I have added a bit of jitter to the points to minimize overlap.

plot(jitter(y, amount = 0.02) ~ x, data = mydata)
yhat <- fitted(z)
lines(yhat[order(x)] ~ x[order(x)])

Approximate confidence bands can also be added to your scatter plot. visreg can do this, but here I show how to calculate and plot them in base R. Confidence intervals are based on a normal approximation and assume a large sample size -- they are not very accurate for small and moderate sample sizes. Replace the 1.96 in the command below with 1 to plot standard error bands instead. First, calculate the upper and lower limits as follows.

zhat <- predict(z, se.fit = TRUE)       # result on logit or log scale
zupper <- zhat$fit + 1.96 * zhat$se.fit
zlower <- zhat$fit - 1.96 * zhat$se.fit

Then convert to the original scale.

yupper <- exp(zupper)/(1 + exp(zupper)) # for logit link
ylower <- exp(zlower)/(1 + exp(zlower))
# or
yupper <- exp(zupper)                   # for log link
ylower <- exp(zlower)

Finally, plot the bands.

lines(yupper[order(x)] ~ x[order(x)], lty = 2)
lines(ylower[order(x)] ~ x[order(x)], lty = 2)

Use summary to view the parameter estimates (slope and intercept) on the logit or log scale, along with standard errors. Note that P-values in the summary() output are based on a normal approximation and are not accurate for small to moderate sample sizes -- use the log likelihood ratio test instead (see below).

summary(z)

To obtain 95% likelihood-based confidence intervals for the parameters, use a command from the MASS library

library(MASS)
confint(z, level = 0.95)

Use the following command to obtain the analysis of deviance table. This provides log-likelihood ratio tests of the most important parameters. Model terms are tested sequentially, and the results will depend on the order in which variables are entered into the model formula.

anova(z, test = "Chisq")

Categorical explanatory variable

Glm models can also be fitted to data when the explanatory variable is categorical rather than continuous. This is analogous to single-factor ANOVA in lm, but here the response variable is binary or has three or more categories.

The glm models are as follows. A is the categorical variable identifying treatment groups.

z <- glm(y ~ A, family = binomial(link="logit"), data = mydata)
z <- glm(y ~ A, family = poisson(link="log"), data = mydata)

As in lm, it may be useful to order the different groups (categories of A) so that a control group is first, followed by treatment groups. This can help interpretation of the model parameters. For example, if there are 4 groups in A and group "c" is the control group, set the order of the group

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

As in lm, you can also change the way that R models the categorical variable. For example, a useful approach is to fit a model in which the parameters are the group means (in glm these means are fitted on the logit or log scale before being back-transformed to the original scale). To accomplish this, refit the linear model with a "-1" in the formula, which leaves out an intercept.

z <- glm(y ~ A - 1, family = binomial(link = "logit"), data = mydata)
z <- glm(y ~ A - 1, family = poisson(link = "log"), data = mydata)

Use lsmeans to generate model-based fitted means for treatments.

library(lsmeans)                   # load lsmeans library
lsmeans(z, c("A"), data = mydata)  # means of "A" treatment on transformed scale

You'll need to transform the means and the lower and upper limits of the confidence intervals to the original scale using the inverse of the link function.

Modeling contingency tables

Binary response data ("success" = 1, "failure" = 0) are often compared among groups using contingency tables. For example, here are fictitious data on the number of successes and failures of independent units belonging to one of 4 groups identified as levels of the variable A. The first row below has the variable names of each column.

A     successes    failures    n
a         5           10      15
b        15           20      35
c        25           30      55
d        35           40      75

There are two approaches to analyzing these data with glm. Both yield identical results when the response variable is binary. Only the second approach can be used when the response variable has more than two categories.

The first approach is to fit a glm model to the proportion of successes, taking care to indicate the number of observations on which each proportion is based using the weights option. Assume that the variables "successes" and n are in your data frame mydata. First calculate the proportion of successes.

mydata$prop <- mydata$successes/mydata$n
z <- glm(prop ~ A, weights = n, family = binomial(link = "logit"), data = mydata)

The second approach is to model the interaction between the grouping variable A and a single outcome variable (success/failure). To use this approach we need to reshape the data as a data frame of frequencies,

A      freq    outcome
a         5    success
b        15    success
c        25    success
d        35    success
a        10    failure
b        20    failure
c        30    failure
d        40    failure

We then fit the frequencies with a glm model whose linear predictor includes the grouping variable A, the outcome variable outcome, and their interaction

z <- glm(freq ~ A * outcome, family = poisson(link = "log"), data = mydata)
summary(z)
anova(z, test="Chisq")

The summary coefficients of interest are those corresponding to the interaction between A and outcome, rather than coefficients for the main effects for A and outcome. The interaction measures the extent to which probability of success or failure depends on A, which is the point of the analysis. Similarly, the interaction is the term of principal interest in the analysis of deviance table.

Checking assumptions of logistic and log-linear regression

Generalized linear models are freed from the assumption that residuals are normally distributed with equal variance, but the method nevertheless makes important assumptions that should be checked.

Unfortunately there is no plot(z) method for glm model objects. If you enter plot(z), the resulting plots will look strange because the data are discrete, and they won't be easy to interpret. (R is using the plot function for lm model objects, which aren't valid for glm objects.) Instead, it is more useful to graph the fitted curves or means onto plots of the data.

To check whether the fitted curve from the logistic or log-linear regression accurately captures the trend in the data, you might compare it to that of a lowess fit. The lowess method calculates a smooth curve that estimates the trend in the data (rather like a running mean). Is it similar to the glm fitted curve?

plot(jitter(y, amount = 0.02) ~ x, data = mydata)
lines(lowess(mydata$x, mydata$y))

To use an analogous smoothing method that also takes the error distribution and link function into account, try fitting a cubic spline instead of using lowess. See the section "Fit a Cubic Spline" elsewhere on this page for help.

Logistic and log-linear regression assume that the "dispersion parameter" is 1. What this means is that the variance of the residuals are as expected from the binomial or Poisson distributions specified by the link function. In the case of count data, for example, a dispersion parameter of 1 means that the variance of the explanatory variable within each group should be approximately equal to the mean for the group. However, count data are notorious for having higher variances than expected (dispersion parameter > 1). One way to check is to tabulate the mean and variance for each group and see how they compare. Another way is to estimate the dispersion parameter directly by refitting the data to a quasi-likelihood model that includes such a parameter. For binary data (logistic regression) the quasi-likelihood model is fitted as follows.

z <- glm(y ~ x, family = quasibinomial(link="logit"), data = mydata)
summary(z)
anova(z, test = "F")

For count data the corresponding method is

z <- glm(y ~ x, family = quasipoisson(link="log"), data = mydata)
summary(z)
anova(z, test = "F")

If the estimated dispersion parameter is much greater than 1, then you should analyze your data using the quasibinomial or quasipoisson link function. This will lead to more reliable standard errors for parameters and more accurate P-values in tests. Note that you should use anova(z, test = "F"), i.e., an F test, when testing hypotheses that use gaussian, quasibinomial or quasipoisson link functions. This is because the quasi-likelihood is not a real likelihood and the generalized log-likelihood ratio test is not accurate.


Model selection

R includes a number of model selection tools to evaluate and compare the fits of alternative models to data. Most of the methods shown below work on lm objects as well as lme, glm, nl, gls and gam objects.

AIC and related quantities

These functions compute log-likelihoods and the Akaike information criterion for fitted models. To use them, begin by fitting the model of interest to the data,

z <- lm(y ~ x1 + x2 + x3, data = mydata)

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
AIC(z, k = log(n)) # Bayesian Information Criterion (n is the sample size)
BIC(z)             # Bayesian Information Criterion

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 + x2 + x3, data = mydata)  # first model
z2 <- lm(y ~ x3 + x4, data = mydata)       # 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)            # relative 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.

BIC and related quantities

The computations are similar for BIC (Bayesian Information Criterion) as for AIC. Under the BIC criterion the model with the smallest BIC is the "best" in the set of models fitted.

To calculate BIC, use one of the following commands on the fitted model object "z".

AIC(z, k = log(n))   # n is the sample size
BIC(z)

What matters is not BIC itself but the difference in BIC between models. The model with the smallest BIC will have a difference of zero. These differences can be used to compute the BIC weights. The model with the highest weight is the "best" model in the set.

x <- c(bic1, bic2)           # stores two BIC values in a vector
delta <- x - min(x)          # BIC differences
L <- exp(-0.5 * delta)       # 
w <- L/sum(L)                # BIC weights

BIC weights can be interpreted as posterior probabilities of the models, assuming that one of the models in the set is the "true" model (whatever that means), if all the models in the set have equal prior probabilities.

All subsets regression

The leaps command in the leaps library carries out an exhaustive search for the best subsets of the explanatory variables for predicting a response variable. The method is only suited to analyzing numeric variables, not categorical variables (a variable having just two categories can be recoded as a numeric variable with indicator 0's and 1's). To begin, load the leaps library,

library(leaps)

To prepare your data for analyzing, you will need to create a new data frame that contains only the explanatory variables you wish to include in the analysis. Don't include the response variable or other variables in that data frame. Any interaction terms to fit will need to be added manually.

x <- mydata[ ,c(2,5,6,7,11)]
x$x1x2 <- x$x1 * x$x2        # interaction between x1 and x2

Then run leaps. Models are compared using Mallow's Cp rather than AIC, but the two quantities are related. We can compute AIC ourselves later.

z <- leaps(x, y, names=names(x))   # note that x comes before y

By default, the resulting object (here called z) will contain the results for the 10 best models for each value of p, the number of predictors, according to the criterion of smallest Cp. To increase the number of models saved for each value of p to 20, for example, modify the nbest option.

z<-leaps(x, y, names = names(x), nbest = 20)

Additional useful commands:

plot(z$size, z$Cp)          # plot Cp vs number of predictors, p
lines(z$size, z$size)       # adds the line for Cp = p
i <- which(z$Cp==min(z$Cp)) # finds the model with smallest Cp
vars <- which(z$which[i,])  # id variables of best model

At the end of this Rtips section on model selection below is a home-made (and rather slow) leaps2aic command that will calculate additional quantities of interest from the leaps output. To use it, copy and paste the text into your R command window (you will need to do this only once in a session). These calculations assume that the the set of saved models from leaps includes that model having the lowest AICc and AIC, which is likely. However, the leaps output might not include all runners up (to be safe, increase the number of models saved using the nbest option in leaps). leaps2aic computes differences in AIC (AIC relative to the best) as well as differences in AICc (AIC adjusted for small sample size) and BIC (Bayesian Information Criterion).

To use the leaps2aic command, make certain that the arguments x, y and z are exactly the same as the ones you passed to the leaps command immediately beforehand. y is the response variable, x is the data frame containing the explanatory variables, and z is the saved output object from the leaps command.

z1 <- leaps2aic(x, y, z)
z1[,c("model","AICc", "AIC")] # prints model and AIC differences

The leaps package also includes the command regsubsets, which has the advantage of a formula method, but the output is not as functional.

z <- regsubsets(y ~ x1 + x2 + x3, data = mydata)
summary(z)
plot(z)

stepAIC

stepAIC is a command in the MASS library that will automatically carry out a restricted search for the "best" model, as measured by either AIC or BIC (Bayesian Information criterion) (but not AICc, unfortunately). stepAIC accepts both categorical and numerical variables. It does not search all possible subsets of variables, but rather it searches in a specific sequence determined by the order of the variables in the formula of the model object. (This can be a problem if a design is severely unbalanced, in which case it may be necessary to try more than one sequence.) The search sequence obeys certain restrictions regarding inclusion of interactions and main effects, such as

  • A quadratic term x2 is not fitted unless x is also present in the model
  • The interaction a:b is not fitted unless both a and b are also present
  • The threeway interaction a:b:c is not fitted unless all two-way interactions of a, b, c, are also present
  • and so on. stepAIC focuses on finding the best model, and it does not yield a list of all the models that are nearly equivalent to the best.

To use stepAIC you will need to load the MASS library

library(MASS)

To begin, fit a "full" model with all the variables of interest included. It can be quite a chore, and take multiple lines of script, to write down all the variables in the model. To simplify this procedure, put the single response variable and all the explanatory variables you want to investigate into a new data frame. Leave out all variables you are not interested in fitting. There's no need to include variables representing interactions. Then you can use one of the following shortcuts

z <- lm(y ~., data = mydata)     # additive model, all variables
z <- lm(y ~(.)^2, data = mydata) # includes 2-way interactions
z <- lm(y ~(.)^3, data = mydata) # includes 3-way interactions

To use stepAIC, include the fitted model output from the previous lm command,

z1 <- stepAIC(z, upper=~., lower=~1, direction="both")

The output will appear on the screen. The output object (here called z1) stores only lm output from fitting the best model identified by stepAIC. To obtain results from fitting the best model, use the usual commands for lm objects, e.g.,

summary(z1)
plot(z1)
anova(z1)

Now, to interpret the output of stepAIC. Start at the top and follow the sequence of steps that the procedure has carried out. First, it fits the full model you gave it and prints the AIC value. Then it lists the AIC values that result when it drops out each term, one at a time, leaving all other variables in the model. (It starts with the highest-order interactions, if you have included any.) It picks the best model of the bunch tested (the one with the lowest AIC) and then starts again. At each iteration it may also add variables one at a time to see if AIC is lowered further. The process continues until neither adding nor removing any single term results in a lower AIC.

To use the BIC criterion with stepAIC, change the option "k" to the log of the sample size

z1 <- stepAIC(z, upper = ~., lower = ~1, direction = "both", k = log(nrow(mydata)))

Here is the leaps2aic function.

--------------
leaps2aic <- function(x,y,z){
  # Calculate AIC and AICc for every model stored in "z", the
  # output of "leaps". This version sorts the output by AICc.
  # Requires the same x and y objects inputted to "leaps"
  rows <- 1:length(z$Cp)
  for(i in 1:nrow(z$which)){
    vars <- which(z$which[i,])
    newdat <- cbind.data.frame(y, x[,vars])
    znew <- lm(y~., data = newdat)
    L <- logLik(znew)
    aL <- attributes(L)
    z$model[i] <- paste(z$label[-1][vars],collapse=" + ")
    z$p[i] <- z$size[i]  # terms in model, including intercept
    z$k[i] <- aL$df      # parameters (including sigma^2)
    z$AIC[i] <- -2*L[1] + 2*aL$df
    z$AICc[i] <- -2*L[1] + 2*aL$df + (2*aL$df*(aL$df+1))/(aL$nobs-aL$df-1)
    z$BIC[i] <- -2*L[1] + log(aL$nobs)*aL$df
    }
  z$AIC <- z$AIC - min(z$AIC)
  z$AICc <- z$AICc - min(z$AICc)
  z$BIC <- z$BIC - min(z$BIC)
  result <- do.call("cbind.data.frame",
             list(z[c("model","p","k","Cp","AIC","AICc","BIC")], 
             stringsAsFactors=FALSE))
  result <- result[order(result$AICc),]
  rownames(result) <- 1:nrow(result)
  newdat <- cbind.data.frame(y,x)
  print(summary( lm( as.formula(paste("y ~ ",result$model[1])), 
             data=cbind.data.frame(y, x) )))
  return(result)
  }

Fit nonlinear models to data

The function nls in R allows fitting of nonlinear relationships between a response variable and one or more explanatory variables. Table 20.1 in Crawley's The R book lists a variety of nonlinear functions useful in biology. Examples of nonlinear functions include the power function,

y = axb                # when y intercept is 0
y = a + bxc            # when y intercept is not zero

the exponential function,

y = a(1 - e-bx)         # basic version
y = a - be-cx           # similar, but more flexible

and the Michaelis-Menten function,

y = ax/(b+x)           # when y intercept is 0
y = a + bx/(c+x)       # when y intercept is not zero

Unlike lm, nls requires that the formula statement includes all parameters, including an intercept if you want one fitted. Whereas in lm you would write the formula for linear regression as

y ~ x

to fit the same model in "nls" you would instead write the formula as

y ~ a + b*x

which includes "a" for the intercept and "b" for the slope. Including them in the formula tells the nls function what to estimate.

To fit a power function to the data, the complete "nls" command would be

z <- nls(y ~ a * x^b, data = mydata, start = list(a=1, b=1))

Note that nls requires that you provide reasonable starting values for the parameters, in this case "a" and "b", so that it knows where to begin searching. The program will then run through a series of iterations that bring it closer and closer to the least squares estimates of the parameter values. Your starting values don't have to be terribly close, but the program is more likely to succeed in finding the least squares estimates the closer you are.

The resulting object (which I've named z) is an nls 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:

summary(z)             # parameter estimates and overall model fit
coef(z)                # model coefficients (means, slopes, intercepts)
confint(z)             # confidence intervals for parameters

resid(z)               # residuals
fitted(z)              # predicted values

predict(z, newdata=x1) # predicted values for new observations
anova(z1, z2)          # compare fits of 2 models, "full" vs "reduced"
logLik(z)              # log-likelihood of the parameters
AIC(z)                 # Akaike Information Criterion
BIC(z                  # Bayesian Information Criterion

If nls has difficulty finding the least squares estimates, it may stop before it has succeeded. For example, it will stop and print an error if it runs out of iterations before it has converged on the solution (by default, the maximum number of iterations allowed is 50). You can help it along by increasing the number of iterations using the "control" option,

z <- nls(y ~ a * x^b, data = mydata, start = list(a=1,b=1),
               control = list(maxiter=200, warnOnly=TRUE)))

The control option warnOnly=TRUE ensures that when the program stops prematurely, it prints a warning rather than an error message. The difference is that with a warning you get to evaluate the fit of the not-yet-converged model stored in the model object (here stored in "z"). Type ?nls.control to see other options.


Fit a cubic spline

The function "gam" allows estimation of a function, f, relating a response variable Y to a continuous explanatory variable X,

Y = f(X) + random error.

For example, data on individual fitness, such as survival or reproductive success (Y), may be obtained to estimate the form of natural selection f on a trait X. In this case f is a "fitness function".

A cubic spline can be used to estimate f, which may take any smooth form that the data warrant. The cubic spline is a type of generalized additive model (gam). As in the case of generalized linear models, the random errors may have a distribution that is normal, binomial, Poisson, or otherwise. The advantage of the spline method is that it is flexible and avoids the need to make prior assumptions about the shape of f. See Schluter (1988, Evolution 42: 849–861; reprint) for a more detailed explanation of its use in estimating a fitness function.

Key to fitting the data with a cubic spline is the choice of smoothing parameter that controls the "roughness" of the fitted curve. A small value for the smoothing parameter (close to 0) will lead to a very wiggly curve that passes close to each of the data points, whereas a large value for the smoothing parameter is equivalent to fitting a straight line through the data. In the procedure used here, the smoothing parameter is chosen as that value minimizing the generalized cross-validation (GCV) score. Minimizing this score is tantamount to maximizing the predictive ability of the fitted model.

To begin, load the built-in mgcv library in R,

library(mgcv)

The following command fits a cubic spline to the relationship between a response variable y and a single continuous variable x, both of which are columns of the data frame, mydata. The method will find and fit the curve with a smoothing parameter that minimizes the GCV score. If the response variable is binary, such as data on survival (0 or 1), then a logit link function should be used.

z <- gam(y ~ s(x), data = mydata, family = binomial(link = "logit"),
                      method = "GCV.Cp")

Use instead family=poisson(link = "log") for Poisson data and family=gaussian(link = "identity") for normal errors.

The results are stored in the fitted "gam" object, here called "z". To extract results use the following commands.

summary(z)   # regression coefficients
fitted(z)    # predicted (fitted) values on original scale
predict(z)   # predicted values on the logit (or log, etc) scale
plot(z)      # the spline on the logit (or log, etc) scale
residuals(z, type = "response")  # residuals

Other useful results stored in the gam object:

z$sp         # "best" value for the smoothing parameter
z$gcv.ubre   # the corresponding (minimum) GCV score
sum(z$hat)   # effective number of parameters of fitted curve

To plot predicted values with Bayesian standard errors (plotted on the logit or log scale if residuals are binomial or Poisson).

plot(z, se = TRUE, seWithMean=TRUE) # plus or minus 2 standard error
plot(z, se = 1, seWithMean=TRUE)    # plus or minus 1 standard error

To plot predicted values with standard errors on the original scale, in the case of non-normal errors, supply the appropriate inverse of the link function. For example,

plot(z, se = 1, seWithMean = TRUE, rug = FALSE, shift = mean(predict(z)),
        trans = exp)                             # Poisson data
plot(z, se = 1, seWithMean = TRUE, rug = FALSE, shift = mean(predict(z)),
        trans = function(x){exp(x)/(1+exp(x))})  # binomial data

To plot predicted values for nicely spaced x-values, not necessarily the original x-values, with 1 standard error for prediction above and below, use this modified procedure instead. The commands below assume that your variables are named x and y and are in a data frame named mydata.

newx <- seq(min(mydata$x), max(mydata$x), length.out = 100)
z1 <- predict(z, newdata=list(x = newx), se.fit = TRUE)
invlogit <- function(x){exp(x)/(exp(x) + 1)}
yhat <- invlogit(z1$fit)
upper <- invlogit(z1$fit + z1$se.fit)
lower <- invlogit(z1$fit - z1$se.fit)
plot(newx, yhat, type="l", ylim = c(0,1))
lines(newx, upper, lty = 2)
lines(newx, lower, lty = 2)

What is the relationship between the choice of smoothing parameter (lambda) and the GCV score? The following commands will plot the relationship for binary data.

lambda <- exp(seq(-4,0, by=.05))        # fit a range of lambdas >0
gcvscore <- sapply(lambda, function(lambda, mydata){
      gam(y ~ s(x), data = mydata, family = binomial, 
      sp = lambda, method="GCV.Cp")$gcv.ubre},mydata)

plot(lambda, gcvscore, type = "l") # or
plot(log(lambda), gcvscore, type = "l")